added Riemann Zeta function (contribution from Jean-Luc Re'my)

git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@2289 280ebfd0-de03-0410-8827-d642c229c3f4

6 files changed, 593 insertions, 124 deletions

diff --git a/Makefile.am b/Makefile.am
index 462e0fe34..977d7fd9d 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -7,7 +7,7 @@ include_HEADERS = mpfr.h mpf2mpfr.h
 
 lib_LIBRARIES = libmpfr.a
 
-libmpfr_a_SOURCES = mpfr.h mpf2mpfr.h mpfr-impl.h mpfr-math.h mpfr-test.h log_b2.h exceptions.c save_expo.c extract.c uceil_exp2.c uceil_log2.c ufloor_log2.c add.c add1.c add_one_ulp.c add_ui.c agm.c clear.c cmp.c cmp_abs.c cmp_si.c cmp_ui.c comparisons.c div_2exp.c div_2si.c div_2ui.c div.c div_ui.c dump.c eq.c exp2.c exp3.c exp.c frac.c get_d.c get_exp.c get_str.c init.c inp_str.c isinteger.c isinf.c isnan.c isnum.c const_log2.c log.c mul_2exp.c mul_2si.c mul_2ui.c mul.c mul_ui.c neg.c next.c out_str.c const_pi.c pow.c pow_si.c pow_ui.c print_raw.c print_rnd_mode.c random2.c random.c reldiff.c round_prec.c set.c setmax.c setmin.c set_d.c set_dfl_prec.c set_exp.c set_rnd.c set_f.c set_prc_raw.c set_prec.c set_q.c set_si.c set_str.c set_str_raw.c set_ui.c set_z.c sqrt.c sqrt_ui.c sub.c sub1.c sub_one_ulp.c sub_ui.c rint.c ui_div.c ui_sub.c urandomb.c get_z_exp.c swap.c  factorial.c cosh.c sinh.c tanh.c acosh.c asinh.c atanh.c atan.c cmp2.c exp_2.c asin.c const_euler.c cos.c sin.c tan.c fma.c hypot.c log1p.c expm1.c log2.c log10.c ui_pow.c ui_pow_ui.c minmax.c dim.c copysign.c gmp_op.c init2.c acos.c sin_cos.c set_nan.c set_inf.c powerof2.c gamma.c set_ld.c get_ld.c cbrt.c volatile.c fits_s.h fits_sshort.c fits_sint.c fits_slong.c fits_u.h fits_ushort.c fits_uint.c fits_ulong.c get_si.c get_ui.c
+libmpfr_a_SOURCES = mpfr.h mpf2mpfr.h mpfr-impl.h mpfr-math.h mpfr-test.h log_b2.h exceptions.c save_expo.c extract.c uceil_exp2.c uceil_log2.c ufloor_log2.c add.c add1.c add_one_ulp.c add_ui.c agm.c clear.c cmp.c cmp_abs.c cmp_si.c cmp_ui.c comparisons.c div_2exp.c div_2si.c div_2ui.c div.c div_ui.c dump.c eq.c exp2.c exp3.c exp.c frac.c get_d.c get_exp.c get_str.c init.c inp_str.c isinteger.c isinf.c isnan.c isnum.c const_log2.c log.c mul_2exp.c mul_2si.c mul_2ui.c mul.c mul_ui.c neg.c next.c out_str.c const_pi.c pow.c pow_si.c pow_ui.c print_raw.c print_rnd_mode.c random2.c random.c reldiff.c round_prec.c set.c setmax.c setmin.c set_d.c set_dfl_prec.c set_exp.c set_rnd.c set_f.c set_prc_raw.c set_prec.c set_q.c set_si.c set_str.c set_str_raw.c set_ui.c set_z.c sqrt.c sqrt_ui.c sub.c sub1.c sub_one_ulp.c sub_ui.c rint.c ui_div.c ui_sub.c urandomb.c get_z_exp.c swap.c  factorial.c cosh.c sinh.c tanh.c acosh.c asinh.c atanh.c atan.c cmp2.c exp_2.c asin.c const_euler.c cos.c sin.c tan.c fma.c hypot.c log1p.c expm1.c log2.c log10.c ui_pow.c ui_pow_ui.c minmax.c dim.c copysign.c gmp_op.c init2.c acos.c sin_cos.c set_nan.c set_inf.c powerof2.c gamma.c set_ld.c get_ld.c cbrt.c volatile.c fits_s.h fits_sshort.c fits_sint.c fits_slong.c fits_u.h fits_ushort.c fits_uint.c fits_ulong.c get_si.c get_ui.c zeta.c
 
 libmpfr_a_LIBADD = @LIBOBJS@
 
diff --git a/mpfr.h b/mpfr.h
index 4d328706e..9a8e23b2f 100644
--- a/mpfr.h
+++ b/mpfr.h
@@ -283,6 +283,7 @@ int mpfr_log1p _PROTO ((mpfr_ptr, mpfr_srcptr, mp_rnd_t));
 int mpfr_expm1 _PROTO ((mpfr_ptr, mpfr_srcptr, mp_rnd_t)); 
 int mpfr_cbrt _PROTO ((mpfr_ptr, mpfr_srcptr, mp_rnd_t)); 
 int mpfr_gamma _PROTO ((mpfr_ptr, mpfr_srcptr, mp_rnd_t));
+void mpfr_zeta _PROTO ((mpfr_ptr, mpfr_srcptr, mp_rnd_t));
 
 int mpfr_min _PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t));
 int mpfr_max _PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t));
