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

move hypgeom code to its own module, +improvements and docs

Browse source
This commit is contained in:
Fredrik Johansson 2013-01-23 15:54:38 +01:00
parent 2a68e4bcb2
commit 3f0901e71a
15 changed files with 510 additions and 198 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
Makefile.in
View file

@ -102,5 +102,5 @@ build/%.lo: %.c
build/%.o: %.c
$(CC) -fPIC $(CFLAGS) $(INCS) -c $< -o $@
BUILD_DIRS = fmpr fmprb fmprb_poly fmprb_mat fmpz_holonomic fmpcb fmpcb_poly fmpcb_mat bernoulli
BUILD_DIRS = fmpr fmprb fmprb_poly fmprb_mat fmpz_holonomic fmpcb fmpcb_poly fmpcb_mat bernoulli hypgeom

2
doc/source/bernoulli.rst
View file

@ -1,4 +1,4 @@
**bernoulli.h** -- support code for Bernoulli numbers
**bernoulli.h** -- support for Bernoulli numbers
===============================================================================
Generation of Bernoulli numbers

10
doc/source/fmpcb.rst
View file

@ -371,11 +371,11 @@ Special functions
.. math ::
\begin{align}
\sum_{k=N}^U f(k) = & \int_N^U f(t) dt + \frac{1}{2} \left(f(N) + f(U)\right) \\
& + \sum_{k=1}^{M} \frac{B_{2k}}{(2k)!} \left( f^{(2k-1)}(U) - f^{(2k-1)}(N) \right) \\
& - \int_N^U \frac{\tilde B_{2M}(t)}{(2M)!} f^{(2M)}(t) dt
\end{align}
\sum_{k=N}^U f(k) = \int_N^U f(t) dt + \frac{1}{2} \left(f(N) + f(U)\right)
+ \sum_{k=1}^{M} \frac{B_{2k}}{(2k)!} \left( f^{(2k-1)}(U) - f^{(2k-1)}(N) \right)
- \int_N^U \frac{\tilde B_{2M}(t)}{(2M)!} f^{(2M)}(t) dt
where `f` is a sufficiently differentiable function (for example,
analytic), `B_n` is a Bernoulli number, and

195
doc/source/hypgeom.rst Normal file
View file

