From 52474120cee3fae344cf7c44b12962d28399d463 Mon Sep 17 00:00:00 2001 From: fredrik Date: Sun, 19 Sep 2021 15:07:07 +0200 Subject: [PATCH] begin make old gamma functions call new code --- acb.h | 6 - acb/digamma.c | 4 +- acb/gamma.c | 249 +------ acb/rising2_ui.c | 23 +- acb/rising2_ui_bs.c | 83 --- acb/rising2_ui_rs.c | 113 --- acb/rising_ui.c | 33 +- acb/rising_ui_bs.c | 80 --- acb/rising_ui_rec.c | 22 - acb/rising_ui_rs.c | 115 ---- acb/test/t-rising2_ui_bs.c | 110 --- acb/test/t-rising2_ui_rs.c | 110 --- acb/test/t-rising_ui_bs.c | 99 --- acb/test/t-rising_ui_rec.c | 99 --- acb/test/t-rising_ui_rs.c | 106 --- acb_hypgeom/gamma.c | 6 +- acb_poly/digamma_series.c | 4 +- acb_poly/gamma_series.c | 4 +- acb_poly/lgamma_series.c | 30 +- acb_poly/rgamma_series.c | 4 +- arb.h | 6 - arb/digamma.c | 4 +- arb/gamma.c | 946 +------------------------- arb/rising2_ui.c | 23 +- arb/rising2_ui_bs.c | 83 --- arb/rising2_ui_rs.c | 162 ----- arb/rising_ui.c | 33 +- arb/rising_ui_bs.c | 80 --- arb/rising_ui_rec.c | 22 - arb/rising_ui_rs.c | 154 ----- arb/test/t-rising2_ui_bs.c | 110 --- arb/test/t-rising2_ui_rs.c | 110 --- arb/test/t-rising_ui_bs.c | 109 --- arb/test/t-rising_ui_rec.c | 109 --- arb/test/t-rising_ui_rs.c | 112 --- arb_hypgeom.h | 2 +- arb_hypgeom/gamma_fmpq.c | 72 +- arb_hypgeom/test/t-gamma_taylor_tab.c | 2 +- arb_poly/digamma_series.c | 4 +- arb_poly/gamma_series.c | 4 +- arb_poly/lgamma_series.c | 16 +- arb_poly/rgamma_series.c | 4 +- doc/source/acb.rst | 33 +- doc/source/arb.rst | 39 +- 44 files changed, 179 insertions(+), 3360 deletions(-) delete mode 100644 acb/rising2_ui_bs.c delete mode 100644 acb/rising2_ui_rs.c delete mode 100644 acb/rising_ui_bs.c delete mode 100644 acb/rising_ui_rec.c delete mode 100644 acb/rising_ui_rs.c delete mode 100644 acb/test/t-rising2_ui_bs.c delete mode 100644 acb/test/t-rising2_ui_rs.c delete mode 100644 acb/test/t-rising_ui_bs.c delete mode 100644 acb/test/t-rising_ui_rec.c delete mode 100644 acb/test/t-rising_ui_rs.c delete mode 100644 arb/rising2_ui_bs.c delete mode 100644 arb/rising2_ui_rs.c delete mode 100644 arb/rising_ui_bs.c delete mode 100644 arb/rising_ui_rec.c delete mode 100644 arb/rising_ui_rs.c delete mode 100644 arb/test/t-rising2_ui_bs.c delete mode 100644 arb/test/t-rising2_ui_rs.c delete mode 100644 arb/test/t-rising_ui_bs.c delete mode 100644 arb/test/t-rising_ui_rec.c delete mode 100644 arb/test/t-rising_ui_rs.c diff --git a/acb.h b/acb.h index d8370952..5d438486 100644 --- a/acb.h +++ b/acb.h @@ -761,14 +761,8 @@ void acb_chebyshev_t2_ui(acb_t a, acb_t b, ulong n, const acb_t x, slong prec); void acb_chebyshev_u_ui(acb_t a, ulong n, const acb_t x, slong prec); void acb_chebyshev_u2_ui(acb_t a, acb_t b, ulong n, const acb_t x, slong prec); -void acb_rising_ui_bs(acb_t y, const acb_t x, ulong n, slong prec); -void acb_rising_ui_rs(acb_t y, const acb_t x, ulong n, ulong m, slong prec); -void acb_rising_ui_rec(acb_t y, const acb_t x, ulong n, slong prec); void acb_rising_ui(acb_t z, const acb_t x, ulong n, slong prec); void acb_rising(acb_t z, const acb_t x, const acb_t n, slong prec); - -void acb_rising2_ui_bs(acb_t u, acb_t v, const acb_t x, ulong n, slong prec); -void acb_rising2_ui_rs(acb_t u, acb_t v, const acb_t x, ulong n, ulong m, slong prec); void acb_rising2_ui(acb_t u, acb_t v, const acb_t x, ulong n, slong prec); void acb_rising_ui_get_mag(mag_t bound, const acb_t s, ulong n); diff --git a/acb/digamma.c b/acb/digamma.c index 7cc66349..359c0af1 100644 --- a/acb/digamma.c +++ b/acb/digamma.c @@ -11,7 +11,7 @@ #include "acb.h" -void acb_gamma_stirling_choose_param(int * reflect, slong * r, slong * n, +void acb_hypgeom_gamma_stirling_choose_param(int * reflect, slong * r, slong * n, const acb_t x, int use_reflect, int digamma, slong prec); void acb_gamma_stirling_eval(acb_t s, const acb_t z, slong nterms, int digamma, slong prec); @@ -32,7 +32,7 @@ acb_digamma(acb_t y, const acb_t x, slong prec) wp = prec + FLINT_BIT_COUNT(prec); - acb_gamma_stirling_choose_param(&reflect, &r, &n, x, 1, 1, wp); + acb_hypgeom_gamma_stirling_choose_param(&reflect, &r, &n, x, 1, 1, wp); acb_init(t); acb_init(u); diff --git a/acb/gamma.c b/acb/gamma.c index 21d62a90..ddfb6b3a 100644 --- a/acb/gamma.c +++ b/acb/gamma.c @@ -11,6 +11,7 @@ #include "bernoulli.h" #include "acb.h" +#include "acb_hypgeom.h" void acb_gamma_stirling_choose_param(int * reflect, slong * r, slong * n, @@ -20,6 +21,7 @@ void acb_gamma_stirling_bound(mag_ptr err, const acb_t z, slong k0, slong knum, void arb_gamma_stirling_coeff(arb_t b, ulong k, int digamma, slong prec); + void acb_gamma_stirling_eval(acb_t s, const acb_t z, slong nterms, int digamma, slong prec) { @@ -110,261 +112,20 @@ acb_gamma_stirling_eval(acb_t s, const acb_t z, slong nterms, int digamma, slong arb_clear(b); } -static void -_acb_gamma(acb_t y, const acb_t x, slong prec, int inverse) -{ - int reflect; - slong r, n, wp; - acb_t t, u, v; - double acc; - - wp = prec + FLINT_BIT_COUNT(prec); - - /* todo: for large x (if exact or accurate enough), increase precision */ - acc = acb_rel_accuracy_bits(x); - acc = FLINT_MAX(acc, 0); - wp = FLINT_MIN(prec, acc + 20); - wp = FLINT_MAX(wp, 2); - wp = wp + FLINT_BIT_COUNT(wp); - - if (acc < 3) /* try to avoid divisions blowing up */ - { - if (arf_cmp_d(arb_midref(acb_realref(x)), -0.5) < 0) - { - reflect = 1; - r = 0; - } - else if (arf_cmp_si(arb_midref(acb_realref(x)), 1) < 0) - { - reflect = 0; - r = 1; - } - else - { - reflect = 0; - r = 0; - } - - n = 1; - } - else - { - acb_gamma_stirling_choose_param(&reflect, &r, &n, x, 1, 0, wp); - } - - acb_init(t); - acb_init(u); - acb_init(v); - - if (reflect) - { - acb_sub_ui(t, x, 1, wp); - acb_neg(t, t); - acb_rising_ui_rec(u, t, r, wp); - arb_const_pi(acb_realref(v), wp); - acb_mul_arb(u, u, acb_realref(v), wp); - acb_add_ui(t, t, r, wp); - acb_gamma_stirling_eval(v, t, n, 0, wp); - - if (inverse) - { - /* rgamma(x) = gamma(1-x+r) sin(pi x) / ((rf(1-x, r) * pi) */ - acb_exp(v, v, wp); - acb_sin_pi(t, x, wp); - acb_mul(v, v, t, wp); - acb_mul(y, u, v, wp); - acb_div(y, v, u, prec); - } - else - { - /* gamma(x) = (rf(1-x, r) * pi) rgamma(1-x+r) csc(pi x) */ - acb_neg(v, v); - acb_exp(v, v, wp); - acb_csc_pi(t, x, wp); - acb_mul(v, v, t, wp); - acb_mul(y, v, u, prec); - } - } - else - { - acb_add_ui(t, x, r, wp); - acb_gamma_stirling_eval(u, t, n, 0, wp); - - if (inverse) - { - /* rgamma(x) = rf(x,r) rgamma(x+r) */ - acb_neg(u, u); - acb_exp(u, u, prec); - acb_rising_ui_rec(v, x, r, wp); - acb_mul(y, v, u, prec); - } - else - { - /* gamma(x) = gamma(x+r) / rf(x,r) */ - acb_exp(u, u, prec); - acb_rising_ui_rec(v, x, r, wp); - acb_div(y, u, v, prec); - } - } - - acb_clear(t); - acb_clear(u); - acb_clear(v); -} - void acb_gamma(acb_t y, const acb_t x, slong prec) { - if (acb_is_real(x)) - { - arb_gamma(acb_realref(y), acb_realref(x), prec); - arb_zero(acb_imagref(y)); - return; - } - - _acb_gamma(y, x, prec, 0); + acb_hypgeom_gamma(y, x, prec); } void acb_rgamma(acb_t y, const acb_t x, slong prec) { - if (acb_is_real(x)) - { - arb_rgamma(acb_realref(y), acb_realref(x), prec); - arb_zero(acb_imagref(y)); - return; - } - - _acb_gamma(y, x, prec, 1); -} - -/* corrects branch cut of sum_{k=0}^{r-1} log(z+k), given the - logarithm of the product */ -void -_acb_log_rising_correct_branch(acb_t t, - const acb_t t_wrong, const acb_t z, ulong r, slong prec) -{ - acb_t f; - arb_t pi, u, v; - fmpz_t pi_mult; - slong i, argprec; - - acb_init(f); - - arb_init(u); - arb_init(pi); - arb_init(v); - - fmpz_init(pi_mult); - - argprec = FLINT_MIN(prec, 40); - - arb_zero(u); - for (i = 0; i < r; i++) - { - acb_add_ui(f, z, i, argprec); - acb_arg(v, f, argprec); - arb_add(u, u, v, argprec); - } - - if (argprec == prec) - { - arb_set(acb_imagref(t), u); - } - else - { - arb_sub(v, u, acb_imagref(t), argprec); - arb_const_pi(pi, argprec); - arb_div(v, v, pi, argprec); - - if (arb_get_unique_fmpz(pi_mult, v)) - { - arb_const_pi(v, prec); - arb_mul_fmpz(v, v, pi_mult, prec); - arb_add(acb_imagref(t), acb_imagref(t), v, prec); - } - else - { - arb_zero(u); - for (i = 0; i < r; i++) - { - acb_add_ui(f, z, i, prec); - acb_arg(v, f, prec); - arb_add(u, u, v, prec); - } - arb_set(acb_imagref(t), u); - } - } - - acb_clear(f); - - arb_clear(u); - arb_clear(v); - arb_clear(pi); - - fmpz_clear(pi_mult); + acb_hypgeom_rgamma(y, x, prec); } void acb_lgamma(acb_t y, const acb_t x, slong prec) { - int reflect; - slong r, n, wp; - acb_t t, u, v; - - if (acb_is_real(x) && arb_is_positive(acb_realref(x))) - { - arb_lgamma(acb_realref(y), acb_realref(x), prec); - arb_zero(acb_imagref(y)); - return; - } - - wp = prec + FLINT_BIT_COUNT(prec); - - acb_gamma_stirling_choose_param(&reflect, &r, &n, x, 1, 0, wp); - - acb_init(t); - acb_init(u); - acb_init(v); - - if (reflect) - { - /* log gamma(x) = log rf(1-x, r) - log gamma(1-x+r) - log sin(pi x) + log(pi) */ - acb_sub_ui(u, x, 1, wp); - acb_neg(u, u); - - acb_rising_ui_rec(t, u, r, prec); - acb_log(t, t, wp); - _acb_log_rising_correct_branch(t, t, u, r, wp); - - acb_add_ui(u, u, r, wp); - acb_gamma_stirling_eval(v, u, n, 0, wp); - acb_sub(t, t, v, wp); - - acb_log_sin_pi(u, x, wp); - acb_sub(t, t, u, wp); - - acb_const_pi(u, wp); - acb_log(u, u, wp); - - acb_add(y, t, u, wp); - } - else - { - /* log gamma(x) = log gamma(x+r) - log rf(x,r) */ - - acb_add_ui(t, x, r, wp); - acb_gamma_stirling_eval(u, t, n, 0, wp); - - acb_rising_ui_rec(t, x, r, prec); - acb_log(t, t, wp); - _acb_log_rising_correct_branch(t, t, x, r, wp); - - acb_sub(y, u, t, prec); - } - - acb_clear(t); - acb_clear(u); - acb_clear(v); + acb_hypgeom_lgamma(y, x, prec); } - diff --git a/acb/rising2_ui.c b/acb/rising2_ui.c index ec45a030..eb95c66d 100644 --- a/acb/rising2_ui.c +++ b/acb/rising2_ui.c @@ -10,13 +10,30 @@ */ #include "acb.h" +#include "acb_hypgeom.h" void acb_rising2_ui(acb_t u, acb_t v, const acb_t x, ulong n, slong prec) { - if (prec < 256 || n < 8 || acb_bits(x) < prec / 8) - acb_rising2_ui_bs(u, v, x, n, prec); + if (x == u || x == v) + { + acb_t t; + acb_init(t); + acb_set(t, x); + acb_rising2_ui(u, v, t, n, prec); + acb_clear(t); + } else - acb_rising2_ui_rs(u, v, x, n, 0, prec); + { + acb_struct tmp[2]; + + tmp[0] = *u; + tmp[1] = *v; + + acb_hypgeom_rising_ui_jet(tmp, x, n, 2, prec); + + *u = tmp[0]; + *v = tmp[1]; + } } diff --git a/acb/rising2_ui_bs.c b/acb/rising2_ui_bs.c deleted file mode 100644 index 65f5b7e1..00000000 --- a/acb/rising2_ui_bs.c +++ /dev/null @@ -1,83 +0,0 @@ -/* - Copyright (C) 2013 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 "acb.h" - -static void -bsplit(acb_t p, acb_t q, const acb_t x, ulong a, ulong b, slong prec) -{ - if (b - a < 8) - { - ulong k; - acb_t t; - - acb_one(p); - acb_add_ui(q, x, a, prec); - - acb_init(t); - - for (k = a + 1; k < b; k++) - { - acb_add_ui(t, x, k, prec); - acb_mul(p, p, t, prec); - acb_add(p, p, q, prec); - acb_mul(q, q, t, prec); - } - - acb_clear(t); - } - else - { - acb_t r, s; - ulong m; - - acb_init(r); - acb_init(s); - - m = a + (b - a) / 2; - bsplit(p, q, x, a, m, prec); - bsplit(r, s, x, m, b, prec); - - acb_mul(p, p, s, prec); - acb_mul(r, r, q, prec); - acb_add(p, p, r, prec); - acb_mul(q, q, s, prec); - - acb_clear(r); - acb_clear(s); - } -} - -void -acb_rising2_ui_bs(acb_t u, acb_t v, const acb_t x, ulong n, slong prec) -{ - if (n == 0) - { - acb_zero(v); - acb_one(u); - } - else if (n == 1) - { - acb_set(u, x); - acb_one(v); - } - else - { - acb_t t; - slong wp = ARF_PREC_ADD(prec, FLINT_BIT_COUNT(n)); - - acb_init(t); /* support aliasing */ - acb_set(t, x); - bsplit(v, u, t, 0, n, wp); - acb_clear(t); - } -} - diff --git a/acb/rising2_ui_rs.c b/acb/rising2_ui_rs.c deleted file mode 100644 index d60cb38f..00000000 --- a/acb/rising2_ui_rs.c +++ /dev/null @@ -1,113 +0,0 @@ -/* - Copyright (C) 2013 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 "flint/fmpz_poly.h" -#include "acb.h" - -void _gamma_rf_bsplit(fmpz * A, ulong a, ulong b); - -void -acb_rising2_ui_rs(acb_t u, acb_t v, - const acb_t x, ulong n, ulong m, slong prec) -{ - if (n == 0) - { - acb_zero(v); - acb_one(u); - } - else if (n == 1) - { - acb_set(u, x); - acb_one(v); - } - else - { - slong wp; - ulong i, j, a, b; - acb_ptr xs; - acb_t S, T, U, V; - fmpz *A, *B; - - wp = ARF_PREC_ADD(prec, FLINT_BIT_COUNT(n)); - - if (m == 0) - { - ulong m1, m2; - m1 = 0.6 * pow(wp, 0.4); - m2 = n_sqrt(n); - m = FLINT_MIN(m1, m2); - } - - m = FLINT_MAX(m, 1); - - xs = _acb_vec_init(m + 1); - A = _fmpz_vec_init(2 * m + 1); - B = A + (m + 1); - - acb_init(S); - acb_init(T); - acb_init(U); - acb_init(V); - _acb_vec_set_powers(xs, x, m + 1, wp); - - for (i = 0; i < n; i += m) - { - a = i; - b = FLINT_MIN(n, a + m); - - if (a == 0 || b != a + m) - { - _gamma_rf_bsplit(A, a, b); - } - else - { - fmpz tt = m; - _fmpz_poly_taylor_shift(A, &tt, m + 1); - } - - _fmpz_poly_derivative(B, A, b - a + 1); - - acb_set_fmpz(S, A); - - for (j = 1; j <= b - a; j++) - acb_addmul_fmpz(S, xs + j, A + j, wp); - - acb_set_fmpz(T, B); - - for (j = 1; j < b - a; j++) - acb_addmul_fmpz(T, xs + j, B + j, wp); - - if (i == 0) - { - acb_set(U, S); - acb_set(V, T); - } - else - { - acb_mul(V, V, S, wp); - acb_addmul(V, U, T, wp); - acb_mul(U, U, S, wp); - } - } - - acb_set(u, U); - acb_set(v, V); - - _acb_vec_clear(xs, m + 1); - _fmpz_vec_clear(A, 2 * m + 1); - - acb_clear(S); - acb_clear(T); - acb_clear(U); - acb_clear(V); - } -} - diff --git a/acb/rising_ui.c b/acb/rising_ui.c index 26a49f91..8e997c8c 100644 --- a/acb/rising_ui.c +++ b/acb/rising_ui.c @@ -10,44 +10,17 @@ */ #include "acb.h" +#include "acb_hypgeom.h" void acb_rising_ui(acb_t y, const acb_t x, ulong n, slong prec) { - if (n < FLINT_MAX(prec, 100)) - { - acb_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); - } + acb_hypgeom_rising_ui(y, x, n, prec); } void acb_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_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(y, x, n, prec); } diff --git a/acb/rising_ui_bs.c b/acb/rising_ui_bs.c deleted file mode 100644 index d617178c..00000000 --- a/acb/rising_ui_bs.c +++ /dev/null @@ -1,80 +0,0 @@ -/* - 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 "acb.h" - -/* assumes y and x are not aliased */ -static void -bsplit(acb_t y, const acb_t x, ulong a, ulong b, slong prec) -{ - if (b - a == 1) - { - acb_set_round(y, x, prec); - } - else if (b - a <= 10) - { - slong i; - acb_t t; - acb_init(t); - - acb_add_ui(t, x, a, prec); - acb_add_ui(y, x, a + 1, prec); - acb_mul(y, y, t, prec); - - for (i = a + 2; i < b; i++) - { - acb_add_ui(t, x, i, prec); - acb_mul(y, y, t, prec); - } - - acb_clear(t); - } - 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_rising_ui_bs(acb_t y, const acb_t x, ulong n, slong prec) -{ - if (n == 0) - { - acb_one(y); - } - else if (n == 1) - { - acb_set_round(y, x, prec); - } - else - { - 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(y, t, prec); - acb_clear(t); - } -} - diff --git a/acb/rising_ui_rec.c b/acb/rising_ui_rec.c deleted file mode 100644 index 0d579fad..00000000 --- a/acb/rising_ui_rec.c +++ /dev/null @@ -1,22 +0,0 @@ -/* - 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 "acb.h" - -void -acb_rising_ui_rec(acb_t y, const acb_t x, ulong n, slong prec) -{ - if (prec < 256 || n < 8 || acb_bits(x) < prec / 8) - acb_rising_ui_bs(y, x, n, prec); - else - acb_rising_ui_rs(y, x, n, 0, prec); -} - diff --git a/acb/rising_ui_rs.c b/acb/rising_ui_rs.c deleted file mode 100644 index dcda2b22..00000000 --- a/acb/rising_ui_rs.c +++ /dev/null @@ -1,115 +0,0 @@ -/* - 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 "flint/fmpz_poly.h" -#include "acb.h" - -void rising_difference_polynomial(fmpz * s, fmpz * c, ulong m); - -void -acb_rising_ui_rs(acb_t y, const acb_t x, ulong n, ulong m, slong prec) -{ - acb_ptr xs; - acb_t t, u, v; - ulong i, k, rem; - fmpz_t c, h; - fmpz *s, *d; - slong wp; - - if (n == 0) - { - acb_one(y); - return; - } - - if (n == 1) - { - acb_set_round(y, x, prec); - return; - } - - wp = ARF_PREC_ADD(prec, FLINT_BIT_COUNT(n)); - - acb_init(t); - acb_init(u); - acb_init(v); - fmpz_init(c); - fmpz_init(h); - - if (m == 0) - { - ulong m1, m2; - m1 = 0.2 * pow(2.0 * wp, 0.4); - m2 = n_sqrt(n); - m = FLINT_MIN(m1, m2); - } - - m = FLINT_MIN(m, n); - m = FLINT_MAX(m, 1); - - xs = _acb_vec_init(m + 1); - d = _fmpz_vec_init(m * m); - s = _fmpz_vec_init(m + 1); - - _acb_vec_set_powers(xs, x, m + 1, wp); - - rising_difference_polynomial(s, d, m); - - /* tail */ - rem = m; - while (rem + m <= n) - rem += m; - acb_one(y); - for (k = rem; k < n; k++) - { - acb_add_ui(t, xs + 1, k, wp); - acb_mul(y, y, t, wp); - } - - /* initial rising factorial */ - acb_zero(t); - for (i = 1; i <= m; i++) - acb_addmul_fmpz(t, xs + i, s + i, wp); - - acb_mul(y, y, t, wp); - - /* the leading coefficient is always the same */ - acb_mul_fmpz(xs + m - 1, xs + m - 1, d + m - 1 + 0, wp); - - for (k = 0; k + 2 * m <= n; k += m) - { - for (i = 0; i < m - 1; i++) - { - fmpz_set_ui(h, k); - _fmpz_poly_evaluate_horner_fmpz(c, d + i * m, m - i, h); - - if (i == 0) - acb_add_fmpz(t, t, c, wp); - else - acb_addmul_fmpz(t, xs + i, c, wp); - } - - acb_add(t, t, xs + m - 1, wp); - acb_mul(y, y, t, wp); - } - - acb_set_round(y, y, prec); - - acb_clear(t); - acb_clear(u); - acb_clear(v); - _acb_vec_clear(xs, m + 1); - _fmpz_vec_clear(d, m * m); - _fmpz_vec_clear(s, m + 1); - fmpz_clear(c); - fmpz_clear(h); -} - diff --git a/acb/test/t-rising2_ui_bs.c b/acb/test/t-rising2_ui_bs.c deleted file mode 100644 index 70331aa0..00000000 --- a/acb/test/t-rising2_ui_bs.c +++ /dev/null @@ -1,110 +0,0 @@ -/* - Copyright (C) 2013 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 "flint/arith.h" -#include "acb_poly.h" - -int main() -{ - slong iter; - flint_rand_t state; - - flint_printf("rising2_ui_bs...."); - fflush(stdout); - - flint_randinit(state); - - for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++) - { - acb_t a, u, v, u2, v2; - fmpz *f; - acb_ptr g; - ulong n; - slong i, prec; - - acb_init(a); - acb_init(u); - acb_init(v); - acb_init(u2); - acb_init(v2); - - acb_randtest(a, state, 1 + n_randint(state, 4000), 10); - acb_randtest(u, state, 1 + n_randint(state, 4000), 10); - acb_randtest(v, state, 1 + n_randint(state, 4000), 10); - n = n_randint(state, 120); - - f = _fmpz_vec_init(n + 1); - g = _acb_vec_init(n + 1); - - prec = 2 + n_randint(state, 4000); - acb_rising2_ui_bs(u, v, a, n, prec); - - arith_stirling_number_1u_vec(f, n, n + 1); - for (i = 0; i <= n; i++) - acb_set_fmpz(g + i, f + i); - _acb_poly_evaluate(u2, g, n + 1, a, prec); - - _acb_poly_derivative(g, g, n + 1, prec); - _acb_poly_evaluate(v2, g, n, a, prec); - - if (!acb_overlaps(u, u2) || !acb_overlaps(v, v2)) - { - flint_printf("FAIL: overlap\n\n"); - flint_printf("n = %wu\n", n); - flint_printf("a = "); acb_printd(a, 15); flint_printf("\n\n"); - flint_printf("u = "); acb_printd(u, 15); flint_printf("\n\n"); - flint_printf("u2 = "); acb_printd(u2, 15); flint_printf("\n\n"); - flint_printf("v = "); acb_printd(v, 15); flint_printf("\n\n"); - flint_printf("v2 = "); acb_printd(v2, 15); flint_printf("\n\n"); - flint_abort(); - } - - acb_set(u2, a); - acb_rising2_ui_bs(u2, v, u2, n, prec); - - if (!acb_equal(u2, u)) - { - flint_printf("FAIL: aliasing 1\n\n"); - flint_printf("a = "); acb_printd(a, 15); flint_printf("\n\n"); - flint_printf("u = "); acb_printd(u, 15); flint_printf("\n\n"); - flint_printf("u2 = "); acb_printd(u2, 15); flint_printf("\n\n"); - flint_printf("n = %wu\n", n); - flint_abort(); - } - - acb_set(v2, a); - acb_rising2_ui_bs(u, v2, v2, n, prec); - - if (!acb_equal(v2, v)) - { - flint_printf("FAIL: aliasing 2\n\n"); - flint_printf("a = "); acb_printd(a, 15); flint_printf("\n\n"); - flint_printf("v = "); acb_printd(v, 15); flint_printf("\n\n"); - flint_printf("v2 = "); acb_printd(v2, 15); flint_printf("\n\n"); - flint_printf("n = %wu\n", n); - flint_abort(); - } - - acb_clear(a); - acb_clear(u); - acb_clear(v); - acb_clear(u2); - acb_clear(v2); - _fmpz_vec_clear(f, n + 1); - _acb_vec_clear(g, n + 1); - } - - flint_randclear(state); - flint_cleanup(); - flint_printf("PASS\n"); - return EXIT_SUCCESS; -} - diff --git a/acb/test/t-rising2_ui_rs.c b/acb/test/t-rising2_ui_rs.c deleted file mode 100644 index 10d7fa13..00000000 --- a/acb/test/t-rising2_ui_rs.c +++ /dev/null @@ -1,110 +0,0 @@ -/* - Copyright (C) 2013 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 "flint/arith.h" -#include "acb_poly.h" - -int main() -{ - slong iter; - flint_rand_t state; - - flint_printf("rising2_ui_rs...."); - fflush(stdout); - - flint_randinit(state); - - for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++) - { - acb_t a, u, v, u2, v2; - fmpz *f; - acb_ptr g; - ulong n; - slong i, prec; - - acb_init(a); - acb_init(u); - acb_init(v); - acb_init(u2); - acb_init(v2); - - acb_randtest(a, state, 1 + n_randint(state, 4000), 10); - acb_randtest(u, state, 1 + n_randint(state, 4000), 10); - acb_randtest(v, state, 1 + n_randint(state, 4000), 10); - n = n_randint(state, 120); - - f = _fmpz_vec_init(n + 1); - g = _acb_vec_init(n + 1); - - prec = 2 + n_randint(state, 4000); - acb_rising2_ui_rs(u, v, a, n, 0, prec); - - arith_stirling_number_1u_vec(f, n, n + 1); - for (i = 0; i <= n; i++) - acb_set_fmpz(g + i, f + i); - _acb_poly_evaluate(u2, g, n + 1, a, prec); - - _acb_poly_derivative(g, g, n + 1, prec); - _acb_poly_evaluate(v2, g, n, a, prec); - - if (!acb_overlaps(u, u2) || !acb_overlaps(v, v2)) - { - flint_printf("FAIL: overlap\n\n"); - flint_printf("n = %wu\n", n); - flint_printf("a = "); acb_printd(a, 15); flint_printf("\n\n"); - flint_printf("u = "); acb_printd(u, 15); flint_printf("\n\n"); - flint_printf("u2 = "); acb_printd(u2, 15); flint_printf("\n\n"); - flint_printf("v = "); acb_printd(v, 15); flint_printf("\n\n"); - flint_printf("v2 = "); acb_printd(v2, 15); flint_printf("\n\n"); - flint_abort(); - } - - acb_set(u2, a); - acb_rising2_ui_rs(u2, v, u2, n, 0, prec); - - if (!acb_equal(u2, u)) - { - flint_printf("FAIL: aliasing 1\n\n"); - flint_printf("a = "); acb_printd(a, 15); flint_printf("\n\n"); - flint_printf("u = "); acb_printd(u, 15); flint_printf("\n\n"); - flint_printf("u2 = "); acb_printd(u2, 15); flint_printf("\n\n"); - flint_printf("n = %wu\n", n); - flint_abort(); - } - - acb_set(v2, a); - acb_rising2_ui_rs(u, v2, v2, n, 0, prec); - - if (!acb_equal(v2, v)) - { - flint_printf("FAIL: aliasing 2\n\n"); - flint_printf("a = "); acb_printd(a, 15); flint_printf("\n\n"); - flint_printf("v = "); acb_printd(v, 15); flint_printf("\n\n"); - flint_printf("v2 = "); acb_printd(v2, 15); flint_printf("\n\n"); - flint_printf("n = %wu\n", n); - flint_abort(); - } - - acb_clear(a); - acb_clear(u); - acb_clear(v); - acb_clear(u2); - acb_clear(v2); - _fmpz_vec_clear(f, n + 1); - _acb_vec_clear(g, n + 1); - } - - flint_randclear(state); - flint_cleanup(); - flint_printf("PASS\n"); - return EXIT_SUCCESS; -} - diff --git a/acb/test/t-rising_ui_bs.c b/acb/test/t-rising_ui_bs.c deleted file mode 100644 index 2eef9d95..00000000 --- a/acb/test/t-rising_ui_bs.c +++ /dev/null @@ -1,99 +0,0 @@ -/* - Copyright (C) 2012 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 "acb.h" - -int main() -{ - slong iter; - flint_rand_t state; - - flint_printf("rising_ui_bs...."); - fflush(stdout); - - flint_randinit(state); - - /* check functional equation */ - for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++) - { - acb_t x, xn, y, z; - ulong n, m; - - acb_init(x); - acb_init(xn); - acb_init(y); - acb_init(z); - - acb_randtest(x, state, 1 + n_randint(state, 4000), 10); - n = n_randint(state, 40); - m = n_randint(state, 40); - acb_add_ui(xn, x, n, 1 + n_randint(state, 4000)); - - acb_rising_ui_bs(y, x, n, 2 + n_randint(state, 4000)); - acb_rising_ui_bs(z, xn, m, 2 + n_randint(state, 4000)); - acb_mul(y, y, z, 2 + n_randint(state, 4000)); - - acb_rising_ui_bs(z, x, n + m, 2 + n_randint(state, 4000)); - - if (!acb_overlaps(y, z)) - { - flint_printf("FAIL: overlap\n\n"); - flint_printf("n = %wu\n", n); - flint_printf("m = %wu\n", m); - flint_printf("x = "); acb_print(x); flint_printf("\n\n"); - flint_printf("xn = "); acb_print(xn); flint_printf("\n\n"); - flint_printf("y = "); acb_print(y); flint_printf("\n\n"); - flint_printf("z = "); acb_print(z); flint_printf("\n\n"); - flint_abort(); - } - - acb_clear(x); - acb_clear(xn); - acb_clear(y); - acb_clear(z); - } - - /* aliasing of y and x */ - for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++) - { - acb_t x, y; - ulong n; - slong prec; - - acb_init(x); - acb_init(y); - - acb_randtest(x, state, 1 + n_randint(state, 200), 10); - acb_randtest(y, state, 1 + n_randint(state, 200), 10); - n = n_randint(state, 100); - - prec = 2 + n_randint(state, 4000); - acb_rising_ui_bs(y, x, n, prec); - acb_rising_ui_bs(x, x, n, prec); - - if (!acb_equal(x, y)) - { - flint_printf("FAIL: aliasing\n\n"); - flint_printf("x = "); acb_print(x); flint_printf("\n\n"); - flint_printf("y = "); acb_print(y); flint_printf("\n\n"); - flint_printf("n = %wu\n", n); - flint_abort(); - } - - acb_clear(x); - acb_clear(y); - } - - flint_randclear(state); - flint_cleanup(); - flint_printf("PASS\n"); - return EXIT_SUCCESS; -} diff --git a/acb/test/t-rising_ui_rec.c b/acb/test/t-rising_ui_rec.c deleted file mode 100644 index 997059f1..00000000 --- a/acb/test/t-rising_ui_rec.c +++ /dev/null @@ -1,99 +0,0 @@ -/* - Copyright (C) 2012 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 "acb.h" - -int main() -{ - slong iter; - flint_rand_t state; - - flint_printf("rising_ui_rec...."); - fflush(stdout); - - flint_randinit(state); - - /* check functional equation */ - for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++) - { - acb_t x, xn, y, z; - ulong n, m; - - acb_init(x); - acb_init(xn); - acb_init(y); - acb_init(z); - - acb_randtest(x, state, 1 + n_randint(state, 4000), 10); - n = n_randint(state, 80); - m = n_randint(state, 40); - acb_add_ui(xn, x, n, 1 + n_randint(state, 4000)); - - acb_rising_ui_rec(y, x, n, 2 + n_randint(state, 4000)); - acb_rising_ui_rec(z, xn, m, 2 + n_randint(state, 4000)); - acb_mul(y, y, z, 2 + n_randint(state, 4000)); - - acb_rising_ui_rec(z, x, n + m, 2 + n_randint(state, 4000)); - - if (!acb_overlaps(y, z)) - { - flint_printf("FAIL: overlap\n\n"); - flint_printf("n = %wu\n", n); - flint_printf("m = %wu\n", m); - flint_printf("x = "); acb_print(x); flint_printf("\n\n"); - flint_printf("xn = "); acb_print(xn); flint_printf("\n\n"); - flint_printf("y = "); acb_print(y); flint_printf("\n\n"); - flint_printf("z = "); acb_print(z); flint_printf("\n\n"); - flint_abort(); - } - - acb_clear(x); - acb_clear(xn); - acb_clear(y); - acb_clear(z); - } - - /* aliasing of y and x */ - for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++) - { - acb_t x, y; - ulong n; - slong prec; - - acb_init(x); - acb_init(y); - - acb_randtest(x, state, 1 + n_randint(state, 200), 10); - acb_randtest(y, state, 1 + n_randint(state, 200), 10); - n = n_randint(state, 100); - - prec = 2 + n_randint(state, 4000); - acb_rising_ui_rec(y, x, n, prec); - acb_rising_ui_rec(x, x, n, prec); - - if (!acb_equal(x, y)) - { - flint_printf("FAIL: aliasing\n\n"); - flint_printf("x = "); acb_print(x); flint_printf("\n\n"); - flint_printf("y = "); acb_print(y); flint_printf("\n\n"); - flint_printf("n = %wu\n", n); - flint_abort(); - } - - acb_clear(x); - acb_clear(y); - } - - flint_randclear(state); - flint_cleanup(); - flint_printf("PASS\n"); - return EXIT_SUCCESS; -} diff --git a/acb/test/t-rising_ui_rs.c b/acb/test/t-rising_ui_rs.c deleted file mode 100644 index 2118c91d..00000000 --- a/acb/test/t-rising_ui_rs.c +++ /dev/null @@ -1,106 +0,0 @@ -/* - Copyright (C) 2012 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 "acb.h" - -int main() -{ - slong iter; - flint_rand_t state; - - flint_printf("rising_ui_rs...."); - fflush(stdout); - - flint_randinit(state); - - /* check functional equation */ - for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++) - { - acb_t x, xn, y, z; - ulong n, m, step1, step2, step3; - - acb_init(x); - acb_init(xn); - acb_init(y); - acb_init(z); - - acb_randtest(x, state, 1 + n_randint(state, 4000), 10); - n = n_randint(state, 80); - m = n_randint(state, 40); - acb_add_ui(xn, x, n, 1 + n_randint(state, 4000)); - - step1 = n_randint(state, 20); - step2 = n_randint(state, 20); - step3 = n_randint(state, 20); - - acb_rising_ui_rs(y, x, n, step1, 2 + n_randint(state, 4000)); - acb_rising_ui_rs(z, xn, m, step2, 2 + n_randint(state, 4000)); - acb_mul(y, y, z, 2 + n_randint(state, 4000)); - - acb_rising_ui_rs(z, x, n + m, step3, 2 + n_randint(state, 4000)); - - if (!acb_overlaps(y, z)) - { - flint_printf("FAIL: overlap\n\n"); - flint_printf("n = %wu\n", n); - flint_printf("m = %wu\n", m); - flint_printf("x = "); acb_print(x); flint_printf("\n\n"); - flint_printf("xn = "); acb_print(xn); flint_printf("\n\n"); - flint_printf("y = "); acb_print(y); flint_printf("\n\n"); - flint_printf("z = "); acb_print(z); flint_printf("\n\n"); - flint_abort(); - } - - acb_clear(x); - acb_clear(xn); - acb_clear(y); - acb_clear(z); - } - - /* aliasing of y and x */ - for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++) - { - acb_t x, y; - ulong n; - slong prec; - ulong step; - - acb_init(x); - acb_init(y); - - acb_randtest(x, state, 1 + n_randint(state, 200), 10); - acb_randtest(y, state, 1 + n_randint(state, 200), 10); - n = n_randint(state, 100); - - step = n_randint(state, 20); - - prec = 2 + n_randint(state, 4000); - acb_rising_ui_rs(y, x, n, step, prec); - acb_rising_ui_rs(x, x, n, step, prec); - - if (!acb_equal(x, y)) - { - flint_printf("FAIL: aliasing\n\n"); - flint_printf("x = "); acb_print(x); flint_printf("\n\n"); - flint_printf("y = "); acb_print(y); flint_printf("\n\n"); - flint_printf("n = %wu\n", n); - flint_abort(); - } - - acb_clear(x); - acb_clear(y); - } - - flint_randclear(state); - flint_cleanup(); - flint_printf("PASS\n"); - return EXIT_SUCCESS; -} diff --git a/acb_hypgeom/gamma.c b/acb_hypgeom/gamma.c index 4d8c3d34..5ae70db4 100644 --- a/acb_hypgeom/gamma.c +++ b/acb_hypgeom/gamma.c @@ -105,7 +105,7 @@ acb_hypgeom_gamma_stirling(acb_t y, const acb_t x, int reciprocal, slong prec) { acb_sub_ui(t, x, 1, wp); acb_neg(t, t); - acb_rising_ui_rec(u, t, r, wp); + acb_hypgeom_rising_ui_rec(u, t, r, wp); arb_const_pi(acb_realref(v), wp); acb_mul_arb(u, u, acb_realref(v), wp); acb_add_ui(t, t, r, wp); @@ -140,14 +140,14 @@ acb_hypgeom_gamma_stirling(acb_t y, const acb_t x, int reciprocal, slong prec) /* rgamma(x) = rf(x,r) rgamma(x+r) */ acb_neg(u, u); acb_exp(u, u, prec); - acb_rising_ui_rec(v, x, r, wp); + acb_hypgeom_rising_ui_rec(v, x, r, wp); acb_mul(y, v, u, prec); } else { /* gamma(x) = gamma(x+r) / rf(x,r) */ acb_exp(u, u, prec); - acb_rising_ui_rec(v, x, r, wp); + acb_hypgeom_rising_ui_rec(v, x, r, wp); acb_div(y, u, v, prec); } } diff --git a/acb_poly/digamma_series.c b/acb_poly/digamma_series.c index 7cc2d427..c1e47dc2 100644 --- a/acb_poly/digamma_series.c +++ b/acb_poly/digamma_series.c @@ -11,7 +11,7 @@ #include "acb_poly.h" -void acb_gamma_stirling_choose_param(int * reflect, slong * r, slong * n, +void acb_hypgeom_gamma_stirling_choose_param(int * reflect, slong * r, slong * n, const acb_t x, int use_reflect, int digamma, slong prec); void @@ -60,7 +60,7 @@ _acb_poly_digamma_series(acb_ptr res, acb_srcptr h, slong hlen, slong len, slong acb_init(zr); /* use Stirling series */ - acb_gamma_stirling_choose_param(&reflect, &r, &n, h, 1, 1, wp); + acb_hypgeom_gamma_stirling_choose_param(&reflect, &r, &n, h, 1, 1, wp); /* psi(x) = psi((1-x)+r) - h(1-x,r) - pi*cot(pi*x) */ if (reflect) diff --git a/acb_poly/gamma_series.c b/acb_poly/gamma_series.c index d03cc6d0..39d533f9 100644 --- a/acb_poly/gamma_series.c +++ b/acb_poly/gamma_series.c @@ -13,7 +13,7 @@ void acb_gamma_stirling_bound(mag_ptr err, const acb_t x, slong k0, slong knum, slong n); -void acb_gamma_stirling_choose_param(int * reflect, slong * r, slong * n, +void acb_hypgeom_gamma_stirling_choose_param(int * reflect, slong * r, slong * n, const acb_t x, int use_reflect, int digamma, slong prec); void arb_gamma_stirling_coeff(arb_t b, ulong k, int digamma, slong prec); @@ -232,7 +232,7 @@ _acb_poly_gamma_series(acb_ptr res, acb_srcptr h, slong hlen, slong len, slong p acb_init(f + 1); /* use Stirling series */ - acb_gamma_stirling_choose_param(&reflect, &r, &n, h, 1, 0, wp); + acb_hypgeom_gamma_stirling_choose_param(&reflect, &r, &n, h, 1, 0, wp); /* gamma(h) = (rf(1-h, r) * pi) / (gamma(1-h+r) sin(pi h)), h = h0 + t*/ if (reflect) diff --git a/acb_poly/lgamma_series.c b/acb_poly/lgamma_series.c index dc3fd525..504b0cb4 100644 --- a/acb_poly/lgamma_series.c +++ b/acb_poly/lgamma_series.c @@ -15,34 +15,12 @@ void _acb_log_rising_correct_branch(acb_t t, const acb_t t_wrong, const acb_t z, ulong r, slong prec); -void acb_gamma_stirling_choose_param(int * reflect, slong * r, slong * n, +void acb_hypgeom_gamma_stirling_choose_param(int * reflect, slong * r, slong * n, const acb_t x, int use_reflect, int digamma, slong prec); void _acb_poly_gamma_stirling_eval(acb_ptr res, const acb_t z, slong n, slong num, slong prec); -static __inline__ void -_log_rising_ui_series(acb_ptr t, const acb_t x, slong r, slong len, slong prec) -{ - acb_struct f[2]; - slong rflen; - - acb_init(f); - acb_init(f + 1); - - acb_set(f, x); - acb_one(f + 1); - - rflen = FLINT_MIN(len, r + 1); - _acb_poly_rising_ui_series(t, f, FLINT_MIN(2, len), r, rflen, prec); - _acb_poly_log_series(t, t, rflen, len, prec); - - _acb_log_rising_correct_branch(t, t, x, r, prec); - - acb_clear(f); - acb_clear(f + 1); -} - void _acb_poly_lgamma_series(acb_ptr res, acb_srcptr h, slong hlen, slong len, slong prec) { @@ -95,7 +73,7 @@ _acb_poly_lgamma_series(acb_ptr res, acb_srcptr h, slong hlen, slong len, slong acb_init(zr); /* use Stirling series */ - acb_gamma_stirling_choose_param(&reflect, &r, &n, h, 1, 0, wp); + acb_hypgeom_gamma_stirling_choose_param(&reflect, &r, &n, h, 1, 0, wp); if (reflect) { @@ -104,7 +82,7 @@ _acb_poly_lgamma_series(acb_ptr res, acb_srcptr h, slong hlen, slong len, slong { acb_sub_ui(u, h, 1, wp); acb_neg(u, u); - _log_rising_ui_series(t, u, r, len, wp); + acb_hypgeom_log_rising_ui_jet(t, u, r, len, wp); for (i = 1; i < len; i += 2) acb_neg(t + i, t + i); } @@ -143,7 +121,7 @@ _acb_poly_lgamma_series(acb_ptr res, acb_srcptr h, slong hlen, slong len, slong if (r != 0) { - _log_rising_ui_series(t, h, r, len, wp); + acb_hypgeom_log_rising_ui_jet(t, h, r, len, wp); _acb_vec_sub(u, u, t, len, wp); } } diff --git a/acb_poly/rgamma_series.c b/acb_poly/rgamma_series.c index 34686434..844d8792 100644 --- a/acb_poly/rgamma_series.c +++ b/acb_poly/rgamma_series.c @@ -11,7 +11,7 @@ #include "acb_poly.h" -void acb_gamma_stirling_choose_param(int * reflect, slong * r, slong * n, +void acb_hypgeom_gamma_stirling_choose_param(int * reflect, slong * r, slong * n, const acb_t x, int use_reflect, int digamma, slong prec); void @@ -56,7 +56,7 @@ _acb_poly_rgamma_series(acb_ptr res, acb_srcptr h, slong hlen, slong len, slong acb_init(f + 1); /* otherwise use Stirling series */ - acb_gamma_stirling_choose_param(&reflect, &r, &n, h, 1, 0, wp); + acb_hypgeom_gamma_stirling_choose_param(&reflect, &r, &n, h, 1, 0, wp); /* rgamma(h) = (gamma(1-h+r) sin(pi h)) / (rf(1-h, r) * pi), h = h0 + t*/ if (reflect) diff --git a/arb.h b/arb.h index e84b2ad6..09abd0e3 100644 --- a/arb.h +++ b/arb.h @@ -567,15 +567,9 @@ void arb_digamma(arb_t y, const arb_t x, slong prec); void arb_zeta(arb_t z, const arb_t s, slong prec); void arb_hurwitz_zeta(arb_t z, const arb_t s, const arb_t a, slong prec); -void arb_rising_ui_bs(arb_t y, const arb_t x, ulong n, slong prec); -void arb_rising_ui_rs(arb_t y, const arb_t x, ulong n, ulong m, slong prec); -void arb_rising_ui_rec(arb_t y, const arb_t x, ulong n, slong prec); void arb_rising_ui(arb_t z, const arb_t x, ulong n, slong prec); void arb_rising_fmpq_ui(arb_t y, const fmpq_t x, ulong n, slong prec); void arb_rising(arb_t z, const arb_t x, const arb_t n, slong prec); - -void arb_rising2_ui_rs(arb_t u, arb_t v, const arb_t x, ulong n, ulong m, slong prec); -void arb_rising2_ui_bs(arb_t u, arb_t v, const arb_t x, ulong n, slong prec); void arb_rising2_ui(arb_t u, arb_t v, const arb_t x, ulong n, slong prec); void arb_log_ui_from_prev(arb_t s, ulong k, arb_t log_prev, ulong prev, slong prec); diff --git a/arb/digamma.c b/arb/digamma.c index 0d0774ce..e44d1626 100644 --- a/arb/digamma.c +++ b/arb/digamma.c @@ -12,7 +12,7 @@ #include "flint/arith.h" #include "arb.h" -void arb_gamma_stirling_choose_param(int * reflect, slong * r, slong * n, +void arb_hypgeom_gamma_stirling_choose_param(int * reflect, slong * r, slong * n, const arb_t x, int use_reflect, int digamma, slong prec); void arb_gamma_stirling_eval(arb_t s, const arb_t z, slong nterms, int digamma, slong prec); @@ -63,7 +63,7 @@ arb_digamma(arb_t y, const arb_t x, slong prec) wp = prec + FLINT_BIT_COUNT(prec); - arb_gamma_stirling_choose_param(&reflect, &r, &n, x, 1, 1, wp); + arb_hypgeom_gamma_stirling_choose_param(&reflect, &r, &n, x, 1, 1, wp); arb_init(t); arb_init(u); diff --git a/arb/gamma.c b/arb/gamma.c index 3391cdf0..467abb46 100644 --- a/arb/gamma.c +++ b/arb/gamma.c @@ -10,182 +10,10 @@ */ #include "arb.h" -#include "acb.h" +#include "arb_hypgeom.h" #include "bernoulli.h" -#include "hypgeom.h" -/* tuning factor */ -#define GAMMA_STIRLING_BETA 0.27 - -#define PI 3.1415926535897932385 - -static slong -choose_n(double log2z, double argz, int digamma, slong prec) -{ - double argf, boundn; - slong n; - - argf = 1.0 / cos(0.5 * argz); - argf = log(argf) * (1. / log(2)); - - for (n = 1; ; n++) - { - if (digamma) - boundn = bernoulli_bound_2exp_si(2*n) - (2*n)*log2z + (2*n+1)*argf; - else - boundn = bernoulli_bound_2exp_si(2*n) - (2*n-1)*log2z + (2*n)*argf; - - /* success */ - if (boundn <= -prec) - return n; - - /* if the term magnitude does not decrease, r is too small */ - if (boundn > 1) - { - flint_printf("exception: gamma_stirling_choose_param failed to converge\n"); - flint_abort(); - } - } -} - -static void -choose_small(int * reflect, slong * r, slong * n, - double x, double y, int use_reflect, int digamma, slong prec) -{ - double w, argz, log2z; - slong rr; - - /* use reflection formula if very negative */ - if (x < -5.0 && use_reflect) - { - *reflect = 1; - x = 1.0 - x; - } - else - { - *reflect = 0; - } - - /* argument reduction until |z| >= w */ - w = FLINT_MAX(1.0, GAMMA_STIRLING_BETA * prec); - - rr = 0; - while (x < 1.0 || x*x + y*y < w*w) - { - x++; - rr++; - } - - log2z = 0.5 * log(x*x + y*y) * 1.44269504088896341; - argz = atan2(y, x); - - *r = rr; - *n = choose_n(log2z, argz, digamma, prec); -} - -static void -choose_large(int * reflect, slong * r, slong * n, - const arf_t a, const arf_t b, int use_reflect, int digamma, slong prec) -{ - if (use_reflect && arf_sgn(a) < 0) - *reflect = 1; - else - *reflect = 0; - - *r = 0; - - /* so big that we will certainly have n = 0 */ - if (arf_cmpabs_2exp_si(a, WORD_MAX / 8) >= 0 || - arf_cmpabs_2exp_si(b, WORD_MAX / 8) >= 0) - { - *n = 0; - } - else - { - slong ab, bb; - double log2z, argz; - - ab = arf_abs_bound_lt_2exp_si(a); - bb = arf_abs_bound_lt_2exp_si(b); - - log2z = FLINT_MAX(ab, bb); - - /* piecewise approximation of the argument */ - if (arf_is_zero(b)) - { - if ((arf_sgn(a) < 0) && !(*reflect)) - argz = PI; - else - argz = 0.0; - } - else - { - if ((arf_sgn(a) < 0) && !(*reflect)) - if (arf_cmpabs(a, b) <= 0) - argz = PI * 0.75; - else - argz = PI; - else - if (arf_cmpabs(a, b) <= 0) - argz = PI * 0.25; - else - argz = PI * 0.5; - } - - if (argz == PI) - *n = 0; - else - *n = choose_n(log2z, argz, digamma, prec); - } -} - - -void -acb_gamma_stirling_choose_param(int * reflect, slong * r, slong * n, - const acb_t z, int use_reflect, int digamma, slong prec) -{ - const arf_struct * a = arb_midref(acb_realref(z)); - const arf_struct * b = arb_midref(acb_imagref(z)); - - if (!arf_is_finite(a) || !arf_is_finite(b)) - { - *reflect = *r = *n = 0; - } - else if (arf_cmpabs_2exp_si(a, 40) > 0 || arf_cmpabs_2exp_si(b, 40) > 0) - { - choose_large(reflect, r, n, a, b, use_reflect, digamma, prec); - } - else - { - choose_small(reflect, r, n, - arf_get_d(a, ARF_RND_UP), - arf_get_d(b, ARF_RND_UP), use_reflect, digamma, prec); - } -} - -void -arb_gamma_stirling_choose_param(int * reflect, slong * r, slong * n, - const arb_t x, int use_reflect, int digamma, slong prec) -{ - const arf_struct * a = arb_midref(x); - - if (arf_is_inf(a) || arf_is_nan(a)) - { - *reflect = *r = *n = 0; - } - else if (arf_cmpabs_2exp_si(a, 40) > 0) - { - arf_t b; - arf_init(b); - choose_large(reflect, r, n, a, b, use_reflect, digamma, prec); - arf_clear(b); - } - else - { - choose_small(reflect, r, n, - arf_get_d(a, ARF_RND_UP), 0.0, use_reflect, digamma, prec); - } -} +/* todo: move/cleanup helper functions */ void acb_gamma_bound_phase(mag_t bound, const acb_t z) @@ -349,6 +177,8 @@ arb_gamma_stirling_coeff(arb_t b, ulong k, int digamma, slong prec) fmpz_clear(d); } + + void arb_gamma_stirling_eval(arb_t s, const arb_t z, slong nterms, int digamma, slong prec) { @@ -436,791 +266,33 @@ arb_gamma_stirling_eval(arb_t s, const arb_t z, slong nterms, int digamma, slong arb_clear(logz); } -void -arb_gamma_fmpq_stirling(arb_t y, const fmpq_t a, slong prec) -{ - int reflect; - slong r, n, wp; - arb_t x, t, u, v; - - wp = prec + FLINT_BIT_COUNT(prec); - - arb_init(x); - arb_init(t); - arb_init(u); - arb_init(v); - - arb_set_fmpq(x, a, wp); - arb_gamma_stirling_choose_param(&reflect, &r, &n, x, 1, 0, wp); - - if (reflect) - { - /* gamma(x) = (rf(1-x, r) * pi) / (gamma(1-x+r) sin(pi x)) */ - fmpq_t b; - fmpq_init(b); - fmpz_sub(fmpq_numref(b), fmpq_denref(a), fmpq_numref(a)); - fmpz_set(fmpq_denref(b), fmpq_denref(a)); - arb_rising_fmpq_ui(u, b, r, wp); - fmpq_clear(b); - arb_const_pi(v, wp); - arb_mul(u, u, v, wp); - arb_sub_ui(t, x, 1, wp); - arb_neg(t, t); - arb_add_ui(t, t, r, wp); - arb_gamma_stirling_eval(v, t, n, 0, wp); - arb_exp(v, v, wp); - arb_sin_pi_fmpq(t, a, wp); - arb_mul(v, v, t, wp); - } - else - { - /* gamma(x) = gamma(x+r) / rf(x,r) */ - arb_add_ui(t, x, r, wp); - arb_gamma_stirling_eval(u, t, n, 0, wp); - arb_exp(u, u, prec); - arb_rising_fmpq_ui(v, a, r, wp); - } - - arb_div(y, u, v, prec); - - arb_clear(t); - arb_clear(u); - arb_clear(v); - arb_clear(x); -} - -void -arb_gamma_const_1_3_eval(arb_t s, slong prec) -{ - hypgeom_t series; - arb_t t, u; - slong wp = prec + 4 + 2 * FLINT_BIT_COUNT(prec); - - arb_init(t); - arb_init(u); - - hypgeom_init(series); - - fmpz_poly_set_str(series->A, "1 1"); - fmpz_poly_set_str(series->B, "1 1"); - fmpz_poly_set_str(series->P, "4 5 -46 108 -72"); - fmpz_poly_set_str(series->Q, "4 0 0 0 512000"); - - prec += FLINT_CLOG2(prec); - - arb_hypgeom_infsum(s, t, series, wp, wp); - - arb_sqrt_ui(u, 10, wp); - arb_mul(t, t, u, wp); - - arb_const_pi(u, wp); - arb_pow_ui(u, u, 4, wp); - arb_mul_ui(u, u, 12, wp); - arb_mul(s, s, u, wp); - - arb_div(s, s, t, wp); - arb_root_ui(s, s, 2, wp); - arb_root_ui(s, s, 3, prec); - - hypgeom_clear(series); - arb_clear(t); - arb_clear(u); -} - -ARB_DEF_CACHED_CONSTANT(arb_gamma_const_1_3, arb_gamma_const_1_3_eval) - -void -arb_gamma_const_1_4_eval(arb_t x, slong prec) -{ - arb_t t, u; - slong wp = prec + 4 + 2 * FLINT_BIT_COUNT(prec); - - arb_init(t); - arb_init(u); - - arb_one(t); - arb_sqrt_ui(u, 2, wp); - arb_agm(x, t, u, wp); - - arb_const_pi(t, wp); - arb_mul_2exp_si(t, t, 1); - arb_sqrt(u, t, wp); - arb_mul(u, u, t, wp); - - arb_div(x, u, x, wp); - arb_sqrt(x, x, wp); - - arb_clear(t); - arb_clear(u); -} - -ARB_DEF_CACHED_CONSTANT(arb_gamma_const_1_4, arb_gamma_const_1_4_eval) - -void -arb_gamma_small_frac(arb_t y, unsigned int p, unsigned int q, slong prec) -{ - slong wp = prec + 4; - - if (q == 1) - { - arb_one(y); - } - else if (q == 2) /* p = 1 */ - { - arb_const_sqrt_pi(y, prec); - } - else if (q == 3) - { - if (p == 1) - { - arb_gamma_const_1_3(y, prec); - } - else /* p = 2 */ - { - arb_t t; - arb_init(t); - arb_gamma_const_1_3(y, wp); - arb_sqrt_ui(t, 3, wp); - arb_mul(y, y, t, wp); - arb_const_pi(t, wp); - arb_div(y, t, y, prec); - arb_mul_2exp_si(y, y, 1); - arb_clear(t); - } - } - else if (q == 4) - { - if (p == 1) - { - arb_gamma_const_1_4(y, prec); - } - else /* p = 3 */ - { - arb_t t; - arb_init(t); - arb_sqrt_ui(y, 2, wp); - arb_const_pi(t, wp); - arb_mul(y, y, t, wp); - arb_gamma_const_1_4(t, wp); - arb_div(y, y, t, prec); - arb_clear(t); - } - } - else if (q == 6) - { - arb_t t; - arb_init(t); - arb_const_pi(t, wp); - arb_div_ui(t, t, 3, wp); - arb_sqrt(t, t, wp); - arb_set_ui(y, 2); - arb_root_ui(y, y, 3, wp); - arb_mul(t, t, y, wp); - arb_gamma_const_1_3(y, wp); - arb_mul(y, y, y, prec); - - if (p == 1) - { - arb_div(y, y, t, prec); - } - else /* p = 5 */ - { - arb_div(y, t, y, wp); - arb_const_pi(t, wp); - arb_mul(y, y, t, prec); - arb_mul_2exp_si(y, y, 1); - } - - arb_clear(t); - } - else - { - flint_printf("small fraction not implemented!\n"); - flint_abort(); - } -} - -slong _arb_compute_bs_exponents(slong * tab, slong n); -slong _arb_get_exp_pos(const slong * tab, slong step); - -static void -bsplit2(arb_t P, arb_t Q, const fmpz_t zp, const fmpz_t zq, - const slong * xexp, arb_srcptr xpow, - ulong N, slong a, slong b, int cont, slong prec) -{ - if (b - a == 1) - { - fmpz_t t; - fmpz_init(t); - fmpz_set(t, zp); - fmpz_addmul_ui(t, zq, a + 1); - arb_set_fmpz(P, t); - arb_set(Q, P); - fmpz_clear(t); - } - else - { - arb_t Pb, Qb; - slong step, i, m; - - arb_init(Pb); - arb_init(Qb); - - step = (b - a) / 2; - m = a + step; - - bsplit2(P, Q, zp, zq, xexp, xpow, N, a, m, 1, prec); - bsplit2(Pb, Qb, zp, zq, xexp, xpow, N, m, b, 1, prec); - - arb_mul(P, P, Pb, prec); - arb_mul(Q, Q, Pb, prec); - - i = (step == 1) ? 0 : _arb_get_exp_pos(xexp, step); - arb_addmul(Q, Qb, xpow + i, prec); - - arb_clear(Pb); - arb_clear(Qb); - } -} - -static void -bsplit3(arb_t P, arb_t Q, const fmpz_t zp, const fmpz_t zq, - const slong * xexp, arb_srcptr xpow, - ulong N, slong a, slong b, int cont, slong prec) -{ - if (b - a == 1) - { - fmpz_t t; - fmpz_init(t); - arb_set(P, xpow + 0); /* N zq */ - fmpz_set(t, zp); - fmpz_submul_ui(t, zq, a + 1); /* zp - (a + 1) zq */ - arb_set_fmpz(Q, t); - fmpz_clear(t); - } - else - { - arb_t Pb, Qb; - slong step, i, m; - - arb_init(Pb); - arb_init(Qb); - - step = (b - a) / 2; - m = a + step; - - bsplit3(P, Q, zp, zq, xexp, xpow, N, a, m, 1, prec); - bsplit3(Pb, Qb, zp, zq, xexp, xpow, N, m, b, 1, prec); - - i = _arb_get_exp_pos(xexp, m - a); - - arb_mul(P, P, xpow + i, prec); - if (b - m != m - a) - arb_mul(P, P, xpow + 0, prec); - - arb_addmul(P, Q, Pb, prec); - - if (cont) - { - arb_mul(Q, Q, Qb, prec); - } - else - { - i = _arb_get_exp_pos(xexp, m - a); - - arb_mul(Q, xpow + i, xpow + i, prec); - if (b - m != m - a) - arb_mul(Q, Q, xpow + 0, prec); - } - - arb_clear(Pb); - arb_clear(Qb); - } -} - -double d_lambertw_branch1(double x); - -static ulong -more_trailing_zeros(ulong N) -{ - ulong bc, N2; - - bc = FLINT_BIT_COUNT(N); - - if (bc >= 8) - { - N2 = (N >> (bc - 5)) << (bc - 5); - N2 += ((N2 != N) << (bc - 5)); - return N2; - } - else - { - return N; - } -} - -#define C_LOG2 0.69314718055994530942 -#define C_EXP1 2.7182818284590452354 - -static void -build_bsplit_power_table(arb_ptr xpow, const slong * xexp, slong len, slong prec) -{ - slong i; - - for (i = 1; i < len; i++) - { - if (xexp[i] == 2 * xexp[i-1]) - { - arb_mul(xpow + i, xpow + i - 1, xpow + i - 1, prec); - } - else if (xexp[i] == 2 * xexp[i-2]) /* prefer squaring if possible */ - { - arb_mul(xpow + i, xpow + i - 2, xpow + i - 2, prec); - } - else if (xexp[i] == 2 * xexp[i-1] + 1) - { - arb_mul(xpow + i, xpow + i - 1, xpow + i - 1, prec); - arb_mul(xpow + i, xpow + i, xpow, prec); - } - else if (xexp[i] == 2 * xexp[i-2] + 1) - { - arb_mul(xpow + i, xpow + i - 2, xpow + i - 2, prec); - arb_mul(xpow + i, xpow + i, xpow, prec); - } - else - { - flint_printf("power table has the wrong structure!\n"); - flint_abort(); - } - } -} - -/* assumes z in [1, 2] */ -static void -arb_gamma_fmpq_general_off1(arb_t res, const fmpq_t z, slong prec) -{ - slong wp, N, n, n2, length, length2, wp2; - double alpha; - arb_t P, Q; - slong *xexp, *xexp2; - arb_ptr xpow; - mag_t err, err2; - - wp = prec + 30; - - alpha = 0.52; /* tuning parameter between 0.5 and 1 */ - - N = alpha * C_LOG2 * wp; - N = more_trailing_zeros(N); - alpha = N / (C_LOG2 * wp); - - /* Terms in convergent series */ - n = (1 - alpha) / d_lambertw((1 - alpha) / (alpha * C_EXP1)) * C_LOG2 * wp; - - /* Precision and terms in asymptotic series */ - wp2 = wp * (1 - alpha); - wp2 = FLINT_MAX(wp2, 30); - n2 = (alpha - 1) / d_lambertw_branch1((alpha - 1) / (alpha * C_EXP1)) * C_LOG2 * wp; - n2 = FLINT_MAX(n2, 2); /* binary splitting correctness */ - - mag_init(err); - mag_init(err2); - arb_init(P); - arb_init(Q); - - /* compute the powers of x = N*zq that will appear (at least x^1) */ - xexp = flint_calloc(2 * FLINT_BITS, sizeof(slong)); - xexp2 = flint_calloc(2 * FLINT_BITS, sizeof(slong)); - - length = _arb_compute_bs_exponents(xexp, n); - length2 = _arb_compute_bs_exponents(xexp2, n2); - - xpow = _arb_vec_init(FLINT_MAX(length, length2)); - - arb_set_fmpz(xpow + 0, fmpq_denref(z)); - arb_mul_ui(xpow + 0, xpow + 0, N, wp); - - build_bsplit_power_table(xpow, xexp, length, wp); - - /* 1F1(1, 1+z, N) */ - bsplit2(P, Q, fmpq_numref(z), fmpq_denref(z), xexp, xpow, N, 0, n, 0, wp); - arb_div(P, Q, P, wp); - - /* Convergent series error bound: N^n / n! (1 + (N/n) + ...) */ - mag_set_ui(err, N); - mag_pow_ui(err, err, n); - mag_rfac_ui(err2, n); - mag_mul(err, err, err2); - mag_set_ui(err2, N); - mag_div_ui(err2, err2, n); - mag_geom_series(err2, err2, 0); - mag_mul(err, err, err2); - arb_add_error_mag(P, err); - - /* divide 1F1 by z */ - arb_mul_fmpz(P, P, fmpq_denref(z), wp); - arb_div_fmpz(P, P, fmpq_numref(z), wp); - arb_swap(res, P); - - build_bsplit_power_table(xpow, xexp2, length2, wp2); - - bsplit3(P, Q, fmpq_numref(z), fmpq_denref(z), xexp2, xpow, N, 0, n2, 0, wp2); - arb_div(P, P, Q, wp2); - - /* 2F0 error bound (bounded by first omitted term) */ - mag_fac_ui(err, n2); - mag_set_ui_lower(err2, N); - mag_pow_ui_lower(err2, err2, n2); - mag_div(err, err, err2); - arb_add_error_mag(P, err); - - /* N^z * exp(-N) * (1F1/z + 2F0/N) */ - arb_div_ui(P, P, N, wp2); - - arb_add(res, res, P, wp); - arb_set_ui(Q, N); - arb_pow_fmpq(Q, Q, z, wp); - arb_mul(res, res, Q, wp); - - arb_set_si(Q, -N); - arb_exp(Q, Q, wp); - arb_mul(res, res, Q, wp); - - _arb_vec_clear(xpow, FLINT_MAX(length, length2)); - flint_free(xexp); - flint_free(xexp2); - - arb_clear(P); - arb_clear(Q); - mag_clear(err); - mag_clear(err2); -} - -/* assumes z in (0, 1] */ -void -arb_gamma_fmpq_hyp(arb_t res, const fmpq_t z, slong prec) -{ - fmpq_t t; - fmpq_init(t); - fmpq_add_ui(t, z, 1); - arb_gamma_fmpq_general_off1(res, t, prec); - arb_mul_fmpz(res, res, fmpq_denref(z), prec + 30); - arb_div_fmpz(res, res, fmpq_numref(z), prec); - fmpq_clear(t); -} - -void -arb_gamma_fmpq_outward(arb_t y, const fmpq_t x, slong prec) -{ - fmpq_t a; - fmpz_t n; - fmpz p, q; - slong m; - arb_t t, u; - - fmpq_init(a); - fmpz_init(n); - arb_init(t); - arb_init(u); - - /* write x = a + n with 0 < a <= 1 */ - if (fmpz_is_one(fmpq_denref(x))) - { - fmpq_one(a); - fmpz_sub_ui(n, fmpq_numref(x), 1); - } - else - { - fmpz_fdiv_qr(n, fmpq_numref(a), fmpq_numref(x), fmpq_denref(x)); - fmpz_set(fmpq_denref(a), fmpq_denref(x)); - } - - if (!fmpz_fits_si(n)) - { - flint_printf("gamma: too large fmpq to reduce to 0!\n"); - flint_abort(); - } - - m = fmpz_get_si(n); - - /* evaluate gamma(a) */ - p = *fmpq_numref(a); - q = *fmpq_denref(a); - - if (q == 1 || q == 2 || q == 3 || q == 4 || q == 6) - { - arb_gamma_small_frac(t, p, q, prec); - } - else - { - arb_gamma_fmpq_hyp(t, a, prec); - } - - /* argument reduction */ - if (m >= 0) - { - arb_rising_fmpq_ui(u, a, m, prec); - arb_mul(y, t, u, prec); - } - else - { - arb_rising_fmpq_ui(u, x, -m, prec); - arb_div(y, t, u, prec); - } - - fmpq_clear(a); - fmpz_clear(n); - arb_clear(t); - arb_clear(u); -} - void arb_gamma_fmpq(arb_t y, const fmpq_t x, slong prec) { - fmpz p, q; - - p = *fmpq_numref(x); - q = *fmpq_denref(x); - - if ((q == 1 || q == 2 || q == 3 || q == 4 || q == 6) && !COEFF_IS_MPZ(p)) - { - if (q == 1) - { - if (p <= 0) - { - arb_indeterminate(y); - return; - } - - if (p < 1200 || 1.44265 * (p*log(p) - p) < 15.0 * prec) - { - fmpz_t t; - fmpz_init(t); - fmpz_fac_ui(t, p - 1); - arb_set_round_fmpz(y, t, prec); - fmpz_clear(t); - return; - } - } - - p = FLINT_ABS(p); - - if (p < q * 500.0 || p < q * (500.0 + 0.1 * prec * sqrt(prec))) - { - arb_gamma_fmpq_outward(y, x, prec); - return; - } - } - - if (q != 1 && prec > 7000 + 300 * fmpz_bits(fmpq_denref(x)) && - (slong) fmpz_bits(fmpq_numref(x)) - (slong) fmpz_bits(fmpq_denref(x)) < FLINT_BITS && - fabs(fmpq_get_d(x)) < 0.03 * prec * sqrt(prec)) - { - arb_gamma_fmpq_outward(y, x, prec); - return; - } - - arb_gamma_fmpq_stirling(y, x, prec); + arb_hypgeom_gamma_fmpq(y, x, prec); } void arb_gamma_fmpz(arb_t y, const fmpz_t x, slong prec) { - fmpq_t t; - *fmpq_numref(t) = *x; - *fmpq_denref(t) = WORD(1); - arb_gamma_fmpq(y, t, prec); -} - -static void -_arb_gamma(arb_t y, const arb_t x, slong prec, int inverse) -{ - int reflect; - slong r, n, wp; - arb_t t, u, v; - double acc; - - if (arb_is_exact(x)) - { - const arf_struct * mid = arb_midref(x); - - if (arf_is_special(mid)) - { - if (!inverse && arf_is_pos_inf(mid)) - arb_set(y, x); - else if (arf_is_nan(mid) || arf_is_neg_inf(mid) || !inverse) - arb_indeterminate(y); - else - arb_zero(y); - return; - } - else if (inverse && arf_is_int(mid) && arf_sgn(mid) < 0) - { - arb_zero(y); - } - else - { - /* fast gamma(n), gamma(n/2) or gamma(n/4), ... */ - if (arf_cmpabs_2exp_si(mid, prec) < 0 && - (arf_is_int_2exp_si(mid, -2) || (prec > 1000 && arf_is_int_2exp_si(mid, -1000)))) - { - fmpq_t a; - fmpq_init(a); - arf_get_fmpq(a, mid); - arb_gamma_fmpq(y, a, prec + 2 * inverse); - if (inverse) - arb_inv(y, y, prec); - fmpq_clear(a); - return; - } - } - } - - /* todo: for large x (if exact or accurate enough), increase precision */ - acc = arb_rel_accuracy_bits(x); - acc = FLINT_MAX(acc, 0); - wp = FLINT_MIN(prec, acc + 20); - wp = FLINT_MAX(wp, 2); - wp = wp + FLINT_BIT_COUNT(wp); - - if (acc < 3) /* try to avoid divisions blowing up */ - { - if (arf_cmp_d(arb_midref(x), -0.5) < 0) - { - reflect = 1; - r = 0; - } - else if (arf_cmp_si(arb_midref(x), 1) < 0) - { - reflect = 0; - r = 1; - } - else - { - reflect = 0; - r = 0; - } - - n = 1; - } - else - { - arb_gamma_stirling_choose_param(&reflect, &r, &n, x, 1, 0, wp); - } - - arb_init(t); - arb_init(u); - arb_init(v); - - if (reflect) - { - arb_sub_ui(t, x, 1, wp); - arb_neg(t, t); - arb_rising_ui_rec(u, t, r, wp); - arb_const_pi(v, wp); - arb_mul(u, u, v, wp); - arb_add_ui(t, t, r, wp); - arb_gamma_stirling_eval(v, t, n, 0, wp); - - if (inverse) - { - /* rgamma(x) = gamma(1-x+r) sin(pi x) / ((rf(1-x, r) * pi) */ - arb_exp(v, v, wp); - arb_sin_pi(t, x, wp); - arb_mul(v, v, t, wp); - arb_mul(y, u, v, wp); - arb_div(y, v, u, prec); - } - else - { - /* gamma(x) = (rf(1-x, r) * pi) rgamma(1-x+r) csc(pi x) */ - arb_neg(v, v); - arb_exp(v, v, wp); - arb_csc_pi(t, x, wp); - arb_mul(v, v, t, wp); - arb_mul(y, v, u, prec); - } - } - else - { - arb_add_ui(t, x, r, wp); - arb_gamma_stirling_eval(u, t, n, 0, wp); - - if (inverse) - { - /* rgamma(x) = rf(x,r) rgamma(x+r) */ - arb_neg(u, u); - arb_exp(u, u, prec); - arb_rising_ui_rec(v, x, r, wp); - arb_mul(y, v, u, prec); - } - else - { - /* gamma(x) = gamma(x+r) / rf(x,r) */ - arb_exp(u, u, prec); - arb_rising_ui_rec(v, x, r, wp); - arb_div(y, u, v, prec); - } - } - - arb_clear(t); - arb_clear(u); - arb_clear(v); + arb_hypgeom_gamma_fmpz(y, x, prec); } void arb_gamma(arb_t y, const arb_t x, slong prec) { - _arb_gamma(y, x, prec, 0); + arb_hypgeom_gamma(y, x, prec); } void arb_rgamma(arb_t y, const arb_t x, slong prec) { - _arb_gamma(y, x, prec, 1); + arb_hypgeom_rgamma(y, x, prec); } void arb_lgamma(arb_t y, const arb_t x, slong prec) { - int reflect; - slong r, n, wp; - arb_t t, u; - - if (!arb_is_positive(x)) - { - arb_indeterminate(y); - return; - } - - /* fast gamma(n), gamma(n/2) or gamma(n/4), ... */ - if (arb_is_exact(x) && arf_cmpabs_2exp_si(arb_midref(x), prec) < 0 && - (arf_is_int_2exp_si(arb_midref(x), -2) || (prec > 10000 && arf_is_int_2exp_si(arb_midref(x), -1000)))) - { - fmpq_t a; - fmpq_init(a); - arf_get_fmpq(a, arb_midref(x)); - arb_gamma_fmpq(y, a, prec); - arb_log(y, y, prec); - fmpq_clear(a); - return; - } - - wp = prec + FLINT_BIT_COUNT(prec); - - arb_gamma_stirling_choose_param(&reflect, &r, &n, x, 0, 0, wp); - - /* log(gamma(x)) = log(gamma(x+r)) - log(rf(x,r)) */ - arb_init(t); - arb_init(u); - - arb_add_ui(t, x, r, wp); - arb_gamma_stirling_eval(u, t, n, 0, wp); - arb_rising_ui_rec(t, x, r, wp); - arb_log(t, t, wp); - arb_sub(y, u, t, prec); - - arb_clear(t); - arb_clear(u); + arb_hypgeom_lgamma(y, x, prec); } diff --git a/arb/rising2_ui.c b/arb/rising2_ui.c index 2dd66fc9..683ecc30 100644 --- a/arb/rising2_ui.c +++ b/arb/rising2_ui.c @@ -10,13 +10,30 @@ */ #include "arb.h" +#include "arb_hypgeom.h" void arb_rising2_ui(arb_t u, arb_t v, const arb_t x, ulong n, slong prec) { - if (prec < 512 || n < 8 || arb_bits(x) < prec / 8) - arb_rising2_ui_bs(u, v, x, n, prec); + if (x == u || x == v) + { + arb_t t; + arb_init(t); + arb_set(t, x); + arb_rising2_ui(u, v, t, n, prec); + arb_clear(t); + } else - arb_rising2_ui_rs(u, v, x, n, 0, prec); + { + arb_struct tmp[2]; + + tmp[0] = *u; + tmp[1] = *v; + + arb_hypgeom_rising_ui_jet(tmp, x, n, 2, prec); + + *u = tmp[0]; + *v = tmp[1]; + } } diff --git a/arb/rising2_ui_bs.c b/arb/rising2_ui_bs.c deleted file mode 100644 index a815c138..00000000 --- a/arb/rising2_ui_bs.c +++ /dev/null @@ -1,83 +0,0 @@ -/* - Copyright (C) 2013 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.h" - -static void -bsplit(arb_t p, arb_t q, const arb_t x, ulong a, ulong b, slong prec) -{ - if (b - a < 8) - { - ulong k; - arb_t t; - - arb_one(p); - arb_add_ui(q, x, a, prec); - - arb_init(t); - - for (k = a + 1; k < b; k++) - { - arb_add_ui(t, x, k, prec); - arb_mul(p, p, t, prec); - arb_add(p, p, q, prec); - arb_mul(q, q, t, prec); - } - - arb_clear(t); - } - else - { - arb_t r, s; - ulong m; - - arb_init(r); - arb_init(s); - - m = a + (b - a) / 2; - bsplit(p, q, x, a, m, prec); - bsplit(r, s, x, m, b, prec); - - arb_mul(p, p, s, prec); - arb_mul(r, r, q, prec); - arb_add(p, p, r, prec); - arb_mul(q, q, s, prec); - - arb_clear(r); - arb_clear(s); - } -} - -void -arb_rising2_ui_bs(arb_t u, arb_t v, const arb_t x, ulong n, slong prec) -{ - if (n == 0) - { - arb_zero(v); - arb_one(u); - } - else if (n == 1) - { - arb_set(u, x); - arb_one(v); - } - else - { - arb_t t; - slong wp = ARF_PREC_ADD(prec, FLINT_BIT_COUNT(n)); - - arb_init(t); /* support aliasing */ - arb_set(t, x); - bsplit(v, u, t, 0, n, wp); - arb_clear(t); - } -} - diff --git a/arb/rising2_ui_rs.c b/arb/rising2_ui_rs.c deleted file mode 100644 index 17fe2d0a..00000000 --- a/arb/rising2_ui_rs.c +++ /dev/null @@ -1,162 +0,0 @@ -/* - Copyright (C) 2013 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 "flint/fmpz_poly.h" -#include "arb.h" - -void -_gamma_rf_bsplit(fmpz * A, ulong a, ulong b) -{ - ulong n = b - a; - - if (n == 0) - { - fmpz_one(A); - } - else if (n < 8) - { - ulong j, k; - - fmpz_set_ui(A, a); - fmpz_one(A + 1); - - for (j = 1; j < n; j++) - { - fmpz_one(A + j + 1); - - for (k = j; k > 0; k--) - { - fmpz_mul_ui(A + k, A + k, a + j); - fmpz_add(A + k, A + k, A + k - 1); - } - - fmpz_mul_ui(A, A, a + j); - } - } - else - { - ulong m = a + (b - a) / 2; - ulong w = m - a; - ulong v = b - m; - - fmpz *t, *A1, *A2; - - t = _fmpz_vec_init(w + v + 2); - - A1 = t; - A2 = A1 + w + 1; - - _gamma_rf_bsplit(A1, a, m); - _gamma_rf_bsplit(A2, m, b); - - _fmpz_poly_mul(A, A2, v + 1, A1, w + 1); - - _fmpz_vec_clear(t, w + v + 2); - } -} - -void -arb_rising2_ui_rs(arb_t u, arb_t v, - const arb_t x, ulong n, ulong m, slong prec) -{ - if (n == 0) - { - arb_zero(v); - arb_one(u); - } - else if (n == 1) - { - arb_set(u, x); - arb_one(v); - } - else - { - slong wp; - ulong i, j, a, b; - arb_ptr xs; - arb_t S, T, U, V; - fmpz *A, *B; - - wp = ARF_PREC_ADD(prec, FLINT_BIT_COUNT(n)); - - if (m == 0) - { - ulong m1, m2; - m1 = 0.6 * pow(wp, 0.4); - m2 = n_sqrt(n); - m = FLINT_MIN(m1, m2); - } - - m = FLINT_MAX(m, 1); - - xs = _arb_vec_init(m + 1); - A = _fmpz_vec_init(2 * m + 1); - B = A + (m + 1); - - arb_init(S); - arb_init(T); - arb_init(U); - arb_init(V); - _arb_vec_set_powers(xs, x, m + 1, wp); - - for (i = 0; i < n; i += m) - { - a = i; - b = FLINT_MIN(n, a + m); - - if (a == 0 || b != a + m) - { - _gamma_rf_bsplit(A, a, b); - } - else - { - fmpz tt = m; - _fmpz_poly_taylor_shift(A, &tt, m + 1); - } - - _fmpz_poly_derivative(B, A, b - a + 1); - - arb_set_fmpz(S, A); - - for (j = 1; j <= b - a; j++) - arb_addmul_fmpz(S, xs + j, A + j, wp); - - arb_set_fmpz(T, B); - - for (j = 1; j < b - a; j++) - arb_addmul_fmpz(T, xs + j, B + j, wp); - - if (i == 0) - { - arb_set(U, S); - arb_set(V, T); - } - else - { - arb_mul(V, V, S, wp); - arb_addmul(V, U, T, wp); - arb_mul(U, U, S, wp); - } - } - - arb_set(u, U); - arb_set(v, V); - - _arb_vec_clear(xs, m + 1); - _fmpz_vec_clear(A, 2 * m + 1); - - arb_clear(S); - arb_clear(T); - arb_clear(U); - arb_clear(V); - } -} - diff --git a/arb/rising_ui.c b/arb/rising_ui.c index 429144d0..7b7c9606 100644 --- a/arb/rising_ui.c +++ b/arb/rising_ui.c @@ -10,44 +10,17 @@ */ #include "arb.h" +#include "arb_hypgeom.h" void arb_rising_ui(arb_t y, const arb_t x, ulong n, slong prec) { - if (n < FLINT_MAX(prec, 100)) - { - arb_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); - } + arb_hypgeom_rising_ui(y, x, n, prec); } void arb_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_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(y, x, n, prec); } diff --git a/arb/rising_ui_bs.c b/arb/rising_ui_bs.c deleted file mode 100644 index 6c1fe90c..00000000 --- a/arb/rising_ui_bs.c +++ /dev/null @@ -1,80 +0,0 @@ -/* - 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.h" - -/* assumes y and x are not aliased */ -static void -bsplit(arb_t y, const arb_t x, ulong a, ulong b, slong prec) -{ - if (b - a == 1) - { - arb_set_round(y, x, prec); - } - else if (b - a <= 10) - { - slong i; - arb_t t; - arb_init(t); - - arb_add_ui(t, x, a, prec); - arb_add_ui(y, x, a + 1, prec); - arb_mul(y, y, t, prec); - - for (i = a + 2; i < b; i++) - { - arb_add_ui(t, x, i, prec); - arb_mul(y, y, t, prec); - } - - arb_clear(t); - } - 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_rising_ui_bs(arb_t y, const arb_t x, ulong n, slong prec) -{ - if (n == 0) - { - arb_one(y); - } - else if (n == 1) - { - arb_set_round(y, x, prec); - } - else - { - 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(y, t, prec); - arb_clear(t); - } -} - diff --git a/arb/rising_ui_rec.c b/arb/rising_ui_rec.c deleted file mode 100644 index 654f2253..00000000 --- a/arb/rising_ui_rec.c +++ /dev/null @@ -1,22 +0,0 @@ -/* - 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.h" - -void -arb_rising_ui_rec(arb_t y, const arb_t x, ulong n, slong prec) -{ - if (prec < 512 || n < 8 || arb_bits(x) < prec / 8) - arb_rising_ui_bs(y, x, n, prec); - else - arb_rising_ui_rs(y, x, n, 0, prec); -} - diff --git a/arb/rising_ui_rs.c b/arb/rising_ui_rs.c deleted file mode 100644 index 35737cf0..00000000 --- a/arb/rising_ui_rs.c +++ /dev/null @@ -1,154 +0,0 @@ -/* - 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 "flint/arith.h" -#include "arb.h" - -void -rising_difference_polynomial(fmpz * s, fmpz * c, ulong m) -{ - slong i, j, v; - - fmpz_t t; - fmpz_init(t); - - arith_stirling_number_1u_vec(s, m, m + 1); - - /* Compute the first row */ - for (i = 0; i < m; i++) - { - fmpz_set_ui(t, m); - fmpz_mul_ui(t, t, (i + 1)); - fmpz_mul(c + i, s + (i + 1), t); - - for (j = i + 2; j < m + 1; j++) - { - fmpz_mul_ui(t, t, m * j); - fmpz_divexact_ui(t, t, j - i); - fmpz_addmul(c + i, s + j, t); - } - } - - /* Extend using recurrence and symmetry */ - for (v = 1; v < m; v++) - { - for (i = v; i < m - v; i++) - { - fmpz_mul_ui(t, c + (v - 1) * m + (i + 1), i + 1); - fmpz_divexact_ui(c + v * m + i, t, v); - } - - for (i = 0; i < v; i++) - fmpz_set(c + v * m + i, c + i * m + v); - } - - fmpz_clear(t); -} - -void -arb_rising_ui_rs(arb_t y, const arb_t x, ulong n, ulong m, slong prec) -{ - arb_ptr xs; - arb_t t, u, v; - ulong i, k, rem; - fmpz_t c, h; - fmpz *s, *d; - slong wp; - - if (n == 0) - { - arb_one(y); - return; - } - - if (n == 1) - { - arb_set_round(y, x, prec); - return; - } - - wp = ARF_PREC_ADD(prec, FLINT_BIT_COUNT(n)); - - arb_init(t); - arb_init(u); - arb_init(v); - fmpz_init(c); - fmpz_init(h); - - if (m == 0) - { - ulong m1, m2; - m1 = 0.2 * pow(wp, 0.4); - m2 = n_sqrt(n); - m = FLINT_MIN(m1, m2); - } - - m = FLINT_MIN(m, n); - m = FLINT_MAX(m, 1); - - xs = _arb_vec_init(m + 1); - d = _fmpz_vec_init(m * m); - s = _fmpz_vec_init(m + 1); - - _arb_vec_set_powers(xs, x, m + 1, wp); - - rising_difference_polynomial(s, d, m); - - /* tail */ - rem = m; - while (rem + m <= n) - rem += m; - arb_one(y); - for (k = rem; k < n; k++) - { - arb_add_ui(t, xs + 1, k, wp); - arb_mul(y, y, t, wp); - } - - /* initial rising factorial */ - arb_zero(t); - for (i = 1; i <= m; i++) - arb_addmul_fmpz(t, xs + i, s + i, wp); - - arb_mul(y, y, t, wp); - - /* the leading coefficient is always the same */ - arb_mul_fmpz(xs + m - 1, xs + m - 1, d + m - 1 + 0, wp); - - for (k = 0; k + 2 * m <= n; k += m) - { - for (i = 0; i < m - 1; i++) - { - fmpz_set_ui(h, k); - _fmpz_poly_evaluate_horner_fmpz(c, d + i * m, m - i, h); - - if (i == 0) - arb_add_fmpz(t, t, c, wp); - else - arb_addmul_fmpz(t, xs + i, c, wp); - } - - arb_add(t, t, xs + m - 1, wp); - arb_mul(y, y, t, wp); - } - - arb_set_round(y, y, prec); - - arb_clear(t); - arb_clear(u); - arb_clear(v); - _arb_vec_clear(xs, m + 1); - _fmpz_vec_clear(d, m * m); - _fmpz_vec_clear(s, m + 1); - fmpz_clear(c); - fmpz_clear(h); -} - diff --git a/arb/test/t-rising2_ui_bs.c b/arb/test/t-rising2_ui_bs.c deleted file mode 100644 index f87bad81..00000000 --- a/arb/test/t-rising2_ui_bs.c +++ /dev/null @@ -1,110 +0,0 @@ -/* - Copyright (C) 2013 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 "flint/arith.h" -#include "arb_poly.h" - -int main() -{ - slong iter; - flint_rand_t state; - - flint_printf("rising2_ui_bs...."); - fflush(stdout); - - flint_randinit(state); - - for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++) - { - arb_t a, u, v, u2, v2; - fmpz *f; - arb_ptr g; - ulong n; - slong i, prec; - - arb_init(a); - arb_init(u); - arb_init(v); - arb_init(u2); - arb_init(v2); - - arb_randtest(a, state, 1 + n_randint(state, 4000), 10); - arb_randtest(u, state, 1 + n_randint(state, 4000), 10); - arb_randtest(v, state, 1 + n_randint(state, 4000), 10); - n = n_randint(state, 120); - - f = _fmpz_vec_init(n + 1); - g = _arb_vec_init(n + 1); - - prec = 2 + n_randint(state, 4000); - arb_rising2_ui_bs(u, v, a, n, prec); - - arith_stirling_number_1u_vec(f, n, n + 1); - for (i = 0; i <= n; i++) - arb_set_fmpz(g + i, f + i); - _arb_poly_evaluate(u2, g, n + 1, a, prec); - - _arb_poly_derivative(g, g, n + 1, prec); - _arb_poly_evaluate(v2, g, n, a, prec); - - if (!arb_overlaps(u, u2) || !arb_overlaps(v, v2)) - { - flint_printf("FAIL: overlap\n\n"); - flint_printf("n = %wu\n", n); - flint_printf("a = "); arb_printd(a, 15); flint_printf("\n\n"); - flint_printf("u = "); arb_printd(u, 15); flint_printf("\n\n"); - flint_printf("u2 = "); arb_printd(u2, 15); flint_printf("\n\n"); - flint_printf("v = "); arb_printd(v, 15); flint_printf("\n\n"); - flint_printf("v2 = "); arb_printd(v2, 15); flint_printf("\n\n"); - flint_abort(); - } - - arb_set(u2, a); - arb_rising2_ui_bs(u2, v, u2, n, prec); - - if (!arb_equal(u2, u)) - { - flint_printf("FAIL: aliasing 1\n\n"); - flint_printf("a = "); arb_printd(a, 15); flint_printf("\n\n"); - flint_printf("u = "); arb_printd(u, 15); flint_printf("\n\n"); - flint_printf("u2 = "); arb_printd(u2, 15); flint_printf("\n\n"); - flint_printf("n = %wu\n", n); - flint_abort(); - } - - arb_set(v2, a); - arb_rising2_ui_bs(u, v2, v2, n, prec); - - if (!arb_equal(v2, v)) - { - flint_printf("FAIL: aliasing 2\n\n"); - flint_printf("a = "); arb_printd(a, 15); flint_printf("\n\n"); - flint_printf("v = "); arb_printd(v, 15); flint_printf("\n\n"); - flint_printf("v2 = "); arb_printd(v2, 15); flint_printf("\n\n"); - flint_printf("n = %wu\n", n); - flint_abort(); - } - - arb_clear(a); - arb_clear(u); - arb_clear(v); - arb_clear(u2); - arb_clear(v2); - _fmpz_vec_clear(f, n + 1); - _arb_vec_clear(g, n + 1); - } - - flint_randclear(state); - flint_cleanup(); - flint_printf("PASS\n"); - return EXIT_SUCCESS; -} - diff --git a/arb/test/t-rising2_ui_rs.c b/arb/test/t-rising2_ui_rs.c deleted file mode 100644 index 8bd1b837..00000000 --- a/arb/test/t-rising2_ui_rs.c +++ /dev/null @@ -1,110 +0,0 @@ -/* - Copyright (C) 2013 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 "flint/arith.h" -#include "arb_poly.h" - -int main() -{ - slong iter; - flint_rand_t state; - - flint_printf("rising2_ui_rs...."); - fflush(stdout); - - flint_randinit(state); - - for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++) - { - arb_t a, u, v, u2, v2; - fmpz *f; - arb_ptr g; - ulong n; - slong i, prec; - - arb_init(a); - arb_init(u); - arb_init(v); - arb_init(u2); - arb_init(v2); - - arb_randtest(a, state, 1 + n_randint(state, 4000), 10); - arb_randtest(u, state, 1 + n_randint(state, 4000), 10); - arb_randtest(v, state, 1 + n_randint(state, 4000), 10); - n = n_randint(state, 120); - - f = _fmpz_vec_init(n + 1); - g = _arb_vec_init(n + 1); - - prec = 2 + n_randint(state, 4000); - arb_rising2_ui_rs(u, v, a, n, 0, prec); - - arith_stirling_number_1u_vec(f, n, n + 1); - for (i = 0; i <= n; i++) - arb_set_fmpz(g + i, f + i); - _arb_poly_evaluate(u2, g, n + 1, a, prec); - - _arb_poly_derivative(g, g, n + 1, prec); - _arb_poly_evaluate(v2, g, n, a, prec); - - if (!arb_overlaps(u, u2) || !arb_overlaps(v, v2)) - { - flint_printf("FAIL: overlap\n\n"); - flint_printf("n = %wu\n", n); - flint_printf("a = "); arb_printd(a, 15); flint_printf("\n\n"); - flint_printf("u = "); arb_printd(u, 15); flint_printf("\n\n"); - flint_printf("u2 = "); arb_printd(u2, 15); flint_printf("\n\n"); - flint_printf("v = "); arb_printd(v, 15); flint_printf("\n\n"); - flint_printf("v2 = "); arb_printd(v2, 15); flint_printf("\n\n"); - flint_abort(); - } - - arb_set(u2, a); - arb_rising2_ui_rs(u2, v, u2, n, 0, prec); - - if (!arb_equal(u2, u)) - { - flint_printf("FAIL: aliasing 1\n\n"); - flint_printf("a = "); arb_printd(a, 15); flint_printf("\n\n"); - flint_printf("u = "); arb_printd(u, 15); flint_printf("\n\n"); - flint_printf("u2 = "); arb_printd(u2, 15); flint_printf("\n\n"); - flint_printf("n = %wu\n", n); - flint_abort(); - } - - arb_set(v2, a); - arb_rising2_ui_rs(u, v2, v2, n, 0, prec); - - if (!arb_equal(v2, v)) - { - flint_printf("FAIL: aliasing 2\n\n"); - flint_printf("a = "); arb_printd(a, 15); flint_printf("\n\n"); - flint_printf("v = "); arb_printd(v, 15); flint_printf("\n\n"); - flint_printf("v2 = "); arb_printd(v2, 15); flint_printf("\n\n"); - flint_printf("n = %wu\n", n); - flint_abort(); - } - - arb_clear(a); - arb_clear(u); - arb_clear(v); - arb_clear(u2); - arb_clear(v2); - _fmpz_vec_clear(f, n + 1); - _arb_vec_clear(g, n + 1); - } - - flint_randclear(state); - flint_cleanup(); - flint_printf("PASS\n"); - return EXIT_SUCCESS; -} - diff --git a/arb/test/t-rising_ui_bs.c b/arb/test/t-rising_ui_bs.c deleted file mode 100644 index 5edf322e..00000000 --- a/arb/test/t-rising_ui_bs.c +++ /dev/null @@ -1,109 +0,0 @@ -/* - Copyright (C) 2012 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.h" - -int main() -{ - slong iter; - flint_rand_t state; - - flint_printf("rising_ui_bs...."); - fflush(stdout); - - flint_randinit(state); - - /* compare with fmpq */ - for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++) - { - arb_t a, b; - fmpq_t x, y, z; - ulong n; - slong i; - - arb_init(a); - arb_init(b); - - fmpq_init(x); - fmpq_init(y); - fmpq_init(z); - - arb_randtest(a, state, 1 + n_randint(state, 1000), 10); - arb_randtest(b, state, 1 + n_randint(state, 1000), 10); - n = n_randint(state, 80); - - arb_get_rand_fmpq(x, state, a, 1 + n_randint(state, 10)); - arb_rising_ui_bs(b, a, n, 2 + n_randint(state, 1000)); - - fmpq_one(y); - for (i = 0; i < n; i++) - { - fmpq_set_si(z, i, 1); - fmpq_add(z, x, z); - fmpq_mul(y, y, z); - } - - if (!arb_contains_fmpq(b, y)) - { - flint_printf("FAIL: containment\n\n"); - flint_printf("n = %wu\n", n); - flint_printf("a = "); arb_print(a); flint_printf("\n\n"); - flint_printf("x = "); fmpq_print(x); flint_printf("\n\n"); - flint_printf("b = "); arb_print(b); flint_printf("\n\n"); - flint_printf("y = "); fmpq_print(y); flint_printf("\n\n"); - flint_abort(); - } - - arb_clear(a); - arb_clear(b); - - fmpq_clear(x); - fmpq_clear(y); - fmpq_clear(z); - } - - /* aliasing of y and x */ - for (iter = 0; iter < 500 * arb_test_multiplier(); iter++) - { - arb_t x, y; - ulong n; - slong prec; - - arb_init(x); - arb_init(y); - - arb_randtest(x, state, 1 + n_randint(state, 200), 10); - arb_randtest(y, state, 1 + n_randint(state, 200), 10); - n = n_randint(state, 100); - - prec = 2 + n_randint(state, 1000); - - arb_rising_ui_bs(y, x, n, prec); - arb_rising_ui_bs(x, x, n, prec); - - if (!arb_equal(x, y)) - { - flint_printf("FAIL: aliasing\n\n"); - flint_printf("x = "); arb_print(x); flint_printf("\n\n"); - flint_printf("y = "); arb_print(y); flint_printf("\n\n"); - flint_printf("n = %wu\n", n); - flint_abort(); - } - - arb_clear(x); - arb_clear(y); - } - - flint_randclear(state); - flint_cleanup(); - flint_printf("PASS\n"); - return EXIT_SUCCESS; -} diff --git a/arb/test/t-rising_ui_rec.c b/arb/test/t-rising_ui_rec.c deleted file mode 100644 index 312649bd..00000000 --- a/arb/test/t-rising_ui_rec.c +++ /dev/null @@ -1,109 +0,0 @@ -/* - Copyright (C) 2012 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.h" - -int main() -{ - slong iter; - flint_rand_t state; - - flint_printf("rising_ui_rec...."); - fflush(stdout); - - flint_randinit(state); - - /* compare with fmpq */ - for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++) - { - arb_t a, b; - fmpq_t x, y, z; - ulong n; - slong i; - - arb_init(a); - arb_init(b); - - fmpq_init(x); - fmpq_init(y); - fmpq_init(z); - - arb_randtest(a, state, 1 + n_randint(state, 1000), 10); - arb_randtest(b, state, 1 + n_randint(state, 1000), 10); - n = n_randint(state, 80); - - arb_get_rand_fmpq(x, state, a, 1 + n_randint(state, 10)); - arb_rising_ui_rec(b, a, n, 2 + n_randint(state, 1000)); - - fmpq_one(y); - for (i = 0; i < n; i++) - { - fmpq_set_si(z, i, 1); - fmpq_add(z, x, z); - fmpq_mul(y, y, z); - } - - if (!arb_contains_fmpq(b, y)) - { - flint_printf("FAIL: containment\n\n"); - flint_printf("n = %wu\n", n); - flint_printf("a = "); arb_print(a); flint_printf("\n\n"); - flint_printf("x = "); fmpq_print(x); flint_printf("\n\n"); - flint_printf("b = "); arb_print(b); flint_printf("\n\n"); - flint_printf("y = "); fmpq_print(y); flint_printf("\n\n"); - flint_abort(); - } - - arb_clear(a); - arb_clear(b); - - fmpq_clear(x); - fmpq_clear(y); - fmpq_clear(z); - } - - /* aliasing of y and x */ - for (iter = 0; iter < 500 * arb_test_multiplier(); iter++) - { - arb_t x, y; - ulong n; - slong prec; - - arb_init(x); - arb_init(y); - - arb_randtest(x, state, 1 + n_randint(state, 200), 10); - arb_randtest(y, state, 1 + n_randint(state, 200), 10); - n = n_randint(state, 100); - - prec = 2 + n_randint(state, 1000); - - arb_rising_ui_rec(y, x, n, prec); - arb_rising_ui_rec(x, x, n, prec); - - if (!arb_equal(x, y)) - { - flint_printf("FAIL: aliasing\n\n"); - flint_printf("x = "); arb_print(x); flint_printf("\n\n"); - flint_printf("y = "); arb_print(y); flint_printf("\n\n"); - flint_printf("n = %wu\n", n); - flint_abort(); - } - - arb_clear(x); - arb_clear(y); - } - - flint_randclear(state); - flint_cleanup(); - flint_printf("PASS\n"); - return EXIT_SUCCESS; -} diff --git a/arb/test/t-rising_ui_rs.c b/arb/test/t-rising_ui_rs.c deleted file mode 100644 index ccbc2b33..00000000 --- a/arb/test/t-rising_ui_rs.c +++ /dev/null @@ -1,112 +0,0 @@ -/* - Copyright (C) 2012 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.h" - -int main() -{ - slong iter; - flint_rand_t state; - - flint_printf("rising_ui_rs...."); - fflush(stdout); - - flint_randinit(state); - - /* compare with fmpq */ - for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++) - { - arb_t a, b; - fmpq_t x, y, z; - ulong n, step; - slong i; - - arb_init(a); - arb_init(b); - - fmpq_init(x); - fmpq_init(y); - fmpq_init(z); - - arb_randtest(a, state, 1 + n_randint(state, 1000), 10); - arb_randtest(b, state, 1 + n_randint(state, 1000), 10); - n = n_randint(state, 80); - - arb_get_rand_fmpq(x, state, a, 1 + n_randint(state, 10)); - - step = n_randint(state, 20); - arb_rising_ui_rs(b, a, n, step, 2 + n_randint(state, 1000)); - - fmpq_one(y); - for (i = 0; i < n; i++) - { - fmpq_set_si(z, i, 1); - fmpq_add(z, x, z); - fmpq_mul(y, y, z); - } - - if (!arb_contains_fmpq(b, y)) - { - flint_printf("FAIL: containment\n\n"); - flint_printf("n = %wu\n", n); - flint_printf("a = "); arb_print(a); flint_printf("\n\n"); - flint_printf("x = "); fmpq_print(x); flint_printf("\n\n"); - flint_printf("b = "); arb_print(b); flint_printf("\n\n"); - flint_printf("y = "); fmpq_print(y); flint_printf("\n\n"); - flint_abort(); - } - - arb_clear(a); - arb_clear(b); - - fmpq_clear(x); - fmpq_clear(y); - fmpq_clear(z); - } - - /* aliasing of y and x */ - for (iter = 0; iter < 500 * arb_test_multiplier(); iter++) - { - arb_t x, y; - ulong n, step; - slong prec; - - arb_init(x); - arb_init(y); - - arb_randtest(x, state, 1 + n_randint(state, 200), 10); - arb_randtest(y, state, 1 + n_randint(state, 200), 10); - n = n_randint(state, 100); - - prec = 2 + n_randint(state, 1000); - - step = n_randint(state, 20); - arb_rising_ui_rs(y, x, n, step, prec); - arb_rising_ui_rs(x, x, n, step, prec); - - if (!arb_equal(x, y)) - { - flint_printf("FAIL: aliasing\n\n"); - flint_printf("x = "); arb_print(x); flint_printf("\n\n"); - flint_printf("y = "); arb_print(y); flint_printf("\n\n"); - flint_printf("n = %wu\n", n); - flint_abort(); - } - - arb_clear(x); - arb_clear(y); - } - - flint_randclear(state); - flint_cleanup(); - flint_printf("PASS\n"); - return EXIT_SUCCESS; -} diff --git a/arb_hypgeom.h b/arb_hypgeom.h index 9b0b6892..12a62fab 100644 --- a/arb_hypgeom.h +++ b/arb_hypgeom.h @@ -50,7 +50,7 @@ typedef struct char negative; } arb_hypgeom_gamma_coeff_t; -extern arb_hypgeom_gamma_coeff_t arb_hypgeom_gamma_coeffs[ARB_HYPGEOM_GAMMA_TAB_NUM]; +ARB_DLL extern arb_hypgeom_gamma_coeff_t arb_hypgeom_gamma_coeffs[ARB_HYPGEOM_GAMMA_TAB_NUM]; int _arb_hypgeom_gamma_coeff_shallow(arf_t c, mag_t err, slong i, slong prec); void arb_hypgeom_gamma_stirling(arb_t res, const arb_t x, int reciprocal, slong prec); diff --git a/arb_hypgeom/gamma_fmpq.c b/arb_hypgeom/gamma_fmpq.c index eac59618..996a87e5 100644 --- a/arb_hypgeom/gamma_fmpq.c +++ b/arb_hypgeom/gamma_fmpq.c @@ -10,6 +10,75 @@ */ #include "arb_hypgeom.h" +#include "hypgeom.h" + +void +arb_gamma_const_1_3_eval(arb_t s, slong prec) +{ + hypgeom_t series; + arb_t t, u; + slong wp = prec + 4 + 2 * FLINT_BIT_COUNT(prec); + + arb_init(t); + arb_init(u); + + hypgeom_init(series); + + fmpz_poly_set_str(series->A, "1 1"); + fmpz_poly_set_str(series->B, "1 1"); + fmpz_poly_set_str(series->P, "4 5 -46 108 -72"); + fmpz_poly_set_str(series->Q, "4 0 0 0 512000"); + + prec += FLINT_CLOG2(prec); + + arb_hypgeom_infsum(s, t, series, wp, wp); + + arb_sqrt_ui(u, 10, wp); + arb_mul(t, t, u, wp); + + arb_const_pi(u, wp); + arb_pow_ui(u, u, 4, wp); + arb_mul_ui(u, u, 12, wp); + arb_mul(s, s, u, wp); + + arb_div(s, s, t, wp); + arb_root_ui(s, s, 2, wp); + arb_root_ui(s, s, 3, prec); + + hypgeom_clear(series); + arb_clear(t); + arb_clear(u); +} + +ARB_DEF_CACHED_CONSTANT(arb_gamma_const_1_3, arb_gamma_const_1_3_eval) + +void +arb_gamma_const_1_4_eval(arb_t x, slong prec) +{ + arb_t t, u; + slong wp = prec + 4 + 2 * FLINT_BIT_COUNT(prec); + + arb_init(t); + arb_init(u); + + arb_one(t); + arb_sqrt_ui(u, 2, wp); + arb_agm(x, t, u, wp); + + arb_const_pi(t, wp); + arb_mul_2exp_si(t, t, 1); + arb_sqrt(u, t, wp); + arb_mul(u, u, t, wp); + + arb_div(x, u, x, wp); + arb_sqrt(x, x, wp); + + arb_clear(t); + arb_clear(u); +} + +ARB_DEF_CACHED_CONSTANT(arb_gamma_const_1_4, arb_gamma_const_1_4_eval) + void arb_hypgeom_gamma_stirling_choose_param(int * reflect, slong * r, slong * n, const arb_t x, int use_reflect, int digamma, slong prec); void arb_hypgeom_gamma_stirling_inner(arb_t s, const arb_t z, slong N, slong prec); @@ -67,9 +136,6 @@ arb_hypgeom_gamma_fmpq_stirling(arb_t y, const fmpq_t a, slong prec) arb_clear(x); } -void arb_gamma_const_1_3(arb_t x, slong prec); -void arb_gamma_const_1_4(arb_t x, slong prec); - void arb_hypgeom_gamma_small_frac(arb_t y, unsigned int p, unsigned int q, slong prec) { diff --git a/arb_hypgeom/test/t-gamma_taylor_tab.c b/arb_hypgeom/test/t-gamma_taylor_tab.c index a123b919..d95de4c0 100644 --- a/arb_hypgeom/test/t-gamma_taylor_tab.c +++ b/arb_hypgeom/test/t-gamma_taylor_tab.c @@ -11,7 +11,7 @@ #include "arb_hypgeom.h" -extern arb_hypgeom_gamma_coeff_t arb_hypgeom_gamma_coeffs[ARB_HYPGEOM_GAMMA_TAB_NUM]; +ARB_DLL extern arb_hypgeom_gamma_coeff_t arb_hypgeom_gamma_coeffs[ARB_HYPGEOM_GAMMA_TAB_NUM]; int main() { diff --git a/arb_poly/digamma_series.c b/arb_poly/digamma_series.c index 511c35b3..14ccc0dc 100644 --- a/arb_poly/digamma_series.c +++ b/arb_poly/digamma_series.c @@ -11,7 +11,7 @@ #include "arb_poly.h" -void arb_gamma_stirling_choose_param(int * reflect, slong * r, slong * n, +void arb_hypgeom_gamma_stirling_choose_param(int * reflect, slong * r, slong * n, const arb_t x, int use_reflect, int digamma, slong prec); void _arb_poly_gamma_stirling_eval2(arb_ptr res, const arb_t z, slong n, slong num, int diff, slong prec); @@ -80,7 +80,7 @@ _arb_poly_digamma_series(arb_ptr res, arb_srcptr h, slong hlen, slong len, slong else { /* use Stirling series */ - arb_gamma_stirling_choose_param(&reflect, &r, &n, h, 1, 1, wp); + arb_hypgeom_gamma_stirling_choose_param(&reflect, &r, &n, h, 1, 1, wp); /* psi(x) = psi((1-x)+r) - h(1-x,r) - pi*cot(pi*x) */ if (reflect) diff --git a/arb_poly/gamma_series.c b/arb_poly/gamma_series.c index 09813142..ac4242c3 100644 --- a/arb_poly/gamma_series.c +++ b/arb_poly/gamma_series.c @@ -13,7 +13,7 @@ void arb_gamma_stirling_bound(mag_ptr err, const arb_t x, slong k0, slong knum, slong n); -void arb_gamma_stirling_choose_param(int * reflect, slong * r, slong * n, +void arb_hypgeom_gamma_stirling_choose_param(int * reflect, slong * r, slong * n, const arb_t x, int use_reflect, int digamma, slong prec); void arb_gamma_stirling_coeff(arb_t b, ulong k, int digamma, slong prec); @@ -277,7 +277,7 @@ _arb_poly_gamma_series(arb_ptr res, arb_srcptr h, slong hlen, slong len, slong p else { /* otherwise use Stirling series */ - arb_gamma_stirling_choose_param(&reflect, &r, &n, h, 1, 0, wp); + arb_hypgeom_gamma_stirling_choose_param(&reflect, &r, &n, h, 1, 0, wp); /* gamma(h) = (rf(1-h, r) * pi) / (gamma(1-h+r) sin(pi h)), h = h0 + t*/ if (reflect) diff --git a/arb_poly/lgamma_series.c b/arb_poly/lgamma_series.c index 95e75361..6ed3faa8 100644 --- a/arb_poly/lgamma_series.c +++ b/arb_poly/lgamma_series.c @@ -10,12 +10,13 @@ */ #include "arb_poly.h" +#include "arb_hypgeom.h" slong arf_get_si(const arf_t x, arf_rnd_t rnd); void _arb_poly_lgamma_series_at_one(arb_ptr u, slong len, slong prec); -void arb_gamma_stirling_choose_param(int * reflect, slong * r, slong * n, +void arb_hypgeom_gamma_stirling_choose_param(int * reflect, slong * r, slong * n, const arb_t x, int use_reflect, int digamma, slong prec); void _arb_poly_gamma_stirling_eval(arb_ptr res, const arb_t z, slong n, slong num, slong prec); @@ -23,20 +24,11 @@ void _arb_poly_gamma_stirling_eval(arb_ptr res, const arb_t z, slong n, slong nu static __inline__ void _log_rising_ui_series(arb_ptr t, const arb_t x, slong r, slong len, slong prec) { - arb_struct f[2]; slong rflen; - arb_init(f); - arb_init(f + 1); - arb_set(f, x); - arb_one(f + 1); - rflen = FLINT_MIN(len, r + 1); - _arb_poly_rising_ui_series(t, f, FLINT_MIN(2, len), r, rflen, prec); + arb_hypgeom_rising_ui_jet(t, x, r, rflen, prec); _arb_poly_log_series(t, t, rflen, len, prec); - - arb_clear(f); - arb_clear(f + 1); } void @@ -91,7 +83,7 @@ _arb_poly_lgamma_series(arb_ptr res, arb_srcptr h, slong hlen, slong len, slong else { /* otherwise use Stirling series */ - arb_gamma_stirling_choose_param(&reflect, &r, &n, h, 0, 0, wp); + arb_hypgeom_gamma_stirling_choose_param(&reflect, &r, &n, h, 0, 0, wp); arb_add_ui(zr, h, r, wp); _arb_poly_gamma_stirling_eval(u, zr, n, len, wp); diff --git a/arb_poly/rgamma_series.c b/arb_poly/rgamma_series.c index a04a8b28..32aa6fc5 100644 --- a/arb_poly/rgamma_series.c +++ b/arb_poly/rgamma_series.c @@ -15,7 +15,7 @@ slong arf_get_si(const arf_t x, arf_rnd_t rnd); void _arb_poly_lgamma_series_at_one(arb_ptr u, slong len, slong prec); -void arb_gamma_stirling_choose_param(int * reflect, slong * r, slong * n, +void arb_hypgeom_gamma_stirling_choose_param(int * reflect, slong * r, slong * n, const arb_t x, int use_reflect, int digamma, slong prec); void _arb_poly_gamma_stirling_eval(arb_ptr res, const arb_t z, slong n, slong num, slong prec); @@ -101,7 +101,7 @@ _arb_poly_rgamma_series(arb_ptr res, arb_srcptr h, slong hlen, slong len, slong else { /* otherwise use Stirling series */ - arb_gamma_stirling_choose_param(&reflect, &r, &n, h, 1, 0, wp); + arb_hypgeom_gamma_stirling_choose_param(&reflect, &r, &n, h, 1, 0, wp); /* rgamma(h) = (gamma(1-h+r) sin(pi h)) / (rf(1-h, r) * pi), h = h0 + t*/ if (reflect) diff --git a/doc/source/acb.rst b/doc/source/acb.rst index 0da54bc9..25fdf53f 100644 --- a/doc/source/acb.rst +++ b/doc/source/acb.rst @@ -864,36 +864,18 @@ Lambert W function Rising factorials ------------------------------------------------------------------------------- -.. function:: void acb_rising_ui_bs(acb_t z, const acb_t x, ulong n, slong prec) - -.. function:: void acb_rising_ui_rs(acb_t z, const acb_t x, ulong n, ulong step, slong prec) - -.. function:: void acb_rising_ui_rec(acb_t z, const acb_t x, ulong n, slong prec) - -.. function:: void acb_rising_ui(acb_t z, const acb_t x, ulong n, slong prec) - -.. function:: void acb_rising(acb_t z, const acb_t x, const acb_t n, slong prec) +.. function:: void acb_rising_ui(acb_t z, const acb_t x, const acb_t n, slong prec) + void acb_rising(acb_t z, const acb_t x, const acb_t n, slong prec) Computes the rising factorial `z = x (x+1) (x+2) \cdots (x+n-1)`. - - The *bs* version uses binary splitting. The *rs* version uses rectangular - splitting. The *rec* version uses either *bs* or *rs* depending - on the input. The default version uses the gamma function unless - *n* is a small integer. - - The *rs* version takes an optional *step* parameter for tuning - purposes (to use the default step length, pass zero). - -.. function :: void acb_rising2_ui_bs(acb_t u, acb_t v, const acb_t x, ulong n, slong prec) - -.. function :: void acb_rising2_ui_rs(acb_t u, acb_t v, const acb_t x, ulong n, ulong step, slong prec) + These functions are aliases for :func:`acb_hypgeom_rising_ui` + and :func:`acb_hypgeom_rising`. .. function :: void acb_rising2_ui(acb_t u, acb_t v, const acb_t x, ulong n, slong prec) Letting `u(x) = x (x+1) (x+2) \cdots (x+n-1)`, simultaneously compute - `u(x)` and `v(x) = u'(x)`, respectively using binary splitting, - rectangular splitting (with optional nonzero step length *step* - to override the default choice), and an automatic algorithm choice. + `u(x)` and `v(x) = u'(x)`. + This function is a wrapper of :func:`acb_hypgeom_rising_ui_jet`. .. function :: void acb_rising_ui_get_mag(mag_t bound, const acb_t x, ulong n) @@ -907,15 +889,18 @@ Gamma function .. function:: void acb_gamma(acb_t y, const acb_t x, slong prec) Computes the gamma function `y = \Gamma(x)`. + This is an alias for :func:`acb_hypgeom_gamma`. .. function:: void acb_rgamma(acb_t y, const acb_t x, slong prec) Computes the reciprocal gamma function `y = 1/\Gamma(x)`, avoiding division by zero at the poles of the gamma function. + This is an alias for :func:`acb_hypgeom_rgamma`. .. function:: void acb_lgamma(acb_t y, const acb_t x, slong prec) Computes the logarithmic gamma function `y = \log \Gamma(x)`. + This is an alias for :func:`acb_hypgeom_lgamma`. The branch cut of the logarithmic gamma function is placed on the negative half-axis, which means that diff --git a/doc/source/arb.rst b/doc/source/arb.rst index 4161d06a..3ab56e53 100644 --- a/doc/source/arb.rst +++ b/doc/source/arb.rst @@ -1318,25 +1318,12 @@ Lambert W function Gamma function and factorials ------------------------------------------------------------------------------- -.. function:: void arb_rising_ui_bs(arb_t z, const arb_t x, ulong n, slong prec) - -.. function:: void arb_rising_ui_rs(arb_t z, const arb_t x, ulong n, ulong step, slong prec) - -.. function:: void arb_rising_ui_rec(arb_t z, const arb_t x, ulong n, slong prec) - .. function:: void arb_rising_ui(arb_t z, const arb_t x, ulong n, slong prec) - -.. function:: void arb_rising(arb_t z, const arb_t x, const arb_t n, slong prec) + void arb_rising(arb_t z, const arb_t x, const arb_t n, slong prec) Computes the rising factorial `z = x (x+1) (x+2) \cdots (x+n-1)`. - - The *bs* version uses binary splitting. The *rs* version uses rectangular - splitting. The *rec* version uses either *bs* or *rs* depending - on the input. The default version uses the gamma function unless - *n* is a small integer. - - The *rs* version takes an optional *step* parameter for tuning - purposes (to use the default step length, pass zero). + These functions are aliases for :func:`arb_hypgeom_rising_ui` + and :func:`arb_hypgeom_rising`. .. function:: void arb_rising_fmpq_ui(arb_t z, const fmpq_t x, ulong n, slong prec) @@ -1345,16 +1332,11 @@ Gamma function and factorials compared to *prec*, it is more efficient to convert *x* to an approximation and use :func:`arb_rising_ui`. -.. function :: void arb_rising2_ui_bs(arb_t u, arb_t v, const arb_t x, ulong n, slong prec) - -.. function :: void arb_rising2_ui_rs(arb_t u, arb_t v, const arb_t x, ulong n, ulong step, slong prec) - .. function :: void arb_rising2_ui(arb_t u, arb_t v, const arb_t x, ulong n, slong prec) Letting `u(x) = x (x+1) (x+2) \cdots (x+n-1)`, simultaneously compute - `u(x)` and `v(x) = u'(x)`, respectively using binary splitting, - rectangular splitting (with optional nonzero step length *step* - to override the default choice), and an automatic algorithm choice. + `u(x)` and `v(x) = u'(x)`. + This function is a wrapper of :func:`arb_hypgeom_rising_ui_jet`. .. function:: void arb_fac_ui(arb_t z, ulong n, slong prec) @@ -1372,23 +1354,26 @@ Gamma function and factorials rising factorial as `{n \choose k} = (n-k+1)_k / k!`. .. function:: void arb_gamma(arb_t z, const arb_t x, slong prec) - -.. function:: void arb_gamma_fmpq(arb_t z, const fmpq_t x, slong prec) - -.. function:: void arb_gamma_fmpz(arb_t z, const fmpz_t x, slong prec) + void arb_gamma_fmpq(arb_t z, const fmpq_t x, slong prec) + void arb_gamma_fmpz(arb_t z, const fmpz_t x, slong prec) Computes the gamma function `z = \Gamma(x)`. + These functions are aliases for :func:`arb_hypgeom_gamma`, + :func:`arb_hypgeom_gamma_fmpq`, :func:`arb_hypgeom_gamma_fmpz`. + .. function:: void arb_lgamma(arb_t z, const arb_t x, slong prec) Computes the logarithmic gamma function `z = \log \Gamma(x)`. The complex branch structure is assumed, so if `x \le 0`, the result is an indeterminate interval. + This function is an alias for :func:`arb_hypgeom_lgamma`. .. function:: void arb_rgamma(arb_t z, const arb_t x, slong prec) Computes the reciprocal gamma function `z = 1/\Gamma(x)`, avoiding division by zero at the poles of the gamma function. + This function is an alias for :func:`arb_hypgeom_rgamma`. .. function:: void arb_digamma(arb_t y, const arb_t x, slong prec)