diff --git a/mpfr.texi b/mpfr.texi
index 37889247f..6e490d1f6 100644
--- a/mpfr.texi
+++ b/mpfr.texi
@@ -1279,6 +1279,11 @@ rounded to the direction @var{rnd} with the precision of @var{rop}.
 Return 0 iff the result is exact.
 @end deftypefun
 
+@deftypefun void mpfr_zeta (mpfr_t @var{rop}, mpfr_t @var{op}, mp_rnd_t @var{rnd})
+Set @var{rop} to the value of the Riemann Zeta function on @var{op},
+rounded to the direction @var{rnd} with the precision of @var{rop}.
+@end deftypefun
+
 @deftypefun int mpfr_fma (mpfr_t @var{rop}, mpfr_t @var{opx},mpfr_t @var{opy},mpfr_t @var{opz}, mp_rnd_t @var{rnd}) 
 Set @var{rop} to @math{@var{opx} @GMPtimes{} @var{opy} + @var{opz}}, rounded to the direction 
 @var{rnd} with the precision of @var{rop}. Return 0 iff the result is exact,
@@ -1579,6 +1584,7 @@ David Daney contributed the hyperbolic and inverse hyperbolic functions,
 the base-2 exponential, and the factorial function. Fabrice Rouillier
 contributed the original version of @file{mul_ui.c}, the @file{gmp_op.c}
 file, and helped to the Windows porting.
+Jean-Luc Rémy contributed the @code{mpfr_zeta} code.
 
 @node References, GNU Free Documentation License, Contributors, Top
 @comment  node-name,  next,  previous,  up
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 3d86093e4..1fa17ecc4 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -1,6 +1,6 @@
 AUTOMAKE_OPTIONS = gnu
 
-check_PROGRAMS = reuse tabs tadd tagm tcan_round tcmp tcmp2 tcmp_ui tdiv tdiv_ui tdump teq texp tfrac tget_str tlog tconst_log2 tmul tmul_2exp tmul_ui tout_str tconst_pi tpow trandom tround_prec tset_d tset_f tset_q tset_si tset_str tset_z tsqrt tsqrt_ui tui_div tui_sub tswap ttrunc trint tisnan tget_d tatan tcosh tsinh ttanh tacosh tasinh tatanh thyperbolic texp2 tfactorial tsub tasin tconst_euler tcos tsin ttan tsub_ui tset tlog1p texpm1 tlog2 tlog10 tui_pow tpow3 tadd_ui texceptions tfma thypot tacos tgamma tset_ld tcbrt tsin_cos mpf_compat mpfr_compat
+check_PROGRAMS = reuse tabs tadd tagm tcan_round tcmp tcmp2 tcmp_ui tdiv tdiv_ui tdump teq texp tfrac tget_str tlog tconst_log2 tmul tmul_2exp tmul_ui tout_str tconst_pi tpow trandom tround_prec tset_d tset_f tset_q tset_si tset_str tset_z tsqrt tsqrt_ui tui_div tui_sub tswap ttrunc trint tisnan tget_d tatan tcosh tsinh ttanh tacosh tasinh tatanh thyperbolic texp2 tfactorial tsub tasin tconst_euler tcos tsin ttan tsub_ui tset tlog1p texpm1 tlog2 tlog10 tui_pow tpow3 tadd_ui texceptions tfma thypot tacos tgamma tset_ld tcbrt tsin_cos mpf_compat mpfr_compat tzeta
 EXTRA_DIST = tgeneric.c mpf_compat.h
 
 LDADD = libfrtests.a $(MPFR_LIBM) $(top_builddir)/libmpfr.a  @LDADD@
diff --git a/tests/tzeta.c b/tests/tzeta.c
index ece65a197..54679198a 100644
--- a/tests/tzeta.c
+++ b/tests/tzeta.c
@@ -1,6 +1,7 @@
-/* Test file for mpfr_zeta.
+/* tzeta -- test file for the Riemann Zeta function
 
-Copyright 1999, 2000, 2001, 2002, 2003 Free Software Foundation.
+Copyright 2003 Free Software Foundation.
+Contributed by Jean-Luc Re'my and the Spaces project, INRIA Lorraine.
 
 This file is part of the MPFR Library.
 
@@ -20,59 +21,213 @@ the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 MA 02111-1307, USA. */
 
 #include <stdio.h>
-#include <math.h>
+#include <stdlib.h>
 #include "gmp.h"
-#include "gmp-impl.h"
-#include "longlong.h"
 #include "mpfr.h"
 
