hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-04 17:01:40 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

Coulomb wave functions

Browse source
This commit is contained in:
fredrik 2019-02-26 15:12:19 +01:00
parent 3f2a79b7d5
commit 8223465324
8 changed files with 1090 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

4
acb_hypgeom.h
View file

@ -137,6 +137,10 @@ void acb_hypgeom_airy_jet(acb_ptr ai, acb_ptr bi, const acb_t z, slong len, slon
void _acb_hypgeom_airy_series(acb_ptr ai, acb_ptr ai_prime, acb_ptr bi, acb_ptr bi_prime, acb_srcptr z, slong zlen, slong len, slong prec);
void acb_hypgeom_airy_series(acb_poly_t ai, acb_poly_t ai_prime, acb_poly_t bi, acb_poly_t bi_prime, const acb_poly_t z, slong len, slong prec);
void acb_hypgeom_coulomb(acb_t F, acb_t G, acb_t Hpos, acb_t Hneg, const acb_t l, const acb_t eta, const acb_t z, slong prec);
void acb_hypgeom_coulomb_jet(acb_ptr F, acb_ptr G, acb_ptr Hpos, acb_ptr Hneg, const acb_t l, const acb_t eta, const acb_t z, slong len, slong prec);
void acb_hypgeom_coulomb_series(acb_poly_t F, acb_poly_t G, acb_poly_t Hpos, acb_poly_t Hneg, const acb_t l, const acb_t eta, const acb_poly_t z, slong len, slong prec);
void acb_hypgeom_gamma_upper_asymp(acb_t res, const acb_t s, const acb_t z, int modified, slong prec);
void acb_hypgeom_gamma_upper_1f1a(acb_t res, const acb_t s, const acb_t z, int modified, slong prec);
void acb_hypgeom_gamma_upper_1f1b(acb_t res, const acb_t s, const acb_t z, int modified, slong prec);

346
acb_hypgeom/coulomb.c Normal file
View file

@ -0,0 +1,346 @@
/*
Copyright (C) 2019 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "acb_hypgeom.h"
static void
acb_hypgeom_coulomb_is_real(int * C, int * F, int * G, const acb_t l1, const acb_t eta, const acb_t z)
{
*C = *F = *G = 0;
if (acb_is_real(l1) && acb_is_real(eta))
{
if (arb_is_positive(acb_realref(l1)) || arb_is_nonzero(acb_realref(eta)))
{
*C = 1;
}
if (acb_is_real(z))
{
if (arb_is_positive(acb_realref(z)))
{
*F = *G = 1;
}
if (acb_is_int(l1))
*F = 1;
}
}
}
void
_acb_hypgeom_coulomb(acb_t F, acb_t G, acb_t Hpos, acb_t Hneg, const acb_t l, const acb_t eta, const acb_t z, int asymp, slong prec)
{
acb_t u, v, lu, lv, z1, z2, m, h, T1, T2, U1, U2, H1, H2, C, theta;
int C_real, F_real, G_real;
int want_U1, want_U2, cut;
acb_init(u); acb_init(v); acb_init(lu); acb_init(lv);
acb_init(z1); acb_init(z2); acb_init(m); acb_init(h);
acb_init(T1); acb_init(T2); acb_init(U1); acb_init(U2);
acb_init(H1); acb_init(H2); acb_init(C); acb_init(theta);
acb_indeterminate(U1);
acb_indeterminate(U2);
/* z1 = 2iz, z2 = -2iz, */
acb_mul_onei(z1, z);
acb_mul_2exp_si(z1, z1, 1);
acb_neg(z2, z1);
if (asymp == -1)
asymp = acb_hypgeom_u_use_asymp(z1, prec);
/* Need the union of both sides of the branch cut for G, H+, H-. */
if (arb_is_nonnegative(acb_imagref(z)) || arb_is_positive(acb_realref(z)))
cut = 0;
else
cut = 1;
want_U1 = want_U2 = 0;
if (asymp)
{
want_U1 = want_U2 = 1;
}
else
{
if (G != NULL || Hpos != NULL || Hneg != NULL)
{
if (arf_sgn(arb_midref(acb_imagref(z))) >= 0)
want_U1 = 1;
else
want_U2 = 1;
if (cut)
want_U1 = want_U2 = 1;
}
}
/* m = l+1 */
acb_add_ui(m, l, 1, prec);
acb_hypgeom_coulomb_is_real(&C_real, &F_real, &G_real, m, eta, z);
/* u = 1+l+i eta, v = 1+l-i eta */
acb_mul_onei(u, eta);
acb_add(u, u, m, prec);
acb_div_onei(v, eta);
acb_add(v, v, m, prec);
/* lu = lgamma(u), v = lgamma(v) */
acb_lgamma(lu, u, prec);
if (C_real)
acb_conj(lv, lu);
else
acb_lgamma(lv, v, prec);
/* m = 2l+2 */
acb_mul_2exp_si(m, m, 1);
if (asymp)
{
if (want_U1 && want_U2 && G_real)
{
acb_hypgeom_u_asymp(U1, u, m, z2, -1, prec);
acb_conj(U2, U1);
}
else
{
if (want_U1) acb_hypgeom_u_asymp(U1, u, m, z2, -1, prec);
if (want_U2) acb_hypgeom_u_asymp(U2, v, m, z1, -1, prec);
}
}
else
{
if (want_U1 && want_U2 && G_real)
{
acb_hypgeom_u(U1, u, m, z2, prec);
acb_pow(h, z2, u, prec);
acb_mul(U1, U1, h, prec);
acb_conj(U2, U1);
}
else
{
if (want_U1)
{
acb_hypgeom_u(U1, u, m, z2, prec);
acb_pow(h, z2, u, prec);
acb_mul(U1, U1, h, prec);
}
if (want_U2)
{
acb_hypgeom_u(U2, v, m, z1, prec);
acb_pow(h, z1, v, prec);
acb_mul(U2, U2, h, prec);
}
}
}
/* h = log(2z) */
acb_mul_2exp_si(h, z, 1);
acb_log(h, h, prec);
/* C = exp(h l + (-pi eta + lu + lv)/2) */
acb_const_pi(C, prec);
acb_mul(C, C, eta, prec);
acb_neg(C, C);
if (C_real)
{
acb_mul_2exp_si(T1, lu, 1);
arb_zero(acb_imagref(T1));
acb_add(C, C, T1, prec);
}
else
{
acb_add(C, C, lu, prec);
acb_add(C, C, lv, prec);
}
acb_mul_2exp_si(C, C, -1);
acb_addmul(C, h, l, prec);
if (asymp)
{
/* T1 = exp(-(-iz + lv + u log(z1)) U1 */
/* T2 = exp(-(+iz + lu + v log(z2)) U2 */
acb_log(T1, z1, prec);
acb_mul(T1, T1, u, prec);
acb_add(T1, T1, lv, prec);
acb_mul_2exp_si(z1, z1, -1);
acb_sub(T1, T1, z1, prec);
acb_mul_2exp_si(z1, z1, 1);
acb_neg(T1, T1);
acb_exp(T1, T1, prec);
acb_mul(T1, T1, U1, prec);
if (F_real)
{
acb_mul_2exp_si(F, T1, 1);
arb_zero(acb_imagref(F));
}
else
{
acb_log(T2, z2, prec);
acb_mul(T2, T2, v, prec);
acb_add(T2, T2, lu, prec);
acb_mul_2exp_si(z2, z2, -1);
acb_sub(T2, T2, z2, prec);
acb_mul_2exp_si(z2, z2, 1);
acb_neg(T2, T2);
acb_exp(T2, T2, prec);
acb_mul(T2, T2, U2, prec);
/* F = (T1 + T2) z C */
acb_add(F, T1, T2, prec);
}
}
else
{
/* C *= exp(-iz) */
acb_div_onei(F, z);
acb_add(C, C, F, prec);
acb_hypgeom_m(F, v, m, z1, 1, prec);
}
acb_exp(C, C, prec);
acb_mul(F, F, C, prec);
if (F_real)
arb_zero(acb_imagref(F));
acb_mul(F, F, z, prec);
if (G != NULL || Hpos != NULL || Hneg != NULL)
{
/* theta = z - eta h - 0.5 l pi + (lu - lv) / (2i) */
acb_sub(theta, lu, lv, prec);
acb_div_onei(theta, theta);
acb_mul_2exp_si(theta, theta, -1);
acb_const_pi(H1, prec);
acb_mul_2exp_si(H1, H1, -1);
acb_submul(theta, H1, l, prec);
acb_submul(theta, eta, h, prec);
acb_add(theta, theta, z, prec);
/* H1 = exp(+i theta) U1, H2 = exp(-i theta) U2 */
acb_mul_onei(H1, theta);
acb_exp_invexp(H1, H2, H1, prec);
acb_mul(H1, H1, U1, prec);
acb_mul(H2, H2, U2, prec);
if (G != NULL)
{
if (asymp && arb_is_positive(acb_realref(z)))
{
if (G_real)
{
if (arf_sgn(arb_midref(acb_imagref(z))) >= 0)
acb_set(G, H1);
else
acb_set(G, H2);
arb_zero(acb_imagref(G));
}
else
{
acb_add(G, H1, H2, prec);
acb_mul_2exp_si(G, G, -1);
}
}
else
{
acb_div_onei(u, F);
acb_add(u, H1, u, prec);
acb_mul_onei(v, F);
acb_add(v, H2, v, prec);
if (cut)
acb_union(G, u, v, prec);
else if (arf_sgn(arb_midref(acb_imagref(z))) >= 0)
acb_set(G, u);
else
acb_set(G, v);
if (G_real)
arb_zero(acb_imagref(G));
}
}
if (Hpos != NULL)
{
acb_set(u, H1);
acb_mul_onei(v, F);
acb_mul_2exp_si(v, v, 1);
acb_add(v, H2, v, prec);
if (cut)
acb_union(Hpos, u, v, prec);
else if (arf_sgn(arb_midref(acb_imagref(z))) >= 0)
acb_set(Hpos, u);
else
acb_set(Hpos, v);
if (G_real)
arb_set(acb_imagref(Hpos), acb_realref(F));
}
if (Hneg != NULL)
{
acb_div_onei(u, F);
acb_mul_2exp_si(u, u, 1);
acb_add(u, H1, u, prec);
acb_set(v, H2);
if (cut)
acb_union(Hneg, u, v, prec);
else if (arf_sgn(arb_midref(acb_imagref(z))) >= 0)
acb_set(Hneg, u);
else
acb_set(Hneg, v);
if (G_real)
arb_neg(acb_imagref(Hneg), acb_realref(F));
}
}
acb_clear(u); acb_clear(v); acb_clear(lu); acb_clear(lv);
acb_clear(z1); acb_clear(z2); acb_clear(m); acb_clear(h);
acb_clear(T1); acb_clear(T2); acb_clear(U1); acb_clear(U2);
acb_clear(H1); acb_clear(H2); acb_clear(C); acb_clear(theta);
}
void
acb_hypgeom_coulomb(acb_t F, acb_t G, acb_t Hpos, acb_t Hneg, const acb_t l, const acb_t eta, const acb_t z, slong prec)
{
/* We always compute F. Also handle aliasing. */
acb_t F2, l2, eta2, z2;
acb_init(F2);
acb_init(l2);
acb_init(eta2);
acb_init(z2);
acb_set(l2, l);
acb_set(eta2, eta);
acb_set(z2, z);
_acb_hypgeom_coulomb(F2, G, Hpos, Hneg, l2, eta2, z2, -1, prec);
if (F != NULL)
acb_set(F, F2);
acb_clear(F2);
acb_clear(l2);
acb_clear(eta2);
acb_clear(z2);
}

