rising factorial algorithms extravaganza

arb.h

@@ -652,6 +652,18 @@ _arb_vec_entry_ptr(arb_ptr vec, slong i)
     return vec + i;
 }
 
+ARB_INLINE void
+_arb_vec_printn(arb_srcptr vec, slong len, slong ndigits, ulong flags)
+{
+    slong i;
+    for (i = 0; i < len; i++)
+    {
+        arb_printn(vec + i, ndigits, flags);
+        if (i < len - 1)
+            flint_printf(", ");
+    }
+}
+
 ARB_INLINE void
 _arb_vec_zero(arb_ptr A, slong n)
 {

arb_hypgeom.h

@@ -23,7 +23,17 @@ void _arb_hypgeom_rising_coeffs_1(ulong * c, ulong k, slong l);
 void _arb_hypgeom_rising_coeffs_2(ulong * c, ulong k, slong l);
 void _arb_hypgeom_rising_coeffs_fmpz(fmpz * c, ulong k, slong l);
 
+void arb_hypgeom_rising_ui_forward(arb_t res, const arb_t x, ulong n, slong prec);
 void arb_hypgeom_rising_ui_rs(arb_t res, const arb_t x, ulong n, ulong m, slong prec);
+void arb_hypgeom_rising_ui_bs(arb_t res, const arb_t x, ulong n, slong prec);
+void arb_hypgeom_rising_ui_rec(arb_t res, const arb_t x, ulong n, slong prec);
+void arb_hypgeom_rising_ui(arb_t y, const arb_t x, ulong n, slong prec);
+void arb_hypgeom_rising(arb_t y, const arb_t x, const arb_t n, slong prec);
+
+void arb_hypgeom_rising_ui_jet_powsum(arb_ptr res, const arb_t x, ulong n, slong len, slong prec);
+void arb_hypgeom_rising_ui_jet_rs(arb_ptr res, const arb_t x, ulong n, ulong m, slong len, slong prec);
+void arb_hypgeom_rising_ui_jet_bs(arb_ptr res, const arb_t x, ulong n, slong len, slong prec);
+void arb_hypgeom_rising_ui_jet(arb_ptr res, const arb_t x, ulong n, slong len, slong prec);
 
 void arb_hypgeom_pfq(arb_t res, arb_srcptr a, slong p, arb_srcptr b, slong q,
     const arb_t z, int regularized, slong prec);

arb_hypgeom/rising_ui.c

@@ -0,0 +1,53 @@
+/*
+    Copyright (C) 2014 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_rising_ui(arb_t y, const arb_t x, ulong n, slong prec)
+{
+    if (n < FLINT_MAX(prec, 100))
+    {
+        arb_hypgeom_rising_ui_rec(y, x, n, prec);
+    }
+    else
+    {
+        arb_t t;
+        arb_init(t);
+        arb_add_ui(t, x, n, prec);
+        arb_gamma(t, t, prec);
+        arb_rgamma(y, x, prec);
+        arb_mul(y, y, t, prec);
+        arb_clear(t);
+    }
+}
+
+void
+arb_hypgeom_rising(arb_t y, const arb_t x, const arb_t n, slong prec)
+{
+    if (arb_is_int(n) && arf_sgn(arb_midref(n)) >= 0 &&
+        arf_cmpabs_ui(arb_midref(n), FLINT_MAX(prec, 100)) < 0)
+    {
+        arb_hypgeom_rising_ui_rec(y, x,
+            arf_get_si(arb_midref(n), ARF_RND_DOWN), prec);
+    }
+    else
+    {
+        arb_t t;
+        arb_init(t);
+        arb_add(t, x, n, prec);
+        arb_gamma(t, t, prec);
+        arb_rgamma(y, x, prec);
+        arb_mul(y, y, t, prec);
+        arb_clear(t);
+    }
+}
+

arb_hypgeom/rising_ui_bs.c

@@ -0,0 +1,69 @@
+/*
+    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"
+
+static void
+bsplit(arb_t y, const arb_t x, ulong a, ulong b, slong prec)
+{
+    if (b - a <= 16)
+    {
+        if (a == 0)
+        {
+            arb_hypgeom_rising_ui_forward(y, x, b, prec);
+        }
+        else
+        {
+            arb_add_ui(y, x, a, prec);
+            arb_hypgeom_rising_ui_forward(y, y, b - a, prec);
+        }
+    }
+    else
+    {
+        arb_t t, u;
+        ulong m = a + (b - a) / 2;
+
+        arb_init(t);
+        arb_init(u);
+
+        bsplit(t, x, a, m, prec);
+        bsplit(u, x, m, b, prec);
+
+        arb_mul(y, t, u, prec);
+
+        arb_clear(t);
+        arb_clear(u);
+    }
+}
+
+void
+arb_hypgeom_rising_ui_bs(arb_t res, const arb_t x, ulong n, slong prec)
+{
+    if (n <= 1)
+    {
+        if (n == 0)
+            arb_one(res);
+        else
+            arb_set_round(res, x, prec);
+        return;
+    }
+
+    {
+        arb_t t;
+        slong wp = ARF_PREC_ADD(prec, FLINT_BIT_COUNT(n));
+
+        arb_init(t);
+        bsplit(t, x, 0, n, wp);
+        arb_set_round(res, t, prec);
+        arb_clear(t);
+    }
+}
+

arb_hypgeom/rising_ui_forward.c

@@ -0,0 +1,83 @@
+/*
+    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
+_arf_increment_fast(arf_t x, slong prec)
+{
+    if (arf_sgn(x) > 0)
+    {
+        mp_limb_t hi, v, cy;
+        mp_ptr xptr;
+        mp_size_t xn;
+        slong xexp;
+
+        xexp = ARF_EXP(x);
+
+        if (xexp >= 1 && xexp <= FLINT_BITS - 1)
+        {
+            ARF_GET_MPN_READONLY(xptr, xn, x);
+
+            hi = xptr[xn - 1];
+            v = hi + (UWORD(1) << (FLINT_BITS - xexp));
+            cy = v < hi;
+
+            if (cy == 0)
+            {
+                xptr[xn - 1] = v;
+                return 0;
+            }
+        }
+    }
+
+    return arf_add_ui(x, x, 1, prec, ARF_RND_DOWN);
+}
+
+void
+_arb_increment_fast(arb_t x, slong prec)
+{
+    if (_arf_increment_fast(arb_midref(x), prec))
+        arf_mag_add_ulp(arb_radref(x), arb_radref(x), arb_midref(x), prec);
+}
+
+void
+arb_hypgeom_rising_ui_forward(arb_t res, const arb_t x, ulong n, slong prec)
+{
+    arb_t t;
+    ulong k;
+    slong wp;
+
+    if (n <= 1)
+    {
+        if (n == 0)
+            arb_one(res);
+        else
+            arb_set_round(res, x, prec);
+        return;
+    }
+
+    wp = prec + FLINT_BIT_COUNT(n);
+
+    arb_init(t);
+
+    arb_add_ui(t, x, 1, wp);
+    arb_mul(res, x, t, (n == 2) ? prec : wp);
+
+    for (k = 2; k < n; k++)
+    {
+        _arb_increment_fast(t, wp);
+        arb_mul(res, res, t, k == (n - 1) ? prec : wp);
+    }
+
+    arb_clear(t);
+}
+

arb_hypgeom/rising_ui_jet.c

@@ -0,0 +1,44 @@
+/*
+    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_rising_ui_jet(arb_ptr res, const arb_t x, ulong n, slong len, slong prec)
+{
+    if (n <= 7)
+    {
+        arb_hypgeom_rising_ui_jet_powsum(res, x, n, len, prec);
+    }
+    else if (len == 2)
+    {
+        if (n <= 30 || arb_bits(x) >= prec / 128)
+            arb_hypgeom_rising_ui_jet_rs(res, x, n, 0, len, prec);
+        else
+            arb_hypgeom_rising_ui_jet_bs(res, x, n, len, prec);
+    }
+    else
+    {
+        if (n <= 20 || (n <= 200 && prec > 400 * n && arb_bits(x) >= prec / 4))
+        {
+            arb_hypgeom_rising_ui_jet_powsum(res, x, n, len, prec);
+        }
+        else if (len >= 64 || (arb_bits(x) + 1 < prec / 1024 && n >= 32))
+        {
+            arb_hypgeom_rising_ui_jet_bs(res, x, n, len, prec);
+        }
+        else
+        {
+            arb_hypgeom_rising_ui_jet_rs(res, x, n, 0, len, prec);
+        }
+    }
+}
+

arb_hypgeom/rising_ui_jet_bs.c

@@ -0,0 +1,84 @@
+/*
+    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"
+
+static void
+bsplit(arb_ptr res, const arb_t x, ulong a, ulong b, slong trunc, slong prec)
+{
+    trunc = FLINT_MIN(trunc, b - a + 1);
+
+    if (b - a <= 12)
+    {
+        if (a == 0)
+        {
+            arb_hypgeom_rising_ui_jet_powsum(res, x, b - a, FLINT_MIN(trunc, b - a + 1), prec);
+        }
+        else
+        {
+            arb_t t;
+            arb_init(t);
+            arb_add_ui(t, x, a, prec);
+            arb_hypgeom_rising_ui_jet_powsum(res, t, b - a, FLINT_MIN(trunc, b - a + 1), prec);
+            arb_clear(t);
+        }
+    }
+    else
+    {
+        arb_ptr L, R;
+        slong len1, len2;
+
+        slong m = a + (b - a) / 2;
+
+        len1 = poly_pow_length(2, m - a, trunc);
+        len2 = poly_pow_length(2, b - m, trunc);
+
+        L = _arb_vec_init(len1 + len2);
+        R = L + len1;
+
+        bsplit(L, x, a, m, trunc, prec);
+        bsplit(R, x, m, b, trunc, prec);
+
+        _arb_poly_mullow(res, L, len1, R, len2,
+            FLINT_MIN(trunc, len1 + len2 - 1), prec);
+
+        _arb_vec_clear(L, len1 + len2);
+    }
+}
+
+void
+arb_hypgeom_rising_ui_jet_bs(arb_ptr res, const arb_t x, ulong n, slong len, slong prec)
+{
+    if (len == 0)
+        return;
+
+    if (len > n + 1)
+    {
+        _arb_vec_zero(res + n + 1, len - n - 1);
+        len = n + 1;
+    }
+
+    if (len == n + 1)
+    {
+        arb_one(res + n);
+        len = n;
+    }
+
+    if (n <= 1)
+    {
+        if (n == 1)
+            arb_set_round(res, x, prec);
+        return;
+    }
+
+    bsplit(res, x, 0, n, len, prec);
+}
+

arb_hypgeom/rising_ui_jet_powsum.c

@@ -0,0 +1,159 @@
+/*
+    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_rising_ui_jet_powsum(arb_ptr res, const arb_t x, ulong n, slong len, slong prec)
+{
+    slong i, j, k, wp;
+    arb_ptr xpow;
+    TMP_INIT;
+
+    if (len == 0)
+        return;
+
+    if (len > n + 1)
+    {
+        _arb_vec_zero(res + n + 1, len - n - 1);
+        len = n + 1;
+    }
+
+    if (len == n + 1)
+    {
+        arb_one(res + n);
+        len = n;
+    }
+
+    if (n <= 1)
+    {
+        if (n == 1)
+            arb_set_round(res, x, prec);
+        return;
+    }
+
+    if (len == 1)
+    {
+        arb_hypgeom_rising_ui_rs(res, x, n, 0, prec);
+        return;
+    }
+
+    if (n == 2)
+    {
+        arb_mul_2exp_si(res + 1, x, 1);
+        arb_add_ui(res + 1, res + 1, 1, prec);
+        arb_mul(res, x, x, prec + 4);
+        arb_add(res, res, x, prec);
+        return;
+    }
+
+    if (n <= 12 || (FLINT_BITS == 64 && n <= 20))
+    {
+        mp_ptr c;
+        TMP_START;
+
+        wp = ARF_PREC_ADD(prec, FLINT_BIT_COUNT(n));
+        c = TMP_ALLOC(sizeof(mp_limb_t) * (n + 1) * len);
+
+        _nmod_vec_zero(c, (n + 1) * len);
+
+        c[0] = 0;
+        c[1] = 1;
+        c[(n + 1) + 0] = 1;
+
+        for (i = 2; i <= n; i++)
+        {
+            for (j = FLINT_MIN(len - 1, i); j >= 0; j--)
+            {
+                slong ln, pos;
+
+                ln = i + 1 - j;
+                pos = (n + 1) * j;
+                if (i == j)
+                {
+                    c[pos] = 1;
+                }
+                else
+                {
+                    c[pos + ln - 1] = c[pos + ln - 2];
+                    for (k = ln - 2; k >= 1; k--)
+                        c[pos + k] = c[pos + k] * (i - 1) + c[pos + k - 1];
+                    c[pos + 0] *= (i - 1);
+                    if (j != 0)
+                        for (k = ln - 1; k >= 0; k--)
+                            c[pos + k] += c[pos - (n + 1) + k];
+                }
+            }
+        }
+
+        xpow = _arb_vec_init(n + 1);
+        _arb_vec_set_powers(xpow, x, n + 1, wp);
+
+        arb_dot_ui(res, NULL, 0, xpow + 1, 1, c + 1, 1, n, prec);
+
+        for (i = 1; i < len; i++)
+            arb_dot_ui(res + i, NULL, 0, xpow, 1, c + (n + 1) * i, 1, n + 1 - i, prec);
+
+        _arb_vec_clear(xpow, n + 1);
+        TMP_END;
+    }
+    else
+    {
+        fmpz * c;
+
+        wp = ARF_PREC_ADD(prec, FLINT_BIT_COUNT(n));
+        c = _fmpz_vec_init((n + 1) * len);
+
+        fmpz_one(c + 1);
+        fmpz_one(c + (n + 1) + 0);
+
+        for (i = 2; i <= n; i++)
+        {
+            for (j = FLINT_MIN(len - 1, i); j >= 0; j--)
+            {
+                slong ln, pos;
+
+                ln = i + 1 - j;
+                pos = (n + 1) * j;
+                if (i == j)
+                {
+                    fmpz_one(c + pos);
+                }
+                else
+                {
+                    fmpz_set(c + pos + ln - 1, c + pos + ln - 2);
+                    for (k = ln - 2; k >= 1; k--)
+                    {
+                        fmpz_mul_ui(c + pos + k, c + pos + k, i - 1);
+                        fmpz_add(c + pos + k, c + pos + k, c + pos + k - 1);
+                    }
+
+                    fmpz_mul_ui(c + pos + 0, c + pos + 0, i - 1);
+                    if (j != 0)
+                        for (k = ln - 1; k >= 0; k--)
+                            fmpz_add(c + pos + k, c + pos + k, c + pos - (n + 1) + k);
+                }
+            }
+        }
+
+        xpow = _arb_vec_init(n + 1);
+        _arb_vec_set_powers(xpow, x, n + 1, wp);
+
+        arb_dot_fmpz(res, NULL, 0, xpow + 1, 1, c + 1, 1, n, prec);
+
+        for (i = 1; i < len; i++)
+            arb_dot_fmpz(res + i, NULL, 0, xpow, 1, c + (n + 1) * i, 1, n + 1 - i, prec);
+
+        _arb_vec_clear(xpow, n + 1);
+        _fmpz_vec_clear(c, (n + 1) * len);
+    }
+}
+

arb_hypgeom/rising_ui_jet_rs.c

@@ -0,0 +1,173 @@
+/*
+    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_rising_ui_jet_rs(arb_ptr res, const arb_t x, ulong n, ulong m, slong len, slong prec)
+{
+    slong i, j, k, l, m0, xmlen, tlen, ulen, climbs, climbs_max, wp;
+    arb_ptr tmp, xpow;
+    arb_ptr t, u;
+    mp_ptr c;
+    TMP_INIT;
+
+    if (len == 0)
+        return;
+
+    if (len > n + 1)
+    {
+        _arb_vec_zero(res + n + 1, len - n - 1);
+        len = n + 1;
+    }
+
+    if (len == n + 1)
+    {
+        arb_one(res + n);
+        len = n;
+    }
+
+    if (n <= 1)
+    {
+        if (n == 1)
+            arb_set_round(res, x, prec);
+        return;
+    }
+
+    if (len == 1)
+    {
+        arb_hypgeom_rising_ui_rs(res, x, n, m, prec);
+        return;
+    }
+
+    TMP_START;
+
+    if (m == 0)
+    {
+        if (n <= 7)
+            m = n;
+        else if (n <= 12)
+            m = (n + 1) / 2;
+        else if (n <= 32)
+            m = (n + 2) / 3;
+        else
+        {
+            m0 = n_sqrt(n);
+            m = 8 + 0.5 * pow(prec, 0.4);
+            m = FLINT_MIN(m, m0);
+            m = FLINT_MIN(m, 64);
+        }
+    }
+
+    wp = ARF_PREC_ADD(prec, FLINT_BIT_COUNT(n));
+
+    climbs_max = FLINT_BIT_COUNT(n - 1) * m;
+    c = TMP_ALLOC(sizeof(mp_limb_t) * climbs_max * m);
+
+    /* length of (x+t)^m */
+    xmlen = FLINT_MIN(len, m + 1);
+
+    tmp = _arb_vec_init(2 * len + (m + 1) * xmlen);
+    t = tmp;
+    u = tmp + len;
+    xpow = tmp + 2 * len;
+
+    _arb_vec_set_powers(xpow, x, m + 1, wp);
+
+    tlen = 1;
+
+    /* First derivatives */
+    for (i = 1; i <= m; i++)
+        arb_mul_ui(xpow + (m + 1) + i, xpow + i - 1, i, wp);
+
+    /* Higher derivatives if we need them */
+    if (len >= 3)
+    {
+        fmpz * f = _fmpz_vec_init(len);
+
+        fmpz_one(f + 0);
+        fmpz_one(f + 1);
+
+        for (i = 2; i <= m; i++)
+        {
+            for (j = FLINT_MIN(xmlen - 1, i + 1); j >= 1; j--)
+                fmpz_add(f + j, f + j, f + j - 1);
+
+            for (j = 2; j < FLINT_MIN(xmlen, i + 1); j++)
+                arb_mul_fmpz(xpow + (m + 1) * j + i, xpow + i - j, f + j, wp);
+        }
+
+        _fmpz_vec_clear(f, len);
+    }
+
+    for (k = 0; k < n; k += m)
+    {
+        l = FLINT_MIN(m, n - k);
+        climbs = FLINT_BIT_COUNT(k + l - 1) * l;
+        climbs = (climbs + FLINT_BITS - 1) / FLINT_BITS;
+
+        ulen = FLINT_MIN(len, l + 1);
+
+        if (l == 1)
+        {
+            arb_add_ui(u, x, k, wp);
+            arb_one(u + 1);
+        }
+        else
+        {
+            if (climbs == 1)
+            {
+                _arb_hypgeom_rising_coeffs_1(c, k, l);
+                for (j = 0; j < ulen; j++)
+                    arb_dot_ui(u + j, xpow + (m + 1) * j + l, 0, xpow + (m + 1) * j + j, 1, c + j, 1, l - j, wp);
+            }
+            else if (climbs == 2)
+            {
+                _arb_hypgeom_rising_coeffs_2(c, k, l);
+                for (j = 0; j < ulen; j++)
+                    arb_dot_uiui(u + j, xpow + (m + 1) * j + l, 0, xpow + (m + 1) * j + j, 1, c + 2 * j, 1, l - j, wp);
+            }
+            else
+            {
+                fmpz * f = (fmpz *) c;
+
+                for (i = 0; i < l; i++)
+                    fmpz_init(f + i);
+
+                _arb_hypgeom_rising_coeffs_fmpz(f, k, l);
+
+                for (j = 0; j < ulen; j++)
+                    arb_dot_fmpz(u + j, xpow + (m + 1) * j + l, 0, xpow + (m + 1) * j + j, 1, f + j, 1, l - j, wp);
+
+                for (i = 0; i < l; i++)
+                    fmpz_clear(f + i);
+            }
+        }
+
+        if (k == 0)
+        {
+            tlen = ulen;
+            _arb_vec_swap(t, u, ulen);
+        }
+        else
+        {
+            _arb_poly_mullow(res, t, tlen, u, ulen, FLINT_MIN(len, tlen + ulen - 1), wp);
+            tlen = FLINT_MIN(len, tlen + ulen - 1);
+            _arb_vec_swap(t, res, tlen);
+        }
+    }
+
+    _arb_vec_set_round(res, t, len, prec);
+
+    _arb_vec_clear(tmp, 2 * len + (m + 1) * xmlen);
+    TMP_END;
+}
+

