hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-05 09:21:38 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

add fmpcb_gamma, fmpcb_rgamma, fmpcb_lgamma and refactored Stirling series code

Browse source
This commit is contained in:
Fredrik Johansson 2013-02-20 11:58:08 +01:00
parent 3bd6ad4cdd
commit 6921bdcee4
19 changed files with 1302 additions and 201 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

27
doc/source/fmpcb.rst
View file

@ -279,6 +279,10 @@ Elementary functions
Sets `s = \sin z`, `c = \cos z`.
.. function:: void fmpcb_sin_pi(fmpcb_t s, const fmpcb_t z, long prec)
Sets `s = \sin \pi z`.
.. function:: void fmpcb_pow_fmpz(fmpcb_t y, const fmpcb_t b, const fmpz_t e, long prec)
.. function:: void fmpcb_pow_ui(fmpcb_t y, const fmpcb_t b, ulong e, long prec)
@ -328,6 +332,29 @@ Elementary functions
Special functions
-------------------------------------------------------------------------------
.. function:: void fmpcb_rfac_ui_bsplit(fmpcb_t y, const fmpcb_t x, ulong n, long prec)
Sets *x* to the rising factorial `x (x+1) (x+2) \cdots (x+n-1)`,
computed using binary splitting.
.. function:: void fmpcb_gamma(fmpcb_t y, const fmpcb_t x, long prec)
.. function:: void fmpcb_rgamma(fmpcb_t y, const fmpcb_t x, long prec)
.. function:: void fmpcb_lgamma(fmpcb_t y, const fmpcb_t x, long prec)
Sets, respectively, `y = \Gamma(x)`, `y = 1/\Gamma(x)`,
`y = \log \Gamma(x)`.
The branch cut of the logarithmic gamma function is placed on the
negative half-axis, which means that
`\log \Gamma(z) + \log z = \log \Gamma(z+1)` holds for all `z`,
whereas `\log \Gamma(z) \ne \log(\Gamma(z))` in general.
Warning: the log gamma function does not currently use the reflection
formula, and gets very slow for `z` far into the left half-plane.
These functions are simple wrappers for the Stirling series code in the *gamma* module.
.. function:: void fmpcb_zeta_series_em_sum(fmpcb_struct * z, const fmpcb_t s, const fmpcb_t a, int deflate, ulong N, ulong M, long d, long prec)
.. function:: void fmpcb_zeta_series(fmpcb_struct * z, const fmpcb_t s, const fmpcb_t a, int deflate, long d, long prec)

30
doc/source/fmprb.rst
View file

@ -884,32 +884,6 @@ Gamma function
.. function:: void fmprb_lgamma(fmprb_t y, const fmprb_t x, long prec)
Sets, respectively, `y = \Gamma(x)`, `y = 1/\Gamma(x)`,
`y = \log \Gamma(x)`.
For large `x`, uses Stirling's expansion
.. math ::
\log \Gamma(x) = \left(x-\frac{1}{2}\right)\log x - x +
\frac{\ln {2 \pi}}{2} + \sum_{k=1}^{n-1} t_k + R(n,x)
where
.. math ::
t_k = \frac{B_{2k}}{2k(2k-1)x^{2k-1}}
and `|R(n,x)| < t_n`.
If `x` is too small for the asymptotic expansion to give sufficient
accuracy directly, we translate to `x + r`
using the formula `\Gamma(x) = \Gamma(x+r) /
(x (x+1) (x+2) \cdots (x+r-1))`.
To obtain a remainder smaller than `2^{-b}`, we must choose an `r` such
that `x + r > \beta b`, where `\beta > \log(2) / (2 \pi) \approx 0.11`.
We use a slightly larger factor `\beta \approx 0.2` to more closely
balance `n` and `r`. A much larger `\beta` (e.g. `\beta = 1`) could be
used to reduce the number of Bernoulli numbers that have to be
precomputed, at the expense of slower repeated evaluation.
`y = \log \Gamma(x)`. These functions are simple wrappers for the
Stirling series code in the *gamma* module.

75
doc/source/gamma.rst
View file

@ -1,6 +1,20 @@
**gamma.h** -- support for the gamma function
===============================================================================
This module implements various algorithms for evaluating the
gamma function and related functions. The functions provided here are mainly
intended for internal use, though they may be useful to call directly in some
applications where the default algorithm choices are suboptimal.
Most applications should use the standard, user-friendly top-level functions:
* :func:`fmprb_gamma`
* :func:`fmprb_rgamma`
* :func:`fmprb_lgamma`
* :func:`fmpcb_gamma`
* :func:`fmpcb_rgamma`
* :func:`fmpcb_lgamma`
Evaluation using Taylor series
--------------------------------------------------------------------------------
@ -33,4 +47,65 @@ Evaluation using Taylor series
in order for the translation to be reasonably efficient (too large
values will result in extreme slowness and eventually an exception).
Evaluation using the Stirling series
--------------------------------------------------------------------------------
.. function :: void gamma_stirling_choose_param(int * reflect, long * r, long * n, double x, double y, double beta, int allow_reflection, long prec)
Uses double precision arithmetic to compute parameters `r`, `n` such that
the remainder in the Stirling series approximately is bounded
by `2^{-\mathrm{prec}}`.
The parameter `n` is the truncation point in the asymptotic
Stirling series. If `|z|` is too small for the Stirling series
to give sufficient accuracy directly, we first translate to `z + r`
using the formula `\Gamma(z) = \Gamma(z+r) /
(z (z+1) (z+2) \cdots (z+r-1))`.
If *allow_reflection* is nonzero, the *reflect* flag is set if `z`
should be replaced with `1-z` using the reflection formula.
Note that this function does not guarantee the error bound rigorously;
a rigorous error bound, which also accounts for the radius of `z`,
is computed a posteriori when evaluating the Stirling series.
However, in practice, this function does estimate the bound
very accurately.
To obtain a remainder smaller than `2^{-b}`, we must choose an `r` such
that `x + r > \beta b`, where `\beta > \log(2) / (2 \pi) \approx 0.11`.
In practice, a slightly larger factor `\beta \approx 0.2` more closely
balances `n` and `r`. A much larger `\beta` (e.g. `\beta = 1`) could be
used to reduce the number of Bernoulli numbers that have to be
precomputed, at the expense of slower repeated evaluation.
.. function :: void gamma_stirling_coeff(fmprb_t b, ulong k, long prec)
Sets `b = B_{2k} / (2k (2k-1))`, rounded to *prec* bits.
Assumes that the Bernoulli number has been precomputed.
.. function :: void gamma_stirling_eval_series_fmprb(fmprb_t s, const fmprb_t z, long n, long prec)
.. function :: void gamma_stirling_eval_series_fmpcb(fmpcb_t s, const fmpcb_t z, long n, long prec)
Evaluates the Stirling series
.. math ::
\log \Gamma(z) - R(n,z) = \left(z-\frac{1}{2}\right)\log z - z +
\frac{\ln {2 \pi}}{2} + \sum_{k=1}^{n-1} t_k
where
.. math ::
t_k = \frac{B_{2k}}{2k(2k-1)z^{2k-1}}.
The bound
.. math ::
|R(n,z)| \le \frac{t_n}{\cos(0.5 \arg(z))^{2n}}
is included in the radius of the output.

