all the rising factorial algorithms for acb

acb.h

@@ -861,6 +861,14 @@ _acb_vec_set_round(acb_ptr res, acb_srcptr vec, slong len, slong prec)
         acb_set_round(res + i, vec + i, prec);
 }
 
+ACB_INLINE void
+_acb_vec_swap(acb_ptr res, acb_ptr vec, slong len)
+{
+    slong i;
+    for (i = 0; i < len; i++)
+        acb_swap(res + i, vec + i);
+}
+
 ACB_INLINE void
 _acb_vec_neg(acb_ptr res, acb_srcptr vec, slong len)
 {
@@ -1016,6 +1024,18 @@ acb_printn(const acb_t x, slong digits, ulong flags)
     acb_fprintn(stdout, x, digits, flags);
 }
 
+ACB_INLINE void
+_acb_vec_printn(acb_srcptr vec, slong len, slong ndigits, ulong flags)
+{
+    slong i;
+    for (i = 0; i < len; i++)
+    {
+        acb_printn(vec + i, ndigits, flags);
+        if (i < len - 1)
+            flint_printf(", ");
+    }
+}
+
 void acb_randtest(acb_t z, flint_rand_t state, slong prec, slong mag_bits);
 
 void acb_randtest_special(acb_t z, flint_rand_t state, slong prec, slong mag_bits);

acb_hypgeom.h

