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

use integration as fallback in some cases (experimental)

Browse source
This commit is contained in:
fredrik 2021-12-10 13:36:27 +01:00
parent 956013d9e2
commit 2aa02b9171
7 changed files with 188 additions and 9 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

53
acb_hypgeom/2f1.c
View file

@ -10,6 +10,7 @@
*/
#include "acb_hypgeom.h"
#include "arb_hypgeom.h"
static void
_acb_hypgeom_2f1r_reduced(acb_t res,
@ -31,7 +32,7 @@ _acb_hypgeom_2f1r_reduced(acb_t res,
}
void
acb_hypgeom_2f1(acb_t res, const acb_t a, const acb_t b,
acb_hypgeom_2f1_nointegration(acb_t res, const acb_t a, const acb_t b,
const acb_t c, const acb_t z, int flags, slong prec)
{
int algorithm, regularized;
@ -190,3 +191,53 @@ acb_hypgeom_2f1(acb_t res, const acb_t a, const acb_t b,
}
}
void
acb_hypgeom_2f1(acb_t res, const acb_t a, const acb_t b,
const acb_t c, const acb_t z, int flags, slong prec)
{
acb_t res2;
slong acc, max, t;
acb_init(res2);
acb_hypgeom_2f1_nointegration(res2, a, b, c, z, flags, 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(a);
max = FLINT_MIN(max, t);
t = acb_rel_accuracy_bits(b);
max = FLINT_MIN(max, t);
t = acb_rel_accuracy_bits(c);
max = FLINT_MIN(max, t);
if (max > 2 && acc < 0.5 * max)
{
if (acb_is_real(a) && acb_is_real(b) && acb_is_real(c) && acb_is_real(z) &&
arf_cmpabs_2exp_si(arb_midref(acb_realref(a)), 60) < 0 &&
arf_cmpabs_2exp_si(arb_midref(acb_realref(b)), 60) < 0 &&
arf_cmpabs_2exp_si(arb_midref(acb_realref(c)), 60) < 0 &&
arf_cmpabs_2exp_si(arb_midref(acb_realref(z)), 60) < 0)
{
arb_hypgeom_2f1_integration(acb_realref(res),
acb_realref(a), acb_realref(b), acb_realref(c), acb_realref(z), flags, 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);
}

50
acb_hypgeom/m.c
View file

@ -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
@ -292,7 +293,7 @@ acb_hypgeom_m_choose(int * asymp, int * kummer, slong * wp,
}
void
acb_hypgeom_m(acb_t res, const acb_t a, const acb_t b, const acb_t z, int regularized, slong prec)
acb_hypgeom_m_nointegration(acb_t res, const acb_t a, const acb_t b, const acb_t z, int regularized, slong prec)
{
int asymp, kummer;
slong wp;
@ -311,9 +312,54 @@ acb_hypgeom_m(acb_t res, const acb_t a, const acb_t b, const acb_t z, int regula
acb_set_round(res, res, prec);
}
void
acb_hypgeom_m(acb_t res, const acb_t a, const acb_t b, const acb_t z, int regularized, slong prec)
{
acb_t res2;
slong acc, max, t;
acb_init(res2);
acb_hypgeom_m_nointegration(res2, a, b, 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(a);
max = FLINT_MIN(max, t);
t = acb_rel_accuracy_bits(b);
max = FLINT_MIN(max, t);
if (max > 2 && acc < 0.5 * max)
{
if (acb_is_real(a) && acb_is_real(b) && acb_is_real(z) &&
arf_cmpabs_2exp_si(arb_midref(acb_realref(a)), 60) < 0 &&
arf_cmpabs_2exp_si(arb_midref(acb_realref(b)), 60) < 0 &&
arf_cmpabs_2exp_si(arb_midref(acb_realref(z)), 60) < 0)
{
arb_hypgeom_1f1_integration(acb_realref(res),
acb_realref(a), acb_realref(b), 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);
}
void
acb_hypgeom_1f1(acb_t res, const acb_t a, const acb_t b, const acb_t z, int regularized, slong prec)
{
acb_hypgeom_m(res, a, b, z, regularized, prec);
}

3
acb_hypgeom/test/t-2f1.c
View file

@ -87,6 +87,9 @@ int main()
acb_randtest(w1, state, 1 + n_randint(state, 400), 1 + n_randint(state, ebits));
acb_randtest(w2, state, 1 + n_randint(state, 400), 1 + n_randint(state, ebits));
if (n_randint(state, 2))
arb_zero(acb_imagref(z));
reg1 = n_randint(state, 2);
reg2 = n_randint(state, 2);

7
acb_hypgeom/test/t-m.c
View file

@ -55,6 +55,13 @@ int main()
acb_randtest(w2, state, 1 + n_randint(state, 1000), 1 + n_randint(state, ebits));
regularized = n_randint(state, 2);
if (prec0 <= 300 && prec1 <= 300 && prec2 <= 300 && n_randint(state, 2))
{
arb_zero(acb_imagref(a0));
arb_zero(acb_imagref(b));
arb_zero(acb_imagref(z));
}
acb_add_ui(a1, a0, 1, prec0);
acb_add_ui(a2, a0, 2, prec0);

7
acb_hypgeom/test/t-u.c
View file

@ -59,6 +59,13 @@ int main()
acb_randtest(w1, state, 1 + n_randint(state, 1000), 1 + n_randint(state, 100));
acb_randtest(w2, state, 1 + n_randint(state, 1000), 1 + n_randint(state, 100));
if (prec0 <= 300 && prec1 <= 300 && prec2 <= 300 && n_randint(state, 2))
{
arb_zero(acb_imagref(a0));
arb_zero(acb_imagref(b));
arb_zero(acb_imagref(z));
}
acb_add_ui(a1, a0, 1, prec0);
acb_add_ui(a2, a0, 2, prec0);

49
acb_hypgeom/u.c
View file

@ -9,6 +9,7 @@
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "arb_hypgeom.h"
#include "acb_hypgeom.h"
static void
@ -364,7 +365,7 @@ acb_hypgeom_u_choose(int * asymp, slong * wp,
}
void
acb_hypgeom_u(acb_t res, const acb_t a, const acb_t b, const acb_t z, slong prec)
acb_hypgeom_u_nointegration(acb_t res, const acb_t a, const acb_t b, const acb_t z, slong prec)
{
acb_t t;
arf_srcptr av, tv;
@ -424,3 +425,49 @@ acb_hypgeom_u(acb_t res, const acb_t a, const acb_t b, const acb_t z, slong prec
acb_clear(t);
}
void
acb_hypgeom_u(acb_t res, const acb_t a, const acb_t b, const acb_t z, slong prec)
{
acb_t res2;
slong acc, max, t;
acb_init(res2);
acb_hypgeom_u_nointegration(res2, a, b, z, 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(a);
max = FLINT_MIN(max, t);
t = acb_rel_accuracy_bits(b);
max = FLINT_MIN(max, t);
if (max > 2 && acc < 0.5 * max)
{
if (acb_is_real(a) && acb_is_real(b) && acb_is_real(z) &&
arb_is_positive(acb_realref(z)) &&
arf_cmpabs_2exp_si(arb_midref(acb_realref(a)), 60) < 0 &&
arf_cmpabs_2exp_si(arb_midref(acb_realref(b)), 60) < 0 &&
arf_cmpabs_2exp_si(arb_midref(acb_realref(z)), 60) < 0)
{
arb_hypgeom_u_integration(acb_realref(res),
acb_realref(a), acb_realref(b), acb_realref(z), 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);
}

20
arb_hypgeom/2f1_integration.c
View file

@ -39,7 +39,12 @@ di_integrand_edge(di_t u, di_t v, di_t b1, di_t cb1, di_t nega, di_t z)
di_t X, Y, Z;
X = di_fast_mul(b1, di_fast_log_nonnegative(di_fast_add(di_fast_sqr(u), di_fast_sqr(v))));
if (cb1.a == 0 && cb1.b == 0)
Y = di_interval(0.0, 0.0);
else
Y = di_fast_mul(cb1, di_fast_log_nonnegative(di_fast_add(di_fast_sqr(di_fast_sub_d(u, 1.0)), di_fast_sqr(v))));
Z = di_fast_mul(nega, di_fast_log_nonnegative(di_fast_add(
di_fast_sqr(di_fast_mul(v, z)),
di_fast_sqr(di_fast_sub_d(di_fast_mul(u, z), 1.0)))));
@ -59,7 +64,12 @@ di_integrand_edge_diff(di_t u, di_t v, di_t b1, di_t cb1, di_t nega, di_t z, int
uz1 = di_fast_sub_d(di_fast_mul(u, z), 1.0);
X = di_fast_div(b1, di_fast_add(di_fast_sqr(u), di_fast_sqr(v)));
if (cb1.a == 0 && cb1.b == 0)
Y = di_interval(0.0, 0.0);
else
Y = di_fast_div(cb1, di_fast_add(di_fast_sqr(di_fast_sub_d(u, 1.0)), di_fast_sqr(v)));
Z = di_fast_div(nega, di_fast_add(di_fast_sqr(di_fast_mul(v, z)), di_fast_sqr(uz1)));
if (which == 0)
@ -236,7 +246,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)) || !arb_is_positive(acb_realref(u)) || !arb_is_positive(acb_realref(v)))
if (!arb_is_positive(acb_realref(t)) || !arb_is_positive(acb_realref(u)) ||
(!(arb_is_positive(acb_realref(v)) || arb_is_zero(cb1))))
acb_indeterminate(out);
else
{
@ -281,8 +292,15 @@ integrand(acb_ptr out, const acb_t t, void * param, slong order, slong prec)
acb_log(v, v, prec);
acb_mul_arb(v, v, cb1, prec);
/* (1-zt)^(-a) */
if (arb_is_zero(nega))
{
acb_zero(u);
}
else
{
acb_log(u, u, prec);
acb_mul_arb(u, u, nega, prec);
}
acb_add(out, s, u, prec);
acb_add(out, out, v, prec);