diff --git a/acb_hypgeom/coulomb.c b/acb_hypgeom/coulomb.c index 394c5ef6..e6497f23 100644 --- a/acb_hypgeom/coulomb.c +++ b/acb_hypgeom/coulomb.c @@ -147,11 +147,7 @@ _acb_hypgeom_coulomb(acb_t F, acb_t G, acb_t Hpos, acb_t Hneg, const acb_t l, co } } - /* h = log(2z) */ - acb_mul_2exp_si(h, z, 1); - acb_log(h, h, prec); - - /* C = exp(h l + (-pi eta + lu + lv)/2) */ + /* C = exp((-pi eta + lu + lv)/2) */ acb_const_pi(C, prec); acb_mul(C, C, eta, prec); acb_neg(C, C); @@ -169,7 +165,6 @@ _acb_hypgeom_coulomb(acb_t F, acb_t G, acb_t Hpos, acb_t Hneg, const acb_t l, co } acb_mul_2exp_si(C, C, -1); - acb_addmul(C, h, l, prec); if (asymp) { @@ -214,7 +209,27 @@ _acb_hypgeom_coulomb(acb_t F, acb_t G, acb_t Hpos, acb_t Hneg, const acb_t l, co acb_hypgeom_m(F, v, m, z1, 1, prec); } - acb_exp(C, C, prec); + if (acb_contains_zero(z)) + { + acb_exp(C, C, prec); + /* (2z)^l without logarithm */ + acb_mul_2exp_si(h, z, 1); + acb_pow(h, h, l, prec); + acb_mul(C, C, h, prec); + + /* h = log(2z) */ + acb_indeterminate(h); + } + else + { + /* h = log(2z) */ + acb_mul_2exp_si(h, z, 1); + acb_log(h, h, prec); + + acb_addmul(C, h, l, prec); + acb_exp(C, C, prec); + } + acb_mul(F, F, C, prec); if (F_real) arb_zero(acb_imagref(F)); diff --git a/acb_hypgeom/test/t-coulomb.c b/acb_hypgeom/test/t-coulomb.c index 82cb57c1..f0353d6f 100644 --- a/acb_hypgeom/test/t-coulomb.c +++ b/acb_hypgeom/test/t-coulomb.c @@ -87,6 +87,13 @@ int main() mask = n_randlimb(state); + if (n_randint(state, 2) == 0) + { + acb_randtest_param(t, state, 1 + n_randint(state, 200), 1 + n_randint(state, 100)); + acb_add(z, z, t, prec2); + acb_sub(z, z, t, prec2); + } + /* also test aliasing */ acb_set(G2, z);