-/* #define DEBUG */
+void test_generic _PROTO ((void));
 
+void
+test_generic ()
+{
+  mp_prec_t prec, yprec;
+  int n, err;
+  mp_rnd_t rnd;
+  mpfr_t x, y, z;
+  mp_exp_t e;
+
+  mpfr_init2 (x, MPFR_PREC_MIN);
+  mpfr_init2 (y, MPFR_PREC_MIN);
+  mpfr_init2 (z, MPFR_PREC_MIN);
+
+  for (prec = 2; prec <= 100; prec++)
+    {
+      mpfr_set_prec (x, prec);
+      mpfr_set_prec (z, prec);
+      yprec = prec + 10;
+
+      for (n = 0; n < 10; n++)
+	{
+	  mpfr_random (x); /* x is in [0, 1[ */
+	  mpfr_add_ui (x, x, 1, GMP_RNDN);
+	  e = random () % 5;
+	  mpfr_div_2exp (x, x, 1, GMP_RNDN); /* now in [1/2, 1[ */
+	  mpfr_mul_2exp (x, x, e, GMP_RNDN); /* now in [2^(e-1), 2^e[ */
+          if (random () % 2)
+            mpfr_ui_sub (x, 1, x, GMP_RNDN); /* now less or equal to 1/2 */
+	  rnd = random () % 4;
+	  mpfr_set_prec (y, yprec);
+	  mpfr_zeta (y, x, rnd);
+	  err = (rnd == GMP_RNDN) ? yprec + 1 : yprec;
+	  if (mpfr_can_round (y, err, rnd, rnd, prec))
+	    {
+	      mpfr_round_prec (y, rnd, prec);
+	      mpfr_zeta (z, x, rnd);
+	      if (mpfr_cmp (y, z))
+		{
+		  printf ("results differ for x=");
+		  mpfr_out_str (stdout, 2, prec, x, GMP_RNDN);
+		  printf (" prec=%lu rnd_mode=%s\n", prec,
+			  mpfr_print_rnd_mode (rnd));
+		  printf ("   got ");
+		  mpfr_out_str (stdout, 2, prec, z, GMP_RNDN);
+		  putchar ('\n');
+		  printf ("   expected ");
+		  mpfr_out_str (stdout, 2, prec, y, GMP_RNDN);
+		  putchar ('\n');
+		  exit (1);
+		}
+	    }
+	}
+    }
+
+  mpfr_clear (x);
+  mpfr_clear (y);
+  mpfr_clear (z);
+}
+
+/* Usage: tzeta - generic tests
+          tzeta s prec rnd_mode - compute zeta(s) with precision 'prec'
+                                  and rounding mode 'mode' */
 int