242
acb_hypgeom/coulomb_jet.c Normal file
View file

@ -0,0 +1,242 @@
/*
Copyright (C) 2019 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "acb_hypgeom.h"
static void
_acb_hypgeom_coulomb_jet(acb_ptr F, acb_ptr G, acb_ptr Hpos, acb_ptr Hneg, const acb_t l, const acb_t eta, const acb_t z, slong len, slong prec)
{
acb_t l1, t, R, S;
if (len <= 0)
return;
if (len == 1)
{
acb_hypgeom_coulomb(F, G, Hpos, Hneg, l, eta, z, prec);
return;
}
acb_init(l1);
acb_init(t);
acb_init(R);
acb_init(S);
acb_add_ui(l1, l, 1, prec);
acb_hypgeom_coulomb(F, G, Hpos, Hneg, l, eta, z, prec);
/* todo: somehow recycle the gamma function values for the two evaluations? */
acb_hypgeom_coulomb((F == NULL) ? NULL : (F + 1),
(G == NULL) ? NULL : (G + 1),
(Hpos == NULL) ? NULL : (Hpos + 1),
(Hneg == NULL) ? NULL : (Hneg + 1),
l1, eta, z, prec);
/* R_l = (2l+1) C_l / C_{l+1} = sqrt(l+i eta) sqrt(l-i eta) / l */
if (acb_is_real(l) && acb_is_real(eta) && arb_is_nonzero(acb_realref(eta)))
{
acb_mul(R, l1, l1, prec);
acb_addmul(R, eta, eta, prec);
acb_sqrt(R, R, prec);
}
else
{
acb_mul_onei(t, eta);
acb_add(t, t, l1, prec);
acb_sqrt(t, t, prec);
acb_div_onei(R, eta);
acb_add(R, R, l1, prec);
acb_sqrt(R, R, prec);
acb_mul(R, R, t, prec);
}
acb_div(R, R, l1, prec);
acb_div(S, l1, z, prec);
acb_div(t, eta, l1, prec);
acb_add(S, S, t, prec);
/* todo: fix regular F at origin */
if (F != NULL)
{
acb_mul(F + 1, F + 1, R, prec);
acb_neg(F + 1, F + 1);
acb_addmul(F + 1, F, S, prec);
}
if (G != NULL)
{
acb_mul(G + 1, G + 1, R, prec);
acb_neg(G + 1, G + 1);
acb_addmul(G + 1, G, S, prec);
}
if (Hpos != NULL)
{
acb_mul(Hpos + 1, Hpos + 1, R, prec);
acb_neg(Hpos + 1, Hpos + 1);
acb_addmul(Hpos + 1, Hpos, S, prec);
}
if (Hneg != NULL)
{
acb_mul(Hneg + 1, Hneg + 1, R, prec);
acb_neg(Hneg + 1, Hneg + 1);
acb_addmul(Hneg + 1, Hneg, S, prec);
}
if (len >= 3)
{
acb_t q, q2, w, w2;
acb_init(q);
acb_init(q2);
acb_init(w);
acb_init(w2);
acb_inv(w, z, prec);
acb_mul(w2, w, w, prec);
/* F''/2 = q F, q = (2eta/z + l(l+1)/z^2 - 1)/2 */
acb_mul(q, l, l1, prec);
acb_mul(q, q, w2, prec);
acb_mul_2exp_si(q2, eta, 1);
acb_addmul(q, q2, w, prec);
acb_sub_ui(q, q, 1, prec);
acb_mul_2exp_si(q, q, -1);
if (F != NULL) acb_mul(F + 2, F, q, prec);
if (G != NULL) acb_mul(G + 2, G, q, prec);
if (Hpos != NULL) acb_mul(Hpos + 2, Hpos, q, prec);
if (Hneg != NULL) acb_mul(Hneg + 2, Hneg, q, prec);
/* F'''/6 = (2qF' - q2 F)/6, q2 = 2(eta + l(l+1)/z)/z^2 */
if (len >= 4)
{
acb_mul_2exp_si(q, q, 1);
acb_mul(q2, l, l1, prec);
acb_mul(q2, q2, w, prec);
acb_add(q2, q2, eta, prec);
acb_mul_2exp_si(q2, q2, 1);
acb_mul(q2, q2, w2, prec);
if (F != NULL)
{
acb_mul(F + 3, F + 1, q, prec);
acb_submul(F + 3, F + 0, q2, prec);
acb_div_ui(F + 3, F + 3, 6, prec);
}
if (G != NULL)
{
acb_mul(G + 3, G + 1, q, prec);
acb_submul(G + 3, G, q2, prec);
acb_div_ui(G + 3, G + 3, 6, prec);
}
if (Hpos != NULL)
{
acb_mul(Hpos + 3, Hpos + 1, q, prec);
acb_submul(Hpos + 3, Hpos, q2, prec);
acb_div_ui(Hpos + 3, Hpos + 3, 6, prec);
}
if (Hneg != NULL)
{
acb_mul(Hneg + 3, Hneg + 1, q, prec);
acb_submul(Hneg + 3, Hneg, q2, prec);
acb_div_ui(Hneg + 3, Hneg + 3, 6, prec);
}
}
if (len >= 5)
{
slong k;
acb_ptr s;
s = _acb_vec_init(4);
/*
s4 = -(k(k+7) + 12)
s3 = 2 (k(k+5) + 6) z
s2 = (k(k+3) + z^2 - 2z*eta - l(l+1) + 2)
s1 = 2(z-eta)
s0 = 1
F[k+4] = (s0*F[k] + s1*F[k+1] + s2*F[k+2] + s3*F[k+3]) / s4 / z^2
*/
acb_sub(s + 1, z, eta, prec);
acb_mul_2exp_si(s + 1, s + 1, 1);
acb_mul(q, z, z, prec);
acb_mul(q2, eta, z, prec);
acb_mul_2exp_si(q2, q2, 1);
acb_sub(q, q, q2, prec);
acb_submul(q, l, l1, prec);
for (k = 0; k + 4 < len; k++)
{
acb_mul_si(s + 3, z, 2 * (k * (k + 5) + 6), prec);
acb_add_si(s + 2, q, k * (k + 3) + 2, prec);
if (F != NULL)
{
acb_dot(F + k + 4, F + k, 0, F + k + 1, 1, s + 1, 1, 3, prec);
acb_div_si(F + k + 4, F + k + 4, -(k * (k + 7) + 12), prec);
acb_mul(F + k + 4, F + k + 4, w2, prec);
}
if (G != NULL)
{
acb_dot(G + k + 4, G + k, 0, G + k + 1, 1, s + 1, 1, 3, prec);
acb_div_si(G + k + 4, G + k + 4, -(k * (k + 7) + 12), prec);
acb_mul(G + k + 4, G + k + 4, w2, prec);
}
if (Hpos != NULL)
{
acb_dot(Hpos + k + 4, Hpos + k, 0, Hpos + k + 1, 1, s + 1, 1, 3, prec);
acb_div_si(Hpos + k + 4, Hpos + k + 4, -(k * (k + 7) + 12), prec);
acb_mul(Hpos + k + 4, Hpos + k + 4, w2, prec);
}
if (Hneg != NULL)
{
acb_dot(Hneg + k + 4, Hneg + k, 0, Hneg + k + 1, 1, s + 1, 1, 3, prec);
acb_div_si(Hneg + k + 4, Hneg + k + 4, -(k * (k + 7) + 12), prec);
acb_mul(Hneg + k + 4, Hneg + k + 4, w2, prec);
}
}
_acb_vec_clear(s, 4);
}
acb_clear(q);
acb_clear(q2);
acb_clear(w);
acb_clear(w2);
}
acb_clear(l1);
acb_clear(t);
acb_clear(R);
acb_clear(S);
}
void
acb_hypgeom_coulomb_jet(acb_ptr F, acb_ptr G, acb_ptr Hpos, acb_ptr Hneg, const acb_t l, const acb_t eta, const acb_t z, slong len, slong prec)
{
acb_t t; /* to allow aliasing with z */
acb_init(t);
acb_set(t, z);
_acb_hypgeom_coulomb_jet(F, G, Hpos, Hneg, l, eta, t, len, prec);
acb_clear(t);
}

