mirror of
https://github.com/vale981/arb
synced 2025-03-06 18:01:39 -05:00
172 lines
7.1 KiB
ReStructuredText
172 lines
7.1 KiB
ReStructuredText
**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
|
|
--------------------------------------------------------------------------------
|
|
|
|
.. function :: long gamma_taylor_bound_mag(long n)
|
|
|
|
Letting `1/\Gamma(x) = \sum_{n=1}^{\infty} c_n x^n`, returns an integer
|
|
`p` such that `|c_n| \le 2^{-p}`.
|
|
Uses the bound `\log |c_n| \le b(n) = ((\pi-1)n + (3-5n) \log(n/6))/6`
|
|
(to be published).
|
|
|
|
.. function :: void gamma_taylor_bound_ratio(fmpr_t r, long n)
|
|
|
|
Sets `r` to a bound for `B(n+1) / B(n)` where `B(n) = \exp b(n)`.
|
|
We use `r(n) = 3n^{-5/6}`.
|
|
|
|
.. function :: void gamma_taylor_bound_remainder(fmpr_t err, const fmpr_t z, long n)
|
|
|
|
Given a nonnegative `z` and `n \ge 1`, computes a bound for
|
|
`\sum_{k=n}^{\infty} |c_{k+1}| z^k`. This is given by
|
|
`B(n+1) z^n / (1 - r(n+1) z)`, provided that the denominator is
|
|
positive (otherwise the bound is infinity).
|
|
|
|
.. function :: void gamma_taylor_fmprb(fmprb_t y, const fmprb_t x, long prec)
|
|
|
|
Sets `y = \Gamma(x)`, computed by translating `x` to the interval
|
|
`[-1/2,1/2]` and then evaluating the Taylor series for the
|
|
reciprocal gamma function.
|
|
|
|
Assumes that the absolute value of the midpoint of `x` is not too large
|
|
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 with `z = x+yi`
|
|
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, in the real case, `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.
|
|
|
|
|
|
Rising factorials
|
|
--------------------------------------------------------------------------------
|
|
|
|
.. function :: void gamma_rising_fmprb_ui_bsplit_simple(fmprb_t y, const fmprb_t x, ulong n, long prec)
|
|
|
|
.. function :: void gamma_rising_fmprb_ui_bsplit_eight(fmprb_t y, const fmprb_t x, ulong n, long prec)
|
|
|
|
.. function :: void gamma_rising_fmprb_ui_bsplit_rectangular(fmprb_t y, const fmprb_t x, ulong n, ulong step, long prec)
|
|
|
|
.. function :: void gamma_rising_fmprb_ui_bsplit(fmprb_t y, const fmprb_t x, ulong n, long prec)
|
|
|
|
.. function :: void gamma_rising_fmpcb_ui_bsplit_simple(fmpcb_t y, const fmpcb_t x, ulong n, long prec)
|
|
|
|
.. function :: void gamma_rising_fmpcb_ui_bsplit_eight(fmpcb_t y, const fmpcb_t x, ulong n, long prec)
|
|
|
|
.. function :: void gamma_rising_fmpcb_ui_bsplit_rectangular(fmpcb_t y, const fmpcb_t x, ulong n, ulong step, long prec)
|
|
|
|
.. function :: void gamma_rising_fmpcb_ui_bsplit(fmpcb_t y, const fmpcb_t x, ulong n, long prec)
|
|
|
|
Sets `y` to the rising factorial `x (x+1) (x+2) \cdots (x+n-1)`,
|
|
computed using binary splitting.
|
|
|
|
The different versions of this function process the basecase differently.
|
|
The *simple* version simply multiplies together several factors
|
|
one after another.
|
|
|
|
The *eight* version processes eight factors at a time using the formula
|
|
|
|
.. math ::
|
|
|
|
x(x+1)\cdots(x+7) = (28 + 98x + 63x^2 + 14x^3 + x^4)^2 - 16 (7+2x)^2,
|
|
|
|
replacing 7 full-precision multiplications with 3 squarings,
|
|
1 multiplication, and several linear operations ([CP2005]_, page 316).
|
|
Empirically, if `x` is a full-precision number, this is about twice as
|
|
fast as the *simple* version at high precision. Numerical stability is
|
|
slightly worse.
|
|
|
|
The *rectangular* version processes *step* factors at a time by
|
|
expanding the polynomial `f(t) = t (t+1) (t+2) \cdots (t+\mathrm{step}-1)`
|
|
and evaluating each factor `f(x + \mathrm{step} \, k)`
|
|
using rectangular splitting. At very high precision, if `x` is a
|
|
full-precision number, this asymptotically reduces the number of
|
|
full-precision multiplications required.
|
|
|
|
The function *gamma_rising_fmprb_ui_bsplit* automatically chooses
|
|
an algorithm depending on the inputs.
|
|
|
|
.. function :: void gamma_rising_fmprb_ui_multipoint(fmprb_t f, const fmprb_t c, ulong n, long prec)
|
|
|
|
Sets `y` to the rising factorial `x (x+1) (x+2) \cdots (x+n-1)`,
|
|
computed using fast multipoint evaluation. This only requires
|
|
`O(n^{1/2+\varepsilon})` multiplications, but has high overhead
|
|
and poor numerical stability (adding `O(n)` guard bits to the input
|
|
might be necessary to achieve full accuracy). It can be expected to
|
|
be faster than the binary splitting algorithm if the input is a
|
|
full-precision number, the precision is at least 100000 bits,
|
|
and *n* is of the same order of magnitude as (perhaps slightly
|
|
smaller than) the number of bits.
|
|
|