From bd4359adeb13055bc5c6b34dbdc41173f776fba0 Mon Sep 17 00:00:00 2001 From: fredrik Date: Wed, 27 Feb 2019 20:29:59 +0100 Subject: [PATCH] allow series expansion of Coulomb F at z = 0 (not currently covered by test code) --- acb_hypgeom/coulomb_jet.c | 128 ++++++++++++++++++++++++++++ acb_hypgeom/test/t-coulomb_series.c | 11 ++- arb_hypgeom/test/t-coulomb_series.c | 11 ++- 3 files changed, 142 insertions(+), 8 deletions(-) diff --git a/acb_hypgeom/coulomb_jet.c b/acb_hypgeom/coulomb_jet.c index 07308147..85216c9e 100644 --- a/acb_hypgeom/coulomb_jet.c +++ b/acb_hypgeom/coulomb_jet.c @@ -11,6 +11,117 @@ #include "acb_hypgeom.h" +/* need special case for integer l and z = 0 since the recurrence relations break down */ +static void +_acb_hypgeom_coulomb_f_int_jet(acb_ptr F, const acb_t l, const acb_t eta, const acb_t z, slong len, slong prec) +{ + acb_poly_struct a[1]; + acb_poly_struct b[2]; + acb_poly_t zx, M, zxi; + acb_t t; + slong k; + int real; + + if (arf_cmp_si(arb_midref(acb_realref(l)), -1) < 0) + { + _acb_vec_indeterminate(F, len); + return; + } + + /* F = C * (z+x)^(l+1) e^(-+ i (z+x)) M(l + 1 -+ i eta, 2l+2, +- 2 i (z+x)) */ + + acb_poly_init(a); + acb_poly_init(b); + acb_poly_init(b + 1); + acb_poly_init(zx); + acb_poly_init(M); + acb_poly_init(zxi); + acb_init(t); + + acb_poly_set_coeff_acb(zx, 0, z); + acb_poly_set_coeff_si(zx, 1, 1); + + acb_div_onei(t, eta); + acb_add(t, t, l, prec); + acb_add_ui(t, t, 1, prec); + acb_poly_set_acb(a, t); + + acb_add_ui(t, l, 1, prec); + acb_mul_2exp_si(t, t, 1); + acb_poly_set_acb(b, t); + + acb_poly_one(b + 1); + + acb_onei(t); + acb_mul_2exp_si(t, t, 1); + acb_poly_scalar_mul(zxi, zx, t, prec); + + acb_hypgeom_pfq_series_direct(M, a, 1, b, 2, zxi, 1, -1, len, prec); + + acb_poly_scalar_mul_2exp_si(zxi, zxi, -1); + acb_poly_neg(zxi, zxi); + acb_poly_exp_series(zxi, zxi, len, prec); + acb_poly_mullow(M, M, zxi, len, prec); + + acb_add_ui(t, l, 1, prec); + acb_poly_pow_acb_series(zxi, zx, t, len, prec); + + acb_poly_mullow(M, M, zxi, len, prec); + + { + /* C = 2^l exp((-pi eta + lu + lv)/2) */ + acb_t C, lu, lv; + + acb_init(C); + acb_init(lu); + acb_init(lv); + + acb_add_ui(lu, l, 1, prec); + acb_mul_onei(t, eta); + acb_add(lu, lu, t, prec); + + acb_add_ui(lv, l, 1, prec); + acb_div_onei(t, eta); + acb_add(lv, lv, t, prec); + + acb_lgamma(lu, lu, prec); + acb_lgamma(lv, lv, prec); + + acb_const_pi(t, prec); + acb_add(C, lu, lv, prec); + acb_submul(C, t, eta, prec); + acb_mul_2exp_si(C, C, -1); + acb_exp(C, C, prec); + + acb_set_ui(t, 2); + acb_pow(t, t, l, prec); + acb_mul(C, C, t, prec); + + acb_poly_scalar_mul(M, M, C, prec); + + acb_clear(C); + acb_clear(lu); + acb_clear(lv); + } + + real = acb_is_real(z) && acb_is_real(eta); + + for (k = 0; k < len; k++) + { + acb_poly_get_coeff_acb(F + k, M, k); + if (real) + arb_zero(acb_imagref(F + k)); + } + + acb_poly_clear(a); + acb_poly_clear(b); + acb_poly_clear(b + 1); + acb_poly_clear(zx); + acb_poly_clear(M); + acb_poly_clear(zxi); + acb_clear(t); +} + 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) { @@ -25,6 +136,23 @@ _acb_hypgeom_coulomb_jet(acb_ptr F, acb_ptr G, acb_ptr Hpos, acb_ptr Hneg, const return; } + if (acb_contains_zero(z)) + { + if (F != NULL) + { + if (acb_is_int(l)) + _acb_hypgeom_coulomb_f_int_jet(F, l, eta, z, len, prec); + else + _acb_vec_indeterminate(F, len); + } + + if (G != NULL) _acb_vec_indeterminate(G, len); + if (Hpos != NULL) _acb_vec_indeterminate(Hpos, len); + if (Hneg != NULL) _acb_vec_indeterminate(Hneg, len); + return; + } + + acb_init(l1); acb_init(t); acb_init(R); diff --git a/acb_hypgeom/test/t-coulomb_series.c b/acb_hypgeom/test/t-coulomb_series.c index 35b1cae1..c3616189 100644 --- a/acb_hypgeom/test/t-coulomb_series.c +++ b/acb_hypgeom/test/t-coulomb_series.c @@ -24,7 +24,7 @@ int main() 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; + acb_t c, l, eta, z0; slong n1, n2, prec1, prec2; unsigned int mask; @@ -35,7 +35,7 @@ int main() 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); + acb_init(c); acb_init(l); acb_init(eta); acb_init(z0); prec1 = 2 + n_randint(state, 200); prec2 = 2 + n_randint(state, 200); @@ -72,7 +72,10 @@ int main() acb_poly_derivative(t, z, prec1); acb_poly_truncate(t, FLINT_MAX(n1 - 1, 0)); - if (!acb_poly_overlaps(w, t)) + /* hack: work around mullow(nan, 0) = 0 */ + acb_poly_get_coeff_acb(z0, z, 0); + + if (!acb_contains_zero(z0) && !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"); @@ -141,7 +144,7 @@ int main() 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); + acb_clear(c); acb_clear(l); acb_clear(eta); acb_clear(z0); } flint_randclear(state); diff --git a/arb_hypgeom/test/t-coulomb_series.c b/arb_hypgeom/test/t-coulomb_series.c index cea99085..cc3b9192 100644 --- a/arb_hypgeom/test/t-coulomb_series.c +++ b/arb_hypgeom/test/t-coulomb_series.c @@ -24,7 +24,7 @@ int main() for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++) { arb_poly_t F, G, F2, G2, z, w, t, u; - arb_t c, l, eta; + arb_t c, l, eta, z0; slong n1, n2, prec1, prec2; unsigned int mask; @@ -33,7 +33,7 @@ int main() arb_poly_init(z); arb_poly_init(w); arb_poly_init(t); arb_poly_init(u); - arb_init(c); arb_init(l); arb_init(eta); + arb_init(c); arb_init(l); arb_init(eta); arb_init(z0); prec1 = 2 + n_randint(state, 200); prec2 = 2 + n_randint(state, 200); @@ -61,7 +61,10 @@ int main() arb_poly_derivative(t, z, prec1); arb_poly_truncate(t, FLINT_MAX(n1 - 1, 0)); - if (!arb_poly_overlaps(w, t)) + /* hack: work around mullow(nan, 0) = 0 */ + arb_poly_get_coeff_arb(z0, z, 0); + + if (!arb_contains_zero(z0) && !arb_poly_overlaps(w, t)) { flint_printf("FAIL: wronskian, n1 = %wd\n\n", n1); flint_printf("l = "); arb_printd(l, 30); flint_printf("\n\n"); @@ -114,7 +117,7 @@ int main() arb_poly_clear(z); arb_poly_clear(w); arb_poly_clear(t); arb_poly_clear(u); - arb_clear(c); arb_clear(l); arb_clear(eta); + arb_clear(c); arb_clear(l); arb_clear(eta); arb_clear(z0); } flint_randclear(state);