@ -0,0 +1,195 @@
**hypgeom.h** -- support for hypergeometric series
===============================================================================
This module provides functions for high-precision evaluation of series
of the form
.. math ::
\sum_{k=0}^{n-1} \frac{A(k)}{B(k)} \prod_{j=1}^k \frac{P(k)}{Q(k)} z^k
where `A, B, P, Q` are polynomials. The present version only supports
`A, B, P, Q \in \mathbb{Z}[k]` (represented using the
FLINT *fmpz_poly_t* type). This module also provides functions
for high-precision evaluation of infinite series (`n \to \infty`),
with automatic, rigorous error bounding.
Note that we can standardize to `A = B = 1` by
setting `\tilde P(k) = P(k) A(k) B(k-1), \tilde Q(k) = Q(k) A(k-1) B(k)`.
However, separating out `A` and `B` is convenient and improves
efficiency during evaluation.
Strategy for error bounding
-------------------------------------------------------------------------------
We wish to evaluate `S(z) = \sum_{k=0}^{\infty} T(k) z^k` where `T(k)`
satisfies `T(0) = 1` and
.. math ::
T(k) = R(k) T(k-1) = \left( \frac{P(k)}{Q(k)} \right) T(k-1)
for given polynomials
.. math ::
P(k) = a_p k^p + a_{p-1} k^{p-1} + \ldots a_0
Q(k) = b_q k^q + b_{q-1} k^{q-1} + \ldots b_0.
For convergence, we require `p < q`, or `p = q` with `|z| |a_p| < |b_q|`.
We also assume that `P(k)` and `Q(k)` have no roots among the positive
integers (if there are positive integer roots, the sum is either finite
or undefined). With these conditions satisfied, our goal is to find a
parameter `n \ge 0` such that
.. math ::
\left\lvert \sum_{k=n}^{\infty} T(k) z^k \right\rvert \le 2^{-d}.
We can rewrite the hypergeometric term ratio as
.. math ::
z R(k) = z \frac{P(k)}{Q(k)} =
z \left( \frac{a_p}{b_q} \right) \frac{1}{k^{q-p}} F(k)
where
.. math ::
F(k) = \frac{
1 + \tilde a_{1} / k + \tilde a_{2} / k^2 + \ldots + \tilde a_q / k^p
}{
1 + \tilde b_{1} / k + \tilde b_{2} / k^2 + \ldots + \tilde b_q / k^q
} = 1 + O(1/k)
and where `\tilde a_i = a_{p-i} / a_p`, `\tilde b_i = b_{q-i} / b_q`.
Next, we define
.. math ::
C = \max_{1 \le i \le p} |\tilde a_i|^{(1/i)},
\quad D = \max_{1 \le i \le q} |\tilde b_i|^{(1/i)}.
Now, if `k > C`, the magnitude of the numerator of `F(k)` is
bounded from above by
.. math ::
1 + \sum_{i=1}^p \left(\frac{C}{k}\right)^i \le 1 + \frac{C}{k-C}
and if `k > 2D`, the magnitude of the denominator of `F(k)` is bounded
from below by
.. math ::
1 - \sum_{i=1}^q \left(\frac{D}{k}\right)^i \ge 1 + \frac{D}{D-k}.
Putting the inequalities together gives the following bound,
valid for `k > K = \max(C, 2D)`:
.. math ::
|F(k)| \le \frac{k (k-D)}{(k-C)(k-2D)} = \left(1 + \frac{C}{k-C} \right)
\left(1 + \frac{D}{k-2D} \right).
Let `r = q-p` and `\tilde z = |z a_p / b_q|`. Assuming
`k > \max(C, 2D, {\tilde z}^{1/r})`, we have
.. math ::
|z R(k)| \le G(k) = \frac{\tilde z F(k)}{k^r}
where `G(k)` is monotonically decreasing. Now we just need to find an
`n` such that `G(n) < 1` and for which `|T(n)| / (1 - G(n)) \le 2^{-d}`.
This can be done by computing a floating-point guess for `n` then
trying successively larger values.
This strategy leaves room for some improvement. For example, if
`\tilde b_1` is positive and large, the bound `B` becomes very pessimistic
(a larger positive `\tilde b_1` causes faster convergence,
not slower convergence).
Types, macros and constants
-------------------------------------------------------------------------------
.. type:: hypgeom_struct
.. type:: hypgeom_t
Stores polynomials *A*, *B*, *P*, *Q* and precomputed bounds,
representing a fixed hypergeometric series.
Memory management
-------------------------------------------------------------------------------
.. function:: void hypgeom_init(hypgeom_t hyp)
.. function:: void hypgeom_clear(hypgeom_t hyp)
Error bounding
-------------------------------------------------------------------------------
.. function:: long hypgeom_estimate_terms(const fmpr_t z, int r, long d)
Computes an approximation of the largest `n` such
that `|z|^n/(n!)^r = 2^{-d}`, giving a first-order estimate of the
number of terms needed to approximate the sum of a hypergeometric
series of weight `r \ge 0` and argument `z` to an absolute
precision of `d \ge 0` bits. If `r = 0`, the direct solution of the
equation is given by `n = (\log(1-z) - d \log 2) / \log z`.
If `r > 0`, using `\log n! \approx n \log n - n` gives an equation
that can be solved in terms of the Lambert *W*-function as
`n = (d \log 2) / (r\,W\!(t))` where
`t = (d \log 2) / (e r z^{1/r})`.
The evaluation is done using double precision arithmetic.
The function aborts if the computed value of `n` is greater
than or equal to LONG_MAX / 2.
.. function:: long hypgeom_bound(fmpr_t error, int r, long C, long D, long K, const fmpr_t TK, const fmpr_t z, long prec)
Computes a truncation parameter sufficient to achieve *prec* bits
of absolute accuracy, according to the strategy described above.
The input consists of `r`, `C`, `D`, `K`, precomputed bound for `T(K)`,
and `\tilde z = z (a_p / b_q)`, such that for `k > K`, the hypergeometric
term ratio is bounded by
.. math ::
\frac{\tilde z}{k^r} \frac{k(k-D)}{(k-C)(k-2D)}.
Given this information, we compute a `\varepsilon` and an
integer `n` such that
`\left| \sum_{k=n}^{\infty} T(k) \right| \le \varepsilon \le 2^{-\mathrm{prec}}`.
The output variable *error* is set to the value of `\varepsilon`,
and `n` is returned.
.. function:: void hypgeom_precompute(hypgeom_t hyp)
Precomputes the bounds data `C`, `D`, `K` and an upper bound for `T(K)`.
Summation
-------------------------------------------------------------------------------
.. function:: void hypgeom_fmprb_sum(fmprb_t P, fmprb_t Q, const hypgeom_t hyp, const long n, long prec)
Computes `P, Q` such that `P / Q = \sum_{k=0}^{n-1} T(k)` where `T(k)`
is defined by *hyp*,
using binary splitting and a working precision of *prec* bits.
.. function:: void hypgeom_fmprb_infsum(fmprb_t P, fmprb_t Q, hypgeom_t hyp, long tol, long prec)
Computes `P, Q` such that `P / Q = \sum_{k=0}^{\infty} T(k)` where `T(k)`
is defined by *hyp*, using binary splitting and
working precision of *prec* bits.
The number of terms is chosen automatically to bound the
truncation error by at most `2^{-\mathrm{tol}}`.
The bound for the truncation error is included in the output
as part of *P*.

