dist/src/zeta.c

/* mpfr_zeta -- compute the Riemann Zeta function

Copyright 2003-2018 Free Software Foundation, Inc.
Contributed by the AriC and Caramba projects, INRIA.

This file is part of the GNU MPFR Library.

The GNU MPFR Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 3 of the License, or (at your
option) any later version.

The GNU MPFR Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
License for more details.

You should have received a copy of the GNU Lesser General Public License
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"

/*
   Parameters:
   s - the input floating-point number
   n, p - parameters from the algorithm
   tc - an array of p floating-point numbers tc[1]..tc[p]
   Output:
   b is the result, i.e.
   sum(tc[i]*product((s+2j)*(s+2j-1)/n^2,j=1..i-1), i=1..p)*s*n^(-s-1)
*/
static void
mpfr_zeta_part_b (mpfr_t b, mpfr_srcptr s, int n, int p, mpfr_t *tc)
{
  mpfr_t s1, d, u;
  unsigned long n2;
  int l, t;
  MPFR_GROUP_DECL (group);

  if (p == 0)
    {
      MPFR_SET_ZERO (b);
      MPFR_SET_POS (b);
      return;
    }

  n2 = n * n;
  MPFR_GROUP_INIT_3 (group, MPFR_PREC (b), s1, d, u);

  /* t equals 2p-2, 2p-3, ... ; s1 equals s+t */
  t = 2 * p - 2;
  mpfr_set (d, tc[p], MPFR_RNDN);
  for (l = 1; l < p; l++)
    {
      mpfr_add_ui (s1, s, t, MPFR_RNDN); /* s + (2p-2l) */
      mpfr_mul (d, d, s1, MPFR_RNDN);
      t = t - 1;
      mpfr_add_ui (s1, s, t, MPFR_RNDN); /* s + (2p-2l-1) */
      mpfr_mul (d, d, s1, MPFR_RNDN);
      t = t - 1;
      mpfr_div_ui (d, d, n2, MPFR_RNDN);
      mpfr_add (d, d, tc[p-l], MPFR_RNDN);
      /* since s is positive and the tc[i] have alternate signs,
         the following is unlikely */
      if (MPFR_UNLIKELY (mpfr_cmpabs (d, tc[p-l]) > 0))
        mpfr_set (d, tc[p-l], MPFR_RNDN);
    }
  mpfr_mul (d, d, s, MPFR_RNDN);
  mpfr_add (s1, s, __gmpfr_one, MPFR_RNDN);
  mpfr_neg (s1, s1, MPFR_RNDN);
  mpfr_ui_pow (u, n, s1, MPFR_RNDN);
  mpfr_mul (b, d, u, MPFR_RNDN);

  MPFR_GROUP_CLEAR (group);
}

/* Input: p - an integer
   Output: fills tc[1..p], tc[i] = bernoulli(2i)/(2i)!
   tc[1]=1/12, tc[2]=-1/720, tc[3]=1/30240, ...
*/
static void
mpfr_zeta_c (int p, mpfr_t *tc)
{
  mpfr_t d;
  int k, l;

  if (p > 0)
    {
      mpfr_init2 (d, MPFR_PREC (tc[1]));
      mpfr_div_ui (tc[1], __gmpfr_one, 12, MPFR_RNDN);
      for (k = 2; k <= p; k++)
        {
          mpfr_set_ui (d, k-1, MPFR_RNDN);
          mpfr_div_ui (d, d, 12*k+6, MPFR_RNDN);
          for (l=2; l < k; l++)
            {
              mpfr_div_ui (d, d, 4*(2*k-2*l+3)*(2*k-2*l+2), MPFR_RNDN);
              mpfr_add (d, d, tc[l], MPFR_RNDN);
            }
          mpfr_div_ui (tc[k], d, 24, MPFR_RNDN);
          MPFR_CHANGE_SIGN (tc[k]);
        }
      mpfr_clear (d);
    }
}

/* Input: s - a floating-point number
          n - an integer
   Output: sum - a floating-point number approximating sum(1/i^s, i=1..n-1) */
static void
mpfr_zeta_part_a (mpfr_t sum, mpfr_srcptr s, int n)
{
  mpfr_t u, s1;
  int i;
  MPFR_GROUP_DECL (group);

  MPFR_GROUP_INIT_2 (group, MPFR_PREC (sum), u, s1);

  mpfr_neg (s1, s, MPFR_RNDN);
  mpfr_ui_pow (u, n, s1, MPFR_RNDN);
  mpfr_div_2ui (u, u, 1, MPFR_RNDN);
  mpfr_set (sum, u, MPFR_RNDN);
  for (i=n-1; i>1; i--)
    {
      mpfr_ui_pow (u, i, s1, MPFR_RNDN);
      mpfr_add (sum, sum, u, MPFR_RNDN);
    }
  mpfr_add (sum, sum, __gmpfr_one, MPFR_RNDN);

  MPFR_GROUP_CLEAR (group);
}

/* Input: s - a floating-point number >= 1/2.
          rnd_mode - a rounding mode.
          Assumes s is neither NaN nor Infinite.
   Output: z - Zeta(s) rounded to the precision of z with direction rnd_mode
*/
static int
mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mpfr_rnd_t rnd_mode)
{
  mpfr_t b, c, z_pre, f, s1;
  double beta, sd, dnep;
  mpfr_t *tc1;
  mpfr_prec_t precz, precs, d, dint;
  int p, n, l, add;
  int inex;
  MPFR_GROUP_DECL (group);
  MPFR_ZIV_DECL (loop);

  MPFR_ASSERTD (MPFR_IS_POS (s) && MPFR_GET_EXP (s) >= 0);

  precz = MPFR_PREC (z);
  precs = MPFR_PREC (s);

  /* Zeta(x) = 1+1/2^x+1/3^x+1/4^x+1/5^x+O(1/6^x)
     so with 2^(EXP(x)-1) <= x < 2^EXP(x)
     So for x > 2^3, k^x > k^8, so 2/k^x < 2/k^8
     Zeta(x) = 1 + 1/2^x*(1+(2/3)^x+(2/4)^x+...)
             = 1 + 1/2^x*(1+sum((2/k)^x,k=3..infinity))
            <= 1 + 1/2^x*(1+sum((2/k)^8,k=3..infinity))
     And sum((2/k)^8,k=3..infinity) = -257+128*Pi^8/4725 ~= 0.0438035
     So Zeta(x) <= 1 + 1/2^x*2 for x >= 8
     The error is < 2^(-x+1) <= 2^(-2^(EXP(x)-1)+1) */
  if (MPFR_GET_EXP (s) > 3)
    {
      mpfr_exp_t err;
      err = MPFR_GET_EXP (s) - 1;
      if (err > (mpfr_exp_t) (sizeof (mpfr_exp_t)*CHAR_BIT-2))
        err = MPFR_EMAX_MAX;
      else
        err = ((mpfr_exp_t)1) << err;
      err = 1 - (-err+1); /* GET_EXP(one) - (-err+1) = err :) */
      MPFR_FAST_COMPUTE_IF_SMALL_INPUT (z, __gmpfr_one, err, 0, 1,
                                        rnd_mode, {});
    }

  d = precz + MPFR_INT_CEIL_LOG2(precz) + 10;

  /* we want that s1 = s-1 is exact, i.e. we should have PREC(s1) >= EXP(s) */
  dint = (mpfr_uexp_t) MPFR_GET_EXP (s);
  mpfr_init2 (s1, MAX (precs, dint));
  inex = mpfr_sub (s1, s, __gmpfr_one, MPFR_RNDN);
  MPFR_ASSERTD (inex == 0);

  /* case s=1 should have already been handled */
  MPFR_ASSERTD (!MPFR_IS_ZERO (s1));

  MPFR_GROUP_INIT_4 (group, MPFR_PREC_MIN, b, c, z_pre, f);

  MPFR_ZIV_INIT (loop, d);
  for (;;)
    {
      /* Principal loop: we compute, in z_pre,
         an approximation of Zeta(s), that we send to can_round */
      if (MPFR_GET_EXP (s1) <= -(mpfr_exp_t) ((mpfr_prec_t) (d-3)/2))
        /* Branch 1: when s-1 is very small, one
           uses the approximation Zeta(s)=1/(s-1)+gamma,
           where gamma is Euler's constant */
        {
          dint = MAX (d + 3, precs);
          /* 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);
          mpfr_add (z_pre, z_pre, f, MPFR_RNDN);
        }
      else /* Branch 2 */
        {
          size_t size;

          /* 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 */
#define LOG6dot2832 1.83787940484160805532
          beta = dnep + 0.61 + sd * (LOG6dot2832 - LOG2 *
                                     __gmpfr_floor_log2 (sd));
          if (beta <= 0.0)
            {
              p = 0;
              /* n = 1 + (int) (exp ((dnep - LOG2) / sd)); */
              n = 1 + (int) __gmpfr_ceil_exp2 ((d - 1.0) / sd);
            }
          else
            {
              p = 1 + (int) beta / 2;
              n = 1 + (int) ((sd + 2.0 * (double) p - 1.0) / 6.2832);
            }
          /* 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 */
          add = 4 + (3 * MPFR_INT_CEIL_LOG2 (d)) / 2;
          MPFR_ASSERTD(add >= 10);
          dint = d + add;
          if (dint < precs)
            dint = precs;

          /* internal precision is dint */

          size = (p + 1) * sizeof(mpfr_t);
          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);

          /* precision of z is precz */

          /* Computation of the coefficients c_k */
          mpfr_zeta_c (p, tc1);
          /* 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_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]);
          mpfr_free_func (tc1, size);
          /* End branch 2 */
        }

      if (MPFR_LIKELY (MPFR_CAN_ROUND (z_pre, d-3, precz, rnd_mode)))
        break;
      MPFR_ZIV_NEXT (loop, d);
    }
  MPFR_ZIV_FREE (loop);

  inex = mpfr_set (z, z_pre, rnd_mode);

  MPFR_GROUP_CLEAR (group);
  mpfr_clear (s1);

  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;
  long add;
  mpfr_prec_t precz, prec1, precs, precs1;
  int inex;
  MPFR_GROUP_DECL (group);
  MPFR_ZIV_DECL (loop);
  MPFR_SAVE_EXPO_DECL (expo);

  MPFR_LOG_FUNC (
    ("s[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (s), mpfr_log_prec, s, rnd_mode),
    ("z[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (z), mpfr_log_prec, z, inex));

  /* Zero, Nan or Inf ? */
  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (s)))
    {
      if (MPFR_IS_NAN (s))
        {
          MPFR_SET_NAN (z);
          MPFR_RET_NAN;
        }
      else if (MPFR_IS_INF (s))
        {
          if (MPFR_IS_POS (s))
            return mpfr_set_ui (z, 1, MPFR_RNDN); /* Zeta(+Inf) = 1 */
          MPFR_SET_NAN (z); /* Zeta(-Inf) = NaN */
          MPFR_RET_NAN;
        }
      else /* s iz zero */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO (s));
          return mpfr_set_si_2exp (z, -1, -1, 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| <= 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);

      MPFR_SAVE_EXPO_MARK (expo);
      mpfr_set_si_2exp (z, -1, -1, rnd_mode); /* -1/2 */
      if (rnd_mode == MPFR_RNDA)
        rnd_mode = MPFR_RNDD; /* the result is around -1/2, thus negative */
      if ((rnd_mode == MPFR_RNDU || rnd_mode == MPFR_RNDZ) && signs < 0)
        {
          mpfr_nextabove (z); /* z = -1/2 + epsilon */
          inex = 1;
        }
      else if (rnd_mode == MPFR_RNDD && signs > 0)
        {
          mpfr_nextbelow (z); /* z = -1/2 - epsilon */
          inex = -1;
        }
      else
        {
          if (rnd_mode == MPFR_RNDU) /* s > 0: z = -1/2 */
            inex = 1;
          else if (rnd_mode == MPFR_RNDD)
            inex = -1;              /* s < 0: z = -1/2 */
          else /* (MPFR_RNDZ and s > 0) or MPFR_RNDN: z = -1/2 */
            inex = (signs > 0) ? 1 : -1;
        }
      MPFR_SAVE_EXPO_FREE (expo);
      return mpfr_check_range (z, inex, rnd_mode);
    }

  /* Check for case s= -2n */
  if (MPFR_IS_NEG (s))
    {
      mpfr_t tmp;
      tmp[0] = *s;
      MPFR_EXP (tmp) = MPFR_GET_EXP (s) - 1;
      if (mpfr_integer_p (tmp))
        {
          MPFR_SET_ZERO (z);
          MPFR_SET_POS (z);
          MPFR_RET (0);
        }
    }

  /* Check for case s=1 before changing the exponent range */
  if (mpfr_cmp (s, __gmpfr_one) == 0)
    {
      MPFR_SET_INF (z);
      MPFR_SET_POS (z);
      MPFR_SET_DIVBY0 ();
      MPFR_RET (0);
    }

  MPFR_SAVE_EXPO_MARK (expo);

  /* Compute Zeta */
  if (MPFR_IS_POS (s) && MPFR_GET_EXP (s) >= 0) /* Case s >= 1/2 */
    inex = mpfr_zeta_pos (z, s, rnd_mode);
  else /* use reflection formula
          zeta(s) = 2^s*Pi^(s-1)*sin(Pi*s/2)*gamma(1-s)*zeta(1-s) */
    {
      int overflow = 0;

      precz = MPFR_PREC (z);
      precs = MPFR_PREC (s);

      /* Precision precs1 needed to represent 1 - s, and s + 2,
         without any truncation */
      precs1 = precs + 2 + MAX (0, - MPFR_GET_EXP (s));
      /* Precision prec1 is the precision on elementary computations;
         it ensures a final precision prec1 - add for zeta(s) */
      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_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 */
            {
              /* 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) */

          /* 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 */
          mpfr_pow (y, y, s1, MPFR_RNDN);         /* (2*Pi)^(s-1) */
          mpfr_mul (z_pre, z_pre, y, MPFR_RNDN);
          mpfr_mul_2ui (z_pre, z_pre, 1, MPFR_RNDN);

          /* 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) */
          ey = MPFR_GET_EXP (y);
          if (ey < 0) /* take account of cancellation in sin(p) */
            {
              mpfr_t t;

              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);
              mpfr_sin (y, t, MPFR_RNDN);
              mpfr_clear (t);
            }
          mpfr_mul (z_pre, z_pre, y, MPFR_RNDN);

          if (MPFR_LIKELY (MPFR_CAN_ROUND (z_pre, prec1 - add, precz,
                                           rnd_mode)))
            break;

        next_loop:
          MPFR_ZIV_NEXT (loop, prec1);
          MPFR_GROUP_REPREC_4 (group, prec1, z_pre, s1, y, p);
        }
      MPFR_ZIV_FREE (loop);
      if (overflow != 0)
        {
          inex = mpfr_overflow (z, rnd_mode, overflow);
          MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW);
        }
      else
        inex = mpfr_set (z, z_pre, rnd_mode);
      MPFR_GROUP_CLEAR (group);
    }

  MPFR_SAVE_EXPO_FREE (expo);
  return mpfr_check_range (z, inex, rnd_mode);
}