more hypergeometric helper functions

acb_hypgeom.h

@@ -35,6 +35,15 @@ extern "C" {
 void acb_hypgeom_pfq_bound_factor(mag_t C,
     acb_srcptr a, long p, acb_srcptr b, long q, const acb_t z, ulong n);
 
+long acb_hypgeom_pfq_choose_n(acb_srcptr a, long p,
+                         acb_srcptr b, long q, const acb_t z, long prec);
+
+void acb_hypgeom_pfq_sum_forward(acb_t s, acb_t t, acb_srcptr a, long p, acb_srcptr b, long q,
+    const acb_t z, long n, long prec);
+
+void acb_hypgeom_pfq_sum_rs(acb_t s, acb_t t, acb_srcptr a, long p, acb_srcptr b, long q,
+    const acb_t z, long n, long prec);
+
 void acb_hypgeom_pfq_sum(acb_t s, acb_t t, acb_srcptr a, long p, acb_srcptr b, long q,
     const acb_t z, long n, long prec);
 

acb_hypgeom/pfq_choose_n.c

@@ -0,0 +1,134 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2014 Fredrik Johansson
+
+******************************************************************************/
+
+#include "acb_hypgeom.h"
+
+double mag_get_log2_d_approx(const mag_t x);
+
+long
+acb_hypgeom_pfq_choose_n(acb_srcptr a, long p,
+                         acb_srcptr b, long q, const acb_t z, long prec)
+{
+    long k, n, minimum_n, maximum_n;
+    mag_t zmag;
+
+    /* todo: first test for z = 0 and nonpositive integers */
+
+    double t, u;
+    double log2_z;
+    double log2_term;
+    double log2_factor;
+    double log2_term_max;
+
+    double * are;
+    double * aim;
+    double * bre;
+    double * bim;
+
+    mag_init(zmag);
+
+    are = flint_malloc(sizeof(double) * 2 * (p + q));
+    aim = are + p;
+    bre = aim + p;
+    bim = bre + q;
+
+    acb_get_mag(zmag, z);
+    log2_z = mag_get_log2_d_approx(zmag);
+
+    for (k = 0; k < p; k++)
+    {
+        are[k] = arf_get_d(arb_midref(acb_realref(a + k)), ARF_RND_DOWN);
+        aim[k] = arf_get_d(arb_midref(acb_imagref(a + k)), ARF_RND_DOWN);
+    }
+
+    minimum_n = 1;
+    maximum_n = 50 + 2 * prec;
+
+    for (k = 0; k < q; k++)
+    {
+        bre[k] = arf_get_d(arb_midref(acb_realref(b + k)), ARF_RND_DOWN);
+        bim[k] = arf_get_d(arb_midref(acb_imagref(b + k)), ARF_RND_DOWN);
+
+        if (bre[k] <= 0.25)
+        {
+            minimum_n = FLINT_MAX(minimum_n, 2 - bre[k]);
+        }
+    }
+
+    n = 1;
+    log2_term = 0.0;
+    log2_term_max = log2_term;
+
+    while (n <= maximum_n && minimum_n < maximum_n)
+    {
+        if (log2_term < log2_term_max - prec - 4 && n >= minimum_n)
+            break;
+
+        t = 1.0;
+
+        for (k = 0; k < FLINT_MAX(p, q); k++)
+        {
+            if (k < p)
+            {
+                u = (are[k] + n) * (are[k] + n) + (aim[k] * aim[k]);
+                u = FLINT_ABS(u);
+
+                if (u < 1e-8 || u > 1e100 || t > 1e100)
+                    goto somethingstrange;
+
+                t *= u;
+            }
+
+            if (k < q)
+            {
+                u = (bre[k] + n) * (bre[k] + n) + (bim[k] * bim[k]);
+                u = FLINT_ABS(u);
+
+                if (u < 1e-8 || u > 1e100 || t > 1e100)
+                    goto somethingstrange;
+
+                t /= u;
+            }
+        }
+
+        log2_factor = 0.5 * log(t) * 1.4426950408889634074 + log2_z;
+
+        /* For asymptotic series, require rapid decay */
+        if (p > q && n >= minimum_n && log2_factor > -0.2)
+            break;
+
+        log2_term += log2_factor;
+        log2_term_max = FLINT_MAX(log2_term_max, log2_term);
+        n++;
+    }
+
+somethingstrange:
+
+    flint_free(are);
+    mag_clear(zmag);
+
+    return n;
+}
+

acb_hypgeom/pfq_direct.c

@@ -37,7 +37,10 @@ acb_hypgeom_pfq_direct(acb_t res, acb_srcptr a, long p, acb_srcptr b, long q,
     mag_init(err);
     mag_init(C);
 
-    acb_hypgeom_pfq_sum(s, t, a, p, b, q, z, n, prec);
+    if (n < 0)
+        n = acb_hypgeom_pfq_choose_n(a, p, b, q, z, prec);
+
+    acb_hypgeom_pfq_sum_rs(s, t, a, p, b, q, z, n, prec);
 
     if (!acb_is_zero(t))
     {

acb_hypgeom/pfq_sum.c

@@ -26,52 +26,17 @@
 #include "acb_hypgeom.h"
 
 void
-acb_hypgeom_pfq_sum(acb_t s, acb_t t,
+acb_hypgeom_pfq_sum(acb_t s, acb_t t, acb_srcptr a, long p,
-    acb_srcptr a, long p, acb_srcptr b, long q, const acb_t z, long n, long prec)
+    acb_srcptr b, long q, const acb_t z, long n, long prec)
 {
-    acb_t u, v;
+    if (n > 4 && prec >= 128
-    long k, i;
+        && _acb_vec_bits(a, p) * p + _acb_vec_bits(b, q) * q + 10 < prec / 2)
-
-    acb_init(u);
-    acb_init(v);
-
-    acb_zero(s);
-    acb_one(t);
-
-    for (k = 0; k < n && !acb_is_zero(t); k++)
     {
-        acb_add(s, s, t, prec);
+        acb_hypgeom_pfq_sum_rs(s, t, a, p, b, q, z, n, prec);
-
+    }
-        if (p > 0)
+    else
     {
-            acb_add_ui(u, a, k, prec);
+        acb_hypgeom_pfq_sum_forward(s, t, a, p, b, q, z, n, prec);
-
+    }
-            for (i = 1; i < p; i++)
-            {
-                acb_add_ui(v, a + i, k, prec);
-                acb_mul(u, u, v, prec);
-            }
-
-            acb_mul(t, t, u, prec);
-        }
-
-        if (q > 0)
-        {
-            acb_add_ui(u, b, k, prec);
-
-            for (i = 1; i < q; i++)
-            {
-                acb_add_ui(v, b + i, k, prec);
-                acb_mul(u, u, v, prec);
-            }
-
-            acb_div(t, t, u, prec);
-        }
-
-        acb_mul(t, t, z, prec);
-    }
-
-    acb_clear(u);
-    acb_clear(v);
 }
 

acb_hypgeom/pfq_sum_forward.c

@@ -0,0 +1,77 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2014 Fredrik Johansson
+
+******************************************************************************/
+
+#include "acb_hypgeom.h"
+
+void
+acb_hypgeom_pfq_sum_forward(acb_t s, acb_t t,
+    acb_srcptr a, long p, acb_srcptr b, long q, const acb_t z, long n, long prec)
+{
+    acb_t u, v;
+    long k, i;
+
+    acb_init(u);
+    acb_init(v);
+
+    acb_zero(s);
+    acb_one(t);
+
+    for (k = 0; k < n && !acb_is_zero(t); k++)
+    {
+        acb_add(s, s, t, prec);
+
+        if (p > 0)
+        {
+            acb_add_ui(u, a, k, prec);
+
+            for (i = 1; i < p; i++)
+            {
+                acb_add_ui(v, a + i, k, prec);
+                acb_mul(u, u, v, prec);
+            }
+
+            acb_mul(t, t, u, prec);
+        }
+
+        if (q > 0)
+        {
+            acb_add_ui(u, b, k, prec);
+
+            for (i = 1; i < q; i++)
+            {
+                acb_add_ui(v, b + i, k, prec);
+                acb_mul(u, u, v, prec);
+            }
+
+            acb_div(t, t, u, prec);
+        }
+
+        acb_mul(t, t, z, prec);
+    }
+
+    acb_clear(u);
+    acb_clear(v);
+}
+

acb_hypgeom/pfq_sum_rs.c

@@ -0,0 +1,130 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2014 Fredrik Johansson
+
+******************************************************************************/
+
+#include "acb_hypgeom.h"
+
+void
+acb_hypgeom_pfq_sum_rs(acb_t res, acb_t term, acb_srcptr a, long p,
+                                              acb_srcptr b, long q, const acb_t z, long n, long prec)
+{
+    acb_ptr zpow;
+    acb_t s, t, u;
+    long i, j, k, m;
+    mag_t B, C;
+
+    if (n == 0)
+    {
+        acb_zero(res);
+        acb_one(term);
+        return;
+    }
+
+    if (n < 0)
+        abort();
+
+    m = n_sqrt(n);
+    m = FLINT_MIN(m, 150);
+
+    mag_init(B);
+    mag_init(C);
+    acb_init(s);
+    acb_init(t);
+    acb_init(u);
+    zpow = _acb_vec_init(m + 1);
+
+    _acb_vec_set_powers(zpow, z, m + 1, prec);
+
+    mag_one(B);
+
+    for (k = n; k >= 0; k--)
+    {
+        j = k % m;
+
+        if (k < n)
+            acb_add(s, s, zpow + j, prec);
+
+        if (k > 0)
+        {
+            if (p > 0)
+            {
+                acb_add_ui(u, a, k - 1, prec);
+
+                for (i = 1; i < p; i++)
+                {
+                    acb_add_ui(t, a + i, k - 1, prec);
+                    acb_mul(u, u, t, prec);
+                }
+
+                if (k < n)
+                    acb_mul(s, s, u, prec);
+
+                acb_get_mag(C, u);
+                mag_mul(B, B, C);
+            }
+
+            if (q > 0)
+            {
+                acb_add_ui(u, b, k - 1, prec);
+
+                for (i = 1; i < q; i++)
+                {
+                    acb_add_ui(t, b + i, k - 1, prec);
+                    acb_mul(u, u, t, prec);
+                }
+
+                if (k < n)
+                    acb_div(s, s, u, prec);
+
+                acb_get_mag_lower(C, u);
+                mag_div(B, B, C);
+            }
+
+            if (j == 0 && k < n)
+            {
+                acb_mul(s, s, zpow + m, prec);
+            }
+        }
+    }
+
+    acb_get_mag(C, z);
+    mag_pow_ui(C, C, n);
+    mag_mul(B, B, C);
+
+    acb_zero(term);
+    if (_acb_vec_is_real(a, p) && _acb_vec_is_real(b, q) && acb_is_real(z))
+        arb_add_error_mag(acb_realref(term), B);
+    else
+        acb_add_error_mag(term, B);
+
+    acb_set(res, s);
+
+    mag_clear(B);
+    mag_clear(C);
+    acb_clear(s);
+    acb_clear(t);
+    acb_clear(u);
+    _acb_vec_clear(zpow, m + 1);
+}
+

acb_hypgeom/u_asymp.c

@@ -178,6 +178,8 @@ void acb_hypgeom_u_asymp(acb_t res, const acb_t a, const acb_t b,
     const acb_t z, long n, long prec)
 {
     mag_t C1, Cn, alpha, nu, sigma, rho, zinv, tmp, err;
+    acb_struct aa[3];
+    acb_t s, t, w;
     int R;
 
     if (!acb_is_finite(a) || !acb_is_finite(b) || !acb_is_finite(z))
@@ -196,6 +198,24 @@ void acb_hypgeom_u_asymp(acb_t res, const acb_t a, const acb_t b,
     mag_init(tmp);
     mag_init(err);
 
+    acb_init(aa);
+    acb_init(aa + 1);
+    acb_init(aa + 2);
+    acb_init(s);
+    acb_init(t);
+    acb_init(w);
+
+    acb_set(aa, a);
+    acb_sub(aa + 1, a, b, prec);
+    acb_add_ui(aa + 1, aa + 1, 1, prec);
+    acb_one(aa + 2);
+
+    acb_neg(w, z);
+    acb_inv(w, w, prec);
+
+    if (n < 0)
+        n = acb_hypgeom_pfq_choose_n(aa, 2, aa + 2, 1, w, prec);
+
     acb_hypgeom_u_asymp_bound_factors(&R, alpha, nu,
         sigma, rho, zinv, a, b, z);
 
@@ -205,6 +225,8 @@ void acb_hypgeom_u_asymp(acb_t res, const acb_t a, const acb_t b,
     }
     else
     {
+        acb_hypgeom_pfq_sum(s, t, aa, 2, aa + 2, 1, w, n, prec);
+
         if (R == 1)
         {
             mag_one(C1);
@@ -239,54 +261,12 @@ void acb_hypgeom_u_asymp(acb_t res, const acb_t a, const acb_t b,
         mag_mul(err, err, alpha);
         mag_mul_2exp_si(err, err, 1);
 
-        /* evaluate the sum, naively for now */
+        /* nth term * factor */
-        {
-            acb_t s, t, u, v, w, ab1;
-            long k;
-
-            acb_init(s);
-            acb_init(t);
-            acb_init(u);
-            acb_init(v);
-            acb_init(w);
-            acb_init(ab1);
-
-            acb_one(t);
-            acb_sub(ab1, a, b, prec);
-            acb_add_ui(ab1, ab1, 1, prec);
-
-            if (n != 0)
-            {
-                acb_neg(w, z);
-                acb_inv(w, w, prec);
-            }
-
-            for (k = 0; k < n && !acb_is_zero(t); k++)
-            {
-                acb_add(s, s, t, prec);
-
-                acb_add_ui(u, a, k, prec);
-                acb_add_ui(v, ab1, k, prec);
-                acb_mul(u, u, v, prec);
-                acb_mul(t, t, u, prec);
-                acb_mul(t, t, w, prec);
-                acb_div_ui(t, t, k + 1, prec);
-            }
-
-            /* nth term */
         acb_get_mag(tmp, t);
         mag_mul(err, err, tmp);
         acb_add_error_mag(s, err);
 
         acb_set(res, s);
-
-            acb_clear(s);
-            acb_clear(t);
-            acb_clear(u);
-            acb_clear(v);
-            acb_clear(w);
-            acb_clear(ab1);
-        }
     }
 
     mag_clear(C1);
@@ -298,5 +278,12 @@ void acb_hypgeom_u_asymp(acb_t res, const acb_t a, const acb_t b,
     mag_clear(zinv);
     mag_clear(tmp);
     mag_clear(err);
+
+    acb_clear(aa);
+    acb_clear(aa + 1);
+    acb_clear(aa + 2);
+    acb_clear(s);
+    acb_clear(t);
+    acb_clear(w);
 }
 

doc/source/acb_hypgeom.rst

@@ -66,13 +66,42 @@ or remove a 1 from the `a_i` parameters if there is one.
     As currently implemented, the bound becomes infinite when `n` is
     too small, even if the series converges.
 
+.. function:: long acb_hypgeom_pfq_choose_n(acb_srcptr a, long p, acb_srcptr b, long q, const acb_t z, long prec)
+
+    Heuristically attempts to choose a number of terms *n* to
+    sum of a hypergeometric series at a working precision of *prec* bits.
+
+    Uses double precision arithmetic internally. As currently implemented,
+    it can fail to produce a good result if the parameters are extremely
+    large or extremely close to nonpositive integers.
+
+    Numerical cancellation is assumed to be significant, so truncation
+    is done when the current term is *prec* bits
+    smaller than the largest encountered term.
+
+    This function will also attempt to pick a reasonable
+    truncation point for divergent series.
+
+.. function:: void acb_hypgeom_pfq_sum_forward(acb_t s, acb_t t, acb_srcptr a, long p, acb_srcptr b, long q, const acb_t z, long n, long prec)
+
+.. function:: void acb_hypgeom_pfq_sum_rs(acb_t s, acb_t t, acb_srcptr a, long p, acb_srcptr b, long q, const acb_t z, long n, long prec)
+
 .. function:: void acb_hypgeom_pfq_sum(acb_t s, acb_t t, acb_srcptr a, long p, acb_srcptr b, long q, const acb_t z, long n, long prec)
 
     Computes `s = \sum_{k=0}^{n-1} T(k)` and `t = T(n)`.
-
     Does not allow aliasing between input and output variables.
     We require `n \ge 0`.
 
+    The *forward* version computes the sum using forward
+    recurrence.
+
+    The *rs* version computes the sum in reverse order
+    using rectangular splitting. It only computes a
+    magnitude bound for the value of *t*.
+
+    The default version automatically chooses an algorithm
+    depending on the inputs.
+
 .. function:: void acb_hypgeom_pfq_direct(acb_t res, acb_srcptr a, long p, acb_srcptr b, long q, const acb_t z, long n, long prec)
 
     Computes
@@ -85,7 +114,9 @@ or remove a 1 from the `a_i` parameters if there is one.
 
     directly from the defining series, including a rigorous bound for
     the truncation error `\varepsilon` in the output.
-    We require `n \ge 0`.
+
+    If  `n < 0`, this function chooses a number of terms automatically
+    using :func:`acb_hypgeom_pfq_choose_n`.
 
 Asymptotic series
 -------------------------------------------------------------------------------