129
acb_hypgeom/coulomb_series.c Normal file
View file

@ -0,0 +1,129 @@
/*
Copyright (C) 2019 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "acb_hypgeom.h"
void
_acb_hypgeom_coulomb_series(acb_ptr F, acb_ptr G,
acb_ptr Hpos, acb_ptr Hneg, const acb_t l, const acb_t eta,
acb_srcptr z, slong zlen, slong len, slong prec)
{
acb_ptr t, v;
if (len == 0)
return;
zlen = FLINT_MIN(zlen, len);
if (zlen == 1)
{
acb_hypgeom_coulomb(F, G, Hpos, Hneg, l, eta, z, prec);
if (F != NULL) _acb_vec_zero(F + 1, len - 1);
if (G != NULL) _acb_vec_zero(G + 1, len - 1);
if (Hpos != NULL) _acb_vec_zero(Hpos + 1, len - 1);
if (Hneg != NULL) _acb_vec_zero(Hneg + 1, len - 1);
return;
}
t = _acb_vec_init(len);
v = _acb_vec_init(zlen);
/* copy nonconstant part first to allow aliasing */
acb_zero(v);
_acb_vec_set(v + 1, z + 1, zlen - 1);
acb_hypgeom_coulomb_jet(F, G, Hpos, Hneg, l, eta, z, len, prec);
if (F != NULL)
{
_acb_vec_set(t, F, len);
_acb_poly_compose_series(F, t, len, v, zlen, len, prec);
}
if (G != NULL)
{
_acb_vec_set(t, G, len);
_acb_poly_compose_series(G, t, len, v, zlen, len, prec);
}
if (Hpos != NULL)
{
_acb_vec_set(t, Hpos, len);
_acb_poly_compose_series(Hpos, t, len, v, zlen, len, prec);
}
if (Hneg != NULL)
{
_acb_vec_set(t, Hneg, len);
_acb_poly_compose_series(Hneg, t, len, v, zlen, len, prec);
}
_acb_vec_clear(t, len);
_acb_vec_clear(v, zlen);
}
void
acb_hypgeom_coulomb_series(acb_poly_t F, acb_poly_t G,
acb_poly_t Hpos, acb_poly_t Hneg, const acb_t l, const acb_t eta,
const acb_poly_t z, slong len, slong prec)
{
acb_srcptr zptr;
slong zlen;
acb_t t;
if (len == 0)
{
if (F != NULL) acb_poly_zero(F);
if (G != NULL) acb_poly_zero(G);
if (Hpos != NULL) acb_poly_zero(Hpos);
if (Hneg != NULL) acb_poly_zero(Hneg);
return;
}
zlen = z->length;
if (zlen <= 1)
len = 1;
if (F != NULL) acb_poly_fit_length(F, len);
if (G != NULL) acb_poly_fit_length(G, len);
if (Hpos != NULL) acb_poly_fit_length(Hpos, len);
if (Hneg != NULL) acb_poly_fit_length(Hneg, len);
if (zlen == 0)
{
acb_init(t);
zptr = t;
zlen = 1;
}
else
{
zptr = z->coeffs;
}
_acb_hypgeom_coulomb_series(
F ? F->coeffs : NULL,
G ? G->coeffs : NULL,
Hpos ? Hpos->coeffs : NULL,
Hneg ? Hneg->coeffs : NULL,
l, eta,
zptr, zlen, len, prec);
if (F != NULL) _acb_poly_set_length(F, len);
if (G != NULL) _acb_poly_set_length(G, len);
if (Hpos != NULL) _acb_poly_set_length(Hpos, len);
if (Hneg != NULL) _acb_poly_set_length(Hneg, len);
if (F != NULL) _acb_poly_normalise(F);
if (G != NULL) _acb_poly_normalise(G);
if (Hpos != NULL) _acb_poly_normalise(Hpos);
if (Hneg != NULL) _acb_poly_normalise(Hneg);
}