7
fmpcb.h
View file

@ -488,14 +488,21 @@ void fmpcb_sin(fmpcb_t r, const fmpcb_t z, long prec);
void fmpcb_cos(fmpcb_t r, const fmpcb_t z, long prec);
void fmpcb_sin_cos(fmpcb_t s, fmpcb_t c, const fmpcb_t z, long prec);
void fmpcb_sin_pi(fmpcb_t r, const fmpcb_t z, long prec);
void fmpcb_pow(fmpcb_t r, const fmpcb_t x, const fmpcb_t y, long prec);
void fmpcb_rfac_ui_bsplit(fmpcb_t y, const fmpcb_t x, ulong n, long prec);
void fmpcb_invroot_newton(fmpcb_t r, const fmpcb_t a, ulong m, const fmpcb_t r0, long startprec, long prec);
void fmpcb_root_exp(fmpcb_t r, const fmpcb_t a, long m, long index, long prec);
void fmpcb_root_newton(fmpcb_t r, const fmpcb_t a, long m, long index, long prec);
void fmpcb_root(fmpcb_t r, const fmpcb_t a, long m, long index, long prec);
void fmpcb_gamma(fmpcb_t y, const fmpcb_t x, long prec);
void fmpcb_rgamma(fmpcb_t y, const fmpcb_t x, long prec);
void fmpcb_lgamma(fmpcb_t y, const fmpcb_t x, long prec);
void fmpcb_zeta_series_em_sum(fmpcb_struct * z, const fmpcb_t s, const fmpcb_t a, int deflate, ulong N, ulong M, long d, long prec);
void fmpcb_zeta_series_em_choose_param(fmpr_t bound, ulong * N, ulong * M, const fmpcb_t s, const fmpcb_t a, long d, long target, long prec);
void fmpcb_zeta_series_em_bound(fmpr_t bound, const fmpcb_t s, const fmpcb_t a, long N, long M, long d, long wp);

152
fmpcb/gamma.c Normal file
View file

@ -0,0 +1,152 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2013 Fredrik Johansson
******************************************************************************/
#include "fmpcb.h"
#include "gamma.h"
/* tuning factor */
#define GAMMA_STIRLING_BETA 0.21
static void
choose(int * reflect, long * r, long * n, const fmpcb_t z, int use_refl, long prec)
{
double x, y;
x = fmpr_get_d(fmprb_midref(fmpcb_realref(z)), FMPR_RND_NEAR);
y = fmpr_get_d(fmprb_midref(fmpcb_imagref(z)), FMPR_RND_NEAR);
gamma_stirling_choose_param(reflect, r, n, x, y,
GAMMA_STIRLING_BETA, use_refl, prec);
}
static void
_fmpcb_gamma(fmpcb_t y, const fmpcb_t x, long prec, int inverse)
{
int reflect;
long r, n, wp;
fmpcb_t t, u, v;
wp = prec + FLINT_BIT_COUNT(prec);
choose(&reflect, &r, &n, x, 1, wp);
fmpcb_init(t);
fmpcb_init(u);
fmpcb_init(v);
if (reflect)
{
/* gamma(x) = (rf(1-x, r) * pi) / (gamma(1-x+r) sin(pi x)) */
fmpcb_sub_ui(t, x, 1, wp);
fmpcb_neg(t, t);
fmpcb_rfac_ui_bsplit(u, t, r, wp);
fmprb_const_pi(fmpcb_realref(v), wp);
fmpcb_mul_fmprb(u, u, fmpcb_realref(v), wp);
fmpcb_add_ui(t, t, r, wp);
gamma_stirling_eval_series_fmpcb(v, t, n, wp);
fmpcb_exp(v, v, wp);
fmpcb_sin_pi(t, x, wp);
fmpcb_mul(v, v, t, wp);
}
else
{
/* gamma(x) = gamma(x+r) / rf(x,r) */
fmpcb_add_ui(t, x, r, wp);
gamma_stirling_eval_series_fmpcb(u, t, n, wp);
fmpcb_exp(u, u, prec);
fmpcb_rfac_ui_bsplit(v, x, r, wp);
}
if (inverse)
fmpcb_div(y, v, u, prec);
else
fmpcb_div(y, u, v, prec);
fmpcb_clear(t);
fmpcb_clear(u);
fmpcb_clear(v);
}
void
fmpcb_gamma(fmpcb_t y, const fmpcb_t x, long prec)
{
_fmpcb_gamma(y, x, prec, 0);
}
void
fmpcb_rgamma(fmpcb_t y, const fmpcb_t x, long prec)
{
_fmpcb_gamma(y, x, prec, 1);
}
void
fmpcb_log_rfac_ui(fmpcb_t s, const fmpcb_t z, ulong r, long wp)
{
ulong i;
fmpcb_t t, u;
fmpcb_init(t);
fmpcb_init(u);
fmpcb_zero(t);
for (i = 0; i < r; i++)
{
fmpcb_add_ui(u, z, i, wp);
fmpcb_log(u, u, wp);
fmpcb_add(t, t, u, wp);
}
fmpcb_set(s, t);
fmpcb_clear(t);
fmpcb_clear(u);
}
void
fmpcb_lgamma(fmpcb_t y, const fmpcb_t x, long prec)
{
int reflect;
long r, n, wp;
fmpcb_t t, u;
wp = prec + FLINT_BIT_COUNT(prec);
choose(&reflect, &r, &n, x, 0, wp);
/* log(gamma(x)) = log(gamma(x+r)) - log(rf(x,r)) */
fmpcb_init(t);
fmpcb_init(u);
fmpcb_add_ui(t, x, r, wp);
gamma_stirling_eval_series_fmpcb(u, t, n, wp);
fmpcb_log_rfac_ui(t, x, r, wp);
fmpcb_sub(y, u, t, prec);
fmpcb_clear(t);
fmpcb_clear(u);
}

