summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src/zeta.c
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/zeta.c')
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/zeta.c271
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);
}