160
acb_hypgeom/test/t-coulomb.c Normal file
View file

@ -0,0 +1,160 @@
/*
Copyright (C) 2019 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "acb_hypgeom.h"
int main()
{
slong iter;
flint_rand_t state;
flint_printf("coulomb....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 2000 * arb_test_multiplier(); iter++)
{
acb_t z, t, u, eta, l;
acb_t F1, G1, Hpos1, Hneg1;
acb_t F2, G2, Hpos2, Hneg2;
slong prec1, prec2;
unsigned int mask;
acb_init(z); acb_init(t); acb_init(u); acb_init(eta); acb_init(l);
acb_init(F1); acb_init(G1); acb_init(Hpos1); acb_init(Hneg1);
acb_init(F2); acb_init(G2); acb_init(Hpos2); acb_init(Hneg2);
prec1 = 2 + n_randint(state, 200);
prec2 = 2 + n_randint(state, 200);
acb_randtest_param(eta, state, 1 + n_randint(state, 200), 1 + n_randint(state, 10));
acb_randtest_param(l, state, 1 + n_randint(state, 200), 1 + n_randint(state, 10));
acb_randtest_param(z, state, 1 + n_randint(state, 200), 1 + n_randint(state, 100));
acb_randtest_param(z, state, 1 + n_randint(state, 200), 1 + n_randint(state, 100));
acb_randtest_param(t, state, 1 + n_randint(state, 200), 1 + n_randint(state, 100));
acb_add(z, z, t, 1000);
acb_sub(z, z, t, 1000);
acb_hypgeom_coulomb(F1, G1, Hpos1, Hneg1, l, eta, z, prec1);
acb_hypgeom_coulomb(F2, G2, Hpos2, Hneg2, l, eta, z, prec2);
if (!acb_overlaps(F1, F2) || !acb_overlaps(G1, G2) ||
!acb_overlaps(Hpos1, Hpos2) || !acb_overlaps(Hneg1, Hneg2))
{
flint_printf("FAIL: consistency\n\n");
flint_printf("l = "); acb_printn(l, 30, 0); flint_printf("\n\n");
flint_printf("eta = "); acb_printn(eta, 30, 0); flint_printf("\n\n");
flint_printf("z = "); acb_printn(z, 30, 0); flint_printf("\n\n");
flint_printf("F1 = "); acb_printn(F1, 30, 0); flint_printf("\n");
flint_printf("F2 = "); acb_printn(F2, 30, 0); flint_printf("\n\n");
flint_printf("G1 = "); acb_printn(G1, 30, 0); flint_printf("\n");
flint_printf("G2 = "); acb_printn(G2, 30, 0); flint_printf("\n\n");
flint_printf("Hpos1 = "); acb_printn(Hpos1, 30, 0); flint_printf("\n");
flint_printf("Hpos2 = "); acb_printn(Hpos2, 30, 0); flint_printf("\n\n");
flint_printf("Hneg1 = "); acb_printn(Hneg1, 30, 0); flint_printf("\n");
flint_printf("Hneg2 = "); acb_printn(Hneg2, 30, 0); flint_printf("\n\n");
flint_abort();
}
acb_mul_onei(t, F1);
acb_add(t, t, G1, prec1);
acb_div_onei(u, F1);
acb_add(u, u, G1, prec1);
if (!acb_overlaps(t, Hpos1) || !acb_overlaps(u, Hneg1))
{
flint_printf("FAIL: complex\n\n");
flint_printf("l = "); acb_printn(l, 30, 0); flint_printf("\n\n");
flint_printf("eta = "); acb_printn(eta, 30, 0); flint_printf("\n\n");
flint_printf("z = "); acb_printn(z, 30, 0); flint_printf("\n\n");
flint_printf("F1 = "); acb_printn(F1, 30, 0); flint_printf("\n");
flint_printf("G1 = "); acb_printn(G1, 30, 0); flint_printf("\n");
flint_printf("Hpos1 = "); acb_printn(Hpos1, 30, 0); flint_printf("\n");
flint_printf("Hneg1 = "); acb_printn(Hneg1, 30, 0); flint_printf("\n");
flint_printf("t = "); acb_printn(t, 30, 0); flint_printf("\n");
flint_printf("u = "); acb_printn(u, 30, 0); flint_printf("\n");
flint_abort();
}
mask = n_randlimb(state);
/* also test aliasing */
acb_set(G2, z);
acb_hypgeom_coulomb((mask & 1) ? F2 : NULL,
(mask & 2) ? G2 : NULL,
(mask & 4) ? Hpos2 : NULL,
(mask & 8) ? Hneg2 : NULL, l, eta, G2, prec2);
if (((mask & 1) && !acb_overlaps(F1, F2)) ||
((mask & 2) && !acb_overlaps(G1, G2)) ||
((mask & 4) && !acb_overlaps(Hpos1, Hpos2)) ||
((mask & 8) && !acb_overlaps(Hneg1, Hneg2)))
{
flint_printf("FAIL: consistency (mask)\n\n");
flint_printf("mask = %u\n\n", mask);
flint_printf("l = "); acb_printn(l, 30, 0); flint_printf("\n\n");
flint_printf("eta = "); acb_printn(eta, 30, 0); flint_printf("\n\n");
flint_printf("z = "); acb_printn(z, 30, 0); flint_printf("\n\n");
flint_printf("F1 = "); acb_printn(F1, 30, 0); flint_printf("\n");
flint_printf("F2 = "); acb_printn(F2, 30, 0); flint_printf("\n\n");
flint_printf("G1 = "); acb_printn(G1, 30, 0); flint_printf("\n");
flint_printf("G2 = "); acb_printn(G2, 30, 0); flint_printf("\n\n");
flint_printf("Hpos1 = "); acb_printn(Hpos1, 30, 0); flint_printf("\n");
flint_printf("Hpos2 = "); acb_printn(Hpos2, 30, 0); flint_printf("\n\n");
flint_printf("Hneg1 = "); acb_printn(Hneg1, 30, 0); flint_printf("\n");
flint_printf("Hneg2 = "); acb_printn(Hneg2, 30, 0); flint_printf("\n\n");
flint_abort();
}
/* Check F_{l-1} G_{l} - F_{l} G_{l-1} = l (l^2 + eta^2)^(1/2) */
acb_sub_ui(t, l, 1, prec2);
acb_hypgeom_coulomb(F2, G2, NULL, NULL, t, eta, z, prec2);
acb_mul_onei(t, eta);
acb_add(t, t, l, prec2);
acb_rsqrt(t, t, prec2);
acb_div_onei(u, eta);
acb_add(u, u, l, prec2);
acb_rsqrt(u, u, prec2);
acb_mul(u, t, u, prec2);
acb_mul(u, u, l, prec2);
acb_mul(t, F2, G1, prec2);
acb_submul(t, F1, G2, prec2);
if (!acb_overlaps(t, u))
{
flint_printf("FAIL: cross-product\n\n");
flint_printf("l = "); acb_printn(l, 30, 0); flint_printf("\n\n");
flint_printf("eta = "); acb_printn(eta, 30, 0); flint_printf("\n\n");
flint_printf("z = "); acb_printn(z, 30, 0); flint_printf("\n\n");
flint_printf("F1 = "); acb_printn(F1, 30, 0); flint_printf("\n");
flint_printf("G1 = "); acb_printn(G1, 30, 0); flint_printf("\n");
flint_printf("F2 = "); acb_printn(F2, 30, 0); flint_printf("\n");
flint_printf("G2 = "); acb_printn(G2, 30, 0); flint_printf("\n");
flint_printf("t = "); acb_printn(t, 30, 0); flint_printf("\n");
flint_printf("u = "); acb_printn(u, 30, 0); flint_printf("\n");
flint_abort();
}
acb_clear(z); acb_clear(t); acb_clear(u); acb_clear(eta); acb_clear(l);
acb_clear(F1); acb_clear(G1); acb_clear(Hpos1); acb_clear(Hneg1);
acb_clear(F2); acb_clear(G2); acb_clear(Hpos2); acb_clear(Hneg2);
}
flint_randclear(state);
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}

