1/* tzeta -- test file for the Riemann Zeta function
2
3Copyright 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
4Contributed by Jean-Luc Re'my and the Spaces project, INRIA Lorraine.
5
6This file is part of the GNU MPFR Library.
7
8The GNU MPFR Library is free software; you can redistribute it and/or modify
9it under the terms of the GNU Lesser General Public License as published by
10the Free Software Foundation; either version 3 of the License, or (at your
11option) any later version.
12
13The GNU MPFR Library is distributed in the hope that it will be useful, but
14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
16License for more details.
17
18You should have received a copy of the GNU Lesser General Public License
19along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22
23#include <stdio.h>
24#include <stdlib.h>
25
26#include "mpfr-test.h"
27
28static void
29test1 (void)
30{
31  mpfr_t x, y;
32
33  mpfr_init2 (x, 32);
34  mpfr_init2 (y, 42);
35
36  mpfr_set_str_binary (x, "1.1111111101000111011010010010100e-1");
37  mpfr_zeta (y, x, MPFR_RNDN); /* shouldn't crash */
38
39  mpfr_set_prec (x, 40);
40  mpfr_set_prec (y, 50);
41  mpfr_set_str_binary (x, "1.001101001101000010011010110100110000101e-1");
42  mpfr_zeta (y, x, MPFR_RNDU);
43  mpfr_set_prec (x, 50);
44  mpfr_set_str_binary (x, "-0.11111100011100111111101111100011110111001111111111E1");
45  if (mpfr_cmp (x, y))
46    {
47      printf ("Error for input on 40 bits, output on 50 bits\n");
48      printf ("Expected "); mpfr_print_binary (x); puts ("");
49      printf ("Got      "); mpfr_print_binary (y); puts ("");
50      mpfr_set_str_binary (x, "1.001101001101000010011010110100110000101e-1");
51      mpfr_zeta (y, x, MPFR_RNDU);
52      mpfr_print_binary (x); puts ("");
53      mpfr_print_binary (y); puts ("");
54      exit (1);
55    }
56
57  mpfr_set_prec (x, 2);
58  mpfr_set_prec (y, 55);
59  mpfr_set_str_binary (x, "0.11e3");
60  mpfr_zeta (y, x, MPFR_RNDN);
61  mpfr_set_prec (x, 55);
62  mpfr_set_str_binary (x, "0.1000001000111000010011000010011000000100100100100010010E1");
63  if (mpfr_cmp (x, y))
64    {
65      printf ("Error in mpfr_zeta (1)\n");
66      printf ("Expected "); mpfr_print_binary (x); puts ("");
67      printf ("Got      "); mpfr_print_binary (y); puts ("");
68      exit (1);
69    }
70
71  mpfr_set_prec (x, 3);
72  mpfr_set_prec (y, 47);
73  mpfr_set_str_binary (x, "0.111e4");
74  mpfr_zeta (y, x, MPFR_RNDN);
75  mpfr_set_prec (x, 47);
76  mpfr_set_str_binary (x, "1.0000000000000100000000111001001010111100101011");
77  if (mpfr_cmp (x, y))
78    {
79      printf ("Error in mpfr_zeta (2)\n");
80      exit (1);
81    }
82
83  /* coverage test */
84  mpfr_set_prec (x, 7);
85  mpfr_set_str_binary (x, "1.000001");
86  mpfr_set_prec (y, 2);
87  mpfr_zeta (y, x, MPFR_RNDN);
88  MPFR_ASSERTN(mpfr_cmp_ui (y, 64) == 0);
89
90  /* another coverage test */
91  mpfr_set_prec (x, 24);
92  mpfr_set_ui (x, 2, MPFR_RNDN);
93  mpfr_set_prec (y, 2);
94  mpfr_zeta (y, x, MPFR_RNDN);
95  MPFR_ASSERTN(mpfr_cmp_ui_2exp (y, 3, -1) == 0);
96
97  mpfr_set_nan (x);
98  mpfr_zeta (y, x, MPFR_RNDN);
99  MPFR_ASSERTN(mpfr_nan_p (y));
100
101  mpfr_set_inf (x, 1);
102  mpfr_zeta (y, x, MPFR_RNDN);
103  MPFR_ASSERTN(mpfr_cmp_ui (y, 1) == 0);
104
105  mpfr_set_inf (x, -1);
106  mpfr_zeta (y, x, MPFR_RNDN);
107  MPFR_ASSERTN(mpfr_nan_p (y));
108
109  mpfr_clear (x);
110  mpfr_clear (y);
111}
112
113static const char *const val[] = {
114  "-2000", "0.0",
115  "-2.0", "0.0",
116  "-1.0", "-0.000101010101010101010101010101010101010101010101010101010101010",
117  "-0.9", "-0.000110011110011111010001010001100010111101001010100110001110110",
118  /*  "-0.8", "-0.000111110011101010001011100011010010000001010011110100010001110",
119  "-0.7", "-0.00100101011011111100110011110011111010111111000110110100010110",
120  "-0.6", "-0.00101100101100100100110111111000110010111010110010111000001100",
121  "-0.5", "-0.00110101001110000000100000011001100100010000111100010001111100",
122  "-0.4", "-0.00111111010001100011110001010010111110010001010101111101110001",
123  "-0.3", "-0.0100101100110111010101010100111011000001001010111010110101010",
124  "-0.2", "-0.0101100110000011101110101011011110101111000010000010110101111",
125  "-0.1", "-0.0110101011001111011101001111011000010001111010110011011111011",
126  "-0.0", "-0.100000000000000000000000000000000000000000000000000000000000",
127  "0.1", "-0.100110100110000010101010101110100000101100100011011001000101",
128  "0.2", "-0.10111011111000100011110111100010010001111010010010010100010110",
129  "0.3", "-0.11100111100100010011001000001011001100110010110101101110110110",
130  "0.4", "-1.0010001010000010000110111000100101001000001011101010110101011",
131  "0.5", "-1.0111010111011001110010110000011111100111001111111110111000110",
132  "0.6", "-1.1111001111100001100111101110010001001000001101100110110000100",
133  "0.7", "-10.110001110100010001110111000101010011110011000110010100101000",
134  "0.8", "-100.01110000000000101000010010000011000000111101100101100011010",
135  "0.9", "-1001.0110111000011011111100111100111011100010001111111010000100",
136  "0.99","-0.11000110110110001101011010110001011010011000110001011100101110E7",
137  "0.997", "-0.10100110011000001100111110011111100011110000111011101110001010E9",
138  "0.9995", "-0.11111001111011011000011110111111010111101001000110001111110010E11",
139  "0.99998", "-0.11000011010011110110110000111011101100001000101101011001110100E16",
140  "1.00001", "0.11000011010100000100100111100010001110100000110101110011111011E17",
141  "1.0002", "0.10011100010001001001111000101010111000011011011111110010110100E13",
142  "1.003","0.10100110111101001001010000000110101101110100001010100000110000E9",
143  "1.04", "11001.100101001000001011000111010110011010000001000010111101101",
144  "1.1", "1010.1001010110011110011010100010001100101001001111111101100001",
145  "1.2", "101.10010111011100011111001001100101101111110000110001101100010",
146  "1.3", "11.111011101001010000111001001110100100000101000101101011010100",
147  "1.4", "11.000110110000010100100101011110110001100001110100100100111111",
148  "1.5", "10.100111001100010010100001011111110111101100010011101011011100",
149  "1.6", "10.010010010010011111110000010011000110101001110011101010100110",
150  "1.7", "10.000011011110010111011110001100110010100010011100011111110010",
151  "1.8", "1.1110000111011001110011001101110101010000011011101100010111001",
152  "1.9", "1.1011111111101111011000011110001100100111100110111101101000101",
153  "2.0", "1.1010010100011010011001100010010100110000011111010011001000110",
154  "42.17", "1.0000000000000000000000000000000000000000001110001110001011001",
155  "-17.42", "-11.101110101010101000000001001000001111111101000100001100101100",
156  "-24.17", "-0.10001111010010011111000010001011111010010111101011000010010011E13"*/
157};
158
159static void
160test2 (void)
161{
162  mpfr_t x, y;
163  int i, n = numberof(val);
164
165  mpfr_inits2 (55, x, y, (mpfr_ptr) 0);
166
167  for(i = 0 ; i < n ; i+=2)
168    {
169      mpfr_set_str1 (x, val[i]);
170      mpfr_zeta(y, x, MPFR_RNDZ);
171      if (mpfr_cmp_str (y, val[i+1] , 2, MPFR_RNDZ))
172        {
173          printf("Wrong result for zeta(%s=", val[i]);
174          mpfr_print_binary (x);
175          printf (").\nGot     : ");
176          mpfr_print_binary(y); putchar('\n');
177          printf("Expected: ");
178          mpfr_set_str (y, val[i+1], 2, MPFR_RNDZ);
179          mpfr_print_binary(y); putchar('\n');
180          mpfr_set_prec(y, 65);
181          mpfr_zeta(y, x, MPFR_RNDZ);
182          printf("+ Prec  : ");
183          mpfr_print_binary(y); putchar('\n');
184          exit(1);
185        }
186    }
187  mpfr_clears (x, y, (mpfr_ptr) 0);
188}
189
190#define TEST_FUNCTION mpfr_zeta
191#define TEST_RANDOM_EMIN -48
192#define TEST_RANDOM_EMAX 31
193#include "tgeneric.c"
194
195/* Usage: tzeta - generic tests
196          tzeta s prec rnd_mode - compute zeta(s) with precision 'prec'
197                                  and rounding mode 'mode' */
198int
199main (int argc, char *argv[])
200{
201  mpfr_t s, y, z;
202  mpfr_prec_t prec;
203  mpfr_rnd_t rnd_mode;
204  int inex;
205
206  tests_start_mpfr ();
207
208  if (argc != 1 && argc != 4)
209    {
210      printf ("Usage: tzeta\n"
211              "    or tzeta s prec rnd_mode\n");
212      exit (1);
213    }
214
215  if (argc == 4)
216    {
217      prec = atoi(argv[2]);
218      mpfr_init2 (s, prec);
219      mpfr_init2 (z, prec);
220      mpfr_set_str (s, argv[1], 10, MPFR_RNDN);
221      rnd_mode = (mpfr_rnd_t) atoi(argv[3]);
222
223      mpfr_zeta (z, s, rnd_mode);
224      mpfr_out_str (stdout, 10, 0, z, MPFR_RNDN);
225      printf ("\n");
226
227      mpfr_clear (s);
228      mpfr_clear (z);
229
230      return 0;
231    }
232
233  test1();
234
235  mpfr_init2 (s, MPFR_PREC_MIN);
236  mpfr_init2 (y, MPFR_PREC_MIN);
237  mpfr_init2 (z, MPFR_PREC_MIN);
238
239
240  /* the following seems to loop */
241  mpfr_set_prec (s, 6);
242  mpfr_set_prec (z, 6);
243  mpfr_set_str_binary (s, "1.10010e4");
244  mpfr_zeta (z, s, MPFR_RNDZ);
245
246
247  mpfr_set_prec (s, 53);
248  mpfr_set_prec (y, 53);
249  mpfr_set_prec (z, 53);
250
251  mpfr_set_ui (s, 1, MPFR_RNDN);
252  mpfr_zeta (z, s, MPFR_RNDN);
253  if (!mpfr_inf_p (z) || MPFR_SIGN (z) < 0)
254    {
255      printf ("Error in mpfr_zeta for s = 1 (should be +inf)\n");
256      exit (1);
257    }
258
259  mpfr_set_str_binary (s, "0.1100011101110111111111111010000110010111001011001011");
260  mpfr_set_str_binary (y, "-0.11111101111011001001001111111000101010000100000100100E2");
261  mpfr_zeta (z, s, MPFR_RNDN);
262  if (mpfr_cmp (z, y) != 0)
263    {
264      printf ("Error in mpfr_zeta (1,RNDN)\n");
265      exit (1);
266    }
267  mpfr_zeta (z, s, MPFR_RNDZ);
268  if (mpfr_cmp (z, y) != 0)
269    {
270      printf ("Error in mpfr_zeta (1,RNDZ)\n");
271      exit (1);
272    }
273  mpfr_zeta (z, s, MPFR_RNDU);
274  if (mpfr_cmp (z, y) != 0)
275    {
276      printf ("Error in mpfr_zeta (1,RNDU)\n");
277      exit (1);
278    }
279  mpfr_zeta (z, s, MPFR_RNDD);
280  mpfr_nexttoinf (y);
281  if (mpfr_cmp (z, y) != 0)
282    {
283      printf ("Error in mpfr_zeta (1,RNDD)\n");
284      exit (1);
285    }
286
287  mpfr_set_str_binary (s, "0.10001011010011100110010001100100001011000010011001011");
288  mpfr_set_str_binary (y, "-0.11010011010010101101110111011010011101111101111010110E1");
289  mpfr_zeta (z, s, MPFR_RNDN);
290  if (mpfr_cmp (z, y) != 0)
291    {
292      printf ("Error in mpfr_zeta (2,RNDN)\n");
293      exit (1);
294    }
295  mpfr_zeta (z, s, MPFR_RNDZ);
296  if (mpfr_cmp (z, y) != 0)
297    {
298      printf ("Error in mpfr_zeta (2,RNDZ)\n");
299      exit (1);
300    }
301  mpfr_zeta (z, s, MPFR_RNDU);
302  if (mpfr_cmp (z, y) != 0)
303    {
304      printf ("Error in mpfr_zeta (2,RNDU)\n");
305      exit (1);
306    }
307  mpfr_zeta (z, s, MPFR_RNDD);
308  mpfr_nexttoinf (y);
309  if (mpfr_cmp (z, y) != 0)
310    {
311      printf ("Error in mpfr_zeta (2,RNDD)\n");
312      exit (1);
313    }
314
315  mpfr_set_str_binary (s, "0.1100111110100001111110111000110101111001011101000101");
316  mpfr_set_str_binary (y, "-0.10010111010110000111011111001101100001111011000001010E3");
317  mpfr_zeta (z, s, MPFR_RNDN);
318  if (mpfr_cmp (z, y) != 0)
319    {
320      printf ("Error in mpfr_zeta (3,RNDN)\n");
321      exit (1);
322    }
323  mpfr_zeta (z, s, MPFR_RNDD);
324  if (mpfr_cmp (z, y) != 0)
325    {
326      printf ("Error in mpfr_zeta (3,RNDD)\n");
327      exit (1);
328    }
329  mpfr_nexttozero (y);
330  mpfr_zeta (z, s, MPFR_RNDZ);
331  if (mpfr_cmp (z, y) != 0)
332    {
333      printf ("Error in mpfr_zeta (3,RNDZ)\n");
334      exit (1);
335    }
336  mpfr_zeta (z, s, MPFR_RNDU);
337  if (mpfr_cmp (z, y) != 0)
338    {
339      printf ("Error in mpfr_zeta (3,RNDU)\n");
340      exit (1);
341    }
342
343  mpfr_set_str (s, "-400000001", 10, MPFR_RNDZ);
344  mpfr_zeta (z, s, MPFR_RNDN);
345  if (!(mpfr_inf_p (z) && MPFR_SIGN(z) < 0))
346    {
347      printf ("Error in mpfr_zeta (-400000001)\n");
348      exit (1);
349    }
350  mpfr_set_str (s, "-400000003", 10, MPFR_RNDZ);
351  mpfr_zeta (z, s, MPFR_RNDN);
352  if (!(mpfr_inf_p (z) && MPFR_SIGN(z) > 0))
353    {
354      printf ("Error in mpfr_zeta (-400000003)\n");
355      exit (1);
356    }
357
358  mpfr_set_prec (s, 34);
359  mpfr_set_prec (z, 34);
360  mpfr_set_str_binary (s, "-1.111111100001011110000010001010000e-35");
361  mpfr_zeta (z, s, MPFR_RNDD);
362  mpfr_set_str_binary (s, "-1.111111111111111111111111111111111e-2");
363  if (mpfr_cmp (s, z))
364    {
365      printf ("Error in mpfr_zeta, prec=34, MPFR_RNDD\n");
366      mpfr_dump (z);
367      exit (1);
368    }
369
370  /* bug found by nightly tests on June 7, 2007 */
371  mpfr_set_prec (s, 23);
372  mpfr_set_prec (z, 25);
373  mpfr_set_str_binary (s, "-1.0110110110001000000000e-27");
374  mpfr_zeta (z, s, MPFR_RNDN);
375  mpfr_set_prec (s, 25);
376  mpfr_set_str_binary (s, "-1.111111111111111111111111e-2");
377  if (mpfr_cmp (s, z))
378    {
379      printf ("Error in mpfr_zeta, prec=25, MPFR_RNDN\n");
380      printf ("expected "); mpfr_dump (s);
381      printf ("got      "); mpfr_dump (z);
382      exit (1);
383    }
384
385  /* bug reported by Kevin Rauch on 26 Oct 2007 */
386  mpfr_set_prec (s, 128);
387  mpfr_set_prec (z, 128);
388  mpfr_set_str_binary (s, "-0.1000000000000000000000000000000000000000000000000000000000000001E64");
389  inex = mpfr_zeta (z, s, MPFR_RNDN);
390  MPFR_ASSERTN (mpfr_inf_p (z) && MPFR_SIGN (z) < 0 && inex < 0);
391  inex = mpfr_zeta (z, s, MPFR_RNDU);
392  mpfr_set_inf (s, -1);
393  mpfr_nextabove (s);
394  MPFR_ASSERTN (mpfr_equal_p (z, s) && inex > 0);
395
396  mpfr_clear (s);
397  mpfr_clear (y);
398  mpfr_clear (z);
399
400  test_generic (2, 70, 5);
401  test2 ();
402
403  tests_end_mpfr ();
404  return 0;
405}
406