-main(void)
+main (int argc, char *argv[])
 {
-  mpfr_t p,result,res_p;
-#ifdef DEBUG
-  size_t t;
-#endif
-
-  tests_start_mpfr ();
-
-  mpfr_init2(result, 53);
-  mpfr_init2(res_p, 53);
-  mpfr_init2(p, 53); mpfr_set_d(p, 2.0, GMP_RNDN);
-
-  mpfr_zeta(result,p,GMP_RNDN);
-#ifdef DEBUG
-  printf("Valeur de zeta(2) avec prec=53 et arrondi au plus pres:\n");
-  mpfr_print_binary(result);printf("\n");
-  t=mpfr_out_str(stdout,10,0,result,GMP_RNDN);printf("\n");
-#endif
-  if (mpfr_get_d1 (result) != 1.64493406684822640607e+00) {
-    fprintf(stderr, "mpfr_zeta fails for s=2 and rnd_mode=GMP_RNDN\n");
-    exit(1);
-  }
-
-  /*Test de comparaison avec pi^2/6*/
-  mpfr_set_prec(p, 67);
-  mpfr_const_pi(p,GMP_RNDN);
-  mpfr_mul(p,p,p,GMP_RNDN);
-  mpfr_set_ui(res_p,6,GMP_RNDN);
-  mpfr_div(p,p,res_p,GMP_RNDN);
-  /*mpfr_print_binary(p);printf("\n");
-    t=mpfr_out_str(stdout,10,0,p,GMP_RNDN);printf("\n");*/
-  mpfr_can_round(p,63,GMP_RNDN,GMP_RNDN,53);
-  mpfr_set(res_p,p,GMP_RNDN);
-#ifdef DEBUG
-  printf("Valeur de pi^2/6 avec prec=53 et arrondi au plus pres:\n");
-  mpfr_print_binary(res_p);printf("\n");
-  t=mpfr_out_str(stdout,10,0,res_p,GMP_RNDN);printf("\n");
-#endif
-
-  mpfr_clear(p);
-  mpfr_clear(result);
-  mpfr_clear(res_p);
-
-  tests_end_mpfr ();
+  mpfr_t s, y, z;
+  mp_prec_t prec;
+  mp_rnd_t rnd_mode;
+
+  if (argc != 1 && argc != 4)
+    {
+      fprintf (stderr, "Usage: tzeta\n");
+      fprintf (stderr, "    or tzeta s prec rnd_mode\n");
+      exit (1);
+    }
+
+  if (argc == 4)
+    {
+      prec = atoi(argv[2]);
+      mpfr_init2 (s, prec);
+      mpfr_init2 (z, prec);
+      mpfr_set_str (s, argv[1], 10, GMP_RNDN);
+      rnd_mode = atoi(argv[3]);
+
+      mpfr_zeta (z, s, rnd_mode);
+      mpfr_out_str (stdout, 10, 0, z, GMP_RNDN);
+      printf ("\n");
+
+      mpfr_clear (s);
+      mpfr_clear (z);
+
+      return 0;
+    }
+
+  mpfr_init2 (s, MPFR_PREC_MIN);
+  mpfr_init2 (y, MPFR_PREC_MIN);
+  mpfr_init2 (z, MPFR_PREC_MIN);
+
+
+  /* the following seems to loop */
+  mpfr_set_prec (s, 6);
+  mpfr_set_prec (z, 6);
+  mpfr_set_str_raw (s, "1.10010e4");
+  mpfr_zeta (z, s, GMP_RNDZ);
+
+
+  mpfr_set_prec (s, 53);
+  mpfr_set_prec (y, 53);
+  mpfr_set_prec (z, 53);
+
+  mpfr_set_str_raw (s, "0.1100011101110111111111111010000110010111001011001011");
+  mpfr_set_str_raw (y, "-0.11111101111011001001001111111000101010000100000100100E2");
+  mpfr_zeta (z, s, GMP_RNDN);
+  if (mpfr_cmp (z, y) != 0)
+    {
+      fprintf (stderr, "Error in mpfr_zeta (1,RNDN)\n");
+      exit (1);
+    }
+  mpfr_zeta (z, s, GMP_RNDZ);
+  if (mpfr_cmp (z, y) != 0)
+    {
+      fprintf (stderr, "Error in mpfr_zeta (1,RNDZ)\n");
+      exit (1);
+    }
+  mpfr_zeta (z, s, GMP_RNDU);
+  if (mpfr_cmp (z, y) != 0)
+    {
+      fprintf (stderr, "Error in mpfr_zeta (1,RNDU)\n");
+      exit (1);
+    }
+  mpfr_zeta (z, s, GMP_RNDD);
+  mpfr_add_one_ulp (y, GMP_RNDD);
+  if (mpfr_cmp (z, y) != 0)
+    {
+      fprintf (stderr, "Error in mpfr_zeta (1,RNDD)\n");
+      exit (1);
+    }
+
+  mpfr_set_str_raw (s, "0.10001011010011100110010001100100001011000010011001011");
+  mpfr_set_str_raw (y, "-0.11010011010010101101110111011010011101111101111010110E1");
+  mpfr_zeta (z, s, GMP_RNDN);
+  if (mpfr_cmp (z, y) != 0)
+    {
+      fprintf (stderr, "Error in mpfr_zeta (2,RNDN)\n");
+      exit (1);
+    }
+  mpfr_zeta (z, s, GMP_RNDZ);
+  if (mpfr_cmp (z, y) != 0)
+    {
+      fprintf (stderr, "Error in mpfr_zeta (2,RNDZ)\n");
+      exit (1);
+    }
+  mpfr_zeta (z, s, GMP_RNDU);
+  if (mpfr_cmp (z, y) != 0)
+    {
+      fprintf (stderr, "Error in mpfr_zeta (2,RNDU)\n");
+      exit (1);
+    }
+  mpfr_zeta (z, s, GMP_RNDD);
+  mpfr_add_one_ulp (y, GMP_RNDD);
+  if (mpfr_cmp (z, y) != 0)
+    {
+      fprintf (stderr, "Error in mpfr_zeta (2,RNDD)\n");
+      exit (1);
+    }
+
+  mpfr_set_str_raw (s, "0.1100111110100001111110111000110101111001011101000101");
+  mpfr_set_str_raw (y, "-0.10010111010110000111011111001101100001111011000001010E3");
+  mpfr_zeta (z, s, GMP_RNDN);
+  if (mpfr_cmp (z, y) != 0)
+    {
+      fprintf (stderr, "Error in mpfr_zeta (3,RNDN)\n");
+      exit (1);
+    }
+  mpfr_zeta (z, s, GMP_RNDD);
+  if (mpfr_cmp (z, y) != 0)
+    {
+      fprintf (stderr, "Error in mpfr_zeta (3,RNDD)\n");
+      exit (1);
+    }
+  mpfr_sub_one_ulp (y, GMP_RNDZ);
+  mpfr_zeta (z, s, GMP_RNDZ);
+  if (mpfr_cmp (z, y) != 0)
+    {
+      fprintf (stderr, "Error in mpfr_zeta (3,RNDZ)\n");
+      exit (1);
+    }
+  mpfr_zeta (z, s, GMP_RNDU);
+  if (mpfr_cmp (z, y) != 0)
+    {
+      fprintf (stderr, "Error in mpfr_zeta (3,RNDU)\n");
+      exit (1);
+    }
+
+  test_generic ();
+
+  mpfr_clear (s);
+  mpfr_clear (y);
+  mpfr_clear (z);
+
   return 0;
 }