152
acb_hypgeom/test/t-coulomb_series.c Normal file
View file

@ -0,0 +1,152 @@
/*
Copyright (C) 2019 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "acb_hypgeom.h"
int main()
{
slong iter;
flint_rand_t state;
flint_printf("coulomb_series....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 2000 * arb_test_multiplier(); iter++)
{
acb_poly_t F, G, Hpos, Hneg, F2, G2, Hpos2, Hneg2, z, w, t, u;
acb_t c, l, eta;
slong n1, n2, prec1, prec2;
unsigned int mask;
acb_poly_init(F); acb_poly_init(G);
acb_poly_init(Hpos); acb_poly_init(Hneg);
acb_poly_init(F2); acb_poly_init(G2);
acb_poly_init(Hpos2); acb_poly_init(Hneg2);
acb_poly_init(z); acb_poly_init(w);
acb_poly_init(t);
acb_poly_init(u);
acb_init(c); acb_init(l); acb_init(eta);
prec1 = 2 + n_randint(state, 200);
prec2 = 2 + n_randint(state, 200);
n1 = n_randint(state, 8);
n2 = n_randint(state, 8);
acb_poly_randtest(F, state, 10, prec1, 10);
acb_poly_randtest(G, state, 10, prec1, 10);
acb_poly_randtest(Hpos, state, 10, prec1, 10);
acb_poly_randtest(Hneg, state, 10, prec1, 10);
acb_randtest_param(l, state, 1 + n_randint(state, 200), 1 + n_randint(state, 10));
acb_randtest_param(eta, state, 1 + n_randint(state, 200), 1 + n_randint(state, 10));
acb_poly_randtest(z, state, 1 + n_randint(state, 10), 1 + n_randint(state, 200), 10);
acb_hypgeom_coulomb_series(F, G, Hpos, Hneg, l, eta, z, n1, prec1);
acb_poly_derivative(t, F, prec1);
if (n_randint(state, 2))
acb_poly_set(u, G);
else if (n_randint(state, 2))
acb_poly_set(u, Hpos);
else
acb_poly_set(u, Hneg);
/* F' G - F G' = 1 */
acb_poly_mullow(w, t, u, FLINT_MAX(n1 - 1, 0), prec1);
acb_poly_derivative(u, u, prec1);
acb_poly_mullow(t, F, u, FLINT_MAX(n1 - 1, 0), prec1);
acb_poly_sub(w, w, t, prec1);
acb_poly_derivative(t, z, prec1);
acb_poly_truncate(t, FLINT_MAX(n1 - 1, 0));
if (!acb_poly_overlaps(w, t))
{
flint_printf("FAIL: wronskian, n1 = %wd\n\n", n1);
flint_printf("l = "); acb_printd(l, 30); flint_printf("\n\n");
flint_printf("eta = "); acb_printd(eta, 30); flint_printf("\n\n");
flint_printf("z = "); acb_poly_printd(z, 30); flint_printf("\n\n");
flint_printf("F = "); acb_poly_printd(F, 30); flint_printf("\n\n");
flint_printf("G = "); acb_poly_printd(G, 30); flint_printf("\n\n");
flint_printf("Hpos = "); acb_poly_printd(Hpos, 30); flint_printf("\n\n");
flint_printf("Hneg = "); acb_poly_printd(Hneg, 30); flint_printf("\n\n");
flint_printf("w = "); acb_poly_printd(w, 30); flint_printf("\n\n");
flint_printf("t = "); acb_poly_printd(t, 30); flint_printf("\n\n");
flint_abort();
}
mask = n_randlimb(state);
acb_poly_set(G2, z); /* for aliasing */
if (n_randint(state, 2))
{
acb_poly_randtest(u, state, 1 + n_randint(state, 10), 1 + n_randint(state, 200), 10);
acb_poly_add(G2, G2, u, prec2);
acb_poly_sub(G2, G2, u, prec2);
}
acb_hypgeom_coulomb_series((mask & 1) ? F2 : NULL,
(mask & 2) ? G2 : NULL,
(mask & 4) ? Hpos2 : NULL,
(mask & 8) ? Hneg2 : NULL, l, eta, G2, n2, prec2);
acb_poly_truncate(F, FLINT_MIN(n1, n2));
acb_poly_truncate(G, FLINT_MIN(n1, n2));
acb_poly_truncate(Hpos, FLINT_MIN(n1, n2));
acb_poly_truncate(Hneg, FLINT_MIN(n1, n2));
acb_poly_truncate(F2, FLINT_MIN(n1, n2));
acb_poly_truncate(G2, FLINT_MIN(n1, n2));
acb_poly_truncate(Hpos2, FLINT_MIN(n1, n2));
acb_poly_truncate(Hneg2, FLINT_MIN(n1, n2));
if (((mask & 1) && (!acb_poly_overlaps(F, F2))) ||
((mask & 2) && (!acb_poly_overlaps(G, G2))) ||
((mask & 4) && (!acb_poly_overlaps(Hpos, Hpos2))) ||
((mask & 8) && (!acb_poly_overlaps(Hneg, Hneg2))))
{
flint_printf("FAIL: consistency (mask)\n\n");
flint_printf("mask = %u\n\n", mask);
flint_printf("len1 = %wd, len2 = %wd\n\n", n1, n2);
flint_printf("l = "); acb_printd(l, 30); flint_printf("\n\n");
flint_printf("eta = "); acb_printd(eta, 30); flint_printf("\n\n");
flint_printf("z = "); acb_poly_printd(z, 30); flint_printf("\n\n");
flint_printf("F = "); acb_poly_printd(F, 30); flint_printf("\n\n");
flint_printf("F2 = "); acb_poly_printd(F2, 30); flint_printf("\n\n");
flint_printf("G = "); acb_poly_printd(G, 30); flint_printf("\n\n");
flint_printf("G2 = "); acb_poly_printd(G2, 30); flint_printf("\n\n");
flint_printf("Hpos = "); acb_poly_printd(Hpos, 30); flint_printf("\n\n");
flint_printf("Hpos2 = "); acb_poly_printd(Hpos2, 30); flint_printf("\n\n");
flint_printf("Hneg = "); acb_poly_printd(Hneg, 30); flint_printf("\n\n");
flint_printf("Hneg2 = "); acb_poly_printd(Hneg2, 30); flint_printf("\n\n");
flint_abort();
}
acb_poly_clear(F); acb_poly_clear(G);
acb_poly_clear(Hpos); acb_poly_clear(Hneg);
acb_poly_clear(F2); acb_poly_clear(G2);
acb_poly_clear(Hpos2); acb_poly_clear(Hneg2);
acb_poly_clear(z); acb_poly_clear(w);
acb_poly_clear(t);
acb_poly_clear(u);
acb_clear(c); acb_clear(l); acb_clear(eta);
}
flint_randclear(state);
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}

