complete the 2F1 implementation

acb_hypgeom.h

@@ -147,11 +147,17 @@ void acb_hypgeom_2f1_continuation(acb_t res0, acb_t res1,
 void acb_hypgeom_2f1_series_direct(acb_poly_t res, const acb_poly_t a, const acb_poly_t b,
     const acb_poly_t c, const acb_poly_t z, int regularized, long len, long prec);
 void acb_hypgeom_2f1_direct(acb_t res, const acb_t a, const acb_t b, const acb_t c, const acb_t z, int regularized, long prec);
-void acb_hypgeom_2f1_pfaff(acb_t res, const acb_t a, const acb_t b, const acb_t c, const acb_t z, int regularized, long prec);
-void acb_hypgeom_2f1_inf(acb_t res, const acb_t a, const acb_t b, const acb_t c, const acb_t z, int regularized, long prec);
-void acb_hypgeom_2f1_inf_limit(acb_t res, const acb_t a, const acb_t b, const acb_t c, const acb_t z, int regularized, long prec);
+
+void acb_hypgeom_2f1_transform(acb_t res, const acb_t a, const acb_t b,
+    const acb_t c, const acb_t z, int regularized, int which, long prec);
+
+void acb_hypgeom_2f1_transform_limit(acb_t res, const acb_t a, const acb_t b,
+    const acb_t c, const acb_t z, int regularized, int which, long prec);
+
 void acb_hypgeom_2f1_corner(acb_t res, const acb_t a, const acb_t b, const acb_t c, const acb_t z, int regularized, long prec);
 
+int acb_hypgeom_2f1_choose(const acb_t z);
+
 void acb_hypgeom_2f1(acb_t res, const acb_t a, const acb_t b, const acb_t c, const acb_t z, int regularized, long prec);
 
 #ifdef __cplusplus

acb_hypgeom/2f1.c

@@ -29,7 +29,7 @@ void
 acb_hypgeom_2f1(acb_t res, const acb_t a, const acb_t b,
         const acb_t c, const acb_t z, int regularized, long prec)
 {
-    double x, y;
+    int algorithm;
 
     if (!acb_is_finite(a) || !acb_is_finite(b) || !acb_is_finite(c) || !acb_is_finite(z))
     {
@@ -103,30 +103,19 @@ acb_hypgeom_2f1(acb_t res, const acb_t a, const acb_t b,
         return;
     }
 
-    x = arf_get_d(arb_midref(acb_realref(z)), ARF_RND_DOWN);
-    y = arf_get_d(arb_midref(acb_imagref(z)), ARF_RND_DOWN);
+    algorithm = acb_hypgeom_2f1_choose(z);
 
-    x = FLINT_MAX(FLINT_MIN(x, 1e10), -1e10);
-    y = FLINT_MAX(FLINT_MIN(y, 1e10), -1e10);
-
-    if (x*x + y*y < 0.7)   /* series in z */
+    if (algorithm == 0)
     {
         acb_hypgeom_2f1_direct(res, a, b, c, z, regularized, prec);
     }
-    else if (x < 0.7 && (x*x + y*y)/((x-1)*(x-1) + y*y) < 0.7)  /* series in z/(z-1) */
-    {
-        acb_hypgeom_2f1_pfaff(res, a, b, c, z, regularized, prec);
-    }
-    else if (x*x + y*y > 1.1)  /* series in 1/z */
+    else if (algorithm >= 1 && algorithm <= 5)
     {
-        acb_hypgeom_2f1_inf(res, a, b, c, z, regularized, prec);
+        acb_hypgeom_2f1_transform(res, a, b, c, z, regularized, algorithm, prec);
     }
     else
     {
-        if (x*x + y*y < 1.0)
-            acb_hypgeom_2f1_direct(res, a, b, c, z, regularized, prec);
-        else
-            acb_hypgeom_2f1_inf(res, a, b, c, z, regularized, prec);
+        acb_hypgeom_2f1_corner(res, a, b, c, z, regularized, prec);
     }
 }
 

acb_hypgeom/2f1_pfaff.c → acb_hypgeom/2f1_choose.c

@@ -25,28 +25,49 @@
 
 #include "acb_hypgeom.h"
 
-void
-acb_hypgeom_2f1_pfaff(acb_t res, const acb_t a, const acb_t b,
-    const acb_t c, const acb_t z, int regularized, long prec)
+#define ALWAYS1 (0.25 * 0.25)
+#define ALWAYS2 (0.75 * 0.75)
+#define LIMIT (0.75 * 0.75)
+
+int
+acb_hypgeom_2f1_choose(const acb_t z)
 {
-    acb_t t, u, v;
+    double x, y;
+    double mag[7];
+    int i, pick;
+
+    x = arf_get_d(arb_midref(acb_realref(z)), ARF_RND_DOWN);
+    y = arf_get_d(arb_midref(acb_imagref(z)), ARF_RND_DOWN);
+
+    x = FLINT_MAX(FLINT_MIN(x, 1e10), -1e10);
+    y = FLINT_MAX(FLINT_MIN(y, 1e10), -1e10);
+
+    mag[0] = x*x + y*y;  /* |z|^2 */
+    mag[4] = (1.0-x)*(1.0-x) + y*y;              /* |1-z|^2 */
+
+    if (mag[0] <= ALWAYS1)   return 0;
+
+    mag[1] = mag[0] / FLINT_MAX(mag[4], 1e-10);  /* |z/(z-1)|^2 */
+
+    if (mag[1] <= ALWAYS1)   return 1;
+
+    if (mag[0] <= ALWAYS2 || mag[1] <= ALWAYS2)
+        return mag[0] <= mag[1] ? 0 : 1;
 
-    acb_init(t);
-    acb_init(u);
-    acb_init(v);
+    mag[2] = 1.0 / mag[0];                    /* |1/z|^2 */
+    mag[3] = 1.0 / FLINT_MAX(mag[4], 1e-10);  /* 1/|1-z|^2 */
+    mag[5] = mag[4] / mag[0];                 /* |1-1/z|^2 = |(1-z)/z|^2 */
 
-    acb_sub_ui(t, z, 1, prec); /* t = z-1 */
-    acb_div(u, z, t, prec); /* u = z/(z-1) */
-    acb_neg(t, t);
-    acb_neg(v, a);
-    acb_pow(t, t, v, prec); /* t = (1-z)^-a */
-    acb_sub(v, c, b, prec); /* v = c-b */
+    pick = 0;
+    for (i = 1; i < 6; i++)
+    {
+        if (mag[i] < mag[pick])
+            pick = i;
+    }
 
-    acb_hypgeom_2f1_direct(res, a, v, c, u, regularized, prec);
-    acb_mul(res, res, t, prec);
+    if (mag[pick] <= LIMIT)
+        return pick;
 
-    acb_clear(t);
-    acb_clear(u);
-    acb_clear(v);
+    return 6;
 }
 

acb_hypgeom/2f1_corner.c

@@ -42,11 +42,17 @@ acb_hypgeom_2f1_corner(acb_t res, const acb_t a, const acb_t b,
     upper = arb_is_positive(acb_imagref(z));
 
     /* 0 -> 0.5 +/- 0.5i -> 0.5 +/- 0.75i -> z */
+#if 0
     acb_set_d_d(z1, 0.5, upper ? 0.5 : -0.5);
     acb_set_d_d(z2, 0.5, upper ? 0.75 : -0.75);
+#else
+    acb_set_d_d(z1, 0.375, upper ? 0.625 : -0.625);
+    acb_set_d_d(z2, 0.5, upper ? 0.8125 : -0.8125);
+#endif
 
-    acb_hypgeom_2f1(f1, a, b, c, z1, regularized, prec);
-    acb_hypgeom_2f1(f2, aa, bb, cc, z1, regularized, prec);
+    acb_hypgeom_2f1_direct(f1, a, b, c, z1, regularized, prec);
+
+    acb_hypgeom_2f1_direct(f2, aa, bb, cc, z1, regularized, prec);
     acb_mul(f2, f2, a, prec);
     acb_mul(f2, f2, b, prec);
     if (!regularized)