@@ -19,7 +19,18 @@
 extern "C" {
 #endif
 
+void acb_hypgeom_rising_ui_forward(acb_t res, const acb_t x, ulong n, slong prec);
 void acb_hypgeom_rising_ui_rs(acb_t res, const acb_t x, ulong n, ulong m, slong prec);
+void acb_hypgeom_rising_ui_bs(acb_t res, const acb_t x, ulong n, slong prec);
+void acb_hypgeom_rising_ui_rec(acb_t res, const acb_t x, ulong n, slong prec);
+void acb_hypgeom_rising_ui(acb_t y, const acb_t x, ulong n, slong prec);
+void acb_hypgeom_rising(acb_t y, const acb_t x, const acb_t n, slong prec);
+
+void acb_hypgeom_rising_ui_jet_powsum(acb_ptr res, const acb_t x, ulong n, slong len, slong prec);
+void acb_hypgeom_rising_ui_jet_rs(acb_ptr res, const acb_t x, ulong n, ulong m, slong len, slong prec);
+void acb_hypgeom_rising_ui_jet_bs(acb_ptr res, const acb_t x, ulong n, slong len, slong prec);
+void acb_hypgeom_rising_ui_jet(acb_ptr res, const acb_t x, ulong n, slong len, slong prec);
+
 
 void acb_hypgeom_pfq_bound_factor(mag_t C,
     acb_srcptr a, slong p, acb_srcptr b, slong q, const acb_t z, ulong n);

acb_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 "acb_hypgeom.h"
+
+void
+acb_hypgeom_rising_ui(acb_t y, const acb_t x, ulong n, slong prec)
+{
+    if (n < FLINT_MAX(prec, 100))
+    {
+        acb_hypgeom_rising_ui_rec(y, x, n, prec);
+    }
+    else
+    {
+        acb_t t;
+        acb_init(t);
+        acb_add_ui(t, x, n, prec);
+        acb_gamma(t, t, prec);
+        acb_rgamma(y, x, prec);
+        acb_mul(y, y, t, prec);
+        acb_clear(t);
+    }
+}
+
+void
+acb_hypgeom_rising(acb_t y, const acb_t x, const acb_t n, slong prec)
+{
+    if (acb_is_int(n) && arf_sgn(arb_midref(acb_realref(n))) >= 0 &&
+        arf_cmpabs_ui(arb_midref(acb_realref(n)), FLINT_MAX(prec, 100)) < 0)
+    {
+        acb_hypgeom_rising_ui_rec(y, x,
+            arf_get_si(arb_midref(acb_realref(n)), ARF_RND_DOWN), prec);
+    }
+    else
+    {
+        acb_t t;
+        acb_init(t);
+        acb_add(t, x, n, prec);
+        acb_gamma(t, t, prec);
+        acb_rgamma(y, x, prec);
+        acb_mul(y, y, t, prec);
+        acb_clear(t);
+    }
+}
+

acb_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 "acb_hypgeom.h"
+
+static void
+bsplit(acb_t y, const acb_t x, ulong a, ulong b, slong prec)
+{
+    if (b - a <= 4)
+    {
+        if (a == 0)
+        {
+            acb_hypgeom_rising_ui_forward(y, x, b, prec);
+        }
+        else
+        {
+            acb_add_ui(y, x, a, prec);
+            acb_hypgeom_rising_ui_forward(y, y, b - a, prec);
+        }
+    }
+    else
+    {
+        acb_t t, u;
+        ulong m = a + (b - a) / 2;
+
+        acb_init(t);
+        acb_init(u);
+
+        bsplit(t, x, a, m, prec);
+        bsplit(u, x, m, b, prec);
+
+        acb_mul(y, t, u, prec);
+
+        acb_clear(t);
+        acb_clear(u);
+    }
+}
+
+void
+acb_hypgeom_rising_ui_bs(acb_t res, const acb_t x, ulong n, slong prec)
+{
+    if (n <= 1)
+    {
+        if (n == 0)
+            acb_one(res);
+        else
+            acb_set_round(res, x, prec);
+        return;
+    }
+
+    {
+        acb_t t;
+        slong wp = ARF_PREC_ADD(prec, FLINT_BIT_COUNT(n));
+
+        acb_init(t);
+        bsplit(t, x, 0, n, wp);
+        acb_set_round(res, t, prec);
+        acb_clear(t);
+    }
+}
+

acb_hypgeom/rising_ui_forward.c

@@ -0,0 +1,47 @@
+/*
+    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 "acb_hypgeom.h"
+
+void _arb_increment_fast(arb_t x, slong prec);
+
+void
+acb_hypgeom_rising_ui_forward(acb_t res, const acb_t x, ulong n, slong prec)
+{
+    acb_t t;
+    ulong k;
+    slong wp;
+
+    if (n <= 1)
+    {
+        if (n == 0)
+            acb_one(res);
+        else
+            acb_set_round(res, x, prec);
+        return;
+    }
+
+    wp = prec + FLINT_BIT_COUNT(n);
+
+    acb_init(t);
+
+    acb_add_ui(t, x, 1, wp);
+    acb_mul(res, x, t, (n == 2) ? prec : wp);
+
+    for (k = 2; k < n; k++)
+    {
+        _arb_increment_fast(acb_realref(t), wp);
+        acb_mul(res, res, t, k == (n - 1) ? prec : wp);
+    }
+
+    acb_clear(t);
+}
+

acb_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 "acb_hypgeom.h"
+
+void
+acb_hypgeom_rising_ui_jet(acb_ptr res, const acb_t x, ulong n, slong len, slong prec)
+{
+    if (n <= 7)
+    {
+        acb_hypgeom_rising_ui_jet_powsum(res, x, n, len, prec);
+    }
+    else if (len == 2)
+    {
+        if (n <= 30 || acb_bits(x) >= prec / 128)
+            acb_hypgeom_rising_ui_jet_rs(res, x, n, 0, len, prec);
+        else
+            acb_hypgeom_rising_ui_jet_bs(res, x, n, len, prec);
+    }
+    else
+    {
+        if (n <= 20 || (n <= 200 && prec > 400 * n && acb_bits(x) >= prec / 4))
+        {
+            acb_hypgeom_rising_ui_jet_powsum(res, x, n, len, prec);
+        }
+        else if (len >= 64 || (acb_bits(x) + 1 < prec / 1024 && n >= 32))
+        {
+            acb_hypgeom_rising_ui_jet_bs(res, x, n, len, prec);
+        }
+        else
+        {
+            acb_hypgeom_rising_ui_jet_rs(res, x, n, 0, len, prec);
+        }
+    }
+}
+

acb_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 "acb_hypgeom.h"
+
+static void
+bsplit(acb_ptr res, const acb_t x, ulong a, ulong b, slong trunc, slong prec)
+{
+    trunc = FLINT_MIN(trunc, b - a + 1);
+
+    if (b - a <= 12)
+    {
+        if (a == 0)
+        {
+            acb_hypgeom_rising_ui_jet_powsum(res, x, b - a, FLINT_MIN(trunc, b - a + 1), prec);
+        }
+        else
+        {
+            acb_t t;
+            acb_init(t);
+            acb_add_ui(t, x, a, prec);
+            acb_hypgeom_rising_ui_jet_powsum(res, t, b - a, FLINT_MIN(trunc, b - a + 1), prec);
+            acb_clear(t);
+        }
+    }
+    else
+    {
+        acb_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 = _acb_vec_init(len1 + len2);
+        R = L + len1;
+
+        bsplit(L, x, a, m, trunc, prec);
+        bsplit(R, x, m, b, trunc, prec);
+
+        _acb_poly_mullow(res, L, len1, R, len2,
+            FLINT_MIN(trunc, len1 + len2 - 1), prec);
+
+        _acb_vec_clear(L, len1 + len2);
+    }
+}
+
+void
+acb_hypgeom_rising_ui_jet_bs(acb_ptr res, const acb_t x, ulong n, slong len, slong prec)
+{
+    if (len == 0)
+        return;
+
+    if (len > n + 1)
+    {
+        _acb_vec_zero(res + n + 1, len - n - 1);
+        len = n + 1;
+    }
+
+    if (len == n + 1)
+    {
+        acb_one(res + n);
+        len = n;
+    }
+
+    if (n <= 1)
+    {
+        if (n == 1)
+            acb_set_round(res, x, prec);
+        return;
+    }
+
+    bsplit(res, x, 0, n, len, prec);
+}
+

acb_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 "acb_hypgeom.h"
+
+void
+acb_hypgeom_rising_ui_jet_powsum(acb_ptr res, const acb_t x, ulong n, slong len, slong prec)
+{
+    slong i, j, k, wp;
+    acb_ptr xpow;
+    TMP_INIT;
+
+    if (len == 0)
+        return;
+
+    if (len > n + 1)
+    {
+        _acb_vec_zero(res + n + 1, len - n - 1);
+        len = n + 1;
+    }
+
+    if (len == n + 1)
+    {
+        acb_one(res + n);
+        len = n;
+    }
+
+    if (n <= 1)
+    {
+        if (n == 1)
+            acb_set_round(res, x, prec);
+        return;
+    }
+
+    if (len == 1)
+    {
+        acb_hypgeom_rising_ui_rs(res, x, n, 0, prec);
+        return;
+    }
+
+    if (n == 2)
+    {
+        acb_mul_2exp_si(res + 1, x, 1);
+        acb_add_ui(res + 1, res + 1, 1, prec);
+        acb_mul(res, x, x, prec + 4);
+        acb_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 = _acb_vec_init(n + 1);
+        _acb_vec_set_powers(xpow, x, n + 1, wp);
+
+        acb_dot_ui(res, NULL, 0, xpow + 1, 1, c + 1, 1, n, prec);
+
+        for (i = 1; i < len; i++)
+            acb_dot_ui(res + i, NULL, 0, xpow, 1, c + (n + 1) * i, 1, n + 1 - i, prec);
+
+        _acb_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 = _acb_vec_init(n + 1);
+        _acb_vec_set_powers(xpow, x, n + 1, wp);
+
+        acb_dot_fmpz(res, NULL, 0, xpow + 1, 1, c + 1, 1, n, prec);
+
+        for (i = 1; i < len; i++)
+            acb_dot_fmpz(res + i, NULL, 0, xpow, 1, c + (n + 1) * i, 1, n + 1 - i, prec);
+
+        _acb_vec_clear(xpow, n + 1);
+        _fmpz_vec_clear(c, (n + 1) * len);
+    }
+}
+

acb_hypgeom/rising_ui_jet_rs.c

@@ -0,0 +1,174 @@
+/*
+    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"
+#include "acb_hypgeom.h"
+
+void
+acb_hypgeom_rising_ui_jet_rs(acb_ptr res, const acb_t x, ulong n, ulong m, slong len, slong prec)
+{
+    slong i, j, k, l, m0, xmlen, tlen, ulen, climbs, climbs_max, wp;
+    acb_ptr tmp, xpow;
+    acb_ptr t, u;
+    mp_ptr c;
+    TMP_INIT;
+
+    if (len == 0)
+        return;
+
+    if (len > n + 1)
+    {
+        _acb_vec_zero(res + n + 1, len - n - 1);
+        len = n + 1;
+    }
+
+    if (len == n + 1)
+    {
+        acb_one(res + n);
+        len = n;
+    }
+
+    if (n <= 1)
+    {
+        if (n == 1)
+            acb_set_round(res, x, prec);
+        return;
+    }
+
+    if (len == 1)
+    {
+        acb_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 = _acb_vec_init(2 * len + (m + 1) * xmlen);
+    t = tmp;
+    u = tmp + len;
+    xpow = tmp + 2 * len;
+
+    _acb_vec_set_powers(xpow, x, m + 1, wp);
+
+    tlen = 1;
+
+    /* First derivatives */
+    for (i = 1; i <= m; i++)
+        acb_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++)
+                acb_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)
+        {
+            acb_add_ui(u, x, k, wp);
+            acb_one(u + 1);
+        }
+        else
+        {
+            if (climbs == 1)
+            {
+                _arb_hypgeom_rising_coeffs_1(c, k, l);
+                for (j = 0; j < ulen; j++)
+                    acb_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++)
+                    acb_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++)
+                    acb_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;
+            _acb_vec_swap(t, u, ulen);
+        }
+        else
+        {
+            _acb_poly_mullow(res, t, tlen, u, ulen, FLINT_MIN(len, tlen + ulen - 1), wp);
+            tlen = FLINT_MIN(len, tlen + ulen - 1);
+            _acb_vec_swap(t, res, tlen);
+        }
+    }
+
+    _acb_vec_set_round(res, t, len, prec);
+
+    _acb_vec_clear(tmp, 2 * len + (m + 1) * xmlen);
+    TMP_END;
+}
+