50
doc/source/acb_hypgeom.rst
View file

@ -582,6 +582,56 @@ simultaneously. Any of the four function values can be omitted by passing
truncated to length *len*. As with the other Airy methods, any of the
outputs can be *NULL*.
Coulomb wave functions
-------------------------------------------------------------------------------
Coulomb wave functions are solutions of the Coulomb wave equation
.. math ::
y'' + \left(1 - \frac{2 \eta}{z} - \frac{\ell(\ell+1)}{z}\right) y = 0
which is the radial Schrödinger equation for a charged particle in a
Coulomb potential `1/z`, where `\ell` is the orbital angular momentum and
`\eta` is the Sommerfeld parameter.
The standard solutions are named `F_{\ell}(\eta,z)` (regular
at the origin `z = 0`) and `G_{\ell}(\eta,z)` (irregular at the origin).
The irregular solutions
`H^{\pm}_{\ell}(\eta,z) = G_{\ell}(\eta,z) \pm i F_{\ell}(\eta,z)`
are also used.
Coulomb wave functions are special cases of confluent hypergeometric functions.
For details about the normalization constants and connection formulas
defining the different solutions, see
[DYF1999]_, [Gas2018]_, [Mic2007]_ or chapter 33 in [NIST2012]_.
In this implementation, we define the analytic continuations of all
the functions so that the branch cut with respect to *z* is placed on the
negative real axis.
The following methods optionally compute
`F_{\ell}(\eta,z), G_{\ell}(\eta,z), H^{+}_{\ell}(\eta,z), H^{-}_{\ell}(\eta,z)`
simultaneously. Any of the four function values can be omitted by passing
*NULL* for the unwanted output variables.
The redundant functions `H^{\pm}` are provided explicitly since taking
the linear combination of *F* and *G* suffers from cancellation in
parts of the complex plane.
.. function:: void acb_hypgeom_coulomb(acb_t F, acb_t G, acb_t Hpos, acb_t Hneg, const acb_t l, const acb_t eta, const acb_t z, slong prec)
Writes to *F*, *G*, *Hpos*, *Hneg* the values of the respective
Coulomb wave functions. Any of the outputs can be *NULL*.
.. function:: void acb_hypgeom_coulomb_jet(acb_ptr F, acb_ptr G, acb_ptr Hpos, acb_ptr Hneg, const acb_t l, const acb_t eta, const acb_t z, slong len, slong prec)
Writes to *F*, *G*, *Hpos*, *Hneg* the respective Taylor expansions of the
Coulomb wave functions at the point *z*, truncated to length *len*.
Any of the outputs can be *NULL*.
.. function:: void acb_hypgeom_coulomb_series(acb_poly_t F, acb_poly_t G, acb_poly_t Hpos, acb_poly_t Hneg, const acb_t l, const acb_t eta, const acb_poly_t z, slong len, slong prec)
Computes the Coulomb wave functions evaluated at the power series *z*,
truncated to length *len*. Any of the outputs can be *NULL*.
Incomplete gamma and beta functions
-------------------------------------------------------------------------------

