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