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