7
doc/source/credits.rst
View file

@ -179,12 +179,16 @@ Bibliography
.. [Dup2006] \R. Dupont. "Moyenne arithmético-géométrique, suites de Borchardt et applications." These de doctorat, École polytechnique, Palaiseau (2006). http://http://www.lix.polytechnique.fr/Labo/Regis.Dupont/these_soutenance.pdf
.. [DYF1999] \A. Dzieciol, S. Yngve and P. O. Fröman, "Coulomb wave functions with complex values of the variable and the parameters", J. Math. Phys. 40, 6145 (1999), https://doi.org/10.1063/1.533083
.. [EHJ2016] \A. Enge, W. Hart and F. Johansson, "Short addition sequences for theta functions", preprint (2016), https://arxiv.org/abs/1608.06810
.. [EM2004] \O. Espinosa and V. Moll, "A generalized polygamma function", Integral Transforms and Special Functions (2004), 101-115.
.. [Fil1992] \S. Fillebrown, "Faster Computation of Bernoulli Numbers", Journal of Algorithms 13 (1992) 431-445
.. [Gas2018] \D. Gaspard, "Connection formulas between Coulomb wave functions" (2018), https://arxiv.org/abs/1804.10976
.. [GG2003] \J. von zur Gathen and J. Gerhard, *Modern Computer Algebra*, second edition, Cambridge University Press (2003)
.. [GVL1996] \G. H. Golub and C. F. Van Loan, *Matrix Computations*, third edition, Johns Hopkins University Press (1996).
@ -235,6 +239,8 @@ Bibliography
.. [Leh1970] \R. S. Lehman, "On the Distribution of Zeros of the Riemann Zeta-Function", Proc. of the London Mathematical Society 20(3) (1970), 303-320, https://doi.org/10.1112/plms/s3-20.2.303
.. [Mic2007] \N. Michel, "Precise Coulomb wave functions for a wide range of complex l, η and z", Computer Physics Communications, Volume 176, Issue 3, (2007), 232-249, https://doi.org/10.1016/j.cpc.2006.10.004
.. [Miy2010] \S. Miyajima, "Fast enclosure for all eigenvalues in generalized eigenvalue problems", Journal of Computational and Applied Mathematics, 233 (2010), 2994-3004, https://dx.doi.org/10.1016/j.cam.2009.11.048
.. [MPFR2012] The MPFR team, "MPFR Algorithms" (2012), http://www.mpfr.org/algo.html
@ -264,3 +270,4 @@ Bibliography
.. [Tru2014] \T. S. Trudgian, "An improved upper bound for the argument of the Riemann zeta-function on the critical line II", Journal of Number Theory 134 (2014), 280-292, https://doi.org/10.1016/j.jnt.2013.07.017
.. [Tur1953] \A. M. Turing, "Some Calculations of the Riemann Zeta-Function", Proc. of the London Mathematical Society 3(3) (1953), 99-117, https://doi.org/10.1112/plms/s3-3.1.99