From f8f1c044a44d9bbc1fb19938707d0330eb66e059 Mon Sep 17 00:00:00 2001 From: fredrik Date: Thu, 9 Dec 2021 17:56:09 +0100 Subject: [PATCH] add arb_hypgeom_u_integration; some cleanup --- arb_hypgeom.h | 1 + arb_hypgeom/1f1_integration.c | 81 +++--- arb_hypgeom/test/t-u_integration.c | 100 +++++++ arb_hypgeom/u_integration.c | 453 +++++++++++++++++++++++++++++ doc/source/arb_hypgeom.rst | 19 +- 5 files changed, 608 insertions(+), 46 deletions(-) create mode 100644 arb_hypgeom/test/t-u_integration.c create mode 100644 arb_hypgeom/u_integration.c diff --git a/arb_hypgeom.h b/arb_hypgeom.h index b1c33bb9..a385fb1e 100644 --- a/arb_hypgeom.h +++ b/arb_hypgeom.h @@ -74,6 +74,7 @@ void arb_hypgeom_u(arb_t res, const arb_t a, const arb_t b, const arb_t z, slong void arb_hypgeom_2f1(arb_t res, const arb_t a, const arb_t b, const arb_t c, const arb_t z, int regularized, slong prec); void arb_hypgeom_1f1_integration(arb_t res, const arb_t a, const arb_t b, const arb_t z, int regularized, slong prec); +void arb_hypgeom_u_integration(arb_t res, const arb_t a, const arb_t b, const arb_t z, slong prec); void arb_hypgeom_erf(arb_t res, const arb_t z, slong prec); void _arb_hypgeom_erf_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec); diff --git a/arb_hypgeom/1f1_integration.c b/arb_hypgeom/1f1_integration.c index 50a9ace9..3ea39277 100644 --- a/arb_hypgeom/1f1_integration.c +++ b/arb_hypgeom/1f1_integration.c @@ -38,13 +38,13 @@ polynomial) but for simplicity we just do interval arithmetic. /* z*u + 0.5 [(a-1) log(u^2+v^2) + (b-a-1) log((u-1)^2+v^2)] */ static di_t -di_integrand_edge(di_t u, di_t v, di_t a, di_t b, di_t z) +di_integrand_edge(di_t u, di_t v, di_t a1, di_t ba1, di_t z) { di_t X, Y, Z; X = di_fast_mul(z, u); - Y = di_fast_mul(di_fast_sub_d(a, 1.0), di_fast_log_nonnegative(di_fast_add(di_fast_sqr(u), di_fast_sqr(v)))); - Z = di_fast_mul(di_fast_sub_d(di_fast_sub(b, a), 1.0), di_fast_log_nonnegative(di_fast_add(di_fast_sqr(di_fast_sub_d(u, 1.0)), di_fast_sqr(v)))); + Y = di_fast_mul(a1, di_fast_log_nonnegative(di_fast_add(di_fast_sqr(u), di_fast_sqr(v)))); + Z = di_fast_mul(ba1, di_fast_log_nonnegative(di_fast_add(di_fast_sqr(di_fast_sub_d(u, 1.0)), di_fast_sqr(v)))); return di_fast_add(X, di_fast_mul_d(di_fast_add(Y, Z), 0.5)); } @@ -54,12 +54,12 @@ which == 0 - d/du g(u,v) = z + u*(a-1)/(u^2+v^2) + (u-1)*(b-a-1)/(v^2+(1-u)^2) which == 1 - d/dv g(u,v) = v*(a-1)/(u^2+v^2) + v*(b-a-1)/(v^2+(1-u)^2) */ static di_t -di_integrand_edge_diff(di_t u, di_t v, di_t a, di_t b, di_t z, int which) +di_integrand_edge_diff(di_t u, di_t v, di_t a1, di_t ba1, di_t z, int which) { di_t Y, Z; - Y = di_fast_div(di_fast_sub_d(a, 1.0), di_fast_add(di_fast_sqr(u), di_fast_sqr(v))); - Z = di_fast_div(di_fast_sub_d(di_fast_sub(b, a), 1.0), di_fast_add(di_fast_sqr(di_fast_sub_d(u, 1.0)), di_fast_sqr(v))); + Y = di_fast_div(a1, di_fast_add(di_fast_sqr(u), di_fast_sqr(v))); + Z = di_fast_div(ba1, di_fast_add(di_fast_sqr(di_fast_sub_d(u, 1.0)), di_fast_sqr(v))); if (which == 0) return di_fast_add(z, di_fast_add(di_fast_mul(u, Y), di_fast_mul(di_fast_sub_d(u, 1.0), Z))); @@ -81,10 +81,10 @@ static di_t di_subinterval(di_t x, slong i, slong N) } static void -integrand_wide_bound5(acb_t res, const acb_t t, const arb_t a, const arb_t b, const arb_t z, slong prec) +integrand_wide_bound5(acb_t res, const acb_t t, const arb_t a1, const arb_t ba1, const arb_t z, slong prec) { slong i, N; - di_t du, dv, da, db, dz, dg, dgprime; + di_t du, dv, da1, dba1, dz, dg, dgprime; double radius, bound; double start, end; int which; @@ -93,8 +93,8 @@ integrand_wide_bound5(acb_t res, const acb_t t, const arb_t a, const arb_t b, co N = 8; bound = -D_INF; - da = arb_get_di(a); - db = arb_get_di(b); + da1 = arb_get_di(a1); + dba1 = arb_get_di(ba1); dz = arb_get_di(z); /* left edge: left(u) + [0, right(v)] */ @@ -118,11 +118,11 @@ integrand_wide_bound5(acb_t res, const acb_t t, const arb_t a, const arb_t b, co /* g(u,mid(v)) + g'(u,v) * [0, radius] */ #if 1 - dg = di_integrand_edge(du, di_fast_mid(dv), da, db, dz); - dgprime = di_integrand_edge_diff(du, dv, da, db, dz, 1); + dg = di_integrand_edge(du, di_fast_mid(dv), da1, dba1, dz); + dgprime = di_integrand_edge_diff(du, dv, da1, dba1, dz, 1); dg = di_fast_add(dg, di_fast_mul(dgprime, di_interval(0.0, radius))); #else - dg = di_integrand_edge(du, dv, da, db, dz); + dg = di_integrand_edge(du, dv, da1, dba1, dz); #endif bound = FLINT_MAX(bound, dg.b); @@ -144,11 +144,11 @@ integrand_wide_bound5(acb_t res, const acb_t t, const arb_t a, const arb_t b, co /* g(mid(u),v) + g'(u,v) * [0, radius] */ #if 1 - dg = di_integrand_edge(di_fast_mid(du), dv, da, db, dz); - dgprime = di_integrand_edge_diff(du, dv, da, db, dz, 0); + dg = di_integrand_edge(di_fast_mid(du), dv, da1, dba1, dz); + dgprime = di_integrand_edge_diff(du, dv, da1, dba1, dz, 0); dg = di_fast_add(dg, di_fast_mul(dgprime, di_interval(0.0, radius))); #else - dg = di_integrand_edge(du, dv, da, db, dz); + dg = di_integrand_edge(du, dv, da1, dba1, dz); #endif bound = FLINT_MAX(bound, dg.b); @@ -165,7 +165,8 @@ integrand_wide_bound5(acb_t res, const acb_t t, const arb_t a, const arb_t b, co arb_clear(abound); } -void +/* todo: fix acb_pow(_arb) */ +static void acb_my_pow_arb(acb_t res, const acb_t a, const arb_t b, slong prec) { if (acb_contains_zero(a) && arb_is_positive(b)) @@ -202,18 +203,16 @@ acb_my_pow_arb(acb_t res, const acb_t a, const arb_t b, slong prec) static int integrand(acb_ptr out, const acb_t t, void * param, slong order, slong prec) { - arb_srcptr a, b, z; + arb_srcptr a1, ba1, z; acb_t s, u, v; - arb_t e; - a = ((arb_srcptr) param) + 0; - b = ((arb_srcptr) param) + 1; + a1 = ((arb_srcptr) param) + 0; + ba1 = ((arb_srcptr) param) + 1; z = ((arb_srcptr) param) + 2; acb_init(s); acb_init(u); acb_init(v); - arb_init(e); acb_sub_ui(v, t, 1, prec); acb_neg(v, v); @@ -223,7 +222,7 @@ integrand(acb_ptr out, const acb_t t, void * param, slong order, slong prec) if (!arb_is_positive(acb_realref(t)) || !arb_is_positive(acb_realref(v))) acb_indeterminate(out); else - integrand_wide_bound5(out, t, a, b, z, prec); + integrand_wide_bound5(out, t, a1, ba1, z, prec); } else { @@ -234,13 +233,10 @@ integrand(acb_ptr out, const acb_t t, void * param, slong order, slong prec) acb_exp(s, s, prec); /* t^(a-1) */ - arb_sub_ui(e, a, 1, prec); - acb_my_pow_arb(u, t, e, prec); + acb_my_pow_arb(u, t, a1, prec); /* (1-t)^(b-a-1) */ - arb_sub(e, b, a, prec); - arb_sub_ui(e, e, 1, prec); - acb_my_pow_arb(v, v, e, prec); + acb_my_pow_arb(v, v, ba1, prec); acb_mul(out, s, u, prec); acb_mul(out, out, v, prec); @@ -250,15 +246,12 @@ integrand(acb_ptr out, const acb_t t, void * param, slong order, slong prec) acb_mul_arb(s, t, z, prec); /* t^(a-1) */ - arb_sub_ui(e, a, 1, prec); acb_log(u, t, prec); - acb_mul_arb(u, u, e, prec); + acb_mul_arb(u, u, a1, prec); /* (1-t)^(b-a-1) */ - arb_sub(e, b, a, prec); - arb_sub_ui(e, e, 1, prec); acb_log(v, v, prec); - acb_mul_arb(v, v, e, prec); + acb_mul_arb(v, v, ba1, prec); acb_add(out, s, u, prec); acb_add(out, out, v, prec); @@ -269,7 +262,6 @@ integrand(acb_ptr out, const acb_t t, void * param, slong order, slong prec) acb_clear(s); acb_clear(u); acb_clear(v); - arb_clear(e); return 0; } @@ -332,19 +324,22 @@ arb_hypgeom_1f1_integration(arb_t res, const arb_t a, const arb_t b, const arb_t { acb_calc_integrate_opt_t opt; arb_struct param[3]; - arb_t t; + arb_t t, a1, ba1; acb_t zero, one, I; mag_t abs_tol; int ok; arb_init(t); + arb_init(a1); + arb_init(ba1); + + arb_sub_ui(a1, a, 1, prec); + arb_sub(ba1, b, a, prec); + arb_sub_ui(ba1, ba1, 1, prec); ok = arb_is_finite(z); - arb_sub_ui(t, a, 1, prec); - ok = ok && arb_is_nonnegative(t); - arb_sub(t, b, a, prec); - arb_sub_ui(t, t, 1, prec); - ok = ok && arb_is_nonnegative(t); + ok = ok && arb_is_nonnegative(a1); + ok = ok && arb_is_nonnegative(ba1); if (!ok) { @@ -357,8 +352,8 @@ arb_hypgeom_1f1_integration(arb_t res, const arb_t a, const arb_t b, const arb_t acb_init(one); acb_init(I); - param[0] = *a; - param[1] = *b; + param[0] = *a1; + param[1] = *ba1; param[2] = *z; acb_calc_integrate_opt_init(opt); @@ -390,4 +385,6 @@ arb_hypgeom_1f1_integration(arb_t res, const arb_t a, const arb_t b, const arb_t } arb_clear(t); + arb_clear(a1); + arb_clear(ba1); } diff --git a/arb_hypgeom/test/t-u_integration.c b/arb_hypgeom/test/t-u_integration.c new file mode 100644 index 00000000..888e6c38 --- /dev/null +++ b/arb_hypgeom/test/t-u_integration.c @@ -0,0 +1,100 @@ +/* + Copyright (C) 2021 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 . +*/ + +#include "arb_hypgeom.h" + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("u_integration...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 100 * arb_test_multiplier(); iter++) + { + arb_t a, b, z, r1, r2; + slong prec1, prec2; + + prec1 = 2 + n_randint(state, 80); + prec2 = 2 + n_randint(state, 150); + + arb_init(a); + arb_init(b); + arb_init(z); + arb_init(r1); + arb_init(r2); + + arb_randtest_precise(a, state, 1 + n_randint(state, 200), 1 + n_randint(state, 8)); + arb_randtest_precise(b, state, 1 + n_randint(state, 200), 1 + n_randint(state, 8)); + arb_randtest_precise(z, state, 1 + n_randint(state, 200), 1 + n_randint(state, 8)); + arb_randtest(r1, state, 1 + n_randint(state, 200), 1 + n_randint(state, 5)); + arb_randtest(r2, state, 1 + n_randint(state, 200), 1 + n_randint(state, 5)); + + arb_add_ui(a, a, n_randint(state, 100), prec1 + 100); + arb_add_ui(b, b, n_randint(state, 200), prec1 + 100); + arb_add_ui(z, z, n_randint(state, 100), prec1 + 100); + + arb_hypgeom_u_integration(r1, a, b, z, prec1); + + if (arb_is_finite(r1)) + { + if (n_randint(state, 2)) + arb_hypgeom_u(r2, a, b, z, prec2); + else + arb_hypgeom_u_integration(r2, a, b, z, prec2); + + if (!arb_overlaps(r1, r2)) + { + flint_printf("FAIL: overlap\n\n"); + flint_printf("a = "); arb_printd(a, 30); flint_printf("\n\n"); + flint_printf("b = "); arb_printd(b, 30); flint_printf("\n\n"); + flint_printf("z = "); arb_printd(z, 30); flint_printf("\n\n"); + flint_printf("r1 = "); arb_printd(r1, 30); flint_printf("\n\n"); + flint_printf("r2 = "); arb_printd(r2, 30); flint_printf("\n\n"); + flint_abort(); + } + } + + if (iter == 0) + { + prec1 = 333; + + arb_set_str(a, "1.2e7", prec1); + arb_set_str(b, "1.3e8", prec1); + arb_set_str(z, "1.4e8", prec1); + + arb_hypgeom_u_integration(r1, a, b, z, prec1); + arb_set_str(r2, "7.423720668435352497136323553980928479792327057967344627374385949694801879462740896706198120e-90723524 +/- 6.47e-90723615", prec1); + + if (!arb_overlaps(r1, r2) || arb_rel_accuracy_bits(r1) < arb_rel_accuracy_bits(r2) - 10) + { + flint_printf("FAIL: overlap (1)\n\n"); + flint_printf("r1 = "); arb_printd(r1, 100); flint_printf("\n\n"); + flint_printf("r2 = "); arb_printd(r2, 100); flint_printf("\n\n"); + flint_abort(); + } + } + + arb_clear(a); + arb_clear(b); + arb_clear(z); + arb_clear(r1); + arb_clear(r2); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/arb_hypgeom/u_integration.c b/arb_hypgeom/u_integration.c new file mode 100644 index 00000000..cfaeffa5 --- /dev/null +++ b/arb_hypgeom/u_integration.c @@ -0,0 +1,453 @@ +/* + Copyright (C) 2021 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 . +*/ + +#include "arb_hypgeom.h" +#include "acb_calc.h" +#include "double_interval.h" + +/* + +Integrand (see comments for 1f1_integration): + + exp(f(t)) where f(z) = -z*t + a1*log(t) + ba1*log(1+t) + + g(u,v) = -z*u + 0.5*[a1*log(u^2+v^2) + ba1*log((1+u)^2+v^2)] + + d/du g(u,v) = -z + u*a1/(u^2+v^2) + (1+u)*ba1/(v^2+(1+u)^2) + d/dv g(u,v) = v*a1/(u^2+v^2) + v*ba1/(v^2+(1+u)^2) +*/ + +static di_t +di_integrand_edge(di_t u, di_t v, di_t a1, di_t ba1, di_t z) +{ + di_t X, Y, Z; + + X = di_neg(di_fast_mul(z, u)); + Y = di_fast_mul(a1, di_fast_log_nonnegative(di_fast_add(di_fast_sqr(u), di_fast_sqr(v)))); + Z = di_fast_mul(ba1, di_fast_log_nonnegative(di_fast_add(di_fast_sqr(di_fast_add_d(u, 1.0)), di_fast_sqr(v)))); + + return di_fast_add(X, di_fast_mul_d(di_fast_add(Y, Z), 0.5)); +} + +/* +which == 0 - d/du g(u,v) = -z + u*(a-1)/(u^2+v^2) + (u+1)*(b-a-1)/(v^2+(1+u)^2) +which == 1 - d/dv g(u,v) = v*(a-1)/(u^2+v^2) + v*(b-a-1)/(v^2+(1+u)^2) +*/ +static di_t +di_integrand_edge_diff(di_t u, di_t v, di_t a1, di_t ba1, di_t z, int which) +{ + di_t Y, Z; + + Y = di_fast_div(a1, di_fast_add(di_fast_sqr(u), di_fast_sqr(v))); + Z = di_fast_div(ba1, di_fast_add(di_fast_sqr(di_fast_add_d(u, 1.0)), di_fast_sqr(v))); + + if (which == 0) + return di_fast_add(di_neg(z), di_fast_add(di_fast_mul(u, Y), di_fast_mul(di_fast_add_d(u, 1.0), Z))); + else + return di_fast_mul(v, di_fast_add(Y, Z)); +} + +static di_t di_subinterval(di_t x, slong i, slong N) +{ + di_t res; + double step; + + step = (x.b - x.a) / N; + + res.a = x.a + step * i; + res.b = (i == N - 1) ? x.b : x.a + step * (i + 1); + + return res; +} + +static void +integrand_wide_bound5(acb_t res, const acb_t t, const arb_t a1, const arb_t ba1, const arb_t z, slong prec) +{ + slong i, N; + di_t du, dv, da1, dba1, dz, dg, dgprime; + double radius, bound; + double start, end; + int which; + arb_t abound; + + N = 8; + bound = -D_INF; + + da1 = arb_get_di(a1); + dba1 = arb_get_di(ba1); + dz = arb_get_di(z); + + /* left edge: left(u) + [0, right(v)] */ + /* right edge: right(u) + [0, right(v)] */ + for (which = 0; which < 2; which++) + { + du = arb_get_di(acb_realref(t)); + if (which == 0) + du.b = du.a; + else + du.a = du.b; + + dv = arb_get_di(acb_imagref(t)); + start = 0.0; + end = dv.b; + + for (i = 0; i < N; i++) + { + dv = di_subinterval(di_interval(start, end), i, N); + radius = di_fast_ubound_radius(dv); + + /* g(u,mid(v)) + g'(u,v) * [0, radius] */ +#if 1 + dg = di_integrand_edge(du, di_fast_mid(dv), da1, dba1, dz); + dgprime = di_integrand_edge_diff(du, dv, da1, dba1, dz, 1); + dg = di_fast_add(dg, di_fast_mul(dgprime, di_interval(0.0, radius))); +#else + dg = di_integrand_edge(du, dv, da1, dba1, dz); +#endif + + bound = FLINT_MAX(bound, dg.b); + } + } + + du = arb_get_di(acb_realref(t)); + start = du.a; + end = du.b; + + dv = arb_get_di(acb_imagref(t)); + dv.a = dv.b; + + /* top edge: [left(u), right(u)] + right(v) */ + for (i = 0; i < N; i++) + { + du = di_subinterval(di_interval(start, end), i, N); + radius = di_fast_ubound_radius(du); + + /* g(mid(u),v) + g'(u,v) * [0, radius] */ +#if 1 + dg = di_integrand_edge(di_fast_mid(du), dv, da1, dba1, dz); + dgprime = di_integrand_edge_diff(du, dv, da1, dba1, dz, 0); + dg = di_fast_add(dg, di_fast_mul(dgprime, di_interval(0.0, radius))); +#else + dg = di_integrand_edge(du, dv, da1, dba1, dz); +#endif + + bound = FLINT_MAX(bound, dg.b); + } + + arb_init(abound); + arb_set_d(abound, bound); + arb_exp(abound, abound, prec); + + acb_zero(res); + arb_add_error(acb_realref(res), abound); + arb_add_error(acb_imagref(res), abound); + + arb_clear(abound); +} + +/* todo: fix acb_pow(_arb) */ +static void +acb_my_pow_arb(acb_t res, const acb_t a, const arb_t b, slong prec) +{ + if (acb_contains_zero(a) && arb_is_positive(b)) + { + /* |a^b| <= |a|^b */ + arb_t t, u; + + arb_init(t); + arb_init(u); + + acb_abs(t, a, prec); + arb_get_abs_ubound_arf(arb_midref(t), t, prec); + mag_zero(arb_radref(t)); + + if (arf_cmpabs_2exp_si(arb_midref(t), 0) < 0) + arb_get_abs_lbound_arf(arb_midref(u), b, prec); + else + arb_get_abs_ubound_arf(arb_midref(u), b, prec); + + arb_pow(t, t, u, prec); + + acb_zero(res); + acb_add_error_arb(res, t); + + arb_clear(t); + arb_clear(u); + } + else + { + acb_pow_arb(res, a, b, prec); + } +} + +static int +integrand(acb_ptr out, const acb_t t, void * param, slong order, slong prec) +{ + arb_srcptr a1, ba1, z; + acb_t s, u, v; + + a1 = ((arb_srcptr) param) + 0; + ba1 = ((arb_srcptr) param) + 1; + z = ((arb_srcptr) param) + 2; + + acb_init(s); + acb_init(u); + acb_init(v); + + acb_add_ui(v, t, 1, prec); + + if (order == 1) + { + if (!arb_is_positive(acb_realref(t))) + acb_indeterminate(out); + else + integrand_wide_bound5(out, t, a1, ba1, z, prec); + } + else + { + if (acb_contains_zero(t) || acb_contains_zero(v)) + { + /* exp(-z t) */ + acb_mul_arb(s, t, z, prec); + acb_neg(s, s); + acb_exp(s, s, prec); + + /* t^(a-1) */ + acb_my_pow_arb(u, t, a1, prec); + /* (1+t)^(b-a-1) */ + acb_pow_arb(v, v, ba1, prec); + + acb_mul(out, s, u, prec); + acb_mul(out, out, v, prec); + } + else + { + acb_mul_arb(s, t, z, prec); + acb_neg(s, s); + + /* t^(a-1) */ + acb_log(u, t, prec); + acb_mul_arb(u, u, a1, prec); + + /* (1+t)^(b-a-1) */ + acb_log(v, v, prec); + acb_mul_arb(v, v, ba1, prec); + + acb_add(out, s, u, prec); + acb_add(out, out, v, prec); + acb_exp(out, out, prec); + } + } + + acb_clear(s); + acb_clear(u); + acb_clear(v); + + return 0; +} + +/* estimate integral by magnitude at peak */ +static void +estimate_magnitude(mag_t res, const arb_t ra, const arb_t rb, const arb_t rz) +{ + double a, b, z, t1, t2, u, m; + fmpz_t e; + + a = arf_get_d(arb_midref(ra), ARF_RND_NEAR); + b = arf_get_d(arb_midref(rb), ARF_RND_NEAR); + z = arf_get_d(arb_midref(rz), ARF_RND_NEAR); + + u = 4 - 4*b + b*b + 4*a*z - 2*b*z + z*z; + + if (u >= 0.0) + { + t1 = (-2 + b - z + sqrt(u)) / (2 * z); + t2 = (-2 + b - z - sqrt(u)) / (2 * z); + } + else + { + t1 = 1e-8; + t2 = 1e-8; + } + + m = -1e300; + + if (t1 > 0.0) + { + t1 = -z * t1 + (a - 1) * log(t1) + (b - a - 1) * log(1 + t1); + m = FLINT_MAX(m, t1); + } + + if (t2 > 0.0) + { + t2 = -z * t2 + (a - 1) * log(t2) + (b - a - 1) * log(1 + t2); + m = FLINT_MAX(m, t2); + } + + m /= log(2); + + if (fabs(m) < 1e300) + { + fmpz_init(e); + fmpz_set_d(e, m); + mag_set_d_2exp_fmpz(res, 1.0, e); + fmpz_clear(e); + } + else + { + mag_zero(res); + } +} + +static void +bound_tail(mag_t bound, const arb_t a1, const arb_t ba1, const arb_t z, const arb_t N, slong prec) +{ + arb_t s, u, v, C; + + arb_init(s); + arb_init(u); + arb_init(v); + arb_init(C); + + /* + Assume N >= 1 and t >= 0. + + -z*(N+t) + (a-1)*log(N+t) + (b-a-1)*log(1+N+t) + <= [-z*N + (a-1)*log(N) + (b-a-1)*log(1+N)] + [-z*t + [max(0, a-1) + max(0, b-a-1)]*log(1+t/N)] + <= [-z*N + (a-1)*log(N) + (b-a-1)*log(1+N)] + [-z*t + [max(0, a-1) + max(0, b-a-1)]*(t/N)] + + Let C = max(0, a-1) + max(0, b-a-1). Then the remainder integral + is bounded by integrand(N) * N / (N*z - C), assuming that N*z > C. + */ + + arb_max(u, u, a1, prec); + arb_max(v, v, ba1, prec); + arb_add(C, u, v, prec); + /* s = N*z - C */ + arb_mul(s, N, z, prec); + arb_sub(s, s, C, prec); + + if (arb_is_positive(s)) + { + arb_div(C, N, s, prec); + + /* exp(-z*N) */ + arb_mul(s, N, z, prec); + arb_neg(s, s); + + /* N^(a-1) */ + arb_log(u, N, prec); + arb_mul(u, u, a1, prec); + + /* (1+N)^(b-a-1) */ + arb_add_ui(v, N, 1, prec); + arb_log(v, v, prec); + arb_mul(v, v, ba1, prec); + + arb_add(s, s, u, prec); + arb_add(s, s, v, prec); + arb_exp(s, s, prec); + + arb_mul(s, s, C, prec); + + arb_get_mag(bound, s); + } + else + { + mag_inf(bound); + } + + arb_clear(s); + arb_clear(u); + arb_clear(v); + arb_clear(C); +} + +void +arb_hypgeom_u_integration(arb_t res, const arb_t a, const arb_t b, const arb_t z, slong prec) +{ + acb_calc_integrate_opt_t opt; + arb_struct param[3]; + arb_t t, a1, ba1; + acb_t zero, N, I; + mag_t abs_tol, tail_bound; + slong i; + fmpz_t n; + int ok; + + arb_init(t); + arb_init(a1); + arb_init(ba1); + + arb_sub_ui(a1, a, 1, prec); + arb_sub(ba1, b, a, prec); + arb_sub_ui(ba1, ba1, 1, prec); + + ok = arb_is_finite(z) && arb_is_positive(z); + ok = ok && arb_is_nonnegative(a1); + ok = ok && arb_is_finite(b); + + if (!ok) + { + arb_indeterminate(res); + } + else + { + mag_init(abs_tol); + mag_init(tail_bound); + acb_init(zero); + acb_init(zero); + acb_init(N); + acb_init(I); + fmpz_init(n); + + param[0] = *a1; + param[1] = *ba1; + param[2] = *z; + + acb_calc_integrate_opt_init(opt); + /* opt->verbose = 2; */ + /* opt->eval_limit = WORD_MAX; */ + + estimate_magnitude(abs_tol, a, b, z); + mag_mul_2exp_si(abs_tol, abs_tol, -prec); + + for (i = 1; i < FLINT_BITS; i++) + { + fmpz_one(n); + fmpz_mul_2exp(n, n, i); + acb_one(N); + arb_mul_2exp_fmpz(acb_realref(N), acb_realref(N), n); + bound_tail(tail_bound, a1, ba1, z, acb_realref(N), 64); + if (mag_cmp(tail_bound, abs_tol) < 0) + break; + } + + acb_calc_integrate(I, integrand, param, zero, N, prec, abs_tol, opt, prec); + arb_add_error_mag(acb_realref(I), tail_bound); + + arb_rgamma(t, a, prec); + arb_mul(acb_realref(I), acb_realref(I), t, prec); + + arb_set(res, acb_realref(I)); + + mag_clear(abs_tol); + mag_clear(tail_bound); + acb_clear(zero); + acb_clear(N); + acb_clear(I); + fmpz_clear(n); + } + + arb_clear(t); + arb_clear(a1); + arb_clear(ba1); +} diff --git a/doc/source/arb_hypgeom.rst b/doc/source/arb_hypgeom.rst index b26b3d29..9de46904 100644 --- a/doc/source/arb_hypgeom.rst +++ b/doc/source/arb_hypgeom.rst @@ -157,14 +157,13 @@ Confluent hypergeometric functions .. function:: void arb_hypgeom_1f1_integration(arb_t res, const arb_t a, const arb_t b, const arb_t z, int regularized, slong prec) - Computes the confluent hypergeometric function using numerical integration. - The formula + Computes the confluent hypergeometric function using numerical integration + of the representation .. math :: - {}_1F_1(a,b,z) = \frac{\Gamma(b)}{\Gamma(a) \Gamma(b-a)} \int_0^1 e^{zt} t^{a-1} (1-t)^{b-a-1} dt + {}_1F_1(a,b,z) = \frac{\Gamma(b)}{\Gamma(a) \Gamma(b-a)} \int_0^1 e^{zt} t^{a-1} (1-t)^{b-a-1} dt. - is used. This algorithm can be useful if the parameters are large. This will currently only return a finite enclosure if `a \ge 1` and `b - a \ge 1`. @@ -172,6 +171,18 @@ Confluent hypergeometric functions Computes the confluent hypergeometric function `U(a,b,z)`. +.. function:: void arb_hypgeom_u_integration(arb_t res, const arb_t a, const arb_t b, const arb_t z, int regularized, slong prec) + + Computes the confluent hypergeometric function `U(a,b,z)` using numerical integration + of the representation + + .. math :: + + U(a,b,z) = \frac{1}{\Gamma(a)} \int_0^{\infty} e^{-zt} t^{a-1} (1+t)^{b-a-1} dt. + + This algorithm can be useful if the parameters are large. This will currently + only return a finite enclosure if `a \ge 1` and `z > 0`. + Gauss hypergeometric function -------------------------------------------------------------------------------