154
fmpcb/rfac_ui_bsplit.c Normal file
View file

@ -0,0 +1,154 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2012, 2013 Fredrik Johansson
******************************************************************************/
#include "fmpcb.h"
/* x(x+1)...(x+7) = (28 + 98x + 63x^2 + 14x^3 + x^4)^2 - 16 (7+2x)^2 */
static void
rfac_eight(fmpcb_t t, const fmpcb_t x, long prec)
{
if (prec < 768)
{
fmpcb_t u;
long i;
fmpcb_init(u);
fmpcb_add_ui(u, x, 1, prec);
fmpcb_mul(t, x, u, prec);
for (i = 2; i < 8; i++)
{
fmpcb_add_ui(u, u, 1, prec);
fmpcb_mul(t, t, u, prec);
}
fmpcb_clear(u);
}
else
{
fmpcb_t u, v;
fmpcb_init(u);
fmpcb_init(v);
/* t = x^2, v = x^3, u = x^4 */
fmpcb_mul(t, x, x, prec);
fmpcb_mul(v, x, t, prec);
fmpcb_mul(u, t, t, prec);
/* u = (28 + ...)^2 */
fmpcb_addmul_ui(u, v, 14UL, prec);
fmpcb_addmul_ui(u, t, 63UL, prec);
fmpcb_addmul_ui(u, x, 98UL, prec);
fmpcb_add_ui(u, u, 28UL, prec);
fmpcb_mul(u, u, u, prec);
/* 16 (7+2x)^2 = 784 + 448x + 64x^2 */
fmpcb_sub_ui(u, u, 784UL, prec);
fmpcb_submul_ui(u, x, 448UL, prec);
fmpcb_mul_2exp_si(t, t, 6);
fmpcb_sub(t, u, t, prec);
fmpcb_clear(u);
fmpcb_clear(v);
}
}
/* assumes that the length is a multiple of 8 */
static void
bsplit_eight(fmpcb_t y, const fmpcb_t x,
ulong a, ulong b, long prec)
{
if (b - a == 8)
{
fmpcb_t t;
fmpcb_init(t);
fmpcb_add_ui(t, x, a, prec);
rfac_eight(y, t, prec);
fmpcb_clear(t);
}
else
{
ulong m = a + ((b - a) / 16) * 8;
fmpcb_t L, R;
fmpcb_init(L);
fmpcb_init(R);
bsplit_eight(L, x, a, m, prec);
bsplit_eight(R, x, m, b, prec);
fmpcb_mul(y, L, R, prec);
fmpcb_clear(L);
fmpcb_clear(R);
}
}
void
fmpcb_rfac_ui_bsplit(fmpcb_t y, const fmpcb_t x, ulong n, long prec)
{
if (n == 0)
{
fmpcb_one(y);
}
else if (n == 1)
{
fmpcb_set(y, x);
}
else
{
long k, a, wp;
fmpcb_t t, u;
wp = FMPR_PREC_ADD(prec, FLINT_BIT_COUNT(n));
fmpcb_init(t);
fmpcb_init(u);
if (n >= 8)
{
bsplit_eight(t, x, 0, (n / 8) * 8, wp);
a = (n / 8) * 8;
}
else
{
fmpcb_set(t, x);
a = 1;
}
for (k = a; k < n; k++)
{
fmpcb_add_ui(u, x, k, wp);
fmpcb_mul(t, t, u, wp);
}
fmpcb_set(y, t);
fmpcb_clear(t);
fmpcb_clear(u);
}
}

57
fmpcb/sin_pi.c Normal file
View file

@ -0,0 +1,57 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2013 Fredrik Johansson
******************************************************************************/
#include "fmpcb.h"
void
fmpcb_sin_pi(fmpcb_t r, const fmpcb_t z, long prec)
{
#define a fmpcb_realref(z)
#define b fmpcb_imagref(z)
fmprb_t sa, ca, sb, cb;
fmprb_init(sa);
fmprb_init(ca);
fmprb_init(sb);
fmprb_init(cb);
fmprb_sin_cos_pi(sa, ca, a, prec);
fmprb_const_pi(cb, prec);
fmprb_mul(cb, cb, b, prec);
fmprb_sinh_cosh(sb, cb, cb, prec);
fmprb_mul(fmpcb_realref(r), sa, cb, prec);
fmprb_mul(fmpcb_imagref(r), sb, ca, prec);
fmprb_clear(sa);
fmprb_clear(ca);
fmprb_clear(sb);
fmprb_clear(cb);
#undef a
#undef b
}

