diff --git a/arb.h b/arb.h index 09cc729a..042c4b67 100644 --- a/arb.h +++ b/arb.h @@ -641,6 +641,7 @@ void arb_rising_ui_bs(arb_t y, const arb_t x, ulong n, long prec); void arb_rising_ui_rs(arb_t y, const arb_t x, ulong n, ulong m, long prec); void arb_rising_ui_rec(arb_t y, const arb_t x, ulong n, long prec); void arb_rising_ui(arb_t z, const arb_t x, ulong n, long prec); +void arb_rising_fmpq_ui(arb_t y, const fmpq_t x, ulong n, long prec); /* vector functions */ diff --git a/arb/rising_fmpq_ui.c b/arb/rising_fmpq_ui.c new file mode 100644 index 00000000..b4cc8d30 --- /dev/null +++ b/arb/rising_fmpq_ui.c @@ -0,0 +1,102 @@ +/*============================================================================= + + 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) 2013 Fredrik Johansson + +******************************************************************************/ + +#include "arb.h" + +static void +bsplit(arb_t y, const fmpz_t p, const fmpz_t q, ulong a, ulong b, long prec) +{ + if (b - a <= 8) + { + fmpz_t t, u; + ulong c; + + fmpz_init(t); + fmpz_init(u); + + fmpz_mul_ui(t, q, a); + fmpz_add(t, t, p); + fmpz_set(u, t); + + for (c = a + 1; c < b; c++) + { + fmpz_add(u, u, q); + fmpz_mul(t, t, u); + } + + arb_set_round_fmpz(y, t, prec); + + fmpz_clear(t); + fmpz_clear(u); + } + else + { + arb_t w; + ulong m = a + (b - a) / 2; + arb_init(w); + + bsplit(y, p, q, a, m, prec); + bsplit(w, p, q, m, b, prec); + + arb_mul(y, y, w, prec); + arb_clear(w); + } +} + +void +arb_rising_fmpq_ui(arb_t y, const fmpq_t x, ulong n, long prec) +{ + if (n == 0) + { + arb_one(y); + } + else if (n == 1) + { + arb_set_fmpq(y, x, prec); + } + else + { + long wp; + + wp = ARF_PREC_ADD(prec, FLINT_BIT_COUNT(n)); + + bsplit(y, fmpq_numref(x), fmpq_denref(x), 0, n, wp); + + if (fmpz_is_one(fmpq_denref(x))) + { + arb_set_round(y, y, prec); + } + else + { + arb_t t; + arb_init(t); + arb_set_fmpz(t, fmpq_denref(x)); + arb_pow_ui(t, t, n, wp); + arb_div(y, y, t, prec); + arb_clear(t); + } + } +} + diff --git a/doc/source/arb.rst b/doc/source/arb.rst index 92312fdf..14003ee5 100644 --- a/doc/source/arb.rst +++ b/doc/source/arb.rst @@ -747,6 +747,12 @@ Rising factorials The *rs* version takes an optional *step* parameter for tuning purposes (to use the default step length, pass zero). +.. function:: void arb_rising_fmpq_ui(arb_t z, const fmpq_t x, ulong n, long prec) + + Computes the rising factorial `z = x (x+1) (x+2) \cdots (x+n-1)` using + binary splitting. If the denominator or numerator of *x* is large + compared to *prec*, it is more efficient to convert *x* to an approximation + and use :func:`arb_rising_ui`. Special functions -------------------------------------------------------------------------------