1 /* tzeta -- test file for the Riemann Zeta function
2
3 Copyright 2003-2023 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba projects, INRIA.
5
6 This file is part of the GNU MPFR Library.
7
8 The GNU MPFR Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12
13 The GNU MPFR Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
17
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
20 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22
23 #include "mpfr-test.h"
24
25 static void
26 test1 (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
124 static 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
170 static void
171 test2 (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. */
205 static void
206 intermediate_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' */
251 int
252 main (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 }