1/* tzeta -- test file for the Riemann Zeta function
2
3Copyright 2003-2023 Free Software Foundation, Inc.
4Contributed by the AriC and Caramba projects, INRIA.
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
20https://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 "mpfr-test.h"
24
25static void
26test1 (void)
27{
28  mpfr_t x, y;
29  int inex;
30
31  mpfr_init2 (x, 32);
32  mpfr_init2 (y, 42);
33
34  mpfr_clear_flags ();
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_dump (x);
49      printf ("Got      "); mpfr_dump (y);
50      mpfr_set_str_binary (x, "1.001101001101000010011010110100110000101e-1");
51      mpfr_zeta (y, x, MPFR_RNDU);
52      mpfr_dump (x);
53      mpfr_dump (y);
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_dump (x);
67      printf ("Got      "); mpfr_dump (y);
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  /* yet another coverage test (case beta <= 0.0) */
98  mpfr_set_prec (x, 10);
99  mpfr_set_ui (x, 23, MPFR_RNDN);
100  mpfr_set_prec (y, 15);
101  inex = mpfr_zeta (y, x, MPFR_RNDN);
102  MPFR_ASSERTN(inex < 0);
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_cmp_ui (y, 1) == 0);
108
109  /* Since some tests don't really check that the result is not NaN... */
110  MPFR_ASSERTN (! mpfr_nanflag_p ());
111
112  mpfr_set_inf (x, -1);
113  mpfr_zeta (y, x, MPFR_RNDN);
114  MPFR_ASSERTN(mpfr_nan_p (y));
115
116  mpfr_set_nan (x);
117  mpfr_zeta (y, x, MPFR_RNDN);
118  MPFR_ASSERTN(mpfr_nan_p (y));
119
120  mpfr_clear (x);
121  mpfr_clear (y);
122}
123
124static const char *const val[] = {
125  "-2000", "0.0",
126  "-2.0", "0.0",
127  "-1.0", "-0.000101010101010101010101010101010101010101010101010101010101010",
128  "-0.9", "-0.000110011110011111010001010001100010111101001010100110001110110",
129  /*  "-0.8", "-0.000111110011101010001011100011010010000001010011110100010001110",
130  "-0.7", "-0.00100101011011111100110011110011111010111111000110110100010110",
131  "-0.6", "-0.00101100101100100100110111111000110010111010110010111000001100",
132  "-0.5", "-0.00110101001110000000100000011001100100010000111100010001111100",
133  "-0.4", "-0.00111111010001100011110001010010111110010001010101111101110001",
134  "-0.3", "-0.0100101100110111010101010100111011000001001010111010110101010",
135  "-0.2", "-0.0101100110000011101110101011011110101111000010000010110101111",
136  "-0.1", "-0.0110101011001111011101001111011000010001111010110011011111011",
137  "-0.0", "-0.100000000000000000000000000000000000000000000000000000000000",
138  "0.1", "-0.100110100110000010101010101110100000101100100011011001000101",
139  "0.2", "-0.10111011111000100011110111100010010001111010010010010100010110",
140  "0.3", "-0.11100111100100010011001000001011001100110010110101101110110110",
141  "0.4", "-1.0010001010000010000110111000100101001000001011101010110101011",
142  "0.5", "-1.0111010111011001110010110000011111100111001111111110111000110",
143  "0.6", "-1.1111001111100001100111101110010001001000001101100110110000100",
144  "0.7", "-10.110001110100010001110111000101010011110011000110010100101000",
145  "0.8", "-100.01110000000000101000010010000011000000111101100101100011010",
146  "0.9", "-1001.0110111000011011111100111100111011100010001111111010000100",
147  "0.99","-0.11000110110110001101011010110001011010011000110001011100101110E7",
148  "0.997", "-0.10100110011000001100111110011111100011110000111011101110001010E9",
149  "0.9995", "-0.11111001111011011000011110111111010111101001000110001111110010E11",
150  "0.99998", "-0.11000011010011110110110000111011101100001000101101011001110100E16",
151  "1.00001", "0.11000011010100000100100111100010001110100000110101110011111011E17",
152  "1.0002", "0.10011100010001001001111000101010111000011011011111110010110100E13",
153  "1.003","0.10100110111101001001010000000110101101110100001010100000110000E9",
154  "1.04", "11001.100101001000001011000111010110011010000001000010111101101",
155  "1.1", "1010.1001010110011110011010100010001100101001001111111101100001",
156  "1.2", "101.10010111011100011111001001100101101111110000110001101100010",
157  "1.3", "11.111011101001010000111001001110100100000101000101101011010100",
158  "1.4", "11.000110110000010100100101011110110001100001110100100100111111",
159  "1.5", "10.100111001100010010100001011111110111101100010011101011011100",
160  "1.6", "10.010010010010011111110000010011000110101001110011101010100110",
161  "1.7", "10.000011011110010111011110001100110010100010011100011111110010",
162  "1.8", "1.1110000111011001110011001101110101010000011011101100010111001",
163  "1.9", "1.1011111111101111011000011110001100100111100110111101101000101",
164  "2.0", "1.1010010100011010011001100010010100110000011111010011001000110",
165  "42.17", "1.0000000000000000000000000000000000000000001110001110001011001",
166  "-17.42", "-11.101110101010101000000001001000001111111101000100001100101100",
167  "-24.17", "-0.10001111010010011111000010001011111010010111101011000010010011E13"*/
168};
169
170static void
171test2 (void)
172{
173  mpfr_t x, y;
174  int i, n = numberof(val);
175
176  mpfr_inits2 (55, x, y, (mpfr_ptr) 0);
177
178  for(i = 0 ; i < n ; i+=2)
179    {
180      mpfr_set_str1 (x, val[i]);
181      mpfr_zeta(y, x, MPFR_RNDZ);
182      if (mpfr_cmp_str (y, val[i+1], 2, MPFR_RNDZ))
183        {
184          printf("Wrong result for zeta(%s=", val[i]);
185          mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN);
186          printf (").\nGot     : ");
187          mpfr_dump (y);
188          printf("Expected: ");
189          mpfr_set_str (y, val[i+1], 2, MPFR_RNDZ);
190          mpfr_dump (y);
191          mpfr_set_prec(y, 65);
192          mpfr_zeta(y, x, MPFR_RNDZ);
193          printf("+ Prec  : ");
194          mpfr_dump (y);
195          exit(1);
196        }
197    }
198  mpfr_clears (x, y, (mpfr_ptr) 0);
199}
200
201/* The following test attempts to trigger an intermediate overflow in
202   Gamma(s1) in the reflection formula with a 32-bit ABI (the example
203   depends on the extended exponent range): r10804 fails when the
204   exponent field is on 32 bits. */
205static void
206intermediate_overflow (void)
207{
208  mpfr_t x, y1, y2;
209  mpfr_flags_t flags1, flags2;
210  int inex1, inex2;
211
212  mpfr_inits2 (64, x, y1, y2, (mpfr_ptr) 0);
213
214  mpfr_set_si (x, -44787928, MPFR_RNDN);
215  mpfr_nextabove (x);
216
217  mpfr_set_str (y1, "0x3.0a6ab0ab281742acp+954986780", 0, MPFR_RNDN);
218  inex1 = -1;
219  flags1 = MPFR_FLAGS_INEXACT;
220
221  mpfr_clear_flags ();
222  inex2 = mpfr_zeta (y2, x, MPFR_RNDN);
223  flags2 = __gmpfr_flags;
224
225  if (!(mpfr_equal_p (y1, y2) &&
226        SAME_SIGN (inex1, inex2) &&
227        flags1 == flags2))
228    {
229      printf ("Error in intermediate_overflow\n");
230      printf ("Expected ");
231      mpfr_dump (y1);
232      printf ("with inex = %d and flags =", inex1);
233      flags_out (flags1);
234      printf ("Got      ");
235      mpfr_dump (y2);
236      printf ("with inex = %d and flags =", inex2);
237      flags_out (flags2);
238      exit (1);
239    }
240  mpfr_clears (x, y1, y2, (mpfr_ptr) 0);
241}
242
243#define TEST_FUNCTION mpfr_zeta
244#define TEST_RANDOM_EMIN -48
245#define TEST_RANDOM_EMAX 31
246#include "tgeneric.c"
247
248/* Usage: tzeta - generic tests
249          tzeta s prec rnd_mode - compute zeta(s) with precision 'prec'
250                                  and rounding mode 'mode' */
251int
252main (int argc, char *argv[])
253{
254  mpfr_t s, y, z;
255  mpfr_prec_t prec;
256  mpfr_rnd_t rnd_mode;
257  mpfr_flags_t flags;
258  int inex;
259
260  tests_start_mpfr ();
261
262  if (argc != 1 && argc != 4)
263    {
264      printf ("Usage: tzeta\n"
265              "    or tzeta s prec rnd_mode\n");
266      exit (1);
267    }
268
269  if (argc == 4)
270    {
271      prec = atoi(argv[2]);
272      mpfr_init2 (s, prec);
273      mpfr_init2 (z, prec);
274      mpfr_set_str (s, argv[1], 10, MPFR_RNDN);
275      rnd_mode = (mpfr_rnd_t) atoi(argv[3]);
276
277      mpfr_zeta (z, s, rnd_mode);
278      mpfr_out_str (stdout, 10, 0, z, MPFR_RNDN);
279      printf ("\n");
280
281      mpfr_clear (s);
282      mpfr_clear (z);
283
284      return 0;
285    }
286
287  test1();
288
289  mpfr_init2 (s, MPFR_PREC_MIN);
290  mpfr_init2 (y, MPFR_PREC_MIN);
291  mpfr_init2 (z, MPFR_PREC_MIN);
292
293
294  /* the following seems to loop */
295  mpfr_set_prec (s, 6);
296  mpfr_set_prec (z, 6);
297  mpfr_set_str_binary (s, "1.10010e4");
298  mpfr_zeta (z, s, MPFR_RNDZ);
299
300  mpfr_set_prec (s, 53);
301  mpfr_set_prec (y, 53);
302  mpfr_set_prec (z, 53);
303
304  mpfr_set_ui (s, 1, MPFR_RNDN);
305  mpfr_clear_divby0();
306  mpfr_zeta (z, s, MPFR_RNDN);
307  if (!mpfr_inf_p (z) || MPFR_IS_NEG (z) || !mpfr_divby0_p())
308    {
309      printf ("Error in mpfr_zeta for s = 1 (should be +inf) with divby0 flag\n");
310      exit (1);
311    }
312
313  mpfr_set_str_binary (s, "0.1100011101110111111111111010000110010111001011001011");
314  mpfr_set_str_binary (y, "-0.11111101111011001001001111111000101010000100000100100E2");
315  mpfr_zeta (z, s, MPFR_RNDN);
316  if (mpfr_cmp (z, y) != 0)
317    {
318      printf ("Error in mpfr_zeta (1,RNDN)\n");
319      exit (1);
320    }
321  mpfr_zeta (z, s, MPFR_RNDZ);
322  if (mpfr_cmp (z, y) != 0)
323    {
324      printf ("Error in mpfr_zeta (1,RNDZ)\n");
325      exit (1);
326    }
327  mpfr_zeta (z, s, MPFR_RNDU);
328  if (mpfr_cmp (z, y) != 0)
329    {
330      printf ("Error in mpfr_zeta (1,RNDU)\n");
331      exit (1);
332    }
333  mpfr_zeta (z, s, MPFR_RNDD);
334  mpfr_nexttoinf (y);
335  if (mpfr_cmp (z, y) != 0)
336    {
337      printf ("Error in mpfr_zeta (1,RNDD)\n");
338      exit (1);
339    }
340
341  mpfr_set_str_binary (s, "0.10001011010011100110010001100100001011000010011001011");
342  mpfr_set_str_binary (y, "-0.11010011010010101101110111011010011101111101111010110E1");
343  mpfr_zeta (z, s, MPFR_RNDN);
344  if (mpfr_cmp (z, y) != 0)
345    {
346      printf ("Error in mpfr_zeta (2,RNDN)\n");
347      exit (1);
348    }
349  mpfr_zeta (z, s, MPFR_RNDZ);
350  if (mpfr_cmp (z, y) != 0)
351    {
352      printf ("Error in mpfr_zeta (2,RNDZ)\n");
353      exit (1);
354    }
355  mpfr_zeta (z, s, MPFR_RNDU);
356  if (mpfr_cmp (z, y) != 0)
357    {
358      printf ("Error in mpfr_zeta (2,RNDU)\n");
359      exit (1);
360    }
361  mpfr_zeta (z, s, MPFR_RNDD);
362  mpfr_nexttoinf (y);
363  if (mpfr_cmp (z, y) != 0)
364    {
365      printf ("Error in mpfr_zeta (2,RNDD)\n");
366      exit (1);
367    }
368
369  mpfr_set_str_binary (s, "0.1100111110100001111110111000110101111001011101000101");
370  mpfr_set_str_binary (y, "-0.10010111010110000111011111001101100001111011000001010E3");
371  mpfr_zeta (z, s, MPFR_RNDN);
372  if (mpfr_cmp (z, y) != 0)
373    {
374      printf ("Error in mpfr_zeta (3,RNDN)\n");
375      exit (1);
376    }
377  mpfr_zeta (z, s, MPFR_RNDD);
378  if (mpfr_cmp (z, y) != 0)
379    {
380      printf ("Error in mpfr_zeta (3,RNDD)\n");
381      exit (1);
382    }
383  mpfr_nexttozero (y);
384  mpfr_zeta (z, s, MPFR_RNDZ);
385  if (mpfr_cmp (z, y) != 0)
386    {
387      printf ("Error in mpfr_zeta (3,RNDZ)\n");
388      exit (1);
389    }
390  mpfr_zeta (z, s, MPFR_RNDU);
391  if (mpfr_cmp (z, y) != 0)
392    {
393      printf ("Error in mpfr_zeta (3,RNDU)\n");
394      exit (1);
395    }
396
397  mpfr_set_str (s, "-400000001", 10, MPFR_RNDZ);
398  mpfr_zeta (z, s, MPFR_RNDN);
399  if (!(mpfr_inf_p (z) && MPFR_IS_NEG (z)))
400    {
401      printf ("Error in mpfr_zeta (-400000001)\n");
402      exit (1);
403    }
404  mpfr_set_str (s, "-400000003", 10, MPFR_RNDZ);
405  mpfr_zeta (z, s, MPFR_RNDN);
406  if (!(mpfr_inf_p (z) && MPFR_IS_POS (z)))
407    {
408      printf ("Error in mpfr_zeta (-400000003)\n");
409      exit (1);
410    }
411
412  mpfr_set_prec (s, 34);
413  mpfr_set_prec (z, 34);
414  mpfr_set_str_binary (s, "-1.111111100001011110000010001010000e-35");
415  mpfr_zeta (z, s, MPFR_RNDD);
416  mpfr_set_str_binary (s, "-1.111111111111111111111111111111111e-2");
417  if (mpfr_cmp (s, z))
418    {
419      printf ("Error in mpfr_zeta, prec=34, MPFR_RNDD\n");
420      mpfr_dump (z);
421      exit (1);
422    }
423
424  /* bug found by nightly tests on June 7, 2007 */
425  mpfr_set_prec (s, 23);
426  mpfr_set_prec (z, 25);
427  mpfr_set_str_binary (s, "-1.0110110110001000000000e-27");
428  mpfr_zeta (z, s, MPFR_RNDN);
429  mpfr_set_prec (s, 25);
430  mpfr_set_str_binary (s, "-1.111111111111111111111111e-2");
431  if (mpfr_cmp (s, z))
432    {
433      printf ("Error in mpfr_zeta, prec=25, MPFR_RNDN\n");
434      printf ("expected "); mpfr_dump (s);
435      printf ("got      "); mpfr_dump (z);
436      exit (1);
437    }
438
439  /* bug reported by Kevin Rauch on 26 Oct 2007 */
440  mpfr_set_prec (s, 128);
441  mpfr_set_prec (z, 128);
442  mpfr_set_str_binary (s, "-0.1000000000000000000000000000000000000000000000000000000000000001E64");
443  inex = mpfr_zeta (z, s, MPFR_RNDN);
444  MPFR_ASSERTN (mpfr_inf_p (z) && MPFR_IS_NEG (z) && inex < 0);
445  inex = mpfr_zeta (z, s, MPFR_RNDU);
446  mpfr_set_inf (s, -1);
447  mpfr_nextabove (s);
448  MPFR_ASSERTN (mpfr_equal_p (z, s) && inex > 0);
449
450  /* bug reported by Fredrik Johansson on 19 Jan 2016 */
451  mpfr_set_prec (s, 536);
452  mpfr_set_ui_2exp (s, 1, -424, MPFR_RNDN);
453  mpfr_sub_ui (s, s, 128, MPFR_RNDN);  /* -128 + 2^(-424) */
454  for (prec = 6; prec <= 536; prec += 8) /* should go through 318 */
455    {
456      mpfr_set_prec (z, prec);
457      mpfr_zeta (z, s, MPFR_RNDD);
458      mpfr_set_prec (y, prec + 10);
459      mpfr_zeta (y, s, MPFR_RNDD);
460      mpfr_prec_round (y, prec, MPFR_RNDD);
461      if (! mpfr_equal_p (z, y))
462        {
463          printf ("mpfr_zeta fails near -128 for inprec=%lu outprec=%lu\n",
464                  (unsigned long) mpfr_get_prec (s), (unsigned long) prec);
465          printf ("expected "); mpfr_dump (y);
466          printf ("got      "); mpfr_dump (z);
467          exit (1);
468        }
469    }
470
471  /* The following test yields an overflow in the error computation.
472     With r10864, this is detected and one gets an assertion failure. */
473  mpfr_set_prec (s, 1025);
474  mpfr_set_si_2exp (s, -1, 1024, MPFR_RNDN);
475  mpfr_nextbelow (s);  /* -(2^1024 + 1) */
476  mpfr_clear_flags ();
477  inex = mpfr_zeta (z, s, MPFR_RNDN);
478  flags = __gmpfr_flags;
479  if (flags != (MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT) ||
480      ! mpfr_inf_p (z) || MPFR_IS_POS (z) || inex >= 0)
481    {
482      printf ("Error in mpfr_zeta for s = -(2^1024 + 1)\nGot ");
483      mpfr_dump (z);
484      printf ("with inex = %d and flags =", inex);
485      flags_out (flags);
486      exit (1);
487    }
488
489  mpfr_clear (s);
490  mpfr_clear (y);
491  mpfr_clear (z);
492
493  /* FIXME: change the last argument back to 5 once the working precision
494     in the mpfr_zeta implementation no longer depends on the precision of
495     the input. */
496  test_generic (MPFR_PREC_MIN, 70, 1);
497  test2 ();
498
499  intermediate_overflow ();
500
501  tests_end_mpfr ();
502  return 0;
503}
504