mirror of
https://github.com/vale981/arb
synced 2025-03-05 09:21:38 -05:00
use mag methods for hypgeom bounding code; add arb_hypgeom_sum methods
This commit is contained in:
Fredrik Johansson
parent
1fa4420ee1
commit
81009e0a30
9 changed files with 439 additions and 132 deletions
23
arb/gamma.c
23
arb/gamma.c
|
|
@ -49,29 +49,6 @@ acb_get_mag_lower(mag_t z, const acb_t x)
|
|||
arf_clear(t);
|
||||
}
|
||||
|
||||
void
|
||||
mag_mul_ui(mag_t z, const mag_t x, ulong y)
|
||||
{
|
||||
arf_t t;
|
||||
arf_init(t);
|
||||
arf_set_mag(t, x);
|
||||
arf_mul_ui(t, t, y, MAG_BITS, ARF_RND_UP);
|
||||
arf_get_mag(z, t);
|
||||
arf_clear(t);
|
||||
}
|
||||
|
||||
void
|
||||
mag_div_ui(mag_t z, const mag_t x, ulong y)
|
||||
{
|
||||
arf_t t;
|
||||
arf_init(t);
|
||||
arf_set_mag(t, x);
|
||||
arf_div_ui(t, t, y, MAG_BITS, ARF_RND_UP);
|
||||
arf_get_mag(z, t);
|
||||
arf_clear(t);
|
||||
}
|
||||
|
||||
|
||||
/* tuning factor */
|
||||
#define GAMMA_STIRLING_BETA 0.27
|
||||
|
||||
|
|
|
|||
16
hypgeom.h
16
hypgeom.h
|
|
@ -27,6 +27,8 @@
|
|||
#define HYPGEOM_H
|
||||
|
||||
#include "fmprb.h"
|
||||
#include "arb.h"
|
||||
#include "mag.h"
|
||||
#include "fmpz_poly.h"
|
||||
|
||||
#ifdef __cplusplus
|
||||
|
|
@ -46,7 +48,7 @@ typedef struct
|
|||
long boundC;
|
||||
long boundD;
|
||||
long boundK;
|
||||
fmpr_t MK;
|
||||
mag_t MK;
|
||||
}
|
||||
hypgeom_struct;
|
||||
|
||||
|
|
@ -58,15 +60,19 @@ 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_estimate_terms(const mag_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);
|
||||
long hypgeom_bound(mag_t error, int r,
|
||||
long C, long D, long K, const mag_t TK, const mag_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_sum(fmprb_t P, fmprb_t Q, const hypgeom_t hyp, long n, long prec);
|
||||
|
||||
void fmprb_hypgeom_infsum(fmprb_t P, fmprb_t Q, hypgeom_t hyp, long target_prec, long prec);
|
||||
|
||||
void arb_hypgeom_sum(arb_t P, arb_t Q, const hypgeom_t hyp, long n, long prec);
|
||||
|
||||
void arb_hypgeom_infsum(arb_t P, arb_t Q, hypgeom_t hyp, long target_prec, long prec);
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
|
|
|||
232
hypgeom/bound.c
232
hypgeom/bound.c
|
|
@ -103,7 +103,8 @@ fmpr_gamma_ui_ubound(fmpr_t x, ulong n, long prec)
|
|||
}
|
||||
}
|
||||
|
||||
long hypgeom_root_bound(const fmpr_t z, int r, long prec)
|
||||
long
|
||||
hypgeom_root_bound(const mag_t z, int r)
|
||||
{
|
||||
if (r == 0)
|
||||
{
|
||||
|
|
@ -114,8 +115,9 @@ long hypgeom_root_bound(const fmpr_t z, int r, long prec)
|
|||
fmpr_t t;
|
||||
long v;
|
||||
fmpr_init(t);
|
||||
fmpr_root(t, z, r, prec, FMPR_RND_UP);
|
||||
fmpr_add_ui(t, t, 1, prec, FMPR_RND_UP);
|
||||
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);
|
||||
return v;
|
||||
|
|
@ -138,129 +140,204 @@ z^n *
|
|||
(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)
|
||||
hypgeom_term_bound(mag_t Tn, const mag_t TK, long K, long A, long B, int r, const mag_t z, long n)
|
||||
{
|
||||
fmpr_t t, u, num, den;
|
||||
mag_t t, u, num;
|
||||
long m;
|
||||
|
||||
fmpr_init(t);
|
||||
fmpr_init(u);
|
||||
fmpr_init(num);
|
||||
fmpr_init(den);
|
||||
mag_init(t);
|
||||
mag_init(u);
|
||||
mag_init(num);
|
||||
|
||||
m = n - K;
|
||||
|
||||
if (m < 0)
|
||||
{
|
||||
printf("hypgeom term bound\n");
|
||||
abort();
|
||||
}
|
||||
|
||||
/* TK * z^n */
|
||||
fmpr_pow_sloppy_ui(t, z, n, wp, FMPR_RND_UP);
|
||||
fmpr_mul(num, TK, t, wp, FMPR_RND_UP);
|
||||
mag_pow_ui(t, z, n);
|
||||
mag_mul(num, TK, t);
|
||||
|
||||
/* (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);
|
||||
/* numerator: (K+A)! (K-2B)! (K-B+m)! */
|
||||
mag_fac_ui(t, K+A);
|
||||
mag_mul(num, num, t);
|
||||
|
||||
/* (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);
|
||||
mag_fac_ui(t, K-2*B);
|
||||
mag_mul(num, num, t);
|
||||
|
||||
mag_fac_ui(t, K-B+m);
|
||||
mag_mul(num, num, t);
|
||||
|
||||
/* denominator: (K-B)! (K-A+m)! (K-2B+m)! */
|
||||
mag_rfac_ui(t, K-B);
|
||||
mag_mul(num, num, t);
|
||||
|
||||
mag_rfac_ui(t, K-A+m);
|
||||
mag_mul(num, num, t);
|
||||
|
||||
mag_rfac_ui(t, K-2*B+m);
|
||||
mag_mul(num, num, t);
|
||||
|
||||
/* ((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);
|
||||
mag_fac_ui(t, K+m);
|
||||
mag_mul(num, num, t);
|
||||
|
||||
mag_rfac_ui(t, K);
|
||||
mag_mul(num, num, t);
|
||||
}
|
||||
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);
|
||||
mag_fac_ui(t, K);
|
||||
mag_rfac_ui(u, K+m);
|
||||
mag_mul(t, t, u);
|
||||
|
||||
mag_pow_ui(t, t, r-1);
|
||||
mag_mul(num, num, t);
|
||||
}
|
||||
|
||||
fmpr_div(Tn, num, den, wp, FMPR_RND_UP);
|
||||
mag_set(Tn, num);
|
||||
|
||||
fmpr_clear(t);
|
||||
fmpr_clear(u);
|
||||
fmpr_clear(num);
|
||||
fmpr_clear(den);
|
||||
mag_clear(t);
|
||||
mag_clear(u);
|
||||
mag_clear(num);
|
||||
}
|
||||
|
||||
|
||||
/* z = max(x-y, 0), rounding down [lower bound] */
|
||||
void
|
||||
mag_sub_lower(mag_t z, const mag_t x, const mag_t y)
|
||||
{
|
||||
if (mag_is_special(x) || mag_is_special(y))
|
||||
{
|
||||
if (mag_is_zero(y))
|
||||
mag_set(z, x);
|
||||
else if (mag_is_zero(x) || mag_is_inf(y))
|
||||
mag_zero(z);
|
||||
else
|
||||
mag_inf(z);
|
||||
}
|
||||
else
|
||||
{
|
||||
fmpr_t t, u;
|
||||
|
||||
fmpr_init(t);
|
||||
fmpr_init(u);
|
||||
|
||||
mag_get_fmpr(t, x);
|
||||
mag_get_fmpr(u, y);
|
||||
|
||||
fmpr_sub(t, t, u, MAG_BITS, FMPR_RND_DOWN);
|
||||
|
||||
if (fmpr_sgn(t) <= 0)
|
||||
mag_zero(z);
|
||||
else /* XXX: exact? */
|
||||
mag_set_fmpz_2exp_fmpz_lower(z, fmpr_manref(t), fmpr_expref(t));
|
||||
|
||||
fmpr_clear(t);
|
||||
fmpr_clear(u);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
mag_set_ui_2exp_si(mag_t z, ulong v, long e)
|
||||
{
|
||||
mag_set_ui(z, v);
|
||||
mag_mul_2exp_si(z, z, e);
|
||||
}
|
||||
|
||||
double
|
||||
mag_get_d(const mag_t z)
|
||||
{
|
||||
long exp;
|
||||
|
||||
if (mag_is_zero(z))
|
||||
return 0.0;
|
||||
else if (mag_is_inf(z))
|
||||
return D_INF;
|
||||
|
||||
/* XXX: overflow/underflow properly */
|
||||
exp = MAG_EXP(z);
|
||||
|
||||
if (exp < -1000)
|
||||
return ldexp(1.0, -1000);
|
||||
else if (exp > 1000)
|
||||
return D_INF;
|
||||
else
|
||||
return ldexp(MAG_MAN(z), exp - MAG_BITS);
|
||||
}
|
||||
|
||||
long
|
||||
hypgeom_bound(fmpr_t error, int r,
|
||||
long A, long B, long K, const fmpr_t TK, const fmpr_t z, long prec)
|
||||
hypgeom_bound(mag_t error, int r,
|
||||
long A, long B, long K, const mag_t TK, const mag_t z, long tol_2exp)
|
||||
{
|
||||
fmpr_t Tn, t, u, one, tol, num, den;
|
||||
long wp = FMPRB_RAD_PREC;
|
||||
mag_t Tn, t, u, one, tol, num, den;
|
||||
long n, m;
|
||||
|
||||
fmpr_init(Tn);
|
||||
fmpr_init(t);
|
||||
fmpr_init(u);
|
||||
fmpr_init(one);
|
||||
fmpr_init(tol);
|
||||
fmpr_init(num);
|
||||
fmpr_init(den);
|
||||
mag_init(Tn);
|
||||
mag_init(t);
|
||||
mag_init(u);
|
||||
mag_init(one);
|
||||
mag_init(tol);
|
||||
mag_init(num);
|
||||
mag_init(den);
|
||||
|
||||
fmpr_one(one);
|
||||
fmpr_set_ui_2exp_si(tol, 1UL, -prec);
|
||||
mag_one(one);
|
||||
mag_set_ui_2exp_si(tol, 1UL, -tol_2exp);
|
||||
|
||||
/* approximate number of needed terms */
|
||||
n = hypgeom_estimate_terms(z, r, prec);
|
||||
n = hypgeom_estimate_terms(z, r, tol_2exp);
|
||||
|
||||
/* 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 */
|
||||
m = hypgeom_root_bound(z, r, wp);
|
||||
m = hypgeom_root_bound(z, r);
|
||||
n = FLINT_MAX(n, m);
|
||||
|
||||
/* 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. */
|
||||
as soon as G(k) < 1. */
|
||||
|
||||
/* bound T(n-1) */
|
||||
hypgeom_term_bound(Tn, TK, K, A, B, r, z, n-1, wp);
|
||||
hypgeom_term_bound(Tn, TK, K, A, B, r, z, n-1);
|
||||
|
||||
while (1)
|
||||
{
|
||||
/* 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);
|
||||
mag_mul_ui(num, z, n);
|
||||
mag_mul_ui(num, num, n - B);
|
||||
|
||||
mag_set_ui_lower(den, n - A);
|
||||
mag_mul_ui_lower(den, den, n - 2*B);
|
||||
|
||||
if (r != 0)
|
||||
{
|
||||
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);
|
||||
mag_set_ui_lower(u, n);
|
||||
mag_pow_ui_lower(u, u, r);
|
||||
mag_mul_lower(den, den, u);
|
||||
}
|
||||
|
||||
fmpr_div(t, num, den, wp, FMPR_RND_UP);
|
||||
mag_div(t, num, den);
|
||||
|
||||
/* multiply bound for T(n-1) by bound for R(n) to bound T(n) */
|
||||
fmpr_mul(Tn, Tn, t, wp, FMPR_RND_UP);
|
||||
mag_mul(Tn, Tn, t);
|
||||
|
||||
/* 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);
|
||||
/* u = max(1-t, 0), rounding down [lower bound] */
|
||||
mag_sub_lower(u, one, t);
|
||||
|
||||
if (fmpr_cmp(u, tol) < 0)
|
||||
if (!mag_is_zero(u))
|
||||
{
|
||||
mag_div(u, Tn, u);
|
||||
|
||||
if (mag_cmp(u, tol) < 0)
|
||||
{
|
||||
fmpr_set(error, u);
|
||||
mag_set(error, u);
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
|
@ -269,13 +346,14 @@ hypgeom_bound(fmpr_t error, int r,
|
|||
n++;
|
||||
}
|
||||
|
||||
fmpr_clear(Tn);
|
||||
fmpr_clear(t);
|
||||
fmpr_clear(u);
|
||||
fmpr_clear(one);
|
||||
fmpr_clear(tol);
|
||||
fmpr_clear(num);
|
||||
fmpr_clear(den);
|
||||
mag_clear(Tn);
|
||||
mag_clear(t);
|
||||
mag_clear(u);
|
||||
mag_clear(one);
|
||||
mag_clear(tol);
|
||||
mag_clear(num);
|
||||
mag_clear(den);
|
||||
|
||||
return n;
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -39,13 +39,14 @@ static __inline__ double d_root(double x, int r)
|
|||
return pow(x, 1. / r);
|
||||
}
|
||||
|
||||
double mag_get_d(const mag_t z);
|
||||
|
||||
long
|
||||
hypgeom_estimate_terms(const fmpr_t z, int r, long prec)
|
||||
hypgeom_estimate_terms(const mag_t z, int r, long prec)
|
||||
{
|
||||
double y, t;
|
||||
|
||||
t = fmpr_get_d(z, FMPR_RND_UP);
|
||||
t = fabs(t);
|
||||
t = mag_get_d(z);
|
||||
|
||||
if (t == 0)
|
||||
return 1;
|
||||
|
|
|
|||
|
|
@ -32,7 +32,7 @@ hypgeom_init(hypgeom_t hyp)
|
|||
fmpz_poly_init(hyp->B);
|
||||
fmpz_poly_init(hyp->P);
|
||||
fmpz_poly_init(hyp->Q);
|
||||
fmpr_init(hyp->MK);
|
||||
mag_init(hyp->MK);
|
||||
hyp->have_precomputed = 0;
|
||||
}
|
||||
|
||||
|
|
@ -43,5 +43,5 @@ hypgeom_clear(hypgeom_t hyp)
|
|||
fmpz_poly_clear(hyp->B);
|
||||
fmpz_poly_clear(hyp->P);
|
||||
fmpz_poly_clear(hyp->Q);
|
||||
fmpr_clear(hyp->MK);
|
||||
mag_clear(hyp->MK);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -126,17 +126,15 @@ _hypgeom_precompute(hypgeom_t hyp, const fmpz_poly_t P, const fmpz_poly_t Q)
|
|||
hyp->boundD = hypgeom_root_norm(Q);
|
||||
hyp->boundK = 1 + FLINT_MAX(hyp->boundC, 2 * hyp->boundD);
|
||||
|
||||
fmpr_one(hyp->MK);
|
||||
mag_one(hyp->MK);
|
||||
|
||||
for (k = 1; k <= hyp->boundK; k++)
|
||||
{
|
||||
fmpz_poly_evaluate_si(t, P, k);
|
||||
fmpz_abs(t, t);
|
||||
fmpr_mul_fmpz(hyp->MK, hyp->MK, t, FMPRB_RAD_PREC, FMPR_RND_UP);
|
||||
mag_mul_fmpz(hyp->MK, hyp->MK, t);
|
||||
|
||||
fmpz_poly_evaluate_si(t, Q, k);
|
||||
fmpz_abs(t, t);
|
||||
fmpr_div_fmpz(hyp->MK, hyp->MK, t, FMPRB_RAD_PREC, FMPR_RND_UP);
|
||||
mag_div_fmpz(hyp->MK, hyp->MK, t);
|
||||
}
|
||||
|
||||
fmpz_clear(t);
|
||||
|
|
@ -163,12 +161,10 @@ hypgeom_precompute(hypgeom_t hyp)
|
|||
fmpz_init(t);
|
||||
|
||||
fmpz_poly_evaluate_si(t, hyp->A, 0);
|
||||
fmpz_abs(t, t);
|
||||
fmpr_mul_fmpz(hyp->MK, hyp->MK, t, FMPRB_RAD_PREC, FMPR_RND_UP);
|
||||
mag_mul_fmpz(hyp->MK, hyp->MK, t);
|
||||
|
||||
fmpz_poly_evaluate_si(t, hyp->B, 0);
|
||||
fmpz_abs(t, t);
|
||||
fmpr_div_fmpz(hyp->MK, hyp->MK, t, FMPRB_RAD_PREC, FMPR_RND_UP);
|
||||
mag_div_fmpz(hyp->MK, hyp->MK, t);
|
||||
|
||||
fmpz_clear(t);
|
||||
}
|
||||
|
|
|
|||
157
hypgeom/sum.c
157
hypgeom/sum.c
|
|
@ -159,7 +159,7 @@ bsplit_recursive_fmprb(fmprb_t P, fmprb_t Q, fmprb_t B, fmprb_t T,
|
|||
}
|
||||
|
||||
void
|
||||
fmprb_hypgeom_sum(fmprb_t P, fmprb_t Q, const hypgeom_t hyp, const long n, long prec)
|
||||
fmprb_hypgeom_sum(fmprb_t P, fmprb_t Q, const hypgeom_t hyp, long n, long prec)
|
||||
{
|
||||
if (n < 1)
|
||||
{
|
||||
|
|
@ -183,17 +183,14 @@ fmprb_hypgeom_sum(fmprb_t P, fmprb_t Q, const hypgeom_t hyp, const long n, long
|
|||
void
|
||||
fmprb_hypgeom_infsum(fmprb_t P, fmprb_t Q, hypgeom_t hyp, long target_prec, long prec)
|
||||
{
|
||||
fmpr_t err;
|
||||
fmpr_t z;
|
||||
mag_t err, z;
|
||||
long n;
|
||||
|
||||
fmpr_init(err);
|
||||
fmpr_init(z);
|
||||
mag_init(err);
|
||||
mag_init(z);
|
||||
|
||||
fmpr_fmpz_div_fmpz(z,
|
||||
hyp->P->coeffs + hyp->P->length - 1,
|
||||
hyp->Q->coeffs + hyp->Q->length - 1, FMPRB_RAD_PREC, FMPR_RND_UP);
|
||||
fmpr_abs(z, 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)
|
||||
{
|
||||
|
|
@ -214,14 +211,148 @@ fmprb_hypgeom_infsum(fmprb_t P, fmprb_t Q, hypgeom_t hyp, long target_prec, long
|
|||
|
||||
/* We have p/q = s + err i.e. (p + q*err)/q = s */
|
||||
{
|
||||
fmpr_t u;
|
||||
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, err, 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);
|
||||
}
|
||||
|
||||
fmpr_clear(z);
|
||||
fmpr_clear(err);
|
||||
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, long a, long b, int cont, long 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);
|
||||
|
||||
arb_set_fmpz(P, PP);
|
||||
arb_set_fmpz(Q, QQ);
|
||||
arb_set_fmpz(B, BB);
|
||||
arb_set_fmpz(T, TT);
|
||||
|
||||
fmpz_clear(PP);
|
||||
fmpz_clear(QQ);
|
||||
fmpz_clear(BB);
|
||||
fmpz_clear(TT);
|
||||
}
|
||||
else
|
||||
{
|
||||
long m;
|
||||
arb_t P2, Q2, B2, T2;
|
||||
|
||||
m = (a + b) / 2;
|
||||
|
||||
arb_init(P2);
|
||||
arb_init(Q2);
|
||||
arb_init(B2);
|
||||
arb_init(T2);
|
||||
|
||||
bsplit_recursive_arb(P, Q, B, T, hyp, a, m, 1, prec);
|
||||
bsplit_recursive_arb(P2, Q2, B2, T2, hyp, m, b, 1, prec);
|
||||
|
||||
if (arb_is_one(B) && arb_is_one(B2))
|
||||
{
|
||||
arb_mul(T, T, Q2, prec);
|
||||
arb_addmul(T, P, T2, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
arb_mul(T, T, B2, prec);
|
||||
arb_mul(T, T, Q2, prec);
|
||||
arb_mul(T2, T2, B, prec);
|
||||
arb_addmul(T, P, T2, prec);
|
||||
}
|
||||
|
||||
arb_mul(B, B, B2, prec);
|
||||
arb_mul(Q, Q, Q2, prec);
|
||||
if (cont)
|
||||
arb_mul(P, P, P2, prec);
|
||||
|
||||
arb_clear(P2);
|
||||
arb_clear(Q2);
|
||||
arb_clear(B2);
|
||||
arb_clear(T2);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
arb_hypgeom_sum(arb_t P, arb_t Q, const hypgeom_t hyp, long n, long prec)
|
||||
{
|
||||
if (n < 1)
|
||||
{
|
||||
arb_zero(P);
|
||||
arb_one(Q);
|
||||
}
|
||||
else
|
||||
{
|
||||
arb_t B, T;
|
||||
arb_init(B);
|
||||
arb_init(T);
|
||||
bsplit_recursive_arb(P, Q, B, T, hyp, 0, n, 0, prec);
|
||||
if (!arb_is_one(B))
|
||||
arb_mul(Q, Q, B, prec);
|
||||
arb_swap(P, T);
|
||||
arb_clear(B);
|
||||
arb_clear(T);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
arb_hypgeom_infsum(arb_t P, arb_t Q, hypgeom_t hyp, long target_prec, long prec)
|
||||
{
|
||||
mag_t err, z;
|
||||
long 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);
|
||||
|
||||
arb_hypgeom_sum(P, Q, hyp, n, prec);
|
||||
|
||||
if (arf_sgn(arb_midref(Q)) < 0)
|
||||
{
|
||||
arb_neg(P, P);
|
||||
arb_neg(Q, Q);
|
||||
}
|
||||
|
||||
/* We have p/q = s + err i.e. (p + q*err)/q = s */
|
||||
{
|
||||
mag_t u;
|
||||
mag_init(u);
|
||||
arb_get_mag(u, Q);
|
||||
mag_mul(u, u, err);
|
||||
mag_add(arb_radref(P), arb_radref(P), u);
|
||||
mag_clear(u);
|
||||
}
|
||||
|
||||
mag_clear(z);
|
||||
mag_clear(err);
|
||||
}
|
||||
|
||||
|
|
|
|||
99
mag.h
99
mag.h
|
|
@ -650,6 +650,105 @@ void mag_rfac_ui(mag_t z, ulong n);
|
|||
/* TODO: document/test */
|
||||
void mag_bernoulli_div_fac_ui(mag_t z, ulong n);
|
||||
|
||||
/* TODO: document/test */
|
||||
void mag_set_fmpz_2exp_fmpz_lower(mag_t z, const fmpz_t man, const fmpz_t exp);
|
||||
|
||||
static __inline__ void
|
||||
mag_set_ui(mag_t z, ulong x)
|
||||
{
|
||||
fmpz_t man, exp;
|
||||
fmpz_init_set_ui(man, x);
|
||||
*exp = 0;
|
||||
mag_set_fmpz_2exp_fmpz(z, man, exp);
|
||||
fmpz_clear(man);
|
||||
}
|
||||
|
||||
static __inline__ void
|
||||
mag_set_ui_lower(mag_t z, ulong x)
|
||||
{
|
||||
fmpz_t man, exp;
|
||||
fmpz_init_set_ui(man, x);
|
||||
*exp = 0;
|
||||
mag_set_fmpz_2exp_fmpz_lower(z, man, exp);
|
||||
fmpz_clear(man);
|
||||
}
|
||||
|
||||
static __inline__ void
|
||||
mag_set_fmpz(mag_t z, const fmpz_t x)
|
||||
{
|
||||
fmpz_t exp;
|
||||
*exp = 0;
|
||||
mag_set_fmpz_2exp_fmpz(z, x, exp);
|
||||
}
|
||||
|
||||
static __inline__ void
|
||||
mag_set_fmpz_lower(mag_t z, const fmpz_t x)
|
||||
{
|
||||
fmpz_t exp;
|
||||
*exp = 0;
|
||||
mag_set_fmpz_2exp_fmpz_lower(z, x, exp);
|
||||
}
|
||||
|
||||
static __inline__ void
|
||||
mag_mul_ui(mag_t z, const mag_t x, ulong y)
|
||||
{
|
||||
mag_t t;
|
||||
mag_init(t);
|
||||
mag_set_ui(t, y);
|
||||
mag_mul(z, x, t);
|
||||
mag_clear(t);
|
||||
}
|
||||
|
||||
static __inline__ void
|
||||
mag_mul_ui_lower(mag_t z, const mag_t x, ulong y)
|
||||
{
|
||||
mag_t t;
|
||||
mag_init(t);
|
||||
mag_set_ui_lower(t, y);
|
||||
mag_mul_lower(z, x, t);
|
||||
mag_clear(t);
|
||||
}
|
||||
|
||||
static __inline__ void
|
||||
mag_mul_fmpz(mag_t z, const mag_t x, const fmpz_t y)
|
||||
{
|
||||
mag_t t;
|
||||
mag_init(t);
|
||||
mag_set_fmpz(t, y);
|
||||
mag_mul(z, x, t);
|
||||
mag_clear(t);
|
||||
}
|
||||
|
||||
static __inline__ void
|
||||
mag_mul_fmpz_lower(mag_t z, const mag_t x, const fmpz_t y)
|
||||
{
|
||||
mag_t t;
|
||||
mag_init(t);
|
||||
mag_set_fmpz_lower(t, y);
|
||||
mag_mul_lower(z, x, t);
|
||||
mag_clear(t);
|
||||
}
|
||||
|
||||
static __inline__ void
|
||||
mag_div_ui(mag_t z, const mag_t x, ulong y)
|
||||
{
|
||||
mag_t t;
|
||||
mag_init(t);
|
||||
mag_set_ui_lower(t, y);
|
||||
mag_div(z, x, t);
|
||||
mag_clear(t);
|
||||
}
|
||||
|
||||
static __inline__ void
|
||||
mag_div_fmpz(mag_t z, const mag_t x, const fmpz_t y)
|
||||
{
|
||||
mag_t t;
|
||||
mag_init(t);
|
||||
mag_set_fmpz_lower(t, y);
|
||||
mag_div(z, x, t);
|
||||
mag_clear(t);
|
||||
}
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
|
|
|||
|
|
@ -42,3 +42,22 @@ mag_set_fmpz_2exp_fmpz(mag_t z, const fmpz_t man, const fmpz_t exp)
|
|||
_fmpz_add_fast(MAG_EXPREF(z), exp, cexp + MAG_BITS);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
mag_set_fmpz_2exp_fmpz_lower(mag_t z, const fmpz_t man, const fmpz_t exp)
|
||||
{
|
||||
if (fmpz_is_zero(man))
|
||||
{
|
||||
mag_zero(z);
|
||||
}
|
||||
else
|
||||
{
|
||||
mp_limb_t m;
|
||||
long cexp;
|
||||
|
||||
m = fmpz_abs_lbound_ui_2exp(&cexp, man, MAG_BITS);
|
||||
MAG_MAN(z) = m;
|
||||
_fmpz_add_fast(MAG_EXPREF(z), exp, cexp + MAG_BITS);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue