mirror of
https://github.com/vale981/arb
synced 2025-03-04 08:51:40 -05:00
begin make old gamma functions call new code
This commit is contained in:
fredrik
parent
9916182253
commit
52474120ce
44 changed files with 179 additions and 3360 deletions
6
acb.h
6
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);
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
|
|
|
|||
249
acb/gamma.c
249
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);
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -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];
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#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);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#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);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -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);
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#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);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#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);
|
||||
}
|
||||
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#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);
|
||||
}
|
||||
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#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;
|
||||
}
|
||||
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#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;
|
||||
}
|
||||
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#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;
|
||||
}
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#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;
|
||||
}
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#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;
|
||||
}
|
||||
|
|
@ -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);
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
6
arb.h
6
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);
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
|
|
|
|||
946
arb/gamma.c
946
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);
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -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];
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#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);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#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);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -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);
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#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);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#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);
|
||||
}
|
||||
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#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);
|
||||
}
|
||||
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#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;
|
||||
}
|
||||
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#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;
|
||||
}
|
||||
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#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;
|
||||
}
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#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;
|
||||
}
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#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;
|
||||
}
|
||||
|
|
@ -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);
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
{
|
||||
|
|
|
|||
|
|
@ -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()
|
||||
{
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue