/* mpfr_zeta -- compute the Riemann Zeta function

Copyright 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
Contributed by Jean-Luc Re'my and the Spaces project, INRIA Lorraine.

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 2.1 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.LIB.  If not, write to
the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
MA 02110-1301, USA. */

#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], GMP_RNDN);
  for (l = 1; l < p; l++)
    {
      mpfr_add_ui (s1, s, t, GMP_RNDN); /* s + (2p-2l) */
      mpfr_mul (d, d, s1, GMP_RNDN);
      t = t - 1;
      mpfr_add_ui (s1, s, t, GMP_RNDN); /* s + (2p-2l-1) */
      mpfr_mul (d, d, s1, GMP_RNDN);
      t = t - 1;
      mpfr_div_ui (d, d, n2, GMP_RNDN);
      mpfr_add (d, d, tc[p-l], GMP_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], GMP_RNDN);
    }
  mpfr_mul (d, d, s, GMP_RNDN);
  mpfr_add (s1, s, __gmpfr_one, GMP_RNDN);
  mpfr_neg (s1, s1, GMP_RNDN);
  mpfr_ui_pow (u, n, s1, GMP_RNDN);
  mpfr_mul (b, d, u, GMP_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, GMP_RNDN);
      for (k = 2; k <= p; k++)
        {
          mpfr_set_ui (d, k-1, GMP_RNDN);
          mpfr_div_ui (d, d, 12*k+6, GMP_RNDN);
          for (l=2; l < k; l++)
            {
              mpfr_div_ui (d, d, 4*(2*k-2*l+3)*(2*k-2*l+2), GMP_RNDN);
              mpfr_add (d, d, tc[l], GMP_RNDN);
            }
          mpfr_div_ui (tc[k], d, 24, GMP_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, GMP_RNDN);
  mpfr_ui_pow (u, n, s1, GMP_RNDN);
  mpfr_div_2ui (u, u, 1, GMP_RNDN);
  mpfr_set (sum, u, GMP_RNDN);
  for (i=n-1; i>1; i--)
    {
      mpfr_ui_pow (u, i, s1, GMP_RNDN);
      mpfr_add (sum, sum, u, GMP_RNDN);
    }
  mpfr_add (sum, sum, __gmpfr_one, GMP_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, mp_rnd_t rnd_mode)
{
  mpfr_t b, c, z_pre, f, s1;
  double beta, sd, dnep;
  mpfr_t *tc1;
  mp_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)
    {
      mp_exp_t err;
      err = MPFR_GET_EXP (s) - 1;
      if (err > (mp_exp_t) (sizeof (mp_exp_t)*CHAR_BIT-2))
        err = MPFR_EMAX_MAX;
      else
        err = ((mp_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, GMP_RNDN);
  MPFR_ASSERTD (inex == 0);

  /* case s=1 */
  if (MPFR_IS_ZERO (s1))
    {
      MPFR_SET_INF (z);
      MPFR_SET_POS (z);
      MPFR_ASSERTD (inex == 0);
      goto clear_and_return;
    }

  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) <= -(mp_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);
          MPFR_TRACE (printf ("branch 1\ninternal precision=%lu\n",
                              (unsigned long) dint));
          MPFR_GROUP_REPREC_4 (group, dint, b, c, z_pre, f);
          mpfr_div (z_pre, __gmpfr_one, s1, GMP_RNDN);
          mpfr_const_euler (f, GMP_RNDN);
          mpfr_add (z_pre, z_pre, f, GMP_RNDN);
        }
      else /* Branch 2 */
        {
          size_t size;

          MPFR_TRACE (printf ("branch 2\n"));
          /* Computation of parameters n, p and working precision */
          dnep = (double) d * LOG2;
          sd = mpfr_get_d (s, GMP_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);
            }
          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 */
          add = 4 + (3 * MPFR_INT_CEIL_LOG2 (d)) / 2;
          MPFR_ASSERTD(add >= 10);
          dint = d + add;
          if (dint < precs)
            dint = precs;

          MPFR_TRACE (printf ("internal precision=%lu\n",
                              (unsigned long) dint));

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

          /* Computation of the coefficients c_k */
          mpfr_zeta_c (p, tc1);
          /* Computation of the 3 parts of the fonction 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, GMP_RNDN);
          mpfr_ui_pow (f, n, s1, GMP_RNDN);
          mpfr_div (c, c, f, GMP_RNDN);
          MPFR_TRACE (MPFR_DUMP (c));
          mpfr_add (z_pre, z_pre, c, GMP_RNDN);
          mpfr_add (z_pre, z_pre, b, GMP_RNDN);
          for (l=1; l<=p; l++)
            mpfr_clear (tc1[l]);
          (*__gmp_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);
    }
  MPFR_ZIV_FREE (loop);

  inex = mpfr_set (z, z_pre, rnd_mode);

  MPFR_GROUP_CLEAR (group);
 clear_and_return:
  mpfr_clear (s1);

  return inex;
}

int
mpfr_zeta (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode)
{
  mpfr_t z_pre, s1, y, p;
  double sd, eps, m1, c;
  long add;
  mp_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[%#R]=%R rnd=%d", s, s, rnd_mode),
                 ("z[%#R]=%R inexact=%d", z, 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, GMP_RNDN); /* Zeta(+Inf) = 1 */
          MPFR_SET_NAN (z); /* Zeta(-Inf) = NaN */
          MPFR_RET_NAN;
        }
      else /* s iz zero */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO (s));
          mpfr_set_ui (z, 1, rnd_mode);
          mpfr_div_2ui (z, z, 1, rnd_mode);
          MPFR_CHANGE_SIGN (z);
          MPFR_RET (0);
        }
    }

  /* 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). */
  if (MPFR_EXP(s) + 1 < - (mp_exp_t) MPFR_PREC(z))
    {
      int signs = MPFR_SIGN(s);
      mpfr_set_si_2exp (z, -1, -1, rnd_mode); /* -1/2 */
      if ((rnd_mode == GMP_RNDU || rnd_mode == GMP_RNDZ) && signs < 0)
        {
          mpfr_nextabove (z); /* z = -1/2 + epsilon */
          inex = 1;
        }
      else if (rnd_mode == GMP_RNDD && signs > 0)
        {
          mpfr_nextbelow (z); /* z = -1/2 - epsilon */
          inex = -1;
        }
      else
        {
          if (rnd_mode == GMP_RNDU) /* s > 0: z = -1/2 */
            inex = 1;
          else if (rnd_mode == GMP_RNDD)
            inex = -1;              /* s < 0: z = -1/2 */
          else /* (GMP_RNDZ and s > 0) or GMP_RNDN: z = -1/2 */
            inex = (signs > 0) ? 1 : -1;
        }
      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_EXP (s) - 1;
      if (mpfr_integer_p (tmp))
        {
          MPFR_SET_ZERO (z);
          MPFR_SET_POS (z);
          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));
      sd = mpfr_get_d (s, GMP_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));
      prec1 = precz + add;
      prec1 = MAX (prec1, precs1) + 10;

      MPFR_GROUP_INIT_4 (group, prec1, z_pre, s1, y, p);
      MPFR_ZIV_INIT (loop, prec1);
      for (;;)
        {
          mpfr_sub (s1, __gmpfr_one, s, GMP_RNDN);/* s1 = 1-s */
          mpfr_zeta_pos (z_pre, s1, GMP_RNDN);   /* zeta(1-s)  */
          mpfr_gamma (y, s1, GMP_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 */
            {
              mpfr_div_2ui (s1, s, 2, GMP_RNDN); /* s/4, exact */
              mpfr_frac (s1, s1, GMP_RNDN); /* exact, -1 < s1 < 0 */
              overflow = (mpfr_cmp_si_2exp (s1, -1, -1) > 0) ? -1 : 1;
              break;
            }
          mpfr_mul (z_pre, z_pre, y, GMP_RNDN);  /* gamma(1-s)*zeta(1-s) */
          mpfr_const_pi (p, GMP_RNDD);
          mpfr_mul (y, s, p, GMP_RNDN);
          mpfr_div_2ui (y, y, 1, GMP_RNDN);      /* s*Pi/2 */
          mpfr_sin (y, y, GMP_RNDN);             /* sin(Pi*s/2) */
          mpfr_mul (z_pre, z_pre, y, GMP_RNDN);
          mpfr_mul_2ui (y, p, 1, GMP_RNDN);      /* 2*Pi */
          mpfr_neg (s1, s1, GMP_RNDN);           /* s-1 */
          mpfr_pow (y, y, s1, GMP_RNDN);         /* (2*Pi)^(s-1) */
          mpfr_mul (z_pre, z_pre, y, GMP_RNDN);
          mpfr_mul_2ui (z_pre, z_pre, 1, GMP_RNDN);

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

          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);
}