computes Riemann function Zeta(s)

first version: only works for s=2 with PREC(result)=53 and rnd_mode=nearest


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

1 files changed, 92 insertions, 0 deletions

diff --git a/zeta.c b/zeta.c
new file mode 100644
index 000000000..6b67cd2d4
--- /dev/null
+++ b/zeta.c
@@ -0,0 +1,92 @@
+#include <stdio.h>
+#include <math.h>
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+#include "mpfr.h"
+
+int
+#if __STDC__
+mpfr_zeta(mpfr_ptr result, mpfr_srcptr op, unsigned char rnd_mode) 
+#else
+mpfr_zeta(result, op, rnd_mode)
+     mpfr_ptr result; 
+     mpfr_srcptr op; 
+     unsigned char rnd_mode; 
+#endif
+{
+	mpfr_t s,s2,x,y,u,b,v,nn,z,z2; 
+	int i,n,succes; 
+
+	/* first version */
+	if (mpfr_get_d(op) != 2.0 || rnd_mode != GMP_RNDN 
+	|| PREC(result) != 53) {
+	   fprintf(stderr, "not yet implemented\n"); exit(1);
+	}
+
+mpfr_set_default_prec(67);
+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);
+mpfr_set_ui(u,1,GMP_RNDN);
+mpfr_set_ui(s,0,GMP_RNDN);
+/*s=Somme des 1/i^2 (i=100...2)*/
+n=100;
+for (i=n; i>1; i--)
+{
+mpfr_set_ui(x,i*i,GMP_RNDN);
+mpfr_div(y,u,x,GMP_RNDN);
+mpfr_add(s,s,y,GMP_RNDN);
+};
+/*mpfr_print_raw(s);printf("\n");
+t=mpfr_out_str(stdout,10,0,s,GMP_RNDN);printf("\n");*/
+/*Application d'Euler-Maclaurin, jusqu'au terme 1/n^7 - n=100)*/
+mpfr_set_ui(nn,n,GMP_RNDN);
+mpfr_div(z,u,nn,GMP_RNDN);
+mpfr_set(s2,z,GMP_RNDN);
+mpfr_mul(z2,z,z,GMP_RNDN);
+mpfr_div_2exp(v,z2,1,GMP_RNDN);
+mpfr_sub(s2,s2,v,GMP_RNDN);
+mpfr_set_ui(b,6,GMP_RNDN);
+mpfr_mul(z,z,z2,GMP_RNDN);
+mpfr_div(v,z,b,GMP_RNDN);
+mpfr_add(s2,s2,v,GMP_RNDN);
+mpfr_set_si(b,-30,GMP_RNDN);
+mpfr_mul(z,z,z2,GMP_RNDN);
+mpfr_div(v,z,b,GMP_RNDN);
+mpfr_add(s2,s2,v,GMP_RNDN);
+mpfr_set_si(b,42,GMP_RNDN);
+mpfr_mul(z,z,z2,GMP_RNDN);
+mpfr_div(v,z,b,GMP_RNDN);
+mpfr_add(s2,s2,v,GMP_RNDN);
+/*mpfr_print_raw(s2);printf("\n");
+t=mpfr_out_str(stdout,10,0,s2,GMP_RNDN);printf("\n");*/
+mpfr_add(s,s,s2,GMP_RNDN);
+/*mpfr_print_raw(s);printf("\n");
+t=mpfr_out_str(stdout,10,0,s,GMP_RNDN);printf("\n");*/
+mpfr_add(s,s,u,GMP_RNDN);
+/*mpfr_print_raw(s);printf("\n");
+t=mpfr_out_str(stdout,10,0,s,GMP_RNDN);printf("\n");*/
+/*Peut-on arrondir ? La reponse est oui*/
+succes=mpfr_can_round(s,57,GMP_RNDN,GMP_RNDN,53);
+mpfr_set(result,s,GMP_RNDN);
+
+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 */
+}