5
doc/source/index.rst
View file

@ -26,7 +26,7 @@ Module documentation
::::::::::::::::::::
.. toctree::
:maxdepth: 2
:maxdepth: 1
fmpr.rst
fmprb.rst
@ -35,8 +35,9 @@ Module documentation
fmpcb.rst
fmpcb_poly.rst
fmpcb_mat.rst
fmpz_holonomic.rst
bernoulli.rst
hypgeom.rst
fmpz_holonomic.rst
Credits and references
::::::::::::::::::::::::

34
fmprb.h
View file

@ -288,40 +288,6 @@ void fmprb_rfac_ui_multipoint(fmprb_t y, const fmprb_t x, ulong n, long prec);
void fmprb_fib_fmpz(fmprb_t f, const fmpz_t n, long prec);
void fmprb_fib_ui(fmprb_t f, ulong n, long prec);
typedef struct
{
fmpz_poly_t A;
fmpz_poly_t B;
fmpz_poly_t P;
fmpz_poly_t Q;
/* precomputation data */
int have_precomputed;
long r;
long boundA;
long boundB;
long boundK;
fmpr_t MK;
}
hypgeom_struct;
typedef hypgeom_struct hypgeom_t[1];
void hypgeom_init(hypgeom_t hyp);
void hypgeom_clear(hypgeom_t hyp);
void hypgeom_precompute(hypgeom_t hyp);
long hypgeom_bound(fmpr_t error, int r,
long A, long B, long K, const fmpr_t TK, const fmpr_t z, long prec);
void fmprb_hypgeom_sum(fmprb_t P, fmprb_t Q, const hypgeom_t hyp, const long n, long prec);
void fmprb_hypgeom_infsum(fmprb_t P, fmprb_t Q, hypgeom_t hyp, long target_prec, long prec);
void fmprb_const_pi_chudnovsky(fmprb_t x, long prec);
void fmprb_const_pi(fmprb_t x, long prec);

1
fmprb/const_pi_chudnovsky.c
View file

@ -24,6 +24,7 @@
******************************************************************************/
#include "fmprb.h"
#include "hypgeom.h"
void
fmprb_const_pi_chudnovsky(fmprb_t s, long prec)

2
fmprb/const_zeta3_bsplit.c
View file

@ -24,7 +24,7 @@
******************************************************************************/
#include "fmprb.h"
#include "hypgeom.h"
void
fmprb_const_zeta3_bsplit(fmprb_t s, long prec)

67
hypgeom.h Normal file
View file

@ -0,0 +1,67 @@
/*=============================================================================
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
******************************************************************************/
#ifndef HYPGEOM_H
#define HYPGEOM_H
#include "fmprb.h"
#include "fmpz_poly.h"
typedef struct
{
fmpz_poly_t A;
fmpz_poly_t B;
fmpz_poly_t P;
fmpz_poly_t Q;
/* precomputation data */
int have_precomputed;
long r;
long boundC;
long boundD;
long boundK;
fmpr_t MK;
}
hypgeom_struct;
typedef hypgeom_struct hypgeom_t[1];
void hypgeom_init(hypgeom_t hyp);
void hypgeom_clear(hypgeom_t hyp);
void hypgeom_precompute(hypgeom_t hyp);
long hypgeom_estimate_terms(const fmpr_t z, int r, long prec);
long hypgeom_bound(fmpr_t error, int r,
long C, long D, long K, const fmpr_t TK, const fmpr_t z, long prec);
void fmprb_hypgeom_sum(fmprb_t P, fmprb_t Q, const hypgeom_t hyp, const long n, long prec);
void fmprb_hypgeom_infsum(fmprb_t P, fmprb_t Q, hypgeom_t hyp, long target_prec, long prec);
#endif

43
hypgeom/Makefile Normal file
View file

@ -0,0 +1,43 @@
SOURCES = $(wildcard *.c)
OBJS = $(patsubst %.c, $(BUILD_DIR)/%.o, $(SOURCES))
LIB_OBJS = $(patsubst %.c, $(BUILD_DIR)/%.lo, $(SOURCES))
TEST_SOURCES = $(wildcard test/*.c)
PROF_SOURCES = $(wildcard profile/*.c)
TUNE_SOURCES = $(wildcard tune/*.c)
TESTS = $(patsubst %.c, %, $(TEST_SOURCES))
PROFS = $(patsubst %.c, %, $(PROF_SOURCES))
TUNE = $(patsubst %.c, %, $(TUNE_SOURCES))
all: $(OBJS)
library: $(LIB_OBJS)
profile:
$(foreach prog, $(PROFS), $(CC) -O2 -std=c99 $(INCS) $(prog).c ../profiler.o -o $(BUILD_DIR)/$(prog) $(LIBS) -lflint;)
tune: $(TUNE_SOURCES)
$(foreach prog, $(TUNE), $(CC) -O2 -std=c99 $(INCS) $(prog).c -o $(BUILD_DIR)/$(prog) $(LIBS) -lflint;)
$(BUILD_DIR)/%.o: %.c
$(CC) $(CFLAGS) -c $(INCS) $< -o $@
$(BUILD_DIR)/%.lo: %.c
$(CC) -fPIC $(CFLAGS) $(INCS) -c $< -o $@
clean:
rm -rf $(BUILD_DIR)
check: library
$(foreach prog, $(TESTS), $(CC) $(CFLAGS) $(INCS) $(prog).c -o $(BUILD_DIR)/$(prog) $(LIBS) -lflint;)
$(foreach prog, $(TESTS), $(BUILD_DIR)/$(prog);)
.PHONY: profile clean check all

254
fmprb/hypgeom_bound.c → hypgeom/bound.c
View file

@ -24,15 +24,12 @@
******************************************************************************/
#include <math.h>
#include "fmprb.h"
#include "double_extras.h"
#include "hypgeom.h"
double fmpr_get_d(const fmpr_t x);
#define LOG2 0.69314718055994530942
#define EXP1 2.7182818284590452354
double d_root(double x, int r)
static __inline__ double d_root(double x, int r)
{
if (r == 1)
return x;
@ -41,46 +38,6 @@ double d_root(double x, int r)
return pow(x, 1. / r);
}
/*
Estimate the truncation point to obtain accuracy 2^(-prec) with the
hypergeometric series |z|^k / (k!)^r.
*/
long
estimate_nterms(double z, int r, long prec)
{
double y;
z = fabs(z);
if (z == 0)
return 1;
if (r == 0)
{
if (z >= 1)
{
printf("z must be smaller than 1\n");
abort();
}
y = (log(1-z) - prec * LOG2) / log(z) + 1;
}
else
{
/* Solve k*log(z) - r*(k*log(k)-k) = -prec*log(2) */
y = prec * LOG2 / (d_root(z, r) * EXP1 * r);
y = prec * LOG2 / (r * d_lambertw(y)) + 1;
}
if (y >= LONG_MAX / 2)
{
printf("error: series will converge too slowly\n");
abort();
}
return y;
}
void
fmpr_gamma_ui_lbound(fmpr_t x, ulong n, long prec)
{
@ -157,151 +114,153 @@ fmpr_gamma_ui_ubound(fmpr_t x, ulong n, long prec)
}
}
/* FIXME: assumes no overflow when computing n + r */
static void
fmpr_rfac_uiui_ubound(fmpr_t x, ulong n, ulong r, long prec)
/* FIXME: use more than doubles for this */
long hypgeom_root_bound(const fmpr_t z, int r)
{
if (r == 0)
{
fmpr_one(x);
}
else if (r == 1)
{
fmpr_set_ui(x, n);
return 0;
}
else
{
fmpr_t t;
fmpr_init(t);
fmpr_gamma_ui_ubound(x, n + r, prec);
fmpr_gamma_ui_lbound(t, n, prec);
fmpr_div(x, x, t, prec, FMPR_RND_UP);
fmpr_clear(t);
double zd;
zd = fmpr_get_d(z);
return d_root(zd, r) + 2;
}
}
static void
fmpr_rfac_uiui_lbound(fmpr_t x, ulong n, ulong r, long prec)
{
if (r == 0)
{
fmpr_one(x);
}
else if (r == 1)
{
fmpr_set_ui(x, n);
}
else
{
fmpr_t t;
fmpr_init(t);
fmpr_gamma_ui_lbound(x, n + r, prec);
fmpr_gamma_ui_ubound(t, n, prec);
fmpr_div(x, x, t, prec, FMPR_RND_DOWN);
fmpr_clear(t);
}
}
/*
Given T(K), compute bound for T(n) z^n.
The general term T(k) is z^k / (k!)^r * R(k) where R(k) = 1 + O(1/k).
We need to multiply by
We have precomputed integers A, B, K such that for all k > K,
|R(k)| <= F(k) = k(k-B)/((k-A)(k-2B)) = (1+A/(k-A))(1+B/(k-2B)).
z^n * 1/rf(K+1,m)^r * (rf(K+1,m)/rf(K+1-A,m)) * (rf(K+1-B,m)/rf(K+1-2B,m))
where m = n - K. This is equal to
z^n *
(K+A)! (K-2B)! (K-B+m)!
----------------------- * ((K+m)! / K!)^(1-r)
(K-B)! (K-A+m)! (K-2B+m)!
*/
void
hypgeom_term_bound(fmpr_t Tn, const fmpr_t TK, long K, long A, long B, int r, const fmpr_t z, long n, long wp)
{
fmpr_t t, u, num, den;
long m;
fmpr_init(t);
fmpr_init(u);
fmpr_init(num);
fmpr_init(den);
m = n - K;
if (m < 0)
abort();
/* TK * z^n */
fmpr_pow_sloppy_ui(t, z, n, wp, FMPR_RND_UP);
fmpr_mul(num, TK, t, wp, FMPR_RND_UP);
/* (K+A)! (K-2B)! (K-B+m)!, upper bounding */
fmpr_gamma_ui_ubound(t, K+A+1, wp);
fmpr_mul(num, num, t, wp, FMPR_RND_UP);
fmpr_gamma_ui_ubound(t, K-2*B+1, wp);
fmpr_mul(num, num, t, wp, FMPR_RND_UP);
fmpr_gamma_ui_ubound(t, K-B+m, wp);
fmpr_mul(num, num, t, wp, FMPR_RND_UP);
/* (K-B)! (K-A+m)! (K-2B+m)!, lower bounding */
fmpr_gamma_ui_lbound(den, K-B+1, wp);
fmpr_gamma_ui_lbound(t, K-A+m+1, wp);
fmpr_mul(den, den, t, wp, FMPR_RND_DOWN);
fmpr_gamma_ui_lbound(t, K-2*B+m+1, wp);
fmpr_mul(den, den, t, wp, FMPR_RND_DOWN);
/* ((K+m)! / K!)^(1-r) */
if (r == 0)
{
fmpr_gamma_ui_ubound(t, K+m+1, wp);
fmpr_mul(num, num, t, wp, FMPR_RND_UP);
fmpr_gamma_ui_lbound(t, K+1, wp);
fmpr_mul(den, den, t, wp, FMPR_RND_DOWN);
}
else if (r != 1)
{
fmpr_gamma_ui_ubound(t, K+1, wp);
fmpr_gamma_ui_lbound(u, K+m+1, wp);
fmpr_div(t, t, u, wp, FMPR_RND_UP);
fmpr_pow_sloppy_ui(t, t, r-1, wp, FMPR_RND_UP);
fmpr_mul(num, num, t, wp, FMPR_RND_UP);
}
fmpr_div(Tn, num, den, wp, FMPR_RND_UP);
fmpr_clear(t);
fmpr_clear(u);
fmpr_clear(num);
fmpr_clear(den);
}
long
hypgeom_bound(fmpr_t error, int r,
long A, long B, long K, const fmpr_t TK, const fmpr_t z, long prec)
{
fmpr_t Tn, t, u, one, tol;
fmpr_t Tn, t, u, one, tol, num, den;
long wp = FMPRB_RAD_PREC;
long n;
double zd;
long n, m;
fmpr_init(Tn);
fmpr_init(t);
fmpr_init(u);
fmpr_init(one);
fmpr_init(tol);
fmpr_init(num);
fmpr_init(den);
fmpr_one(one);
fmpr_set_ui_2exp_si(tol, 1UL, -prec);
zd = fmpr_get_d(z);
n = estimate_nterms(zd, r, prec);
/* approximate number of needed terms */
n = hypgeom_estimate_terms(z, r, prec);
/* required for 1 + O(1/k) part to be decreasing */
n = FLINT_MAX(n, K + 1);
/* required for z^k / (k!)^r to be decreasing
(TODO: don't use doubles for this) */
if (r > 0)
{
long nbd = d_root(zd, r) + 2;
n = FLINT_MAX(n, nbd);
}
/* required for z^k / (k!)^r to be decreasing */
m = hypgeom_root_bound(z, r);
n = FLINT_MAX(n, m);
/* We are now sure that |R(k)| is either decreasing or strictly
smaller than 1 for k >= n, which means that we can bound the tail
using a geometric series as soon as soon as |R(k)| < 1 */
/* We now have |R(k)| <= G(k) where G(k) is monotonically decreasing,
and can bound the tail using a geometric series as soon
as soon as G(k) < 1. */
/* Compute an upper bound for T(n) */
/* (z^n) / (n!)^r * TK * [(K+1)(K+2)...(n)] * [(K-B+1)(K-B+2)...(n-B)]
---------------------------------------------
[(K-A+1)(K-A+2)...(n)] * [(K-2B+1)(K-2B+2)...(n-2B)]
*/
/* z^n * TK */
fmpr_pow_sloppy_ui(Tn, z, n, wp, FMPR_RND_UP);
fmpr_mul(Tn, Tn, TK, wp, FMPR_RND_UP);
/* divide by (n!)^r */
if (r != 0)
{
fmpr_gamma_ui_lbound(t, n + 1, wp);
fmpr_ui_div(t, 1UL, t, wp, FMPR_RND_UP);
fmpr_pow_sloppy_ui(t, t, r, wp, FMPR_RND_UP);
fmpr_mul(Tn, Tn, t, wp, FMPR_RND_UP);
}
fmpr_rfac_uiui_ubound(t, K+1, n-K, wp);
fmpr_mul(Tn, Tn, t, wp, FMPR_RND_UP);
fmpr_rfac_uiui_ubound(t, K-B+1, n-K, wp);
fmpr_mul(Tn, Tn, t, wp, FMPR_RND_UP);
fmpr_rfac_uiui_lbound(t, K-A+1, n-K, wp);
fmpr_div(Tn, Tn, t, wp, FMPR_RND_UP);
fmpr_rfac_uiui_lbound(t, K-2*B+1, n-K, wp);
fmpr_div(Tn, Tn, t, wp, FMPR_RND_UP);
/* bound T(n-1) */
hypgeom_term_bound(Tn, TK, K, A, B, r, z, n-1, wp);
while (1)
{
/* bound for term ratio: z * F(n) / n^r */
/* F(n) <= n (n-B) / ((n-A) (n-2B)) */
fmpr_set_ui(t, n);
fmpr_mul_ui(t, t, n - B, wp, FMPR_RND_UP);
fmpr_div_ui(t, t, n - A, wp, FMPR_RND_UP);
fmpr_div_ui(t, t, n - 2*B, wp, FMPR_RND_UP);
fmpr_mul(t, t, z, wp, FMPR_RND_UP);
/* bound R(n) */
fmpr_mul_ui(num, z, n, wp, FMPR_RND_UP);
fmpr_mul_ui(num, num, n - B, wp, FMPR_RND_UP);
fmpr_set_ui(den, n - A);
fmpr_mul_ui(den, den, n - 2*B, wp, FMPR_RND_DOWN);
if (r != 0)
{
fmpr_div_ui(u, one, n, wp, FMPR_RND_UP);
fmpr_pow_sloppy_ui(u, u, r, wp, FMPR_RND_UP);
fmpr_mul(t, t, u, wp, FMPR_RND_UP);
fmpr_set_ui(u, n);
fmpr_pow_sloppy_ui(u, u, r, wp, FMPR_RND_DOWN);
fmpr_mul(den, den, u, wp, FMPR_RND_DOWN);
}
/* bound by geometric series: Tn / (1 - t) */
/* where the term ratio must be < 1 */
fmpr_sub(u, one, t, wp, FMPR_RND_DOWN);
fmpr_div(t, num, den, wp, FMPR_RND_UP);
/* multiply bound for T(n-1) by bound for R(n) to bound T(n) */
fmpr_mul(Tn, Tn, t, wp, FMPR_RND_UP);
/* geometric series termination check */
fmpr_sub(u, one, t, wp, FMPR_RND_DOWN);
if (fmpr_sgn(u) > 0)
{
fmpr_div(u, Tn, u, wp, FMPR_RND_UP);
@ -314,7 +273,6 @@ hypgeom_bound(fmpr_t error, int r,
}
/* move on to next term */
fmpr_mul(Tn, Tn, t, wp, FMPR_RND_UP);
n++;
}
@ -323,6 +281,8 @@ hypgeom_bound(fmpr_t error, int r,
fmpr_clear(u);
fmpr_clear(one);
fmpr_clear(tol);
fmpr_clear(num);
fmpr_clear(den);
return n;
}

79
hypgeom/estimate_terms_d.c Normal file
View file

@ -0,0 +1,79 @@
/*=============================================================================
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 "double_extras.h"
#include "hypgeom.h"
#define LOG2 0.69314718055994530942
#define EXP1 2.7182818284590452354
static __inline__ double d_root(double x, int r)
{
if (r == 1)
return x;
if (r == 2)
return sqrt(x);
return pow(x, 1. / r);
}
double fmpr_get_d(const fmpr_t x);
long
hypgeom_estimate_terms(const fmpr_t z, int r, long prec)
{
double y, t;
t = fmpr_get_d(z);
t = fabs(t);
if (t == 0)
return 1;
if (r == 0)
{
if (t >= 1)
{
printf("z must be smaller than 1\n");
abort();
}
y = (log(1-t) - prec * LOG2) / log(t) + 1;
}
else
{
y = (prec * LOG2) / (d_root(t, r) * EXP1 * r);
y = (prec * LOG2) / (r * d_lambertw(y)) + 1;
}
if (y >= LONG_MAX / 2)
{
printf("error: series will converge too slowly\n");
abort();
}
return y;
}

2
fmprb/hypgeom_init.c → hypgeom/init.c
View file

@ -23,7 +23,7 @@
******************************************************************************/
#include "fmprb.h"
#include "hypgeom.h"
void
hypgeom_init(hypgeom_t hyp)

8
fmprb/hypgeom_precompute.c → hypgeom/precompute.c
View file

@ -24,7 +24,7 @@
******************************************************************************/
#include <math.h>
#include "fmprb.h"
#include "hypgeom.h"
/* Compute a pure ratio P2(k)/Q2(k) for the term
A(k)/B(k) * [P(1)P(2)...P(k)] / [Q(1)Q(2)...Q(k)] */
@ -122,9 +122,9 @@ _hypgeom_precompute(hypgeom_t hyp, const fmpz_poly_t P, const fmpz_poly_t Q)
fmpz_init(t);
hyp->r = fmpz_poly_degree(Q) - fmpz_poly_degree(P);
hyp->boundA = hypgeom_root_norm(P);
hyp->boundB = hypgeom_root_norm(Q);
hyp->boundK = 1 + FLINT_MAX(hyp->boundA, 2 * hyp->boundB);
hyp->boundC = hypgeom_root_norm(P);
hyp->boundD = hypgeom_root_norm(Q);
hyp->boundK = 1 + FLINT_MAX(hyp->boundC, 2 * hyp->boundD);
fmpr_one(hyp->MK);

4
fmprb/hypgeom_sum.c → hypgeom/sum.c
View file

@ -23,7 +23,7 @@
******************************************************************************/
#include "fmprb.h"
#include "hypgeom.h"
static __inline__ void
fmpz_poly_evaluate_si(fmpz_t y, const fmpz_poly_t poly, long x)
@ -201,7 +201,7 @@ fmprb_hypgeom_infsum(fmprb_t P, fmprb_t Q, hypgeom_t hyp, long target_prec, long
hyp->have_precomputed = 1;
}
n = hypgeom_bound(err, hyp->r, hyp->boundA, hyp->boundB,
n = hypgeom_bound(err, hyp->r, hyp->boundC, hyp->boundD,
hyp->boundK, hyp->MK, z, target_prec);
fmprb_hypgeom_sum(P, Q, hyp, n, prec);