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

add arb_hypgeom_u_integration; some cleanup

Browse source
This commit is contained in:
fredrik 2021-12-09 17:56:09 +01:00
parent 555bc0403e
commit f8f1c044a4
5 changed files with 608 additions and 46 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

1
arb_hypgeom.h
View file

@ -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);

81
arb_hypgeom/1f1_integration.c
View file

@ -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);
}

100
arb_hypgeom/test/t-u_integration.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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;
}

453
arb_hypgeom/u_integration.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}

19
doc/source/arb_hypgeom.rst
View file

@ -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
-------------------------------------------------------------------------------