diff options
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/lngamma.c')
-rw-r--r-- | Build/source/libs/mpfr/mpfr-src/src/lngamma.c | 233 |
1 files changed, 179 insertions, 54 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/lngamma.c b/Build/source/libs/mpfr/mpfr-src/src/lngamma.c index 84bd7d1fae7..5511fd1dccd 100644 --- a/Build/source/libs/mpfr/mpfr-src/src/lngamma.c +++ b/Build/source/libs/mpfr/mpfr-src/src/lngamma.c @@ -33,6 +33,10 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., static void mpfr_gamma_alpha (mpfr_t s, mpfr_prec_t p) { + MPFR_LOG_FUNC + (("p=%Pu", p), + ("s[%Pu]=%.*Rg", mpfr_get_prec (s), mpfr_log_prec, s)); + if (p <= 100) mpfr_set_ui_2exp (s, 614, -10, MPFR_RNDN); /* about 0.6 */ else if (p <= 500) @@ -72,7 +76,7 @@ mpfr_explgamma (mpfr_ptr y, mpfr_srcptr x, mpfr_save_expo_t *pexpo, /* s1 = RNDD(lngamma(x)), inexact */ if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags1))) { - if (MPFR_SIGN (s1) > 0) + if (MPFR_IS_POS (s1)) { MPFR_SAVE_EXPO_UPDATE_FLAGS (*pexpo, MPFR_FLAGS_OVERFLOW); return mpfr_overflow (y, rnd, sign); @@ -141,6 +145,12 @@ unit_bit (mpfr_srcptr x) #endif +/* FIXME: There is an internal overflow when z is very large. + Simple overflow detection with possible false negatives? + For the particular cases near the overflow boundary, + scaling by a power of two? +*/ + /* lngamma(x) = log(gamma(x)). We use formula [6.1.40] from Abramowitz&Stegun: lngamma(z) = (z-1/2)*log(z) - z + 1/2*log(2*Pi) @@ -161,17 +171,18 @@ GAMMA_FUNC (mpfr_ptr y, mpfr_srcptr z0, mpfr_rnd_t rnd) { mpfr_prec_t precy, w; /* working precision */ mpfr_t s, t, u, v, z; - unsigned long m, k, maxm; - mpz_t *INITIALIZED(B); /* variable B declared as initialized */ - int compared; - int inexact = 0; /* 0 means: result y not set yet */ + unsigned long m, k, maxm, l; + int compared, inexact; mpfr_exp_t err_s, err_t; - unsigned long Bm = 0; /* number of allocated B[] */ - unsigned long oldBm; double d; MPFR_SAVE_EXPO_DECL (expo); MPFR_ZIV_DECL (loop); + MPFR_LOG_FUNC + (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (z0), mpfr_log_prec, z0, rnd), + ("y[%Pu]=%.*Rg inexact=%d", + mpfr_get_prec (y), mpfr_log_prec, y, inexact)); + compared = mpfr_cmp_ui (z0, 1); MPFR_SAVE_EXPO_MARK (expo); @@ -183,9 +194,52 @@ GAMMA_FUNC (mpfr_ptr y, mpfr_srcptr z0, mpfr_rnd_t rnd) return mpfr_set_ui (y, 0, MPFR_RNDN); /* lngamma(1 or 2) = +0 */ } + /* Deal with very large inputs: according to [6.1.42], if we denote + R_n(z) = lngamma(z) - (z-1/2)*log(z) + z - 1/2*log(2*Pi), we have + |R_n(z)| <= B_2/2/z, thus for z >= 2 we have + |lngamma(z) - [z*(log(z) - 1)]| < 1/2*log(z) + 1. */ + if (compared > 0 && MPFR_GET_EXP (z0) >= (mpfr_exp_t) MPFR_PREC(y) + 2) + { + /* Since PREC(y) >= 2, this ensures EXP(z0) >= 4, thus |z0| >= 8, + thus 1/2*log(z0) + 1 < log(z0). + Since the largest possible z is < 2^(2^62) on a 64-bit machine, + the largest value of log(z) is 2^62*log(2.) < 3.2e18 < 2^62, + thus if we use at least 62 bits of precision, then log(t)-1 will + be exact */ + mpfr_init2 (t, MPFR_PREC(y) >= 52 ? MPFR_PREC(y) + 10 : 62); + mpfr_log (t, z0, MPFR_RNDU); /* error < 1 ulp */ + inexact = mpfr_sub_ui (t, t, 1, MPFR_RNDU); /* err < 2 ulps, since the + exponent of t might have + decreased */ + MPFR_ASSERTD(inexact == 0); + mpfr_mul (t, z0, t, MPFR_RNDU); /* err < 1+2*2=5 ulps according to + "Generic error on multiplication" + in algorithms.tex */ + if (MPFR_IS_INF(t)) + { + mpfr_clear (t); + MPFR_SAVE_EXPO_FREE (expo); + inexact = mpfr_overflow (y, rnd, 1); + return inexact; + } + if (MPFR_GET_EXP(t) - MPFR_PREC(t) >= 62) + { + /* then ulp(t) >= 2^62 > log(z0) thus the total error is bounded + by 6 ulp(t) */ + if (MPFR_CAN_ROUND (t, MPFR_PREC(t) - 3, MPFR_PREC(y), rnd)) + { + inexact = mpfr_set (y, t, rnd); + mpfr_clear (t); + MPFR_SAVE_EXPO_FREE (expo); + return mpfr_check_range (y, inexact, rnd); + } + } + mpfr_clear (t); + } + /* Deal here with tiny inputs. We have for -0.3 <= x <= 0.3: - log|x| - gamma*x <= log|gamma(x)| <= - log|x| - gamma*x + x^2 */ - if (MPFR_EXP(z0) <= - (mpfr_exp_t) MPFR_PREC(y)) + if (MPFR_GET_EXP (z0) <= - (mpfr_exp_t) MPFR_PREC(y)) { mpfr_t l, h, g; int ok, inex1, inex2; @@ -250,7 +304,7 @@ GAMMA_FUNC (mpfr_ptr y, mpfr_srcptr z0, mpfr_rnd_t rnd) which would need precision n. */ MPFR_ZIV_NEXT (loop, prec); } - while (prec <= -MPFR_EXP(z0)); + while (prec <= - MPFR_GET_EXP (z0)); MPFR_ZIV_FREE (loop); } #endif @@ -263,6 +317,8 @@ GAMMA_FUNC (mpfr_ptr y, mpfr_srcptr z0, mpfr_rnd_t rnd) mpfr_init2 (v, MPFR_PREC_MIN); mpfr_init2 (z, MPFR_PREC_MIN); + inexact = 0; /* 0 means: result y not set yet */ + if (compared < 0) { mpfr_exp_t err_u; @@ -294,7 +350,10 @@ GAMMA_FUNC (mpfr_ptr y, mpfr_srcptr z0, mpfr_rnd_t rnd) ulp(u)/2 + (2-z0)*max(1,log(2-z0))*2^(1-w) = (1/2 + (2-z0)*max(1,log(2-z0))*2^(1-E(u))) ulp(u) */ d = (double) MPFR_GET_EXP(s) * 0.694; /* upper bound for log(2-z0) */ - err_u = MPFR_GET_EXP(s) + __gmpfr_ceil_log2 (d) + 1 - MPFR_GET_EXP(u); + if (MPFR_IS_ZERO(u)) /* in that case the error on u is zero */ + err_u = 0; + else + err_u = MPFR_GET_EXP(s) + __gmpfr_ceil_log2 (d) + 1 - MPFR_GET_EXP(u); err_u = (err_u >= 0) ? err_u + 1 : 0; /* now the error on u is bounded by 2^err_u ulps */ @@ -341,19 +400,22 @@ GAMMA_FUNC (mpfr_ptr y, mpfr_srcptr z0, mpfr_rnd_t rnd) } else { - err_s += 1 - MPFR_GET_EXP(v); + /* if v = 0 here, it was 1 before the call to mpfr_log, + thus the error on v was zero */ + if (!MPFR_IS_ZERO(v)) + err_s += 1 - MPFR_GET_EXP(v); err_s = (err_s >= 0) ? err_s + 1 : 0; /* the error on v is bounded by 2^err_s ulps */ err_u += MPFR_GET_EXP(u); /* absolute error on u */ - err_s += MPFR_GET_EXP(v); /* absolute error on v */ + if (!MPFR_IS_ZERO(v)) /* same as above */ + err_s += MPFR_GET_EXP(v); /* absolute error on v */ mpfr_sub (s, v, u, MPFR_RNDN); /* the total error on s is bounded by ulp(s)/2 + 2^(err_u-w) + 2^(err_s-w) <= ulp(s)/2 + 2^(max(err_u,err_s)+1-w) */ err_s = (err_s >= err_u) ? err_s : err_u; err_s += 1 - MPFR_GET_EXP(s); /* error is 2^err_s ulp(s) */ err_s = (err_s >= 0) ? err_s + 1 : 0; - if (mpfr_can_round (s, w - err_s, MPFR_RNDN, MPFR_RNDZ, precy - + (rnd == MPFR_RNDN))) + if (MPFR_CAN_ROUND (s, w - err_s, precy, rnd)) goto end; } MPFR_ZIV_NEXT (loop, w); @@ -379,11 +441,14 @@ GAMMA_FUNC (mpfr_ptr y, mpfr_srcptr z0, mpfr_rnd_t rnd) and we need k steps of argument reconstruction. Assuming k is large with respect to z0, and k = n, we get 1/(Pi*e)^(2n) ~ 2^(-w), i.e., k ~ w*log(2)/2/log(Pi*e) ~ 0.1616 * w. - However, since the series is more expensive to compute, the optimal - value seems to be k ~ 4.5 * w experimentally. */ + However, since the series is slightly more expensive to compute, + the optimal value seems to be k ~ 0.25 * w experimentally (with + caching of Bernoulli numbers). + For only one computation of gamma with large precision, it is better + to set k to a larger value, say k ~ w. */ mpfr_set_prec (s, 53); mpfr_gamma_alpha (s, w); - mpfr_set_ui_2exp (s, 9, -1, MPFR_RNDU); + mpfr_set_ui_2exp (s, 4, -4, MPFR_RNDU); mpfr_mul_ui (s, s, w, MPFR_RNDU); if (mpfr_cmp (z0, s) < 0) { @@ -431,15 +496,8 @@ GAMMA_FUNC (mpfr_ptr y, mpfr_srcptr z0, mpfr_rnd_t rnd) mpfr_mul (u, u, u, MPFR_RNDN); /* 1/z^2 * (1+u)^3 */ - if (Bm == 0) - { - B = mpfr_bernoulli_internal ((mpz_t *) 0, 0); - B = mpfr_bernoulli_internal (B, 1); - Bm = 2; - } - /* m <= maxm ensures that 2*m*(2*m+1) <= ULONG_MAX */ - maxm = 1UL << (GMP_NUMB_BITS / 2 - 1); + maxm = 1UL << (sizeof(unsigned long) * CHAR_BIT / 2 - 1); /* s:(1+u)^15, t:(1+u)^2, t <= 3/128 */ @@ -463,12 +521,7 @@ GAMMA_FUNC (mpfr_ptr y, mpfr_srcptr z0, mpfr_rnd_t rnd) } /* (1+u)^(10m-8) */ /* invariant: t=1/(2m)/(2m-1)/z^(2m-1)/(2m+1)! */ - if (Bm <= m) - { - B = mpfr_bernoulli_internal (B, m); /* B[2m]*(2m+1)!, exact */ - Bm ++; - } - mpfr_mul_z (v, t, B[m], MPFR_RNDN); /* (1+u)^(10m-7) */ + mpfr_mul_z (v, t, mpfr_bernoulli_cache(m), MPFR_RNDN); /* (1+u)^(10m-7) */ MPFR_ASSERTD(MPFR_GET_EXP(v) <= - (2 * m + 3)); mpfr_add (s, s, v, MPFR_RNDN); } @@ -492,22 +545,92 @@ GAMMA_FUNC (mpfr_ptr y, mpfr_srcptr z0, mpfr_rnd_t rnd) /* add 1/2*log(2*Pi) and subtract log(z0*(z0+1)*...*(z0+k-1)) */ mpfr_const_pi (v, MPFR_RNDN); /* v = Pi*(1+u) */ mpfr_mul_2ui (v, v, 1, MPFR_RNDN); /* v = 2*Pi * (1+u) */ - if (k) + /* k >= 3 */ + mpfr_set (t, z0, MPFR_RNDN); /* t = z0*(1+u) */ + l = 1; + +/* replace #if 1 by #if 0 for the naive argument reconstruction */ +#if 1 + + /* We multiply by (z0+1)*(z0+2)*...*(z0+k-1) by blocks of j consecutive + terms where j ~ sqrt(k). + If we multiply naively by z0+1, then by z0+2, ..., then by z0+j, + the multiplicative term for the rounding error is (1+u)^(2j). + The multiplicative term is not larger when we multiply by + Z[j] + c[j-1]*Z[j-1] + ... + c[2]*Z[2] + c[1]*z0 + c[0] + with c[p] integers, and Z[p] = z0^p * (1+u)^(p-1). + Note that all terms are positive. + Indeed, since c[1] is exact, c[1]*z0 corresponds to (1+u), + then c[1]*z0 + c[0] corresponds to (1+u)^2, + c[2]*Z[2] + c[1]*z0 + c[0] to (1+u)^3, ..., + c[j-1]*Z[j-1] + ... + c[0] to (1+u)^j, + and Z[j] + c[j-1]*Z[j-1] + ... + c[1]*z0 + c[0] to (1+u)^(j+1). + With the accumulation in t, we get (1+u)^(j+2) and j+2 <= 2j. */ + { + unsigned long j, i, p; + mpfr_t *Z; + mpz_t *c; + for (j = 2; (j + 1) * (j + 1) < k; j++); + /* Z[i] stores z0^i for i <= j */ + Z = (mpfr_t *) mpfr_allocate_func ((j + 1) * sizeof (mpfr_t)); + for (i = 2; i <= j; i++) + mpfr_init2 (Z[i], w); + mpfr_sqr (Z[2], z0, MPFR_RNDN); + for (i = 3; i <= j; i++) + if ((i & 1) == 0) + mpfr_sqr (Z[i], Z[i >> 1], MPFR_RNDN); + else + mpfr_mul (Z[i], Z[i-1], z0, MPFR_RNDN); + c = (mpz_t *) mpfr_allocate_func ((j + 1) * sizeof (mpz_t)); + for (i = 0; i <= j; i++) + mpz_init (c[i]); + for (; l + j <= k; l += j) + { + /* c[i] is the coefficient of x^i in (x+l)*...*(x+l+j-1) */ + mpz_set_ui (c[0], 1); + for (i = 0; i < j; i++) + /* multiply (x+l)*(x+l+1)*...*(x+l+i-1) by x+l+i: + (b[i]*x^i + b[i-1]*x^(i-1) + ... + b[0])*(x+l+i) = + b[i]*x^(i+1) + (b[i-1]+(l+i)*b[i])*x^i + ... + + (b[0]+(l+i)*b[1])*x + i*b[0] */ + { + mpz_set (c[i+1], c[i]); /* b[i]*x^(i+1) */ + for (p = i; p > 0; p--) + { + mpz_mul_ui (c[p], c[p], l + i); + mpz_add (c[p], c[p], c[p-1]); /* b[p-1]+(l+i)*b[p] */ + } + mpz_mul_ui (c[0], c[0], l+i); /* i*b[0] */ + } + /* now compute z0^j + c[j-1]*z0^(j-1) + ... + c[1]*z0 + c[0] */ + mpfr_set_z (u, c[0], MPFR_RNDN); + for (i = 0; i < j; i++) + { + mpfr_mul_z (z, (i == 0) ? z0 : Z[i+1], c[i+1], MPFR_RNDN); + mpfr_add (u, u, z, MPFR_RNDN); + } + mpfr_mul (t, t, u, MPFR_RNDN); + } + for (i = 0; i <= j; i++) + mpz_clear (c[i]); + mpfr_free_func (c, (j + 1) * sizeof (mpz_t)); + for (i = 2; i <= j; i++) + mpfr_clear (Z[i]); + mpfr_free_func (Z, (j + 1) * sizeof (mpfr_t)); + } +#endif /* end of fast argument reconstruction */ + + for (; l < k; l++) { - unsigned long l; - mpfr_set (t, z0, MPFR_RNDN); /* t = z0*(1+u) */ - for (l = 1; l < k; l++) - { - mpfr_add_ui (u, z0, l, MPFR_RNDN); /* u = (z0+l)*(1+u) */ - mpfr_mul (t, t, u, MPFR_RNDN); /* (1+u)^(2l+1) */ - } - /* now t: (1+u)^(2k-1) */ - /* instead of computing log(sqrt(2*Pi)/t), we compute - 1/2*log(2*Pi/t^2), which trades a square root for a square */ - mpfr_mul (t, t, t, MPFR_RNDN); /* (z0*...*(z0+k-1))^2, (1+u)^(4k-1) */ - mpfr_div (v, v, t, MPFR_RNDN); - /* 2*Pi/(z0*...*(z0+k-1))^2 (1+u)^(4k+1) */ + mpfr_add_ui (u, z0, l, MPFR_RNDN); /* u = (z0+l)*(1+u) */ + mpfr_mul (t, t, u, MPFR_RNDN); /* (1+u)^(2l+1) */ } + /* now t: (1+u)^(2k-1) */ + /* instead of computing log(sqrt(2*Pi)/t), we compute + 1/2*log(2*Pi/t^2), which trades a square root for a square */ + mpfr_mul (t, t, t, MPFR_RNDN); /* (z0*...*(z0+k-1))^2, (1+u)^(4k-1) */ + mpfr_div (v, v, t, MPFR_RNDN); + /* 2*Pi/(z0*...*(z0+k-1))^2 (1+u)^(4k+1) */ #ifdef IS_GAMMA err_s = MPFR_GET_EXP(s); mpfr_exp (s, s, MPFR_RNDN); @@ -524,7 +647,7 @@ GAMMA_FUNC (mpfr_ptr y, mpfr_srcptr z0, mpfr_rnd_t rnd) |h| <= (2m+48)*ulp(s), thus exp(s0) = exp(s) * exp(-h). For |h| <= 1/4, we have |exp(h)-1| <= 1.2*|h| thus |exp(s) - exp(s0)| <= 1.2 * exp(s) * (2m+48)* 2^(EXP(s)-w). */ - d = 1.2 * (2.0 * (double) m + 48.0); + /* d = 1.2 * (2.0 * (double) m + 48.0); */ /* the error on s is bounded by d*2^err_s * 2^(-w) */ mpfr_sqrt (t, v, MPFR_RNDN); /* let v0 be the exact value of v. We have v = v0*(1+u)^(4k+1), @@ -533,7 +656,13 @@ GAMMA_FUNC (mpfr_ptr y, mpfr_srcptr z0, mpfr_rnd_t rnd) /* the error on input s is bounded by (1+u)^(d*2^err_s), and that on t is (1+u)^(2k+3/2), thus the total error is (1+u)^(d*2^err_s+2k+5/2) */ - err_s += __gmpfr_ceil_log2 (d); + /* err_s += __gmpfr_ceil_log2 (d); */ + /* since d = 1.2 * (2m+48), ceil(log2(d)) = 2 + ceil(log2(0.6*m+14.4)) + <= 2 + ceil(log2(0.6*m+15)) */ + { + unsigned long mm = (1 + m / 5) * 3; /* 0.6*m <= mm */ + err_s += 2 + __gmpfr_int_ceil_log2 (mm + 15); + } err_t = __gmpfr_ceil_log2 (2.0 * (double) k + 2.5); err_s = (err_s >= err_t) ? err_s + 1 : err_t + 1; #else @@ -570,10 +699,6 @@ GAMMA_FUNC (mpfr_ptr y, mpfr_srcptr z0, mpfr_rnd_t rnd) #ifdef IS_GAMMA end0: #endif - oldBm = Bm; - while (Bm--) - mpz_clear (B[Bm]); - (*__gmp_free_func) (B, oldBm * sizeof (mpz_t)); end: if (inexact == 0) @@ -614,7 +739,7 @@ mpfr_lngamma (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd) else /* lngamma(+/-Inf) = lngamma(nonpositive integer) = +Inf */ { if (!MPFR_IS_INF (x)) - mpfr_set_divby0 (); + MPFR_SET_DIVBY0 (); MPFR_SET_INF (y); MPFR_SET_POS (y); MPFR_RET (0); /* exact */ @@ -654,7 +779,7 @@ mpfr_lgamma (mpfr_ptr y, int *signp, mpfr_srcptr x, mpfr_rnd_t rnd) else { if (MPFR_IS_ZERO (x)) - mpfr_set_divby0 (); + MPFR_SET_DIVBY0 (); *signp = MPFR_INT_SIGN (x); MPFR_SET_INF (y); MPFR_SET_POS (y); @@ -668,7 +793,7 @@ mpfr_lgamma (mpfr_ptr y, int *signp, mpfr_srcptr x, mpfr_rnd_t rnd) { MPFR_SET_INF (y); MPFR_SET_POS (y); - mpfr_set_divby0 (); + MPFR_SET_DIVBY0 (); MPFR_RET (0); } |