mirror of
https://github.com/vale981/arb
synced 2025-03-04 17:01:40 -05:00
more numerical integration implementations
This commit is contained in:
fredrik
parent
75df786325
commit
e79a80340f
15 changed files with 679 additions and 34 deletions
|
|
@ -1,5 +1,5 @@
|
|||
/*
|
||||
Copyright (C) 2014-2015 Fredrik Johansson
|
||||
Copyright (C) 2014-2015, 2021 Fredrik Johansson
|
||||
|
||||
This file is part of Arb.
|
||||
|
||||
|
|
@ -9,6 +9,7 @@
|
|||
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "arb_hypgeom.h"
|
||||
#include "acb_hypgeom.h"
|
||||
|
||||
void
|
||||
|
|
@ -200,7 +201,7 @@ acb_hypgeom_bessel_i_0f1(acb_t res, const acb_t nu, const acb_t z, int scaled, s
|
|||
}
|
||||
|
||||
void
|
||||
acb_hypgeom_bessel_i(acb_t res, const acb_t nu, const acb_t z, slong prec)
|
||||
acb_hypgeom_bessel_i_nointegration(acb_t res, const acb_t nu, const acb_t z, int scaled, slong prec)
|
||||
{
|
||||
mag_t zmag;
|
||||
|
||||
|
|
@ -209,27 +210,65 @@ acb_hypgeom_bessel_i(acb_t res, const acb_t nu, const acb_t z, slong prec)
|
|||
|
||||
if (mag_cmp_2exp_si(zmag, 4) < 0 ||
|
||||
(mag_cmp_2exp_si(zmag, 64) < 0 && 2 * mag_get_d(zmag) < prec))
|
||||
acb_hypgeom_bessel_i_0f1(res, nu, z, 0, prec);
|
||||
acb_hypgeom_bessel_i_0f1(res, nu, z, scaled, prec);
|
||||
else
|
||||
acb_hypgeom_bessel_i_asymp(res, nu, z, 0, prec);
|
||||
acb_hypgeom_bessel_i_asymp(res, nu, z, scaled, prec);
|
||||
|
||||
mag_clear(zmag);
|
||||
}
|
||||
|
||||
void
|
||||
acb_hypgeom_bessel_i_scaled(acb_t res, const acb_t nu, const acb_t z, slong prec)
|
||||
_acb_hypgeom_bessel_i(acb_t res, const acb_t nu, const acb_t z, int scaled, slong prec)
|
||||
{
|
||||
mag_t zmag;
|
||||
acb_t res2;
|
||||
slong acc, max, t;
|
||||
|
||||
mag_init(zmag);
|
||||
acb_get_mag(zmag, z);
|
||||
acb_init(res2);
|
||||
|
||||
if (mag_cmp_2exp_si(zmag, 4) < 0 ||
|
||||
(mag_cmp_2exp_si(zmag, 64) < 0 && 2 * mag_get_d(zmag) < prec))
|
||||
acb_hypgeom_bessel_i_0f1(res, nu, z, 1, prec);
|
||||
else
|
||||
acb_hypgeom_bessel_i_asymp(res, nu, z, 1, prec);
|
||||
acb_hypgeom_bessel_i_nointegration(res2, nu, z, scaled, prec);
|
||||
|
||||
mag_clear(zmag);
|
||||
acc = acb_rel_accuracy_bits(res2);
|
||||
|
||||
if (acc < 0.5 * prec)
|
||||
{
|
||||
max = prec;
|
||||
t = acb_rel_accuracy_bits(z);
|
||||
max = FLINT_MIN(max, t);
|
||||
t = acb_rel_accuracy_bits(nu);
|
||||
max = FLINT_MIN(max, t);
|
||||
|
||||
if (max > 2 && acc < 0.5 * max)
|
||||
{
|
||||
if (acb_is_real(nu) && acb_is_real(z) && arf_cmp_2exp_si(arb_midref(acb_realref(nu)), -1) > 0 &&
|
||||
arf_cmpabs_2exp_si(arb_midref(acb_realref(nu)), 60) < 0 &&
|
||||
arf_cmpabs_2exp_si(arb_midref(acb_realref(z)), 60) < 0)
|
||||
{
|
||||
|
||||
arb_hypgeom_bessel_i_integration(acb_realref(res),
|
||||
acb_realref(nu), acb_realref(z), scaled, prec);
|
||||
arb_zero(acb_imagref(res));
|
||||
|
||||
if (acb_rel_accuracy_bits(res) > acb_rel_accuracy_bits(res2) ||
|
||||
(acb_is_finite(res) && !acb_is_finite(res2)))
|
||||
{
|
||||
acb_swap(res, res2);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
acb_swap(res, res2);
|
||||
acb_clear(res2);
|
||||
}
|
||||
|
||||
void
|
||||
acb_hypgeom_bessel_i(acb_t res, const acb_t nu, const acb_t z, slong prec)
|
||||
{
|
||||
_acb_hypgeom_bessel_i(res, nu, z, 0, prec);
|
||||
}
|
||||
|
||||
void
|
||||
acb_hypgeom_bessel_i_scaled(acb_t res, const acb_t nu, const acb_t z, slong prec)
|
||||
{
|
||||
_acb_hypgeom_bessel_i(res, nu, z, 1, prec);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -9,6 +9,7 @@
|
|||
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "arb_hypgeom.h"
|
||||
#include "acb_hypgeom.h"
|
||||
|
||||
void
|
||||
|
|
@ -206,7 +207,7 @@ acb_hypgeom_bessel_k_0f1(acb_t res, const acb_t nu, const acb_t z, int scaled, s
|
|||
}
|
||||
|
||||
void
|
||||
acb_hypgeom_bessel_k(acb_t res, const acb_t nu, const acb_t z, slong prec)
|
||||
acb_hypgeom_bessel_k_nointegration(acb_t res, const acb_t nu, const acb_t z, int scaled, slong prec)
|
||||
{
|
||||
mag_t zmag;
|
||||
|
||||
|
|
@ -215,27 +216,65 @@ acb_hypgeom_bessel_k(acb_t res, const acb_t nu, const acb_t z, slong prec)
|
|||
|
||||
if (mag_cmp_2exp_si(zmag, 4) < 0 ||
|
||||
(mag_cmp_2exp_si(zmag, 64) < 0 && 2 * mag_get_d(zmag) < prec))
|
||||
acb_hypgeom_bessel_k_0f1(res, nu, z, 0, prec);
|
||||
acb_hypgeom_bessel_k_0f1(res, nu, z, scaled, prec);
|
||||
else
|
||||
acb_hypgeom_bessel_k_asymp(res, nu, z, 0, prec);
|
||||
acb_hypgeom_bessel_k_asymp(res, nu, z, scaled, prec);
|
||||
|
||||
mag_clear(zmag);
|
||||
}
|
||||
|
||||
void
|
||||
acb_hypgeom_bessel_k_scaled(acb_t res, const acb_t nu, const acb_t z, slong prec)
|
||||
_acb_hypgeom_bessel_k(acb_t res, const acb_t nu, const acb_t z, int scaled, slong prec)
|
||||
{
|
||||
mag_t zmag;
|
||||
acb_t res2;
|
||||
slong acc, max, t;
|
||||
|
||||
mag_init(zmag);
|
||||
acb_get_mag(zmag, z);
|
||||
acb_init(res2);
|
||||
|
||||
if (mag_cmp_2exp_si(zmag, 4) < 0 ||
|
||||
(mag_cmp_2exp_si(zmag, 64) < 0 && 2 * mag_get_d(zmag) < prec))
|
||||
acb_hypgeom_bessel_k_0f1(res, nu, z, 1, prec);
|
||||
else
|
||||
acb_hypgeom_bessel_k_asymp(res, nu, z, 1, prec);
|
||||
acb_hypgeom_bessel_k_nointegration(res2, nu, z, scaled, prec);
|
||||
|
||||
mag_clear(zmag);
|
||||
acc = acb_rel_accuracy_bits(res2);
|
||||
|
||||
if (acc < 0.5 * prec)
|
||||
{
|
||||
max = prec;
|
||||
t = acb_rel_accuracy_bits(z);
|
||||
max = FLINT_MIN(max, t);
|
||||
t = acb_rel_accuracy_bits(nu);
|
||||
max = FLINT_MIN(max, t);
|
||||
|
||||
if (max > 2 && acc < 0.5 * max)
|
||||
{
|
||||
if (acb_is_real(nu) && acb_is_real(z) && arf_cmp_d(arb_midref(acb_realref(nu)), -0.5) > 0 &&
|
||||
arf_cmp_2exp_si(arb_midref(acb_realref(z)), -16) > 0 &&
|
||||
arf_cmpabs_2exp_si(arb_midref(acb_realref(nu)), 60) < 0 &&
|
||||
arf_cmpabs_2exp_si(arb_midref(acb_realref(z)), 60) < 0)
|
||||
{
|
||||
arb_hypgeom_bessel_k_integration(acb_realref(res),
|
||||
acb_realref(nu), acb_realref(z), scaled, prec);
|
||||
arb_zero(acb_imagref(res));
|
||||
|
||||
if (acb_rel_accuracy_bits(res) > acb_rel_accuracy_bits(res2) ||
|
||||
(acb_is_finite(res) && !acb_is_finite(res2)))
|
||||
{
|
||||
acb_swap(res, res2);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
acb_swap(res, res2);
|
||||
acb_clear(res2);
|
||||
}
|
||||
|
||||
void
|
||||
acb_hypgeom_bessel_k(acb_t res, const acb_t nu, const acb_t z, slong prec)
|
||||
{
|
||||
_acb_hypgeom_bessel_k(res, nu, z, 0, prec);
|
||||
}
|
||||
|
||||
void
|
||||
acb_hypgeom_bessel_k_scaled(acb_t res, const acb_t nu, const acb_t z, slong prec)
|
||||
{
|
||||
_acb_hypgeom_bessel_k(res, nu, z, 1, prec);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -9,6 +9,7 @@
|
|||
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "arb_hypgeom.h"
|
||||
#include "acb_hypgeom.h"
|
||||
|
||||
int
|
||||
|
|
@ -390,7 +391,7 @@ _acb_is_nonnegative_real(const acb_t z)
|
|||
}
|
||||
|
||||
void
|
||||
acb_hypgeom_gamma_upper(acb_t res, const acb_t s, const acb_t z, int regularized, slong prec)
|
||||
acb_hypgeom_gamma_upper_nointegration(acb_t res, const acb_t s, const acb_t z, int regularized, slong prec)
|
||||
{
|
||||
if (!acb_is_finite(s) || !acb_is_finite(z))
|
||||
{
|
||||
|
|
@ -525,3 +526,46 @@ acb_hypgeom_gamma_upper(acb_t res, const acb_t s, const acb_t z, int regularized
|
|||
}
|
||||
}
|
||||
|
||||
void
|
||||
acb_hypgeom_gamma_upper(acb_t res, const acb_t s, const acb_t z, int regularized, slong prec)
|
||||
{
|
||||
acb_t res2;
|
||||
slong acc, max, t;
|
||||
|
||||
acb_init(res2);
|
||||
|
||||
acb_hypgeom_gamma_upper_nointegration(res2, s, z, regularized, prec);
|
||||
|
||||
acc = acb_rel_accuracy_bits(res2);
|
||||
|
||||
if (acc < 0.5 * prec)
|
||||
{
|
||||
max = prec;
|
||||
t = acb_rel_accuracy_bits(z);
|
||||
max = FLINT_MIN(max, t);
|
||||
t = acb_rel_accuracy_bits(s);
|
||||
max = FLINT_MIN(max, t);
|
||||
|
||||
if (max > 2 && acc < 0.5 * max)
|
||||
{
|
||||
if (acb_is_real(s) && acb_is_real(z) &&
|
||||
arf_cmp_2exp_si(arb_midref(acb_realref(z)), -16) > 0 &&
|
||||
arf_cmpabs_2exp_si(arb_midref(acb_realref(s)), 60) < 0 &&
|
||||
arf_cmpabs_2exp_si(arb_midref(acb_realref(z)), 60) < 0)
|
||||
{
|
||||
arb_hypgeom_gamma_upper_integration(acb_realref(res),
|
||||
acb_realref(s), acb_realref(z), regularized, prec);
|
||||
arb_zero(acb_imagref(res));
|
||||
|
||||
if (acb_rel_accuracy_bits(res) > acb_rel_accuracy_bits(res2) ||
|
||||
(acb_is_finite(res) && !acb_is_finite(res2)))
|
||||
{
|
||||
acb_swap(res, res2);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
acb_swap(res, res2);
|
||||
acb_clear(res2);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -36,9 +36,16 @@ int main()
|
|||
prec = 2 + n_randint(state, 500);
|
||||
scaled = n_randint(state, 2);
|
||||
|
||||
acb_randtest_param(nu, state, 1 + n_randint(state, 1000), 1 + n_randint(state, 10));
|
||||
acb_randtest_param(nu, state, 1 + n_randint(state, 1000), 1 + n_randint(state, 20));
|
||||
acb_randtest(z, state, 1 + n_randint(state, 1000), 1 + n_randint(state, 100));
|
||||
|
||||
if (n_randint(state, 2) && prec <= 200)
|
||||
{
|
||||
arb_zero(acb_imagref(nu));
|
||||
arb_zero(acb_imagref(z));
|
||||
arb_abs(acb_realref(nu), acb_realref(nu));
|
||||
}
|
||||
|
||||
switch (n_randint(state, 3))
|
||||
{
|
||||
case 0:
|
||||
|
|
|
|||
|
|
@ -51,6 +51,13 @@ int main()
|
|||
if (n_randint(state, 4) == 0)
|
||||
arb_zero(acb_imagref(z));
|
||||
|
||||
if (n_randint(state, 2) && prec0 <= 200 && prec1 <= 200 && prec2 <= 200)
|
||||
{
|
||||
arb_zero(acb_imagref(nu0));
|
||||
arb_zero(acb_imagref(z));
|
||||
arb_abs(acb_realref(nu0), acb_realref(nu0));
|
||||
}
|
||||
|
||||
acb_sub_ui(nu1, nu0, 1, prec0);
|
||||
acb_sub_ui(nu2, nu0, 2, prec0);
|
||||
|
||||
|
|
|
|||
|
|
@ -129,6 +129,9 @@ void arb_hypgeom_bessel_k(arb_t res, const arb_t nu, const arb_t z, slong prec);
|
|||
void arb_hypgeom_bessel_i_scaled(arb_t res, const arb_t nu, const arb_t z, slong prec);
|
||||
void arb_hypgeom_bessel_k_scaled(arb_t res, const arb_t nu, const arb_t z, slong prec);
|
||||
|
||||
void arb_hypgeom_bessel_i_integration(arb_t res, const arb_t nu, const arb_t z, int scaled, slong prec);
|
||||
void arb_hypgeom_bessel_k_integration(arb_t res, const arb_t nu, const arb_t z, int scaled, slong prec);
|
||||
|
||||
void arb_hypgeom_airy(arb_t ai, arb_t aip, arb_t bi, arb_t bip, const arb_t z, slong prec);
|
||||
void arb_hypgeom_airy_jet(arb_ptr ai, arb_ptr bi, const arb_t z, slong len, slong prec);
|
||||
void arb_hypgeom_airy_series(arb_poly_t ai, arb_poly_t ai_prime,
|
||||
|
|
@ -153,6 +156,8 @@ void arb_hypgeom_gamma_upper(arb_t res, const arb_t s, const arb_t z, int regula
|
|||
void _arb_hypgeom_gamma_upper_series(arb_ptr g, const arb_t s, arb_srcptr h, slong hlen, int regularized, slong n, slong prec);
|
||||
void arb_hypgeom_gamma_upper_series(arb_poly_t g, const arb_t s, const arb_poly_t h, int regularized, slong n, slong prec);
|
||||
|
||||
void arb_hypgeom_gamma_upper_integration(arb_t res, const arb_t s, const arb_t z, int regularized, slong prec);
|
||||
|
||||
void arb_hypgeom_beta_lower(arb_t res, const arb_t a, const arb_t c, const arb_t z, int regularized, slong prec);
|
||||
void arb_hypgeom_beta_lower_series(arb_poly_t res, const arb_t a, const arb_t b, const arb_poly_t z, int regularized, slong len, slong prec);
|
||||
void _arb_hypgeom_beta_lower_series(arb_ptr res, const arb_t a, const arb_t b, arb_srcptr z, slong zlen, int regularized, slong len, slong prec);
|
||||
|
|
|
|||
68
arb_hypgeom/bessel_i_integration.c
Normal file
68
arb_hypgeom/bessel_i_integration.c
Normal file
|
|
@ -0,0 +1,68 @@
|
|||
/*
|
||||
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"
|
||||
|
||||
void
|
||||
arb_hypgeom_bessel_i_integration(arb_t res, const arb_t nu, const arb_t z, int scaled, slong prec)
|
||||
{
|
||||
arb_t t, a, b, w;
|
||||
|
||||
arb_init(t);
|
||||
arb_init(a);
|
||||
arb_init(b);
|
||||
arb_init(w);
|
||||
|
||||
arb_one(a);
|
||||
arb_mul_2exp_si(a, a, -1);
|
||||
arb_add(a, a, nu, prec);
|
||||
|
||||
arb_mul_2exp_si(b, nu, 1);
|
||||
arb_add_ui(b, b, 1, prec);
|
||||
|
||||
arb_mul_2exp_si(w, z, 1);
|
||||
|
||||
arb_hypgeom_1f1_integration(t, a, b, w, 0, prec);
|
||||
|
||||
if (arb_is_finite(t))
|
||||
{
|
||||
if (!scaled)
|
||||
{
|
||||
arb_neg(a, z);
|
||||
arb_exp(a, a, prec);
|
||||
arb_mul(t, t, a, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
arb_neg(a, z);
|
||||
arb_mul_2exp_si(a, a, 1);
|
||||
arb_exp(a, a, prec);
|
||||
arb_mul(t, t, a, prec);
|
||||
}
|
||||
|
||||
arb_mul_2exp_si(w, z, -1);
|
||||
arb_pow(w, w, nu, prec);
|
||||
arb_mul(t, t, w, prec);
|
||||
|
||||
arb_add_ui(w, nu, 1, prec);
|
||||
arb_rgamma(w, w, prec);
|
||||
arb_mul(res, t, w, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
arb_indeterminate(res);
|
||||
}
|
||||
|
||||
arb_clear(t);
|
||||
arb_clear(a);
|
||||
arb_clear(b);
|
||||
arb_clear(w);
|
||||
}
|
||||
60
arb_hypgeom/bessel_k_integration.c
Normal file
60
arb_hypgeom/bessel_k_integration.c
Normal file
|
|
@ -0,0 +1,60 @@
|
|||
/*
|
||||
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"
|
||||
|
||||
void
|
||||
arb_hypgeom_bessel_k_integration(arb_t res, const arb_t nu, const arb_t z, int scaled, slong prec)
|
||||
{
|
||||
arb_t t, a, b, w;
|
||||
|
||||
arb_init(t);
|
||||
arb_init(a);
|
||||
arb_init(b);
|
||||
arb_init(w);
|
||||
|
||||
arb_one(a);
|
||||
arb_mul_2exp_si(a, a, -1);
|
||||
arb_add(a, a, nu, prec);
|
||||
|
||||
arb_mul_2exp_si(b, nu, 1);
|
||||
arb_add_ui(b, b, 1, prec);
|
||||
|
||||
arb_mul_2exp_si(w, z, 1);
|
||||
|
||||
arb_hypgeom_u_integration(t, a, b, w, prec);
|
||||
|
||||
if (arb_is_finite(t))
|
||||
{
|
||||
if (!scaled)
|
||||
{
|
||||
arb_neg(a, z);
|
||||
arb_exp(a, a, prec);
|
||||
arb_mul(t, t, a, prec);
|
||||
}
|
||||
|
||||
arb_mul_2exp_si(w, z, 1);
|
||||
arb_pow(w, w, nu, prec);
|
||||
arb_mul(res, t, w, prec);
|
||||
|
||||
arb_const_sqrt_pi(w, prec);
|
||||
arb_mul(res, res, w, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
arb_indeterminate(res);
|
||||
}
|
||||
|
||||
arb_clear(t);
|
||||
arb_clear(a);
|
||||
arb_clear(b);
|
||||
arb_clear(w);
|
||||
}
|
||||
53
arb_hypgeom/gamma_upper_integration.c
Normal file
53
arb_hypgeom/gamma_upper_integration.c
Normal file
|
|
@ -0,0 +1,53 @@
|
|||
/*
|
||||
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"
|
||||
|
||||
void
|
||||
arb_hypgeom_gamma_upper_integration(arb_t res, const arb_t s,
|
||||
const arb_t z, int regularized, slong prec)
|
||||
{
|
||||
arb_t t, u;
|
||||
|
||||
arb_init(t);
|
||||
arb_init(u);
|
||||
|
||||
arb_one(t);
|
||||
arb_add_ui(u, s, 1, prec);
|
||||
|
||||
arb_hypgeom_u_integration(u, t, u, z, prec);
|
||||
|
||||
if (arb_is_finite(u))
|
||||
{
|
||||
if (regularized != 2)
|
||||
{
|
||||
arb_pow(t, z, s, prec);
|
||||
arb_mul(u, u, t, prec);
|
||||
|
||||
if (regularized == 1)
|
||||
{
|
||||
arb_rgamma(t, s, prec);
|
||||
arb_mul(u, u, t, prec);
|
||||
}
|
||||
}
|
||||
|
||||
arb_neg(t, z);
|
||||
arb_exp(t, t, prec);
|
||||
arb_mul(res, t, u, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
arb_indeterminate(res);
|
||||
}
|
||||
|
||||
arb_clear(t);
|
||||
arb_clear(u);
|
||||
}
|
||||
82
arb_hypgeom/test/t-bessel_i_integration.c
Normal file
82
arb_hypgeom/test/t-bessel_i_integration.c
Normal file
|
|
@ -0,0 +1,82 @@
|
|||
/*
|
||||
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("bessel_i_integration....");
|
||||
fflush(stdout);
|
||||
|
||||
flint_randinit(state);
|
||||
|
||||
for (iter = 0; iter < 100 * arb_test_multiplier(); iter++)
|
||||
{
|
||||
arb_t a, z, r1, r2;
|
||||
slong prec1, prec2;
|
||||
int scaled;
|
||||
|
||||
prec1 = 2 + n_randint(state, 80);
|
||||
prec2 = 2 + n_randint(state, 150);
|
||||
|
||||
arb_init(a);
|
||||
arb_init(z);
|
||||
arb_init(r1);
|
||||
arb_init(r2);
|
||||
|
||||
scaled = n_randint(state, 2);
|
||||
arb_randtest_precise(a, 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(z, z, n_randint(state, 100), prec1 + 100);
|
||||
|
||||
arb_hypgeom_bessel_i_integration(r1, a, z, scaled, prec1);
|
||||
|
||||
if (arb_is_finite(r1))
|
||||
{
|
||||
if (n_randint(state, 2))
|
||||
arb_hypgeom_bessel_i_integration(r2, a, z, scaled, prec2);
|
||||
else
|
||||
{
|
||||
if (scaled)
|
||||
arb_hypgeom_bessel_i_scaled(r2, a, z, prec2);
|
||||
else
|
||||
arb_hypgeom_bessel_i(r2, a, 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("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();
|
||||
}
|
||||
}
|
||||
|
||||
arb_clear(a);
|
||||
arb_clear(z);
|
||||
arb_clear(r1);
|
||||
arb_clear(r2);
|
||||
}
|
||||
|
||||
flint_randclear(state);
|
||||
flint_cleanup();
|
||||
flint_printf("PASS\n");
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
82
arb_hypgeom/test/t-bessel_k_integration.c
Normal file
82
arb_hypgeom/test/t-bessel_k_integration.c
Normal file
|
|
@ -0,0 +1,82 @@
|
|||
/*
|
||||
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("bessel_k_integration....");
|
||||
fflush(stdout);
|
||||
|
||||
flint_randinit(state);
|
||||
|
||||
for (iter = 0; iter < 100 * arb_test_multiplier(); iter++)
|
||||
{
|
||||
arb_t a, z, r1, r2;
|
||||
slong prec1, prec2;
|
||||
int scaled;
|
||||
|
||||
prec1 = 2 + n_randint(state, 80);
|
||||
prec2 = 2 + n_randint(state, 150);
|
||||
|
||||
arb_init(a);
|
||||
arb_init(z);
|
||||
arb_init(r1);
|
||||
arb_init(r2);
|
||||
|
||||
scaled = n_randint(state, 2);
|
||||
arb_randtest_precise(a, 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(z, z, n_randint(state, 100), prec1 + 100);
|
||||
|
||||
arb_hypgeom_bessel_k_integration(r1, a, z, scaled, prec1);
|
||||
|
||||
if (arb_is_finite(r1))
|
||||
{
|
||||
if (n_randint(state, 2))
|
||||
arb_hypgeom_bessel_k_integration(r2, a, z, scaled, prec2);
|
||||
else
|
||||
{
|
||||
if (scaled)
|
||||
arb_hypgeom_bessel_k_scaled(r2, a, z, prec2);
|
||||
else
|
||||
arb_hypgeom_bessel_k(r2, a, 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("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();
|
||||
}
|
||||
}
|
||||
|
||||
arb_clear(a);
|
||||
arb_clear(z);
|
||||
arb_clear(r1);
|
||||
arb_clear(r2);
|
||||
}
|
||||
|
||||
flint_randclear(state);
|
||||
flint_cleanup();
|
||||
flint_printf("PASS\n");
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
111
arb_hypgeom/test/t-gamma_upper_integration.c
Normal file
111
arb_hypgeom/test/t-gamma_upper_integration.c
Normal file
|
|
@ -0,0 +1,111 @@
|
|||
/*
|
||||
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("gamma_upper_integration....");
|
||||
fflush(stdout);
|
||||
|
||||
flint_randinit(state);
|
||||
|
||||
for (iter = 0; iter < 100 * arb_test_multiplier(); iter++)
|
||||
{
|
||||
arb_t a, z, r1, r2;
|
||||
slong prec1, prec2;
|
||||
int regularized;
|
||||
|
||||
prec1 = 2 + n_randint(state, 80);
|
||||
prec2 = 2 + n_randint(state, 150);
|
||||
|
||||
arb_init(a);
|
||||
arb_init(z);
|
||||
arb_init(r1);
|
||||
arb_init(r2);
|
||||
|
||||
regularized = n_randint(state, 3);
|
||||
arb_randtest_precise(a, 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(z, z, n_randint(state, 100), prec1 + 100);
|
||||
|
||||
arb_hypgeom_gamma_upper_integration(r1, a, z, regularized, prec1);
|
||||
|
||||
if (arb_is_finite(r1))
|
||||
{
|
||||
if (n_randint(state, 2))
|
||||
arb_hypgeom_gamma_upper_integration(r2, a, z, regularized, prec2);
|
||||
else
|
||||
arb_hypgeom_gamma_upper(r2, a, z, regularized, 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("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, "5.25", prec1);
|
||||
arb_set_str(z, "1.125", prec1);
|
||||
|
||||
arb_hypgeom_gamma_upper_integration(r1, a, z, 0, prec1);
|
||||
arb_set_str(r2, "35.072486673952527119466762267780775062420879980812938235868074922070540736636294202721367124421387 +/- 5.01e-97", 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_set_str(a, "-23515.25", prec1);
|
||||
arb_set_str(z, "5118.125", prec1);
|
||||
|
||||
arb_hypgeom_gamma_upper_integration(r1, a, z, 0, prec1);
|
||||
arb_set_str(r2, "1.2585189324946269645868452847608400602717904014578585865177198314407671242201776909963228597279e-89448 +/- 8.86e-89543", 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(z);
|
||||
arb_clear(r1);
|
||||
arb_clear(r2);
|
||||
}
|
||||
|
||||
flint_randclear(state);
|
||||
flint_cleanup();
|
||||
flint_printf("PASS\n");
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
|
|
@ -86,6 +86,21 @@ int main()
|
|||
flint_printf("r2 = "); arb_printd(r2, 100); flint_printf("\n\n");
|
||||
flint_abort();
|
||||
}
|
||||
|
||||
arb_set_str(a, "1.0", prec1);
|
||||
arb_set_str(b, "1351.25", prec1);
|
||||
arb_set_str(z, "5118.125", prec1);
|
||||
|
||||
arb_hypgeom_u_integration(r1, a, b, z, prec1);
|
||||
arb_set_str(r2, "0.0002653059837025369361090415756761601958818507141969057844945316037867908521084665825645386229229474 +/- 5.74e-101", prec1);
|
||||
|
||||
if (!arb_overlaps(r1, r2) || arb_rel_accuracy_bits(r1) < arb_rel_accuracy_bits(r2) - 10)
|
||||
{
|
||||
flint_printf("FAIL: overlap (2)\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);
|
||||
|
|
|
|||
|
|
@ -31,7 +31,12 @@ 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))));
|
||||
|
||||
if (a1.a == 0.0 && a1.b == 0.0)
|
||||
Y = di_interval(0.0, 0.0);
|
||||
else
|
||||
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));
|
||||
|
|
@ -46,7 +51,11 @@ 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)));
|
||||
if (a1.a == 0.0 && a1.b == 0.0)
|
||||
Y = di_interval(0.0, 0.0);
|
||||
else
|
||||
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)
|
||||
|
|
@ -206,7 +215,8 @@ integrand(acb_ptr out, const acb_t t, void * param, slong order, slong prec)
|
|||
|
||||
if (order == 1)
|
||||
{
|
||||
if (!arb_is_positive(acb_realref(t)))
|
||||
if (!(arb_is_positive(acb_realref(t)) || arb_is_zero(a1)) ||
|
||||
!arb_is_positive(acb_realref(v)))
|
||||
acb_indeterminate(out);
|
||||
else
|
||||
integrand_wide_bound5(out, t, a1, ba1, z, prec);
|
||||
|
|
@ -222,6 +232,7 @@ integrand(acb_ptr out, const acb_t t, void * param, slong order, slong prec)
|
|||
|
||||
/* t^(a-1) */
|
||||
acb_my_pow_arb(u, t, a1, prec);
|
||||
|
||||
/* (1+t)^(b-a-1) */
|
||||
acb_pow_arb(v, v, ba1, prec);
|
||||
|
||||
|
|
@ -234,8 +245,15 @@ integrand(acb_ptr out, const acb_t t, void * param, slong order, slong prec)
|
|||
acb_neg(s, s);
|
||||
|
||||
/* t^(a-1) */
|
||||
acb_log(u, t, prec);
|
||||
acb_mul_arb(u, u, a1, prec);
|
||||
if (arb_is_zero(a1))
|
||||
{
|
||||
acb_zero(u);
|
||||
}
|
||||
else
|
||||
{
|
||||
acb_log(u, t, prec);
|
||||
acb_mul_arb(u, u, a1, prec);
|
||||
}
|
||||
|
||||
/* (1+t)^(b-a-1) */
|
||||
acb_log(v, v, prec);
|
||||
|
|
@ -278,6 +296,10 @@ estimate_magnitude(mag_t res, const arb_t ra, const arb_t rb, const arb_t rz)
|
|||
t2 = 1e-8;
|
||||
}
|
||||
|
||||
/* todo: better estimate when peak is at 0 */
|
||||
t1 = FLINT_MAX(t1, 1e-8);
|
||||
t2 = FLINT_MAX(t2, 1e-8);
|
||||
|
||||
m = -1e300;
|
||||
|
||||
if (t1 > 0.0)
|
||||
|
|
|
|||
|
|
@ -286,6 +286,11 @@ Incomplete gamma and beta functions
|
|||
intended for internal use; :func:`arb_hypgeom_expint` is the intended
|
||||
interface for computing the exponential integral).
|
||||
|
||||
.. function:: void arb_hypgeom_gamma_upper_integration(arb_t res, const arb_t s, const arb_t z, int regularized, slong prec)
|
||||
|
||||
Computes the upper incomplete gamma function using numerical
|
||||
integration.
|
||||
|
||||
.. function:: void _arb_hypgeom_gamma_upper_series(arb_ptr res, const arb_t s, arb_srcptr z, slong zlen, int regularized, slong n, slong prec)
|
||||
void arb_hypgeom_gamma_upper_series(arb_poly_t res, const arb_t s, const arb_poly_t z, int regularized, slong n, slong prec)
|
||||
|
||||
|
|
@ -329,6 +334,7 @@ Incomplete gamma and beta functions
|
|||
The underscore method requires positive lengths and does not support
|
||||
aliasing.
|
||||
|
||||
|
||||
Internal evaluation functions
|
||||
................................................................................
|
||||
|
||||
|
|
@ -495,6 +501,11 @@ Bessel functions
|
|||
|
||||
Computes the function `e^{z} K_{\nu}(z)`.
|
||||
|
||||
.. function:: void arb_hypgeom_bessel_i_integration(arb_t res, const arb_t nu, const arb_t z, int scaled, slong prec)
|
||||
void arb_hypgeom_bessel_k_integration(arb_t res, const arb_t nu, const arb_t z, int scaled, slong prec)
|
||||
|
||||
Computes the modified Bessel functions using numerical integration.
|
||||
|
||||
Airy functions
|
||||
-------------------------------------------------------------------------------
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue