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

exorcise fmpr and fmprb from the hypgeom module

Browse source
This commit is contained in:
Fredrik Johansson 2016-02-26 17:56:09 +01:00
parent 493234ef86
commit 89ed0aab05
4 changed files with 7 additions and 161 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

16
doc/source/hypgeom.rst
View file

@ -179,22 +179,6 @@ Error bounding
Summation
-------------------------------------------------------------------------------
.. function:: void fmprb_hypgeom_sum(fmprb_t P, fmprb_t Q, const hypgeom_t hyp, const slong n, slong 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 fmprb_hypgeom_infsum(fmprb_t P, fmprb_t Q, hypgeom_t hyp, slong tol, slong 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*.
.. function:: void arb_hypgeom_sum(arb_t P, arb_t Q, const hypgeom_t hyp, const slong n, slong prec)
Computes `P, Q` such that `P / Q = \sum_{k=0}^{n-1} T(k)` where `T(k)`

5
hypgeom.h
View file

@ -26,7 +26,6 @@
#ifndef HYPGEOM_H
#define HYPGEOM_H
#include "fmprb.h"
#include "arb.h"
#include "mag.h"
#include "fmpz_poly.h"
@ -65,10 +64,6 @@ slong hypgeom_estimate_terms(const mag_t z, int r, slong prec);
slong hypgeom_bound(mag_t error, int r,
slong C, slong D, slong K, const mag_t TK, const mag_t z, slong prec);
void fmprb_hypgeom_sum(fmprb_t P, fmprb_t Q, const hypgeom_t hyp, slong n, slong prec);
void fmprb_hypgeom_infsum(fmprb_t P, fmprb_t Q, hypgeom_t hyp, slong target_prec, slong prec);
void arb_hypgeom_sum(arb_t P, arb_t Q, const hypgeom_t hyp, slong n, slong prec);
void arb_hypgeom_infsum(arb_t P, arb_t Q, hypgeom_t hyp, slong target_prec, slong prec);

14
hypgeom/bound.c
View file

@ -36,14 +36,14 @@ hypgeom_root_bound(const mag_t z, int r)
}
else
{
fmpr_t t;
arf_t t;
slong v;
fmpr_init(t);
mag_get_fmpr(t, z);
fmpr_root(t, t, r, MAG_BITS, FMPR_RND_UP);
fmpr_add_ui(t, t, 1, MAG_BITS, FMPR_RND_UP);
v = fmpr_get_si(t, FMPR_RND_UP);
fmpr_clear(t);
arf_init(t);
arf_set_mag(t, z);
arf_root(t, t, r, MAG_BITS, ARF_RND_UP);
arf_add_ui(t, t, 1, MAG_BITS, ARF_RND_UP);
v = arf_get_si(t, ARF_RND_UP);
arf_clear(t);
return v;
}
}

133
hypgeom/sum.c
View file

@ -93,139 +93,6 @@ bsplit_recursive_fmpz(fmpz_t P, fmpz_t Q, fmpz_t B, fmpz_t T,
}
}
static void
bsplit_recursive_fmprb(fmprb_t P, fmprb_t Q, fmprb_t B, fmprb_t T,
const hypgeom_t hyp, slong a, slong b, int cont, slong prec)
{
if (b - a < 4)
{
fmpz_t PP, QQ, BB, TT;
fmpz_init(PP);
fmpz_init(QQ);
fmpz_init(BB);
fmpz_init(TT);
bsplit_recursive_fmpz(PP, QQ, BB, TT, hyp, a, b, cont);
fmprb_set_fmpz(P, PP);
fmprb_set_fmpz(Q, QQ);
fmprb_set_fmpz(B, BB);
fmprb_set_fmpz(T, TT);
fmpz_clear(PP);
fmpz_clear(QQ);
fmpz_clear(BB);
fmpz_clear(TT);
}
else
{
slong m;
fmprb_t P2, Q2, B2, T2;
m = (a + b) / 2;
fmprb_init(P2);
fmprb_init(Q2);
fmprb_init(B2);
fmprb_init(T2);
bsplit_recursive_fmprb(P, Q, B, T, hyp, a, m, 1, prec);
bsplit_recursive_fmprb(P2, Q2, B2, T2, hyp, m, b, 1, prec);
if (fmprb_is_one(B) && fmprb_is_one(B2))
{
fmprb_mul(T, T, Q2, prec);
fmprb_addmul(T, P, T2, prec);
}
else
{
fmprb_mul(T, T, B2, prec);
fmprb_mul(T, T, Q2, prec);
fmprb_mul(T2, T2, B, prec);
fmprb_addmul(T, P, T2, prec);
}
fmprb_mul(B, B, B2, prec);
fmprb_mul(Q, Q, Q2, prec);
if (cont)
fmprb_mul(P, P, P2, prec);
fmprb_clear(P2);
fmprb_clear(Q2);
fmprb_clear(B2);
fmprb_clear(T2);
}
}
void
fmprb_hypgeom_sum(fmprb_t P, fmprb_t Q, const hypgeom_t hyp, slong n, slong prec)
{
if (n < 1)
{
fmprb_zero(P);
fmprb_one(Q);
}
else
{
fmprb_t B, T;
fmprb_init(B);
fmprb_init(T);
bsplit_recursive_fmprb(P, Q, B, T, hyp, 0, n, 0, prec);
if (!fmprb_is_one(B))
fmprb_mul(Q, Q, B, prec);
fmprb_swap(P, T);
fmprb_clear(B);
fmprb_clear(T);
}
}
void
fmprb_hypgeom_infsum(fmprb_t P, fmprb_t Q, hypgeom_t hyp, slong target_prec, slong prec)
{
mag_t err, z;
slong n;
mag_init(err);
mag_init(z);
mag_set_fmpz(z, hyp->P->coeffs + hyp->P->length - 1);
mag_div_fmpz(z, z, hyp->Q->coeffs + hyp->Q->length - 1);
if (!hyp->have_precomputed)
{
hypgeom_precompute(hyp);
hyp->have_precomputed = 1;
}
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);
if (fmpr_sgn(fmprb_midref(Q)) < 0)
{
fmprb_neg(P, P);
fmprb_neg(Q, Q);
}
/* We have p/q = s + err i.e. (p + q*err)/q = s */
{
fmpr_t u, v;
fmpr_init(u);
fmpr_init(v);
mag_get_fmpr(v, err);
fmpr_add(u, fmprb_midref(Q), fmprb_radref(Q), FMPRB_RAD_PREC, FMPR_RND_UP);
fmpr_mul(u, u, v, FMPRB_RAD_PREC, FMPR_RND_UP);
fmprb_add_error_fmpr(P, u);
fmpr_clear(u);
fmpr_clear(v);
}
mag_clear(z);
mag_clear(err);
}
static void
bsplit_recursive_arb(arb_t P, arb_t Q, arb_t B, arb_t T,
const hypgeom_t hyp, slong a, slong b, int cont, slong prec)