+
+
diff --git a/zeta.c b/zeta.c
index 1020afaf0..2cb7d7395 100644
--- a/zeta.c
+++ b/zeta.c
@@ -1,6 +1,7 @@
-/* mpfr_zeta -- Riemann Zeta function at a floating-point number
+/* mpfr_zeta -- compute the Riemann Zeta function
 
-Copyright 1999, 2000, 2001, 2002, 2003 Free Software Foundation.
+Copyright 2003 Free Software Foundation.
+Contributed by Jean-Luc Re'my and the Spaces project, INRIA Lorraine.
 
 This file is part of the MPFR Library.
 
@@ -19,103 +20,409 @@ along with the MPFR Library; see the file COPYING.LIB.  If not, write to
 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 MA 02111-1307, USA. */
 
+/* #define DEBUG */
+
 #include <stdio.h>
+#include <stdlib.h>
 #include <math.h>
 #include "gmp.h"
-#include "mpfr.h"
 #include "gmp-impl.h"
-#include "longlong.h"
+#include "mpfr.h"
+#include "mpfr-impl.h"
+
+void mpfr_zeta_part_b _PROTO ((mpfr_t, mpfr_srcptr, int, int, mpfr_t *));
+void mpfr_zeta_c      _PROTO ((int, mpfr_t *));
+void mpfr_zeta_part_a _PROTO ((mpfr_t, mpfr_srcptr, int));
+void mpfr_zeta_pos    _PROTO ((mpfr_t, mpfr_srcptr, mp_rnd_t));
 
