diff options
author | Akira Kakuto <kakuto@fuk.kindai.ac.jp> | 2017-12-26 11:50:03 +0000 |
---|---|---|
committer | Akira Kakuto <kakuto@fuk.kindai.ac.jp> | 2017-12-26 11:50:03 +0000 |
commit | 6ddc203e5368cf1d4108e64c18b0bf6c7d53d176 (patch) | |
tree | dda1d9b65325699c60418cfcd02fa5d800cde1f8 /Build/source/libs/mpfr/mpfr-src/src/zeta.c | |
parent | bbfe0bab9260ac7a06a69f710b053dcc73cecad2 (diff) |
mpfr-4.0.0
git-svn-id: svn://tug.org/texlive/trunk@46140 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/zeta.c')
-rw-r--r-- | Build/source/libs/mpfr/mpfr-src/src/zeta.c | 271 |
1 files changed, 231 insertions, 40 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/zeta.c b/Build/source/libs/mpfr/mpfr-src/src/zeta.c index d9d073bf7ad..7931ce2584b 100644 --- a/Build/source/libs/mpfr/mpfr-src/src/zeta.c +++ b/Build/source/libs/mpfr/mpfr-src/src/zeta.c @@ -20,6 +20,8 @@ along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ +#include <float.h> /* for DBL_MAX */ + #define MPFR_NEED_LONGLONG_H #include "mpfr-impl.h" @@ -201,8 +203,7 @@ mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mpfr_rnd_t rnd_mode) where gamma is Euler's constant */ { dint = MAX (d + 3, precs); - MPFR_TRACE (printf ("branch 1\ninternal precision=%lu\n", - (unsigned long) dint)); + /* branch 1, with internal precision dint */ MPFR_GROUP_REPREC_4 (group, dint, b, c, z_pre, f); mpfr_div (z_pre, __gmpfr_one, s1, MPFR_RNDN); mpfr_const_euler (f, MPFR_RNDN); @@ -212,12 +213,12 @@ mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mpfr_rnd_t rnd_mode) { size_t size; - MPFR_TRACE (printf ("branch 2\n")); + /* branch 2 */ /* Computation of parameters n, p and working precision */ dnep = (double) d * LOG2; sd = mpfr_get_d (s, MPFR_RNDN); /* beta = dnep + 0.61 + sd * log (6.2832 / sd); - but a larger value is ok */ + but a larger value is OK */ #define LOG6dot2832 1.83787940484160805532 beta = dnep + 0.61 + sd * (LOG6dot2832 - LOG2 * __gmpfr_floor_log2 (sd)); @@ -232,7 +233,6 @@ mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mpfr_rnd_t rnd_mode) p = 1 + (int) beta / 2; n = 1 + (int) ((sd + 2.0 * (double) p - 1.0) / 6.2832); } - MPFR_TRACE (printf ("\nn=%d\np=%d\n",n,p)); /* add = 4 + floor(1.5 * log(d) / log (2)). We should have add >= 10, which is always fulfilled since d = precz + 11 >= 12, thus ceil(log2(d)) >= 4 */ @@ -242,37 +242,33 @@ mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mpfr_rnd_t rnd_mode) if (dint < precs) dint = precs; - MPFR_TRACE (printf ("internal precision=%lu\n", - (unsigned long) dint)); + /* internal precision is dint */ size = (p + 1) * sizeof(mpfr_t); - tc1 = (mpfr_t*) (*__gmp_allocate_func) (size); + tc1 = (mpfr_t*) mpfr_allocate_func (size); for (l=1; l<=p; l++) mpfr_init2 (tc1[l], dint); MPFR_GROUP_REPREC_4 (group, dint, b, c, z_pre, f); - MPFR_TRACE (printf ("precision of z = %lu\n", - (unsigned long) precz)); + /* precision of z is precz */ /* Computation of the coefficients c_k */ mpfr_zeta_c (p, tc1); - /* Computation of the 3 parts of the fonction Zeta. */ + /* Computation of the 3 parts of the function Zeta. */ mpfr_zeta_part_a (z_pre, s, n); mpfr_zeta_part_b (b, s, n, p, tc1); /* s1 = s-1 is already computed above */ mpfr_div (c, __gmpfr_one, s1, MPFR_RNDN); mpfr_ui_pow (f, n, s1, MPFR_RNDN); mpfr_div (c, c, f, MPFR_RNDN); - MPFR_TRACE (MPFR_DUMP (c)); mpfr_add (z_pre, z_pre, c, MPFR_RNDN); mpfr_add (z_pre, z_pre, b, MPFR_RNDN); for (l=1; l<=p; l++) mpfr_clear (tc1[l]); - (*__gmp_free_func) (tc1, size); + mpfr_free_func (tc1, size); /* End branch 2 */ } - MPFR_TRACE (MPFR_DUMP (z_pre)); if (MPFR_LIKELY (MPFR_CAN_ROUND (z_pre, d-3, precz, rnd_mode))) break; MPFR_ZIV_NEXT (loop, d); @@ -287,11 +283,141 @@ mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mpfr_rnd_t rnd_mode) return inex; } +/* return add = 1 + floor(log(c^3*(13+m1))/log(2)) + where c = (1+eps)*(1+eps*max(8,m1)), + m1 = 1 + max(1/eps,2*sd)*(1+eps), + eps = 2^(-precz-14) + sd = abs(s-1) + */ +static long +compute_add (mpfr_srcptr s, mpfr_prec_t precz) +{ + mpfr_t t, u, m1; + long add; + + mpfr_inits2 (64, t, u, m1, (mpfr_ptr) 0); + if (mpfr_cmp_ui (s, 1) >= 0) + mpfr_sub_ui (t, s, 1, MPFR_RNDU); + else + mpfr_ui_sub (t, 1, s, MPFR_RNDU); + /* now t = sd = abs(s-1), rounded up */ + mpfr_set_ui_2exp (u, 1, - precz - 14, MPFR_RNDU); + /* u = eps */ + /* since 1/eps = 2^(precz+14), if EXP(sd) >= precz+14, then + sd >= 1/2*2^(precz+14) thus 2*sd >= 2^(precz+14) >= 1/eps */ + if (mpfr_get_exp (t) >= precz + 14) + mpfr_mul_2exp (t, t, 1, MPFR_RNDU); + else + mpfr_set_ui_2exp (t, 1, precz + 14, MPFR_RNDU); + /* now t = max(1/eps,2*sd) */ + mpfr_add_ui (u, u, 1, MPFR_RNDU); /* u = 1+eps, rounded up */ + mpfr_mul (t, t, u, MPFR_RNDU); /* t = max(1/eps,2*sd)*(1+eps) */ + mpfr_add_ui (m1, t, 1, MPFR_RNDU); + if (mpfr_get_exp (m1) <= 3) + mpfr_set_ui (t, 8, MPFR_RNDU); + else + mpfr_set (t, m1, MPFR_RNDU); + /* now t = max(8,m1) */ + mpfr_div_2exp (t, t, precz + 14, MPFR_RNDU); /* eps*max(8,m1) */ + mpfr_add_ui (t, t, 1, MPFR_RNDU); /* 1+eps*max(8,m1) */ + mpfr_mul (t, t, u, MPFR_RNDU); /* t = c */ + mpfr_add_ui (u, m1, 13, MPFR_RNDU); /* 13+m1 */ + mpfr_mul (u, u, t, MPFR_RNDU); /* c*(13+m1) */ + mpfr_sqr (t, t, MPFR_RNDU); /* c^2 */ + mpfr_mul (u, u, t, MPFR_RNDU); /* c^3*(13+m1) */ + add = mpfr_get_exp (u); + mpfr_clears (t, u, m1, (mpfr_ptr) 0); + return add; +} + +/* return in z a lower bound (for rnd = RNDD) or upper bound (for rnd = RNDU) + of |zeta(s)|/2, using: + log(|zeta(s)|/2) = (s-1)*log(2*Pi) + lngamma(1-s) + + log(|sin(Pi*s/2)| * zeta(1-s)). + Assumes s < 1/2 and s1 = 1-s exactly, thus s1 > 1/2. + y and p are temporary variables. + At input, p is Pi rounded down. + The comments in the code are for rnd = RNDD. */ +static void +mpfr_reflection_overflow (mpfr_t z, mpfr_t s1, const mpfr_t s, mpfr_t y, + mpfr_t p, mpfr_rnd_t rnd) +{ + mpz_t sint; + + MPFR_ASSERTD (rnd == MPFR_RNDD || rnd == MPFR_RNDU); + + /* Since log is increasing, we want lower bounds on |sin(Pi*s/2)| and + zeta(1-s). */ + mpz_init (sint); + mpfr_get_z (sint, s, MPFR_RNDD); /* sint = floor(s) */ + /* We first compute a lower bound of |sin(Pi*s/2)|, which is a periodic + function of period 2. Thus: + if 2k < s < 2k+1, then |sin(Pi*s/2)| is increasing; + if 2k-1 < s < 2k, then |sin(Pi*s/2)| is decreasing. + These cases are distinguished by testing bit 0 of floor(s) as if + represented in two's complement (or equivalently, as an unsigned + integer mod 2): + 0: sint = 0 mod 2, thus 2k < s < 2k+1 and |sin(Pi*s/2)| is increasing; + 1: sint = 1 mod 2, thus 2k-1 < s < 2k and |sin(Pi*s/2)| is decreasing. + Let's recall that the comments are for rnd = RNDD. */ + if (mpz_tstbit (sint, 0) == 0) /* |sin(Pi*s/2)| is increasing: round down + Pi*s to get a lower bound. */ + { + mpfr_mul (y, p, s, rnd); + if (rnd == MPFR_RNDD) + mpfr_nextabove (p); /* we will need p rounded above afterwards */ + } + else /* |sin(Pi*s/2)| is decreasing: round up Pi*s to get a lower bound. */ + { + if (rnd == MPFR_RNDD) + mpfr_nextabove (p); + mpfr_mul (y, p, s, MPFR_INVERT_RND(rnd)); + } + mpfr_div_2ui (y, y, 1, MPFR_RNDN); /* exact, rounding mode doesn't matter */ + /* The rounding direction of sin depends on its sign. We have: + if -4k-2 < s < -4k, then -2k-1 < s/2 < -2k, thus sin(Pi*s/2) < 0; + if -4k < s < -4k+2, then -2k < s/2 < -2k+1, thus sin(Pi*s/2) > 0. + These cases are distinguished by testing bit 1 of floor(s) as if + represented in two's complement (or equivalently, as an unsigned + integer mod 4): + 0: sint = {0,1} mod 4, thus -2k < s/2 < -2k+1 and sin(Pi*s/2) > 0; + 1: sint = {2,3} mod 4, thus -2k-1 < s/2 < -2k and sin(Pi*s/2) < 0. + Let's recall that the comments are for rnd = RNDD. */ + if (mpz_tstbit (sint, 1) == 0) /* -2k < s/2 < -2k+1; sin(Pi*s/2) > 0 */ + { + /* Round sin down to get a lower bound of |sin(Pi*s/2)|. */ + mpfr_sin (y, y, rnd); + } + else /* -2k-1 < s/2 < -2k; sin(Pi*s/2) < 0 */ + { + /* Round sin up to get a lower bound of |sin(Pi*s/2)|. */ + mpfr_sin (y, y, MPFR_INVERT_RND(rnd)); + mpfr_abs (y, y, MPFR_RNDN); /* exact, rounding mode doesn't matter */ + } + mpz_clear (sint); + /* now y <= |sin(Pi*s/2)| when rnd=RNDD, y >= |sin(Pi*s/2)| when rnd=RNDU */ + mpfr_zeta_pos (z, s1, rnd); /* zeta(1-s) */ + mpfr_mul (z, z, y, rnd); + /* now z <= |sin(Pi*s/2)|*zeta(1-s) */ + mpfr_log (z, z, rnd); + /* now z <= log(|sin(Pi*s/2)|*zeta(1-s)) */ + mpfr_lngamma (y, s1, rnd); + mpfr_add (z, z, y, rnd); + /* z <= lngamma(1-s) + log(|sin(Pi*s/2)|*zeta(1-s)) */ + /* since s-1 < 0, we want to round log(2*pi) upwards */ + mpfr_mul_2ui (y, p, 1, MPFR_INVERT_RND(rnd)); + mpfr_log (y, y, MPFR_INVERT_RND(rnd)); + mpfr_mul (y, y, s1, MPFR_INVERT_RND(rnd)); + mpfr_sub (z, z, y, rnd); + mpfr_exp (z, z, rnd); + if (rnd == MPFR_RNDD) + mpfr_nextbelow (p); /* restore original p */ +} + int mpfr_zeta (mpfr_t z, mpfr_srcptr s, mpfr_rnd_t rnd_mode) { mpfr_t z_pre, s1, y, p; - double sd, eps, m1, c; long add; mpfr_prec_t precz, prec1, precs, precs1; int inex; @@ -328,10 +454,9 @@ mpfr_zeta (mpfr_t z, mpfr_srcptr s, mpfr_rnd_t rnd_mode) /* s is neither Nan, nor Inf, nor Zero */ /* check tiny s: we have zeta(s) = -1/2 - 1/2 log(2 Pi) s + ... around s=0, - and for |s| <= 0.074, we have |zeta(s) + 1/2| <= |s|. - Thus if |s| <= 1/4*ulp(1/2), we can deduce the correct rounding - (the 1/4 covers the case where |zeta(s)| < 1/2 and rounding to nearest). - A sufficient condition is that EXP(s) + 1 < -PREC(z). */ + and for |s| <= 2^(-4), we have |zeta(s) + 1/2| <= |s|. + EXP(s) + 1 < -PREC(z) is a sufficient condition to be able to round + correctly, for any PREC(z) >= 1 (see algorithms.tex for details). */ if (MPFR_GET_EXP (s) + 1 < - (mpfr_exp_t) MPFR_PREC(z)) { int signs = MPFR_SIGN(s); @@ -382,7 +507,7 @@ mpfr_zeta (mpfr_t z, mpfr_srcptr s, mpfr_rnd_t rnd_mode) { MPFR_SET_INF (z); MPFR_SET_POS (z); - mpfr_set_divby0 (); + MPFR_SET_DIVBY0 (); MPFR_RET (0); } @@ -402,39 +527,96 @@ mpfr_zeta (mpfr_t z, mpfr_srcptr s, mpfr_rnd_t rnd_mode) /* Precision precs1 needed to represent 1 - s, and s + 2, without any truncation */ precs1 = precs + 2 + MAX (0, - MPFR_GET_EXP (s)); - sd = mpfr_get_d (s, MPFR_RNDN) - 1.0; - if (sd < 0.0) - sd = -sd; /* now sd = abs(s-1.0) */ /* Precision prec1 is the precision on elementary computations; it ensures a final precision prec1 - add for zeta(s) */ - /* eps = pow (2.0, - (double) precz - 14.0); */ - eps = __gmpfr_ceil_exp2 (- (double) precz - 14.0); - m1 = 1.0 + MAX(1.0 / eps, 2.0 * sd) * (1.0 + eps); - c = (1.0 + eps) * (1.0 + eps * MAX(8.0, m1)); - /* add = 1 + floor(log(c*c*c*(13 + m1))/log(2)); */ - add = __gmpfr_ceil_log2 (c * c * c * (13.0 + m1)); + add = compute_add (s, precz); prec1 = precz + add; + /* FIXME: To avoid that the working precision (prec1) depends on the + input precision, one would need to take into account the error made + when s1 is not exactly 1-s when computing zeta(s1) and gamma(s1) + below, and also in the case y=Inf (i.e. when gamma(s1) overflows). + Make sure that underflows do not occur in intermediate computations. + Due to the limited precision, they are probably not possible + in practice; add some MPFR_ASSERTN's to be sure that problems + do not remain undetected? */ prec1 = MAX (prec1, precs1) + 10; MPFR_GROUP_INIT_4 (group, prec1, z_pre, s1, y, p); MPFR_ZIV_INIT (loop, prec1); for (;;) { + mpfr_exp_t ey; + mpfr_t z_up; + + mpfr_const_pi (p, MPFR_RNDD); /* p is Pi */ + mpfr_sub (s1, __gmpfr_one, s, MPFR_RNDN); /* s1 = 1-s */ - mpfr_zeta_pos (z_pre, s1, MPFR_RNDN); /* zeta(1-s) */ mpfr_gamma (y, s1, MPFR_RNDN); /* gamma(1-s) */ - if (MPFR_IS_INF (y)) /* Zeta(s) < 0 for -4k-2 < s < -4k, - Zeta(s) > 0 for -4k < s < -4k+2 */ + if (MPFR_IS_INF (y)) /* zeta(s) < 0 for -4k-2 < s < -4k, + zeta(s) > 0 for -4k < s < -4k+2 */ { - mpfr_div_2ui (s1, s, 2, MPFR_RNDN); /* s/4, exact */ - mpfr_frac (s1, s1, MPFR_RNDN); /* exact, -1 < s1 < 0 */ - overflow = (mpfr_cmp_si_2exp (s1, -1, -1) > 0) ? -1 : 1; - break; + /* FIXME: An overflow in gamma(s1) does not imply that + zeta(s) will overflow. A solution: + 1. Compute + log(|zeta(s)|/2) = (s-1)*log(2*pi) + lngamma(1-s) + + log(abs(sin(Pi*s/2)) * zeta(1-s)) + (possibly sharing computations with the normal case) + with a rather good accuracy (see (2)). + Memorize the sign of sin(...) for the final sign. + 2. Take the exponential, ~= |zeta(s)|/2. If there is an + overflow, then this means an overflow on the final result + (due to the multiplication by 2, which has not been done + yet). + 3. Ziv test. + 4. Correct the sign from the sign of sin(...). + 5. Round then multiply by 2. Here, an overflow in either + operation means a real overflow. */ + mpfr_reflection_overflow (z_pre, s1, s, y, p, MPFR_RNDD); + /* z_pre is a lower bound of |zeta(s)|/2, thus if it overflows, + or has exponent emax, then |zeta(s)| overflows too. */ + if (MPFR_IS_INF (z_pre) || MPFR_GET_EXP(z_pre) == __gmpfr_emax) + { /* determine the sign of overflow */ + mpfr_div_2ui (s1, s, 2, MPFR_RNDN); /* s/4, exact */ + mpfr_frac (s1, s1, MPFR_RNDN); /* exact, -1 < s1 < 0 */ + overflow = (mpfr_cmp_si_2exp (s1, -1, -1) > 0) ? -1 : 1; + break; + } + else /* EXP(z_pre) < __gmpfr_emax */ + { + int ok = 0; + mpfr_t z_down; + mpfr_init2 (z_up, mpfr_get_prec (z_pre)); + mpfr_reflection_overflow (z_up, s1, s, y, p, MPFR_RNDU); + /* if the lower approximation z_pre does not overflow, but + z_up does, we need more precision */ + if (MPFR_IS_INF (z_up) || MPFR_GET_EXP(z_up) == __gmpfr_emax) + goto next_loop; + /* check if z_pre and z_up round to the same number */ + mpfr_init2 (z_down, precz); + mpfr_set (z_down, z_pre, rnd_mode); + /* Note: it might be that EXP(z_down) = emax here, in that + case we will have overflow below when we multiply by 2 */ + mpfr_prec_round (z_up, precz, rnd_mode); + ok = mpfr_cmp (z_down, z_up) == 0; + mpfr_clear (z_up); + mpfr_clear (z_down); + if (ok) + { + /* get correct sign and multiply by 2 */ + mpfr_div_2ui (s1, s, 2, MPFR_RNDN); /* s/4, exact */ + mpfr_frac (s1, s1, MPFR_RNDN); /* exact, -1 < s1 < 0 */ + if (mpfr_cmp_si_2exp (s1, -1, -1) > 0) + mpfr_neg (z_pre, z_pre, rnd_mode); + mpfr_mul_2ui (z_pre, z_pre, 1, rnd_mode); + break; + } + else + goto next_loop; + } } + mpfr_zeta_pos (z_pre, s1, MPFR_RNDN); /* zeta(1-s) */ mpfr_mul (z_pre, z_pre, y, MPFR_RNDN); /* gamma(1-s)*zeta(1-s) */ - mpfr_const_pi (p, MPFR_RNDD); /* p is Pi */ - /* multiply z_pre by 2^s*Pi^(s-1) where p=Pi, s1=1-s */ mpfr_mul_2ui (y, p, 1, MPFR_RNDN); /* 2*Pi */ mpfr_neg (s1, s1, MPFR_RNDN); /* s-1 */ @@ -445,11 +627,19 @@ mpfr_zeta (mpfr_t z, mpfr_srcptr s, mpfr_rnd_t rnd_mode) /* multiply z_pre by sin(Pi*s/2) */ mpfr_mul (y, s, p, MPFR_RNDN); mpfr_div_2ui (p, y, 1, MPFR_RNDN); /* p = s*Pi/2 */ + /* FIXME: sinpi will be available, we should replace the mpfr_sin + call below by mpfr_sinpi(s/2), where s/2 will be exact. + Can mpfr_sin underflow? Moreover, the code below should be + improved so that the "if" condition becomes unlikely, e.g. + by taking a slightly larger working precision. */ mpfr_sin (y, p, MPFR_RNDN); /* y = sin(Pi*s/2) */ - if (MPFR_GET_EXP(y) < 0) /* take account of cancellation in sin(p) */ + ey = MPFR_GET_EXP (y); + if (ey < 0) /* take account of cancellation in sin(p) */ { mpfr_t t; - mpfr_init2 (t, prec1 - MPFR_GET_EXP(y)); + + MPFR_ASSERTN (- ey < MPFR_PREC_MAX - prec1); + mpfr_init2 (t, prec1 - ey); mpfr_const_pi (t, MPFR_RNDD); mpfr_mul (t, s, t, MPFR_RNDN); mpfr_div_2ui (t, t, 1, MPFR_RNDN); @@ -462,6 +652,7 @@ mpfr_zeta (mpfr_t z, mpfr_srcptr s, mpfr_rnd_t rnd_mode) rnd_mode))) break; + next_loop: MPFR_ZIV_NEXT (loop, prec1); MPFR_GROUP_REPREC_4 (group, prec1, z_pre, s1, y, p); } |