acb_hypgeom/rising_ui_rec.c

@@ -0,0 +1,54 @@
+/*
+    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"
+#include "acb_hypgeom.h"
+
+void
+acb_hypgeom_rising_ui_rec(acb_t res, const acb_t x, ulong n, slong prec)
+{
+    if (n <= 1)
+    {
+        if (n == 0)
+            acb_one(res);
+        else
+            acb_set_round(res, x, prec);
+        return;
+    }
+
+    if (arb_is_zero(acb_imagref(x)))
+    {
+        arb_hypgeom_rising_ui_rec(acb_realref(res), acb_realref(x), n, prec);
+        arb_zero(acb_imagref(res));
+        return;
+    }
+
+    if (n == 2 && prec <= 1024)
+    {
+        if (res != x)
+            acb_set(res, x);
+        acb_addmul(res, x, x, prec);
+        return;
+    }
+
+    if (n <= 5 && prec <= 512)
+    {
+        acb_hypgeom_rising_ui_forward(res, x, n, prec);
+    }
+    else
+    {
+        if (n >= 20 && acb_bits(x) < prec / 8)
+            acb_hypgeom_rising_ui_bs(res, x, n, prec);
+        else
+            acb_hypgeom_rising_ui_rs(res, x, n, 0, prec);
+    }
+}
+

acb_hypgeom/rising_ui_rs.c

@@ -43,7 +43,7 @@ acb_hypgeom_rising_ui_rs(acb_t res, const acb_t x, ulong n, ulong m, slong prec)
         else
         {
             m0 = n_sqrt(n);
-            m = 8 + 0.27 * pow(FLINT_MAX(0, prec - 1024), 0.4);
+            m = 8 + 0.27 * pow(1.5 * FLINT_MAX(0, prec - 1024), 0.4);
             m = FLINT_MIN(m, m0);
             m = FLINT_MIN(m, 64);
         }

acb_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 "acb_hypgeom.h"
 
+void
+rising_algorithm(acb_t res, const acb_t x, ulong n, ulong m, slong prec, int alg, int alias)
+{
+    if (alias)
+    {
+        acb_set(res, x);
+        rising_algorithm(res, res, n, m, prec, alg, 0);
+        return;
+    }
+
+    if (alg == 0)
+        acb_hypgeom_rising_ui_rs(res, x, n, m, prec);
+    else if (alg == 1)
+        acb_hypgeom_rising_ui_forward(res, x, n, prec);
+    else if (alg == 2)
+        acb_hypgeom_rising_ui_bs(res, x, n, prec);
+    else if (alg == 3)
+        acb_hypgeom_rising_ui_rec(res, x, n, prec);
+    else
+        acb_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++)
     {
         acb_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;
 
         acb_init(x);
         acb_init(xk);
@@ -41,12 +73,12 @@ int main()
         acb_init(yb);
         acb_init(yayb);
 
-        acb_randtest(x, state, prec, 10);
+        acb_randtest(x, state, prec, 10 + n_randint(state, 200));
         acb_add_ui(xk, x, k, prec);
 
-        acb_hypgeom_rising_ui_rs(y, x, n + k, m1, prec);
-        acb_hypgeom_rising_ui_rs(ya, x, k, m2, prec);
-        acb_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);
         acb_mul(yayb, ya, yb, prec);
 
         if (!acb_overlaps(y, yayb))

acb_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 "acb_hypgeom.h"
+
+void
+rising_algorithm(acb_ptr res, const acb_t x, ulong n, ulong m, slong len, slong prec, int alg)
+{
+    if (alg == 0)
+        acb_hypgeom_rising_ui_jet_powsum(res, x, n, len, prec);
+    else if (alg == 1)
+        acb_hypgeom_rising_ui_jet_rs(res, x, n, m, len, prec);
+    else if (alg == 2)
+        acb_hypgeom_rising_ui_jet_bs(res, x, n, len, prec);
+    else
+        acb_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++)
+    {
+        acb_t x, xk;
+        acb_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);
+
+        acb_init(x);
+        acb_init(xk);
+        y = _acb_vec_init(len);
+        ya = _acb_vec_init(len);
+        yb = _acb_vec_init(len);
+        yayb = _acb_vec_init(len);
+
+        acb_randtest(x, state, prec, 10);
+        acb_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);
+        _acb_poly_mullow(yayb, ya, len, yb, len, len, prec);
+
+        if (!_acb_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 = "); acb_printn(x, 100, 0); flint_printf("\n\n");
+            flint_printf("y = "); _acb_vec_printn(y, len, 100, 0); flint_printf("\n\n");
+            flint_printf("ya = "); _acb_vec_printn(ya, len, 100, 0); flint_printf("\n\n");
+            flint_printf("yb = "); _acb_vec_printn(yb, len, 100, 0); flint_printf("\n\n");
+            flint_printf("yayb = "); _acb_vec_printn(yayb, len, 100, 0); flint_printf("\n\n");
+            flint_abort();
+        }
+
+        acb_clear(x);
+        acb_clear(xk);
+        _acb_vec_clear(y, len);
+        _acb_vec_clear(ya, len);
+        _acb_vec_clear(yb, len);
+        _acb_vec_clear(yayb, len);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

doc/source/acb.rst

@@ -1181,6 +1181,10 @@ Vector functions
 
     Sets *res* to a copy of *vec*, rounding each entry to *prec* bits.
 
+.. function:: void _acb_vec_swap(acb_ptr vec1, acb_ptr vec2, slong len)
+
+    Swaps the entries of *vec1* and *vec2*.
+
 .. function:: void _acb_vec_neg(acb_ptr res, acb_srcptr vec, slong len)
 
 .. function:: void _acb_vec_add(acb_ptr res, acb_srcptr vec1, acb_srcptr vec2, slong len, slong prec)

doc/source/acb_hypgeom.rst

@@ -19,10 +19,42 @@ with prefactors that are products of exponentials, powers, and gamma functions.
 Rising factorials
 -------------------------------------------------------------------------------
 
-.. function:: void acb_hypgeom_rising_ui_rs(acb_t res, const acb_t x, ulong n, ulong m, slong prec)
+.. function:: void acb_hypgeom_rising_ui_forward(acb_t res, const acb_t x, ulong n, slong prec)
+              void acb_hypgeom_rising_ui_bs(acb_t res, const acb_t x, ulong n, slong prec)
+              void acb_hypgeom_rising_ui_rs(acb_t res, const acb_t x, ulong n, ulong m, slong prec)
+              void acb_hypgeom_rising_ui_rec(acb_t res, const acb_t x, ulong n, slong prec)
+              void acb_hypgeom_rising_ui(acb_t y, const acb_t x, ulong n, slong prec)
+              void acb_hypgeom_rising(acb_t y, const acb_t x, const acb_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 acb_hypgeom_rising_ui_jet_powsum(acb_ptr res, const acb_t x, ulong n, slong len, slong prec)
+              void acb_hypgeom_rising_ui_jet_bs(acb_ptr res, const acb_t x, ulong n, slong len, slong prec)
+              void acb_hypgeom_rising_ui_jet_rs(acb_ptr res, const acb_t x, ulong n, ulong m, slong len, slong prec)
+              void acb_hypgeom_rising_ui_jet(acb_ptr res, const acb_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.
 
 Convergent series
 -------------------------------------------------------------------------------