-int
-mpfr_zeta (mpfr_ptr result, mpfr_srcptr op, mp_rnd_t rnd_mode)
+/* 
+   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
+*/
+void
+mpfr_zeta_part_b (mpfr_t b, mpfr_srcptr s, int n, int p, mpfr_t *tc)
 {
-  mpfr_t s,s2,x,y,u,b,v,nn,z,z2; 
-  int i, n, succes;
-  int cmp1;
+  int n2, l, t, precb;
+  mpfr_t s1, d, u;
 
-  if (MPFR_IS_NAN(op) || MPFR_SIGN(op) < 0)
+  if (p == 0)
     {
-      MPFR_SET_NAN(result);
-      MPFR_RET_NAN;
+      mpfr_set_ui (b, 0, GMP_RNDN);
+      return;
     }
 
-  if (MPFR_IS_INF(op))  /* +infinity */
-    return mpfr_set_ui(result, 1, rnd_mode);
-
-  cmp1 = mpfr_cmp_ui(op, 1);
-  if (cmp1 < 0)
-    {
-      MPFR_SET_NAN(result);
-      MPFR_RET_NAN;
-    }
-  if (cmp1 == 0)
+  n2 = n * n;
+  precb = mpfr_get_prec (b);
+  mpfr_init2 (s1, precb);
+  mpfr_init2 (d, precb);
+  mpfr_init2 (u, precb);
+  /* 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_CLEAR_NAN(result);
-      MPFR_SET_INF(result);
-      MPFR_SET_POS(result);
-      return 0;
+      mpfr_add_ui (s1, s, t, GMP_RNDN);
+      mpfr_mul (d, d, s1, GMP_RNDN);
+      t = t - 1;
+      mpfr_add_ui (s1, s, t, GMP_RNDN);
+      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);
+      if (mpfr_cmpabs (d, tc[p-l]) > 0)
+	mpfr_set (d, tc[p-l], GMP_RNDN); 
     }
+  mpfr_mul (d, d, s, GMP_RNDN);
+  mpfr_add_ui (s1, s, 1, GMP_RNDN);
+  mpfr_neg (s1, s1, GMP_RNDN);
+  mpfr_ui_pow (u, n, s1, GMP_RNDN);
+  mpfr_mul (d, d, u, GMP_RNDN);
+  mpfr_set (b, d, GMP_RNDN);
+#ifdef DEBUG
+  printf ("\npart b=");
+  mpfr_out_str (stdout, 10, 0, b, GMP_RNDN);
+  printf ("\n");
+#endif
+  mpfr_clear (s1);
+  mpfr_clear (d);
+  mpfr_clear (u);
+}
 
-  /* 1 < op < +infinity */
+/* Input: p - an integer
+   Output: fills tc[1..p]
+*/
+void
+mpfr_zeta_c (int p, mpfr_t *tc)
+{
+  mpfr_t d;
+  int k, l;
 
-  /* first version */
-  if (mpfr_get_d1 (op) != 2.0 || rnd_mode != GMP_RNDN
-      || MPFR_PREC(result) != 53) {
-    fprintf(stderr, "not yet implemented\n"); exit(1);
-  }
+  if (p > 0)
+    {
+      mpfr_init2 (d, mpfr_get_prec (tc[1]));
+      mpfr_set_ui (tc[1], 1, GMP_RNDN);
+      mpfr_div_ui (tc[1], tc[1], 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-1; 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_neg (tc[k], tc[k], GMP_RNDN);
+	}
+      mpfr_clear(d);
+    }
+}
 
-  mpfr_set_default_prec(67);
+/* Input: s - a floating-point number
+          n - an integer
+   Output: sum - a floating-point number */
+void
+mpfr_zeta_part_a (mpfr_t sum, mpfr_srcptr s, int n)
+{
+  int i, preca;
+  mpfr_t u, s1;
 
-  mpfr_init(x); mpfr_init(y); mpfr_init(s); mpfr_init(s2);
-  mpfr_init(u); mpfr_init(b); mpfr_init(v); mpfr_init(nn);
-  mpfr_init(z); mpfr_init(z2);
+  preca = mpfr_get_prec (sum);
+  mpfr_init2 (u, preca);
+  mpfr_init2 (s1, preca);
+  mpfr_neg (s1, s, GMP_RNDN);
+  mpfr_ui_pow (u, n, s1, GMP_RNDN);
+  mpfr_div_2exp (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_ui (sum, sum, 1, GMP_RNDN);
+#ifdef DEBUG
+  printf ("\npart a = ");
+  mpfr_out_str (stdout, 10, 0, sum, GMP_RNDN);
+  printf ("\n");
+#endif
+  mpfr_clear(s1);
+  mpfr_clear(u);
+}
 
-  mpfr_set_ui(u,1,GMP_RNDN);
-  mpfr_set_ui(s,0,GMP_RNDN);
+/* Input: s - a floating-point number >= 1/2.
+          rnd_mode - a rounding mode
+   Output: z - Zeta(s) rounded to the precision of z with direction rnd_mode
+*/
+void
+mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode)
+{
+  int p, n, l, dint, add, d, can_round;
+  double beta, sd, dnep;
+  mpfr_t a,b,c,z_pre,f,g,s1;
+  mpfr_t *tc1;
+  mp_prec_t precz, precs;
 
-  /* s=Somme des 1/i^2 (i=100...2) */
-  n=100;
-  for (i=n; i>1; i--)
+  if (mpfr_nan_p (s))
+    {
+      mpfr_set_nan (z); /* Zeta(NaN) = NaN */
+      return;
+    }
+  
+  if (mpfr_inf_p (s))
     {
-      mpfr_div_ui(y,u,i*i,GMP_RNDN);
-      mpfr_add(s,s,y,GMP_RNDN);
+      if (MPFR_SIGN(s) > 0)
+	{
+	  mpfr_set_ui (z, 1, GMP_RNDN); /* Zeta(+Inf) = 1 */
+	  return;
+	}
+      else
+	{
+	  mpfr_set_nan (z); /* Zeta(-Inf) = NaN */
+	  return;
+	}
     }
 
-  /* Application d'Euler-Maclaurin, jusqu'au terme 1/n^7 - n=100) */
-  /*   mpfr_set_ui(nn,n,GMP_RNDN); */
-  mpfr_div_ui(z,u,n,GMP_RNDN);
-  mpfr_mul(z2,z,z,GMP_RNDN);
-  mpfr_div_2ui(v,z2,1,GMP_RNDN);
+  precz = mpfr_get_prec (z);
+  precs = mpfr_get_prec (s);
 
-  mpfr_set(s2,z,GMP_RNDN);
-  mpfr_sub(s2,s2,v,GMP_RNDN);
+  mpfr_init2 (a, MPFR_PREC_MIN);
+  mpfr_init2 (b, MPFR_PREC_MIN);
+  mpfr_init2 (c, MPFR_PREC_MIN);
+  mpfr_init2 (z_pre, MPFR_PREC_MIN);
+  mpfr_init2 (f, MPFR_PREC_MIN);
+  mpfr_init2 (g, MPFR_PREC_MIN);
+#ifdef DEBUG
+  printf ("Warning: mpfr_zeta assumes that s and Zeta(s) are not both representable in mpfr\n");
+  printf ("s=");
+  mpfr_print_binary (s);
+  printf ("\n");
+#endif
+  d = precz + 11;
+  mpfr_init2 (s1, precs);
+  do
+    {
+      /* Principal loop: we compute, in z_pre,
+	 an approximation of Zeta(s), that we send to mpfr_can_round */
+      mpfr_sub_ui (s1, s, 1, GMP_RNDN);
+      if (MPFR_EXP(s1) <= -(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);
+#ifdef DEBUG
+	  printf ("branch 1\n");
+	  printf ("internal precision=%d\n", dint);
+#endif
+	  mpfr_set_prec (z_pre, dint);
+	  mpfr_set_prec (g, dint);
+	  mpfr_ui_div (z_pre, 1, s1, GMP_RNDN);
+	  mpfr_const_euler (g, GMP_RNDN);
+	  mpfr_add (z_pre, z_pre, g, GMP_RNDN);
+	}
+    else /* Branch 2 */
+      {
+#ifdef DEBUG	
+	printf ("branch 2\n");
+#endif
+	/* Computation of parameters n, p and working precision */
+	dnep = (double) d * log (2.0);
+	sd = mpfr_get_d (s, GMP_RNDN);
+	beta = dnep + 0.61 + sd * log (6.2832 / sd);
+	if (beta <= 0.0)
+	  {
+	    p = 0;
+	    n = 1 + (int) (exp ((dnep - log(2.0)) / sd));
+	  }
+	else
+	  {
+	    p = 1 + (int) beta / 2;
+	    n = 1 + (int) floor((sd + 2.0 * (double) p - 1.0) / 6.2832);
+	  }
+#ifdef DEBUG
+	printf ("\nn=%d\np=%d\n",n,p);
+#endif
+	add = 4 + (int) floor (1.5 * log ((double) d) / log(2.0));
+	if (add < 10)
+	  add = 10;
+	dint = d + add;
+	if (dint < precs)
+	  dint = precs;
+#ifdef DEBUG
+	printf("internal precision=%d\n",dint);
+#endif
+	tc1 = (mpfr_t*) malloc ((p + 1) * sizeof(mpfr_t));
+	for (l=1; l<=p; l++)
+	  mpfr_init2 (tc1[l], dint);
+	mpfr_set_prec (a, dint);
+	mpfr_set_prec (b, dint);
+	mpfr_set_prec (c, dint);
+	mpfr_set_prec (z_pre, dint);
+	mpfr_set_prec (f, dint);
+#ifdef DEBUG
+	printf ("precision of z =%d\n", precz);
+#endif
+	/* Computation of the coefficients c_k */
+	mpfr_zeta_c (p, tc1);
+	/* Computation of the 3 parts of the fonction Zeta. */
+	mpfr_zeta_part_a (a, s, n);
+	mpfr_zeta_part_b (b, s, n, p, tc1);
+	mpfr_sub_ui (s1, s, 1, GMP_RNDN);
+	mpfr_ui_div (c, 1, s1, GMP_RNDN);
+	mpfr_ui_pow (f, n, s1, GMP_RNDN);
+	mpfr_div (c, c, f, GMP_RNDN);
+#ifdef DEBUG
+	printf ("\npart c=");
+	mpfr_out_str (stdout, 10, 0, c, GMP_RNDN);
+	printf ("\n");
+#endif
+	mpfr_add (z_pre, a, c, GMP_RNDN);
+	mpfr_add (z_pre, z_pre, b, GMP_RNDN);
+	for (l=1; l<=p; l++)
+	  mpfr_clear (tc1[l]);
+	free(tc1);
+	/* End branch 2 */
+      }
+#ifdef DEBUG
+      printf ("\nZeta(s) before rounding=");
+      mpfr_print_binary (z_pre);
+#endif
+      can_round = mpfr_can_round (z_pre, d - 3, GMP_RNDN, rnd_mode, precz);
+      if (can_round)
+	{
+#ifdef DEBUG
+	  printf ("\nwe can round");
+#endif
+	}
+      else
+	{
+#ifdef DEBUG	  
+	  printf ("\nwe cannot round");
+#endif
+	  d = d + 5;
+	}
+    }
+  while (can_round == 0);
 
-  mpfr_mul(z,z,z2,GMP_RNDN);
-  mpfr_div_ui(v,z,6,GMP_RNDN);
-  mpfr_add(s2,s2,v,GMP_RNDN);
+  mpfr_set (z, z_pre, rnd_mode);
+#ifdef DEBUG
+  printf ("\nZeta(s) after rounding=");
+  mpfr_print_binary (z);
+  printf ("\n");
+#endif
+  mpfr_clear (a);
+  mpfr_clear (b);
+  mpfr_clear (c);
+  mpfr_clear (z_pre);
+  mpfr_clear (f);
+  mpfr_clear (g);
+  mpfr_clear (s1);
+}
 
-  mpfr_mul(z,z,z2,GMP_RNDN);
-  mpfr_div_ui(v,z,30,GMP_RNDN);
-  mpfr_sub(s2,s2,v,GMP_RNDN);
+void
+mpfr_zeta (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode)
+{
+  double sd, eps, m1, c;
+  int can_round;
+  long add;
+  mpfr_t z_pre, s1, s2, y, sfrac, sint, p;
+  mp_prec_t precz, prec1, precs, precs1;
 
-  mpfr_mul(z,z,z2,GMP_RNDN);
-  mpfr_div_ui(v,z,42,GMP_RNDN);
-  mpfr_add(s2,s2,v,GMP_RNDN);
+  if (mpfr_cmp_ui(s,0) == 0) /* Case s = 0 */
+    {
+      mpfr_set_ui(z, 1, rnd_mode);
+      mpfr_div_2ui(z, z, 1, rnd_mode);
+      mpfr_neg(z,z,rnd_mode);
+      return;
+    }
+  
+  mpfr_init2(s2, mpfr_get_prec(s));
+  mpfr_div_2ui(s2, s, 1, rnd_mode);
+  if ((MPFR_SIGN(s) == -1) && (mpfr_floor(s2, s2) == 0)) /* Case s = -2n */
+    {
+      mpfr_set_ui(z, 0, rnd_mode);
+      mpfr_clear(s2);
+      return;
+    };
+  mpfr_clear(s2);
+  precz = mpfr_get_prec (z);
+  precs = mpfr_get_prec (s);
+  /* Precision precs1 needed to represent 1 - s, and s + 2,
+     without any truncation */
+  precs1 = precs + 2 + MAX(0, - abs(MPFR_EXP(s)));
+  sd = mpfr_get_d (s, GMP_RNDN);
+  /* 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);
+  m1 = 1.0 + MAX(1.0 / eps,  2.0 * fabs (sd - 1.0)) * (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; /* Note that prec1 is still incremented by  10 at the first  entry of loop below */
+  prec1 = MAX(prec1, precs1);
+  if (MPFR_SIGN(s) != -1 && MPFR_EXP(s) >= 0) /* Case s >= 1/2 */
+    mpfr_zeta_pos (z, s, rnd_mode);
+  else
+    {
+      mpfr_init2(z_pre, MPFR_PREC_MIN);
+      mpfr_init2(s1, MPFR_PREC_MIN);
+      mpfr_init2(y, MPFR_PREC_MIN);
+      mpfr_init2(p, MPFR_PREC_MIN);
+      mpfr_init2(sint, MPFR_PREC_MIN);
+      mpfr_init2(sfrac, MPFR_PREC_MIN);
+      do
+	{
+	  prec1 = prec1 + 10;
+      mpfr_set_prec(z_pre, prec1);
+      mpfr_set_prec(s1, prec1);
+      mpfr_set_prec(y, prec1);
+      mpfr_set_prec(p, prec1);
+      mpfr_set_prec(sint, prec1);
+      mpfr_set_prec(sfrac, prec1);
+      mpfr_ui_sub(s1, 1, s, GMP_RNDN);
+      mpfr_add_ui(sfrac, s, 2, GMP_RNDN);
+      mpfr_div_2ui(sfrac, sfrac, 2, GMP_RNDN);
+      mpfr_floor(sint, sfrac);
+      mpfr_sub(sfrac, sfrac, sint, GMP_RNDN);
+      mpfr_mul_2ui(sfrac, sfrac, 2, GMP_RNDN);
+      mpfr_sub_ui(sfrac, sfrac, 2, GMP_RNDN);
+      if (mpfr_cmp_ui(sfrac, 1) > 0)
+        mpfr_ui_sub (sfrac, 2, sfrac, GMP_RNDN);
+      else if (mpfr_cmp_si(sfrac, -1) < 0)
+        {
+          mpfr_neg(sfrac, sfrac, GMP_RNDN);
+          mpfr_sub_ui(sfrac, sfrac, 2, GMP_RNDN);
+        }
+      mpfr_div_2ui(sfrac, sfrac, 1, GMP_RNDN);
+      mpfr_zeta_pos(z_pre, s1, GMP_RNDN);
+      mpfr_gamma(y, s1, GMP_RNDN);
+      mpfr_mul(z_pre, z_pre, y, GMP_RNDN);
+      mpfr_const_pi(p, GMP_RNDD);
+      mpfr_mul(y, sfrac, p, GMP_RNDN);
+      mpfr_sin(y, y, GMP_RNDN);
+      mpfr_mul(z_pre, z_pre, y, GMP_RNDN);
+      mpfr_mul_2ui(y, p, 1, GMP_RNDN);
+      mpfr_neg(s1, s1, GMP_RNDN);
+      mpfr_pow(y, y, s1, GMP_RNDN);
+      mpfr_mul(z_pre, z_pre, y, GMP_RNDN);
+      mpfr_mul_2ui(z_pre, z_pre, 1, GMP_RNDN);
+      can_round = mpfr_can_round (z_pre, prec1 - add, GMP_RNDN, rnd_mode, precz);}
+      while (can_round == 0);
+      mpfr_set (z, z_pre, rnd_mode);
+      mpfr_clear(z_pre);
+      mpfr_clear(s1);
+      mpfr_clear(y);
+      mpfr_clear(p);
+      mpfr_clear(sint);
+      mpfr_clear(sfrac);
+	}
+}
 
-  mpfr_add(s,s,s2,GMP_RNDN);
-  mpfr_add(s,s,u,GMP_RNDN);
 
-  /*Peut-on arrondir ? La reponse est oui*/
-  succes=mpfr_can_round(s, 57, GMP_RNDN,GMP_RNDN, 53);
 
-  if (succes) mpfr_set(result,s,GMP_RNDN);
-  else {
-    fprintf(stderr, "can't round in mpfr_zeta\n"); exit(1);
-  }
 
-  mpfr_clear(x); mpfr_clear(y); mpfr_clear(s); mpfr_clear(s2);
-  mpfr_clear(u); mpfr_clear(b); mpfr_clear(v); mpfr_clear(nn);
-  mpfr_clear(z); mpfr_clear(z2);
 
-  return 1; /* result is inexact */
-}