89
fmpcb/test/t-gamma.c Normal file
View file

@ -0,0 +1,89 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2012, 2013 Fredrik Johansson
******************************************************************************/
#include "fmpcb.h"
int main()
{
long iter;
flint_rand_t state;
printf("gamma....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 3000; iter++)
{
fmpcb_t a, b, c;
long prec1, prec2;
prec1 = 2 + n_randint(state, 1000);
prec2 = prec1 + 30;
fmpcb_init(a);
fmpcb_init(b);
fmpcb_init(c);
fmprb_randtest_precise(fmpcb_realref(a), state, 1 + n_randint(state, 1000), 3);
fmprb_randtest_precise(fmpcb_imagref(a), state, 1 + n_randint(state, 1000), 3);
fmpcb_gamma(b, a, prec1);
fmpcb_gamma(c, a, prec2);
if (!fmpcb_overlaps(b, c))
{
printf("FAIL: overlap\n\n");
printf("a = "); fmpcb_print(a); printf("\n\n");
printf("b = "); fmpcb_print(b); printf("\n\n");
printf("c = "); fmpcb_print(c); printf("\n\n");
abort();
}
/* check gamma(z+1) = z gamma(z) */
fmpcb_mul(b, b, a, prec1);
fmpcb_add_ui(c, a, 1, prec1);
fmpcb_gamma(c, c, prec1);
if (!fmpcb_overlaps(b, c))
{
printf("FAIL: functional equation\n\n");
printf("a = "); fmpcb_print(a); printf("\n\n");
printf("b = "); fmpcb_print(b); printf("\n\n");
printf("c = "); fmpcb_print(c); printf("\n\n");
abort();
}
fmpcb_clear(a);
fmpcb_clear(b);
fmpcb_clear(c);
}
flint_randclear(state);
_fmpz_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

91
fmpcb/test/t-lgamma.c Normal file
View file

@ -0,0 +1,91 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2012, 2013 Fredrik Johansson
******************************************************************************/
#include "fmpcb.h"
int main()
{
long iter;
flint_rand_t state;
printf("lgamma....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 1000; iter++)
{
fmpcb_t a, b, c;
long prec1, prec2;
prec1 = 2 + n_randint(state, 500);
prec2 = prec1 + 30;
fmpcb_init(a);
fmpcb_init(b);
fmpcb_init(c);
fmprb_randtest_precise(fmpcb_realref(a), state, 1 + n_randint(state, 1000), 3);
fmprb_randtest_precise(fmpcb_imagref(a), state, 1 + n_randint(state, 1000), 3);
fmpcb_lgamma(b, a, prec1);
fmpcb_lgamma(c, a, prec2);
if (!fmpcb_overlaps(b, c))
{
printf("FAIL: overlap\n\n");
printf("a = "); fmpcb_print(a); printf("\n\n");
printf("b = "); fmpcb_print(b); printf("\n\n");
printf("c = "); fmpcb_print(c); printf("\n\n");
abort();
}
/* check lgamma(z+1) = lgamma(z) + log(z) */
fmpcb_log(c, a, prec1);
fmpcb_add(b, b, c, prec1);
fmpcb_add_ui(c, a, 1, prec1);
fmpcb_lgamma(c, c, prec1);
if (!fmpcb_overlaps(b, c))
{
printf("FAIL: functional equation\n\n");
printf("a = "); fmpcb_print(a); printf("\n\n");
printf("b = "); fmpcb_print(b); printf("\n\n");
printf("c = "); fmpcb_print(c); printf("\n\n");
abort();
}
fmpcb_clear(a);
fmpcb_clear(b);
fmpcb_clear(c);
}
flint_randclear(state);
_fmpz_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

89
fmpcb/test/t-rgamma.c Normal file
View file

@ -0,0 +1,89 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2012 Fredrik Johansson
******************************************************************************/
#include "fmpcb.h"
int main()
{
long iter;
flint_rand_t state;
printf("rgamma....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 3000; iter++)
{
fmpcb_t a, b, c;
long prec1, prec2;
prec1 = 2 + n_randint(state, 1000);
prec2 = prec1 + 30;
fmpcb_init(a);
fmpcb_init(b);
fmpcb_init(c);
fmprb_randtest_precise(fmpcb_realref(a), state, 1 + n_randint(state, 1000), 3);
fmprb_randtest_precise(fmpcb_imagref(a), state, 1 + n_randint(state, 1000), 3);
fmpcb_rgamma(b, a, prec1);
fmpcb_rgamma(c, a, prec2);
if (!fmpcb_overlaps(b, c))
{
printf("FAIL: overlap\n\n");
printf("a = "); fmpcb_print(a); printf("\n\n");
printf("b = "); fmpcb_print(b); printf("\n\n");
printf("c = "); fmpcb_print(c); printf("\n\n");
abort();
}
/* check 1/gamma(z+1) = 1/gamma(z)/z */
fmpcb_div(b, b, a, prec1);
fmpcb_add_ui(c, a, 1, prec1);
fmpcb_rgamma(c, c, prec1);
if (!fmpcb_overlaps(b, c))
{
printf("FAIL: functional equation\n\n");
printf("a = "); fmpcb_print(a); printf("\n\n");
printf("b = "); fmpcb_print(b); printf("\n\n");
printf("c = "); fmpcb_print(c); printf("\n\n");
abort();
}
fmpcb_clear(a);
fmpcb_clear(b);
fmpcb_clear(c);
}
flint_randclear(state);
_fmpz_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

222
fmprb/gamma.c
View file

@ -19,161 +19,105 @@
=============================================================================*/
/******************************************************************************
Copyright (C) 2012 Fredrik Johansson
Copyright (C) 2012, 2013 Fredrik Johansson
******************************************************************************/
#include <math.h>
#include "arith.h"
#include "fmprb.h"
#include "bernoulli.h"
void
fmprb_stirling_series_coeff(fmprb_t b, ulong k, long prec)
{
fmpz_t d;
fmpz_init(d);
fmprb_set_round_fmpz(b, fmpq_numref(bernoulli_cache + 2 * k), prec);
fmpz_mul2_uiui(d, fmpq_denref(bernoulli_cache + 2 * k), 2 * k, 2 * k - 1);
fmprb_div_fmpz(b, b, d, prec);
fmpz_clear(d);
}
/*
Heuristic: we use Stirling's series if abs(x) > beta * prec.
For convergence, beta must be greater than log(2)/(2*pi) ~= 0.11.
A larger beta gives faster convergence at the expense of extra
argument reduction.
*/
#include "gamma.h"
/* tuning factor */
#define GAMMA_STIRLING_BETA 0.21
long stirling_choose_r(const fmprb_t x, long wp)
static void
choose(int * reflect, long * r, long * n, const fmprb_t z,
int use_reflect, long prec)
{
double t = fmpr_get_d(fmprb_midref(x), FMPR_RND_FLOOR);
double want = FLINT_MAX(5, GAMMA_STIRLING_BETA * wp);
double x;
return (long) FLINT_MAX(0, want - t + 1);
x = fmpr_get_d(fmprb_midref(z), FMPR_RND_NEAR);
gamma_stirling_choose_param(reflect, r, n, x, 0.0,
GAMMA_STIRLING_BETA, use_reflect, prec);
}
/* TODO: speed up by caching Bernoulli number sizes */
long
stirling_choose_nterms(const fmprb_t x, long r, double bits)
static void
_fmprb_gamma(fmprb_t y, const fmprb_t x, long prec, int inverse)
{
long i;
double t, logt;
double mag;
int reflect;
long r, n, wp;
fmprb_t t, u, v;
t = fmpr_get_d(fmprb_midref(x), FMPR_RND_FLOOR) + r;
logt = log(t);
wp = prec + FLINT_BIT_COUNT(prec);
for (i = 1; ; i++)
{
mag = bernoulli_bound_2exp_si(2 * i) - (logt * 1.4426950408889634) * (2 * i - 1);
if (mag < -bits)
return i;
}
}
void
fmprb_gamma_log_stirling(fmprb_t s, const fmprb_t z, long nterms, long prec)
{
fmprb_t t, u, b, w, v;
long k, term_prec;
double z_mag, term_mag;
choose(&reflect, &r, &n, x, 1, wp);
fmprb_init(t);
fmprb_init(u);
fmprb_init(b);
fmprb_init(w);
fmprb_init(v);
fmprb_log(w, z, prec);
bernoulli_cache_compute(2 * (nterms + 1));
nterms = FLINT_MAX(nterms, 1);
fmprb_zero(s);
if (nterms > 1)
if (reflect)
{
fmprb_ui_div(t, 1UL, z, prec);
fmprb_mul(u, t, t, prec);
z_mag = fmpr_get_d(fmprb_midref(w), FMPR_RND_UP) * 1.44269504088896;
for (k = nterms - 1; k >= 1; k--)
{
term_mag = bernoulli_bound_2exp_si(2 * k);
term_mag -= (2 * k - 1) * z_mag;
term_prec = prec + term_mag;
term_prec = FLINT_MIN(term_prec, prec);
term_prec = FLINT_MAX(term_prec, 10);
fmprb_stirling_series_coeff(b, k, term_prec);
if (prec > 2000)
{
fmprb_set_round(v, u, term_prec);
fmprb_mul(s, s, v, term_prec);
}
else
{
fmprb_mul(s, s, u, term_prec);
}
fmprb_add(s, s, b, term_prec);
}
fmprb_mul(s, s, t, prec);
/* gamma(x) = (rf(1-x, r) * pi) / (gamma(1-x+r) sin(pi x)) */
fmprb_sub_ui(t, x, 1, wp);
fmprb_neg(t, t);
fmprb_rfac_ui_bsplit(u, t, r, wp);
fmprb_const_pi(v, wp);
fmprb_mul(u, u, v, wp);
fmprb_add_ui(t, t, r, wp);
gamma_stirling_eval_series_fmprb(v, t, n, wp);
fmprb_exp(v, v, wp);
fmprb_sin_pi(t, x, wp);
fmprb_mul(v, v, t, wp);
}
else
{
/* gamma(x) = gamma(x+r) / rf(x,r) */
fmprb_add_ui(t, x, r, wp);
gamma_stirling_eval_series_fmprb(u, t, n, wp);
fmprb_exp(u, u, prec);
fmprb_rfac_ui_bsplit(v, x, r, wp);
}
/* finally, we add the remainder error bound */
/* B_2n / (2n (2n-1)) */
fmprb_stirling_series_coeff(b, nterms, FMPRB_RAD_PREC);
/* 1/z^(2n-1) */
fmprb_ui_div(t, 1UL, z, FMPRB_RAD_PREC);
fmprb_pow_ui(t, t, 2 * nterms - 1, FMPRB_RAD_PREC);
fmprb_mul(t, t, b, FMPRB_RAD_PREC);
fmprb_add_error(s, t);
/* (z-0.5)*log(z) - z + log(2*pi)/2 */
fmprb_set_ui(t, 1);
fmprb_mul_2exp_si(t, t, -1);
fmprb_sub(t, z, t, prec);
fmprb_mul(t, w, t, prec);
fmprb_add(s, s, t, prec);
fmprb_sub(s, s, z, prec);
fmprb_const_log_sqrt2pi(t, prec);
fmprb_add(s, s, t, prec);
if (inverse)
fmprb_div(y, v, u, prec);
else
fmprb_div(y, u, v, prec);
fmprb_clear(t);
fmprb_clear(u);
fmprb_clear(b);
fmprb_clear(w);
fmprb_clear(v);
}
void
fmprb_gamma(fmprb_t y, const fmprb_t x, long prec)
{
_fmprb_gamma(y, x, prec, 0);
}
void
fmprb_rgamma(fmprb_t y, const fmprb_t x, long prec)
{
_fmprb_gamma(y, x, prec, 1);
}
void
fmprb_lgamma(fmprb_t y, const fmprb_t x, long prec)
{
int reflect;
long r, n, wp;
fmprb_t t, u;
wp = prec + FLINT_BIT_COUNT(prec);
r = stirling_choose_r(x, wp);
n = stirling_choose_nterms(x, r, wp);
choose(&reflect, &r, &n, x, 0, wp);
/* log(gamma(x)) = log(gamma(x+r)) - log(rf(x,r)) */
fmprb_init(t);
fmprb_init(u);
fmprb_add_ui(t, x, r, wp);
fmprb_gamma_log_stirling(u, t, n, wp);
gamma_stirling_eval_series_fmprb(u, t, n, wp);
fmprb_rfac_ui_bsplit(t, x, r, wp);
fmprb_log(t, t, wp);
@ -183,55 +127,3 @@ fmprb_lgamma(fmprb_t y, const fmprb_t x, long prec)
fmprb_clear(u);
}
void
fmprb_gamma(fmprb_t y, const fmprb_t x, long prec)
{
long r, n, wp;
fmprb_t t, u;
wp = prec + FLINT_BIT_COUNT(prec);
r = stirling_choose_r(x, wp);
n = stirling_choose_nterms(x, r, wp);
/* gamma(x) = gamma(x+r) / rf(x,r) */
fmprb_init(t);
fmprb_init(u);
fmprb_add_ui(t, x, r, wp);
fmprb_gamma_log_stirling(u, t, n, wp);
fmprb_exp(u, u, prec);
fmprb_rfac_ui_bsplit(t, x, r, wp);
fmprb_div(y, u, t, prec);
fmprb_clear(t);
fmprb_clear(u);
}
void
fmprb_rgamma(fmprb_t y, const fmprb_t x, long prec)
{
long r, n, wp;
fmprb_t t, u;
wp = prec + FLINT_BIT_COUNT(prec);
r = stirling_choose_r(x, wp);
n = stirling_choose_nterms(x, r, wp);
/* 1/gamma(x) = rf(x,r) / gamma(x+r) */
fmprb_init(t);
fmprb_init(u);
fmprb_add_ui(t, x, r, wp);
fmprb_gamma_log_stirling(u, t, n, wp);
fmprb_neg(u, u);
fmprb_exp(u, u, prec);
fmprb_rfac_ui_bsplit(t, x, r, wp);
fmprb_mul(y, u, t, prec);
fmprb_clear(t);
fmprb_clear(u);
}

14
fmprb/test/t-gamma.c
View file

@ -19,7 +19,7 @@
=============================================================================*/
/******************************************************************************
Copyright (C) 2012 Fredrik Johansson
Copyright (C) 2012, 2013 Fredrik Johansson
******************************************************************************/
@ -61,11 +61,17 @@ int main()
abort();
}
fmprb_gamma(a, a, prec2);
/* check gamma(z+1) = z gamma(z) */
fmprb_mul(b, b, a, prec1);
fmprb_add_ui(c, a, 1, prec1);
fmprb_gamma(c, c, prec1);
if (!fmprb_equal(a, c))
if (!fmprb_overlaps(b, c))
{
printf("FAIL: aliasing\n\n");
printf("FAIL: functional equation\n\n");
printf("a = "); fmprb_print(a); printf("\n\n");
printf("b = "); fmprb_print(b); printf("\n\n");
printf("c = "); fmprb_print(c); printf("\n\n");
abort();
}

15
fmprb/test/t-rgamma.c
View file

@ -19,7 +19,7 @@
=============================================================================*/
/******************************************************************************
Copyright (C) 2012 Fredrik Johansson
Copyright (C) 2012, 2013 Fredrik Johansson
******************************************************************************/
@ -61,11 +61,17 @@ int main()
abort();
}
fmprb_rgamma(a, a, prec2);
/* check 1/gamma(z+1) = 1/gamma(z)/z */
fmprb_div(b, b, a, prec1);
fmprb_add_ui(c, a, 1, prec1);
fmprb_rgamma(c, c, prec1);
if (!fmprb_equal(a, c))
if (!fmprb_overlaps(b, c))
{
printf("FAIL: aliasing\n\n");
printf("FAIL: functional equation\n\n");
printf("a = "); fmprb_print(a); printf("\n\n");
printf("b = "); fmprb_print(b); printf("\n\n");
printf("c = "); fmprb_print(c); printf("\n\n");
abort();
}
@ -79,3 +85,4 @@ int main()
printf("PASS\n");
return EXIT_SUCCESS;
}

12
gamma.h
View file

@ -29,6 +29,7 @@
#include <math.h>
#include "flint.h"
#include "fmprb.h"
#include "fmpcb.h"
extern __thread long * gamma_taylor_bound_mag_cache;
extern __thread fmpr_struct * gamma_taylor_bound_ratio_cache;
@ -70,5 +71,16 @@ void gamma_taylor_eval_series_fmprb(fmprb_t y, const fmprb_t x, long prec);
void gamma_taylor_fmprb(fmprb_t y, const fmprb_t x, long prec);
void gamma_stirling_choose_param(int * reflect, long * r, long * n, double x, double y, double beta, int allow_reflection, long prec);
void gamma_stirling_coeff(fmprb_t b, ulong k, long prec);
void gamma_stirling_bound_remainder(fmpr_t err, const fmprb_t z, long n);
void gamma_stirling_eval_series_fmprb(fmprb_t s, const fmprb_t z, long nterms, long prec);
void gamma_stirling_eval_series_fmpcb(fmpcb_t s, const fmpcb_t z, long nterms, long prec);
#endif

118
gamma/stirling_choose_param.c Normal file
View file

@ -0,0 +1,118 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2013 Fredrik Johansson
******************************************************************************/
#include "gamma.h"
#include "bernoulli.h"
double
gamma_stirling_arg_factor(double x, double y)
{
double t;
if (y == 0)
{
if (x <= 0)
return 1e30;
else
return 0.0;
}
else
{
t = 1.0 / cos(0.5 * atan2(y, x));
return 2.0 * log(t) * (1. / log(2));
}
}
void
gamma_stirling_choose_param(int * reflect, long * r, long * n,
double x, double y, double beta, int allow_reflection, long prec)
{
double w, logz, argz, boundn, prevboundn;
int i;
long rr, nn;
/* use reflection formula if very negative */
if (x < -5.0 && allow_reflection)
{
*reflect = 1;
x = 1.0 - x;
}
else
{
*reflect = 0;
}
/* argument reduction until |z| >= w */
w = FLINT_MAX(1.0, beta * prec);
rr = 0;
while (x < 1.0 || x*x + y*y < w*w)
{
x++;
rr++;
}
/* this should succeed on the first try if the selection strategy is sound
(i.e. if GAMMA_STIRLING_BETA is large enough) */
for (i = 0; i < 1; i++)
{
/* log2(z) */
logz = 0.5 * log(x*x + y*y) * 1.44269504088896341;
/* log2(1/cos(0.5*arg(z))^2) */
argz = gamma_stirling_arg_factor(x, y);
nn = 1;
prevboundn = 1e300;
while (1)
{
boundn = bernoulli_bound_2exp_si(2 * nn) - (2 * nn-1) * logz + argz * nn;
/* success */
if (boundn <= -prec)
{
*r = rr;
*n = nn;
return;
}
/* if the term magnitude does not decrease, try larger r */
if (boundn > 1)
{
int r2 = 4 + rr * 0.5;
x += r2;
rr += r2;
break;
}
prevboundn = boundn;
nn++;
}
}
printf("exception: gamma_stirling_choose_param failed to converge\n");
abort();
}

39
gamma/stirling_coeff.c Normal file
View file

@ -0,0 +1,39 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2013 Fredrik Johansson
******************************************************************************/
#include "gamma.h"
#include "bernoulli.h"
void
gamma_stirling_coeff(fmprb_t b, ulong k, long prec)
{
fmpz_t d;
fmpz_init(d);
fmprb_set_round_fmpz(b, fmpq_numref(bernoulli_cache + 2 * k), prec);
fmpz_mul2_uiui(d, fmpq_denref(bernoulli_cache + 2 * k), 2 * k, 2 * k - 1);
fmprb_div_fmpz(b, b, d, prec);
fmpz_clear(d);
}

161
gamma/stirling_eval_series_fmpcb.c Normal file
View file

@ -0,0 +1,161 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2012 Fredrik Johansson
******************************************************************************/
#include <math.h>
#include "gamma.h"
#include "bernoulli.h"
/*
B_(2n) / (2n (2n-1)) / |z|^(2n-1) * (1/cos(0.5*arg(z))^(2n))
*/
void
gamma_stirling_bound_remainder_fmpcb(fmpr_t err, const fmpcb_t z, long n)
{
fmpr_t t;
fmprb_t b;
if (fmprb_contains_zero(fmpcb_imagref(z)) &&
fmprb_contains_nonpositive(fmpcb_realref(z)))
{
fmpr_pos_inf(err);
return;
}
/* bound 1 / |z|^(2n-1) */
fmpcb_get_abs_lbound_fmpr(err, z, FMPRB_RAD_PREC);
if (fmpr_is_zero(err))
{
fmpr_pos_inf(err);
return;
}
fmpr_ui_div(err, 1, err, FMPRB_RAD_PREC, FMPR_RND_UP);
fmpr_pow_sloppy_ui(err, err, 2 * n - 1, FMPRB_RAD_PREC, FMPR_RND_UP);
/* bound coefficient */
fmprb_init(b);
fmpr_init(t);
gamma_stirling_coeff(b, n, FMPRB_RAD_PREC);
fmprb_get_abs_ubound_fmpr(t, b, FMPRB_RAD_PREC);
fmpr_mul(err, err, t, FMPRB_RAD_PREC, FMPR_RND_UP);
/* bound 1/cos(0.5*arg(z))^(2n) */
fmpcb_arg(b, z, FMPRB_RAD_PREC);
fmprb_mul_2exp_si(b, b, -1);
fmprb_cos(b, b, FMPRB_RAD_PREC);
fmprb_get_abs_lbound_fmpr(t, b, FMPRB_RAD_PREC);
fmpr_ui_div(t, 1, t, FMPRB_RAD_PREC, FMPR_RND_UP);
fmpr_pow_sloppy_ui(t, t, 2 * n, FMPRB_RAD_PREC, FMPR_RND_UP);
fmpr_mul(err, err, t, FMPRB_RAD_PREC, FMPR_RND_UP);
fmprb_clear(b);
fmpr_clear(t);
}
void
gamma_stirling_eval_series_fmpcb(fmpcb_t s, const fmpcb_t z, long nterms, long prec)
{
fmpcb_t t, u, w, v;
fmprb_t b;
fmpr_t err;
long k, term_prec;
double z_mag, term_mag;
fmpcb_init(t);
fmpcb_init(u);
fmpcb_init(w);
fmpcb_init(v);
fmprb_init(b);
fmpcb_log(w, z, prec);
bernoulli_cache_compute(2 * (nterms + 1));
nterms = FLINT_MAX(nterms, 1);
fmpcb_zero(s);
if (nterms > 1)
{
fmpcb_inv(t, z, prec);
fmpcb_mul(u, t, t, prec);
z_mag = fmpr_get_d(fmprb_midref(fmpcb_realref(w)), FMPR_RND_UP) * 1.44269504088896;
for (k = nterms - 1; k >= 1; k--)
{
term_mag = bernoulli_bound_2exp_si(2 * k);
term_mag -= (2 * k - 1) * z_mag;
term_prec = prec + term_mag;
term_prec = FLINT_MIN(term_prec, prec);
term_prec = FLINT_MAX(term_prec, 10);
gamma_stirling_coeff(b, k, term_prec);
if (prec > 2000)
{
fmpcb_set_round(v, u, term_prec);
fmpcb_mul(s, s, v, term_prec);
}
else
{
fmpcb_mul(s, s, u, term_prec);
}
fmprb_add(fmpcb_realref(s), fmpcb_realref(s), b, term_prec);
}
fmpcb_mul(s, s, t, prec);
}
/* remainder bound */
fmpr_init(err);
gamma_stirling_bound_remainder_fmpcb(err, z, nterms);
fmprb_add_error_fmpr(fmpcb_realref(s), err);
fmprb_add_error_fmpr(fmpcb_imagref(s), err);
fmpr_clear(err);
/* (z-0.5)*log(z) - z + log(2*pi)/2 */
fmprb_one(b);
fmprb_mul_2exp_si(b, b, -1);
fmprb_set(fmpcb_imagref(t), fmpcb_imagref(z));
fmprb_sub(fmpcb_realref(t), fmpcb_realref(z), b, prec);
fmpcb_mul(t, w, t, prec);
fmpcb_add(s, s, t, prec);
fmpcb_sub(s, s, z, prec);
fmprb_const_log_sqrt2pi(b, prec);
fmprb_add(fmpcb_realref(s), fmpcb_realref(s), b, prec);
fmpcb_clear(t);
fmpcb_clear(u);
fmpcb_clear(w);
fmpcb_clear(v);
fmprb_clear(b);
}

145
gamma/stirling_eval_series_fmprb.c Normal file
View file

@ -0,0 +1,145 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2012 Fredrik Johansson
******************************************************************************/
#include <math.h>
#include "gamma.h"
#include "bernoulli.h"
void
gamma_stirling_bound_remainder(fmpr_t err, const fmprb_t z, long n)
{
if (fmprb_contains_nonpositive(z))
{
fmpr_pos_inf(err);
}
else
{
fmpr_t t;
fmprb_t b;
fmpr_init(t);
fmprb_init(b);
fmpr_sub(t, fmprb_midref(z), fmprb_radref(z), FMPRB_RAD_PREC, FMPR_RND_FLOOR);
if (fmpr_sgn(t) <= 0)
{
fmpr_pos_inf(err);
}
else
{
fmpr_ui_div(t, 1, t, FMPRB_RAD_PREC, FMPR_RND_UP);
fmpr_pow_sloppy_ui(t, t, 2 * n - 1, FMPRB_RAD_PREC, FMPR_RND_UP);
gamma_stirling_coeff(b, n, FMPRB_RAD_PREC);
fmprb_get_abs_ubound_fmpr(err, b, FMPRB_RAD_PREC);
fmpr_mul(err, err, t, FMPRB_RAD_PREC, FMPR_RND_UP);
}
fmpr_clear(t);
fmprb_clear(b);
}
}
void
gamma_stirling_eval_series_fmprb(fmprb_t s, const fmprb_t z, long nterms, long prec)
{
fmprb_t t, u, b, w, v;
fmpr_t err;
long k, term_prec;
double z_mag, term_mag;
fmprb_init(t);
fmprb_init(u);
fmprb_init(b);
fmprb_init(w);
fmprb_init(v);
fmprb_log(w, z, prec);
bernoulli_cache_compute(2 * (nterms + 1));
nterms = FLINT_MAX(nterms, 1);
fmprb_zero(s);
if (nterms > 1)
{
fmprb_ui_div(t, 1UL, z, prec);
fmprb_mul(u, t, t, prec);
z_mag = fmpr_get_d(fmprb_midref(w), FMPR_RND_UP) * 1.44269504088896;
for (k = nterms - 1; k >= 1; k--)
{
term_mag = bernoulli_bound_2exp_si(2 * k);
term_mag -= (2 * k - 1) * z_mag;
term_prec = prec + term_mag;
term_prec = FLINT_MIN(term_prec, prec);
term_prec = FLINT_MAX(term_prec, 10);
gamma_stirling_coeff(b, k, term_prec);
if (prec > 2000)
{
fmprb_set_round(v, u, term_prec);
fmprb_mul(s, s, v, term_prec);
}
else
{
fmprb_mul(s, s, u, term_prec);
}
fmprb_add(s, s, b, term_prec);
}
fmprb_mul(s, s, t, prec);
}
/* remainder bound */
fmpr_init(err);
gamma_stirling_bound_remainder(err, z, nterms);
fmprb_add_error_fmpr(s, err);
fmpr_clear(err);
/* (z-0.5)*log(z) - z + log(2*pi)/2 */
fmprb_set_ui(t, 1);
fmprb_mul_2exp_si(t, t, -1);
fmprb_sub(t, z, t, prec);
fmprb_mul(t, w, t, prec);
fmprb_add(s, s, t, prec);
fmprb_sub(s, s, z, prec);
fmprb_const_log_sqrt2pi(t, prec);
fmprb_add(s, s, t, prec);
fmprb_clear(t);
fmprb_clear(u);
fmprb_clear(b);
fmprb_clear(w);
fmprb_clear(v);
}

6
gamma/taylor_precompute.c
View file

@ -50,10 +50,16 @@ gamma_taylor_precompute(long num, long prec)
fmprb_poly_t A;
long i;
_fmprb_vec_clear(gamma_taylor_coeffs, gamma_taylor_num);
printf("precomputing: num %ld, prec %ld bits\n", num, prec);
num = FLINT_MAX(gamma_taylor_num * 1.5, num);
prec = FLINT_MAX(gamma_taylor_prec * 1.5, prec);
prec = FLINT_MAX(prec, 64);
printf("precomputing (2): num %ld, prec %ld bits\n", num, prec);
wp = prec * 1.01 + 32;
/* TODO: cleanup */