arb_hypgeom/rising_ui_rec.c

@@ -0,0 +1,46 @@
+/*
+    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_rising_ui_rec(arb_t res, const arb_t x, ulong n, slong prec)
+{
+    if (n <= 1)
+    {
+        if (n == 0)
+            arb_one(res);
+        else
+            arb_set_round(res, x, prec);
+        return;
+    }
+
+    if (n == 2 && prec <= 1024)
+    {
+        if (res != x)
+            arb_set(res, x);
+        arb_addmul(res, x, x, prec);
+        return;
+    }
+
+    if ((prec < 512 && n <= 20) || (n <= FLINT_MIN(80, 6000 / prec)))
+    {
+        arb_hypgeom_rising_ui_forward(res, x, n, prec);
+    }
+    else
+    {
+        if (n >= 20 && arb_bits(x) < prec / 8)
+            arb_hypgeom_rising_ui_bs(res, x, n, prec);
+        else
+            arb_hypgeom_rising_ui_rs(res, x, n, 0, prec);
+    }
+}
+

arb_hypgeom/test/t-rising_ui.c

@@ -1,5 +1,5 @@
 /*
-    Copyright (C) 2018 Fredrik Johansson
+    Copyright (C) 2021 Fredrik Johansson
 
     This file is part of Arb.
 
@@ -11,21 +11,44 @@
 
 #include "arb_hypgeom.h"
 
+void
+rising_algorithm(arb_t res, const arb_t x, ulong n, ulong m, slong prec, int alg, int alias)
+{
+    if (alias)
+    {
+        arb_set(res, x);
+        rising_algorithm(res, res, n, m, prec, alg, 0);
+        return;
+    }
+
+    if (alg == 0)
+        arb_hypgeom_rising_ui_rs(res, x, n, m, prec);
+    else if (alg == 1)
+        arb_hypgeom_rising_ui_forward(res, x, n, prec);
+    else if (alg == 2)
+        arb_hypgeom_rising_ui_bs(res, x, n, prec);
+    else if (alg == 3)
+        arb_hypgeom_rising_ui_rec(res, x, n, prec);
+    else
+        arb_hypgeom_rising_ui(res, x, n, prec);
+}
+
 int main()
 {
     slong iter;
     flint_rand_t state;
 
-    flint_printf("rising_ui_rs....");
+    flint_printf("rising_ui....");
     fflush(stdout);
 
     flint_randinit(state);
 
-    for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
+    for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++)
     {
         arb_t x, xk, y, ya, yb, yayb;
         ulong k, n, m1, m2, m3;
         slong prec;
+        int alg1, alg2, alg3, alias1, alias2, alias3;
 
         prec = 2 + n_randint(state, 200);
         k = n_randint(state, 10);
@@ -33,6 +56,15 @@ int main()
         m1 = n_randint(state, 2) ? 0 : 1 + n_randint(state, FLINT_MAX(n + k, 1));
         m2 = n_randint(state, 2) ? 0 : 1 + n_randint(state, FLINT_MAX(k, 1));
         m3 = n_randint(state, 2) ? 0 : 1 + n_randint(state, FLINT_MAX(n, 1));
+        alg1 = n_randint(state, 5);
+        alg2 = n_randint(state, 5);
+        alg3 = n_randint(state, 5);
+        alias1 = n_randint(state, 2);
+        alias2 = n_randint(state, 2);
+        alias3 = n_randint(state, 2);
+
+        if (n_randint(state, 100) == 0)
+            n += 100;
 
         arb_init(x);
         arb_init(xk);
@@ -41,12 +73,12 @@ int main()
         arb_init(yb);
         arb_init(yayb);
 
-        arb_randtest(x, state, prec, 10);
+        arb_randtest(x, state, prec, 10 + n_randint(state, 200));
         arb_add_ui(xk, x, k, prec);
 
-        arb_hypgeom_rising_ui_rs(y, x, n + k, m1, prec);
-        arb_hypgeom_rising_ui_rs(ya, x, k, m2, prec);
-        arb_hypgeom_rising_ui_rs(yb, xk, n, m3, prec);
+        rising_algorithm(y, x, n + k, m1, prec, alg1, alias1);
+        rising_algorithm(ya, x, k, m2, prec, alg2, alias2);
+        rising_algorithm(yb, xk, n, m3, prec, alg3, alias3);
         arb_mul(yayb, ya, yb, prec);
 
         if (!arb_overlaps(y, yayb))

arb_hypgeom/test/t-rising_ui_jet.c

@@ -0,0 +1,95 @@
+/*
+    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
+rising_algorithm(arb_ptr res, const arb_t x, ulong n, ulong m, slong len, slong prec, int alg)
+{
+    if (alg == 0)
+        arb_hypgeom_rising_ui_jet_powsum(res, x, n, len, prec);
+    else if (alg == 1)
+        arb_hypgeom_rising_ui_jet_rs(res, x, n, m, len, prec);
+    else if (alg == 2)
+        arb_hypgeom_rising_ui_jet_bs(res, x, n, len, prec);
+    else
+        arb_hypgeom_rising_ui_jet(res, x, n, len, prec);
+}
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("rising_ui_jet....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++)
+    {
+        arb_t x, xk;
+        arb_ptr y, ya, yb, yayb;
+        ulong k, n, m1, m2, m3, len;
+        slong prec;
+        int alg1, alg2, alg3;
+
+        prec = 2 + n_randint(state, 200);
+        len = 1 + n_randint(state, 6);
+        k = n_randint(state, 10);
+        n = n_randint(state, 50);
+        m1 = n_randint(state, 2) ? 0 : 1 + n_randint(state, FLINT_MAX(n + k, 1));
+        m2 = n_randint(state, 2) ? 0 : 1 + n_randint(state, FLINT_MAX(k, 1));
+        m3 = n_randint(state, 2) ? 0 : 1 + n_randint(state, FLINT_MAX(n, 1));
+        alg1 = n_randint(state, 4);
+        alg2 = n_randint(state, 4);
+        alg3 = n_randint(state, 4);
+
+        arb_init(x);
+        arb_init(xk);
+        y = _arb_vec_init(len);
+        ya = _arb_vec_init(len);
+        yb = _arb_vec_init(len);
+        yayb = _arb_vec_init(len);
+
+        arb_randtest(x, state, prec, 10);
+        arb_add_ui(xk, x, k, prec);
+
+        rising_algorithm(y, x, n + k, m1, len, prec, alg1);
+        rising_algorithm(ya, x, k, m2, len, prec, alg2);
+        rising_algorithm(yb, xk, n, m3, len, prec, alg3);
+        _arb_poly_mullow(yayb, ya, len, yb, len, len, prec);
+
+        if (!_arb_poly_overlaps(y, len, yayb, len))
+        {
+            flint_printf("FAIL\n\n");
+            flint_printf("len = %wd, k = %wu, n = %wu, m1 = %wu, m2 = %wu, m3 = %wu\n\n", len, k, n, m1, m2, m3);
+            flint_printf("x = "); arb_printn(x, 100, 0); flint_printf("\n\n");
+            flint_printf("y = "); _arb_vec_printn(y, len, 100, 0); flint_printf("\n\n");
+            flint_printf("ya = "); _arb_vec_printn(ya, len, 100, 0); flint_printf("\n\n");
+            flint_printf("yb = "); _arb_vec_printn(yb, len, 100, 0); flint_printf("\n\n");
+            flint_printf("yayb = "); _arb_vec_printn(yayb, len, 100, 0); flint_printf("\n\n");
+            flint_abort();
+        }
+
+        arb_clear(x);
+        arb_clear(xk);
+        _arb_vec_clear(y, len);
+        _arb_vec_clear(ya, len);
+        _arb_vec_clear(yb, len);
+        _arb_vec_clear(yayb, len);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

doc/source/arb_hypgeom.rst

@@ -31,10 +31,42 @@ Rising factorials
     does not use an asymptotically fast algorithm.
     The degree *n* must be at least 2.
 
-.. function:: void arb_hypgeom_rising_ui_rs(arb_t res, const arb_t x, ulong n, ulong m, slong prec)
+.. function:: void arb_hypgeom_rising_ui_forward(arb_t res, const arb_t x, ulong n, slong prec)
+              void arb_hypgeom_rising_ui_bs(arb_t res, const arb_t x, ulong n, slong prec)
+              void arb_hypgeom_rising_ui_rs(arb_t res, const arb_t x, ulong n, ulong m, slong prec)
+              void arb_hypgeom_rising_ui_rec(arb_t res, const arb_t x, ulong n, slong prec)
+              void arb_hypgeom_rising_ui(arb_t y, const arb_t x, ulong n, slong prec)
+              void arb_hypgeom_rising(arb_t y, const arb_t x, const arb_t n, slong prec)
+
+    Computes the rising factorial `(x)_n`.
+
+    The *forward* version uses the forward recurrence.
+    The *bs* version uses binary splitting.
+    The *rs* version uses rectangular splitting. It takes an extra tuning
+    parameter *m* which can be set to zero to choose automatically.
+    The *rec* version chooses an algorithm automatically, avoiding
+    use of the gamma function (so that it can be used in the computation
+    of the gamma function).
+    The default versions (*rising_ui* and *rising_ui*) choose an algorithm
+    automatically and may additionally fall back on the gamma function.
+
+.. function:: void arb_hypgeom_rising_ui_jet_powsum(arb_ptr res, const arb_t x, ulong n, slong len, slong prec)
+              void arb_hypgeom_rising_ui_jet_bs(arb_ptr res, const arb_t x, ulong n, slong len, slong prec)
+              void arb_hypgeom_rising_ui_jet_rs(arb_ptr res, const arb_t x, ulong n, ulong m, slong len, slong prec)
+              void arb_hypgeom_rising_ui_jet(arb_ptr res, const arb_t x, ulong n, slong len, slong prec)
+
+    Computes the jet of the rising factorial `(x)_n`, truncated to length *len*.
+    In other words, constructs the polynomial `(X + x)_n \in \mathbb{R}[X]`,
+    truncated if `\operatorname{len} < n + 1` (and zero-extended
+    if `\operatorname{len} > n + 1`).
+
+    The *powsum* version computes the sequence of powers of *x* and forms integral
+    linear combinations of these.
+    The *bs* version uses binary splitting.
+    The *rs* version uses rectangular splitting. It takes an extra tuning
+    parameter *m* which can be set to zero to choose automatically.
+    The default version chooses an algorithm automatically.
 
-    Computes the rising factorial `(x)_n` using rectangular splitting.
-    The splitting parameter *m* can be set to zero to choose automatically.
 
 Binomial coefficients
 -------------------------------------------------------------------------------