diff --git a/arb.h b/arb.h index 10d93c4a..558b4800 100644 --- a/arb.h +++ b/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) { diff --git a/arb_hypgeom.h b/arb_hypgeom.h index 3e633153..0d0d6e46 100644 --- a/arb_hypgeom.h +++ b/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); diff --git a/arb_hypgeom/rising_ui.c b/arb_hypgeom/rising_ui.c new file mode 100644 index 00000000..778257f2 --- /dev/null +++ b/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 . +*/ + +#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); + } +} + diff --git a/arb_hypgeom/rising_ui_bs.c b/arb_hypgeom/rising_ui_bs.c new file mode 100644 index 00000000..fb9363a1 --- /dev/null +++ b/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 . +*/ + +#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); + } +} + diff --git a/arb_hypgeom/rising_ui_forward.c b/arb_hypgeom/rising_ui_forward.c new file mode 100644 index 00000000..1884d5eb --- /dev/null +++ b/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 . +*/ + +#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); +} + diff --git a/arb_hypgeom/rising_ui_jet.c b/arb_hypgeom/rising_ui_jet.c new file mode 100644 index 00000000..1ed2420b --- /dev/null +++ b/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 . +*/ + +#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); + } + } +} + diff --git a/arb_hypgeom/rising_ui_jet_bs.c b/arb_hypgeom/rising_ui_jet_bs.c new file mode 100644 index 00000000..9f9adc7d --- /dev/null +++ b/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 . +*/ + +#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); +} + diff --git a/arb_hypgeom/rising_ui_jet_powsum.c b/arb_hypgeom/rising_ui_jet_powsum.c new file mode 100644 index 00000000..2c5e5e88 --- /dev/null +++ b/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 . +*/ + +#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); + } +} + diff --git a/arb_hypgeom/rising_ui_jet_rs.c b/arb_hypgeom/rising_ui_jet_rs.c new file mode 100644 index 00000000..dfba87c8 --- /dev/null +++ b/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 . +*/ + +#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; +} + diff --git a/arb_hypgeom/rising_ui_rec.c b/arb_hypgeom/rising_ui_rec.c new file mode 100644 index 00000000..9cba3590 --- /dev/null +++ b/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 . +*/ + +#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); + } +} + diff --git a/arb_hypgeom/test/t-rising_ui_rs.c b/arb_hypgeom/test/t-rising_ui.c similarity index 60% rename from arb_hypgeom/test/t-rising_ui_rs.c rename to arb_hypgeom/test/t-rising_ui.c index 53046a93..37db8de9 100644 --- a/arb_hypgeom/test/t-rising_ui_rs.c +++ b/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)) diff --git a/arb_hypgeom/test/t-rising_ui_jet.c b/arb_hypgeom/test/t-rising_ui_jet.c new file mode 100644 index 00000000..5379c60d --- /dev/null +++ b/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 . +*/ + +#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; +} diff --git a/doc/source/arb_hypgeom.rst b/doc/source/arb_hypgeom.rst index d1172aa4..3b3db277 100644 --- a/doc/source/arb_hypgeom.rst +++ b/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 -------------------------------------------------------------------------------