diff --git a/acb_hypgeom/2f1.c b/acb_hypgeom/2f1.c
index 2754520e..ea28531e 100644
--- a/acb_hypgeom/2f1.c
+++ b/acb_hypgeom/2f1.c
@@ -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);
+}
diff --git a/acb_hypgeom/m.c b/acb_hypgeom/m.c
index bae4f2d3..fdcb059c 100644
--- a/acb_hypgeom/m.c
+++ b/acb_hypgeom/m.c
@@ -9,6 +9,7 @@
(at your option) any later version. See .
*/
+#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);
}
-
diff --git a/acb_hypgeom/test/t-2f1.c b/acb_hypgeom/test/t-2f1.c
index 03bd76ce..ccb7dfec 100644
--- a/acb_hypgeom/test/t-2f1.c
+++ b/acb_hypgeom/test/t-2f1.c
@@ -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);
diff --git a/acb_hypgeom/test/t-m.c b/acb_hypgeom/test/t-m.c
index ddef8aef..d2423d0c 100644
--- a/acb_hypgeom/test/t-m.c
+++ b/acb_hypgeom/test/t-m.c
@@ -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);
diff --git a/acb_hypgeom/test/t-u.c b/acb_hypgeom/test/t-u.c
index 81b642dd..1dd5d46f 100644
--- a/acb_hypgeom/test/t-u.c
+++ b/acb_hypgeom/test/t-u.c
@@ -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);
diff --git a/acb_hypgeom/u.c b/acb_hypgeom/u.c
index a0c5266e..c13ac669 100644
--- a/acb_hypgeom/u.c
+++ b/acb_hypgeom/u.c
@@ -9,6 +9,7 @@
(at your option) any later version. See .
*/
+#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);
+}
diff --git a/arb_hypgeom/2f1_integration.c b/arb_hypgeom/2f1_integration.c
index 4790259f..813aba3d 100644
--- a/arb_hypgeom/2f1_integration.c
+++ b/arb_hypgeom/2f1_integration.c
@@ -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))));
- 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))));
+
+ 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)));
- Y = di_fast_div(cb1, di_fast_add(di_fast_sqr(di_fast_sub_d(u, 1.0)), 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) */
- acb_log(u, u, prec);
- acb_mul_arb(u, u, nega, prec);
+ 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);