1 /* tzeta -- test file for the Riemann Zeta function 2 3 Copyright 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4 Contributed by the AriC and Caramel 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 <stdio.h> 24 #include <stdlib.h> 25 26 #include "mpfr-test.h" 27 28 static void 29 test1 (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 113 static 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 159 static void 160 test2 (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' */ 198 int 199 main (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_clear_divby0(); 253 mpfr_zeta (z, s, MPFR_RNDN); 254 if (!mpfr_inf_p (z) || MPFR_SIGN (z) < 0 || !mpfr_divby0_p()) 255 { 256 printf ("Error in mpfr_zeta for s = 1 (should be +inf) with divby0 flag\n"); 257 exit (1); 258 } 259 260 mpfr_set_str_binary (s, "0.1100011101110111111111111010000110010111001011001011"); 261 mpfr_set_str_binary (y, "-0.11111101111011001001001111111000101010000100000100100E2"); 262 mpfr_zeta (z, s, MPFR_RNDN); 263 if (mpfr_cmp (z, y) != 0) 264 { 265 printf ("Error in mpfr_zeta (1,RNDN)\n"); 266 exit (1); 267 } 268 mpfr_zeta (z, s, MPFR_RNDZ); 269 if (mpfr_cmp (z, y) != 0) 270 { 271 printf ("Error in mpfr_zeta (1,RNDZ)\n"); 272 exit (1); 273 } 274 mpfr_zeta (z, s, MPFR_RNDU); 275 if (mpfr_cmp (z, y) != 0) 276 { 277 printf ("Error in mpfr_zeta (1,RNDU)\n"); 278 exit (1); 279 } 280 mpfr_zeta (z, s, MPFR_RNDD); 281 mpfr_nexttoinf (y); 282 if (mpfr_cmp (z, y) != 0) 283 { 284 printf ("Error in mpfr_zeta (1,RNDD)\n"); 285 exit (1); 286 } 287 288 mpfr_set_str_binary (s, "0.10001011010011100110010001100100001011000010011001011"); 289 mpfr_set_str_binary (y, "-0.11010011010010101101110111011010011101111101111010110E1"); 290 mpfr_zeta (z, s, MPFR_RNDN); 291 if (mpfr_cmp (z, y) != 0) 292 { 293 printf ("Error in mpfr_zeta (2,RNDN)\n"); 294 exit (1); 295 } 296 mpfr_zeta (z, s, MPFR_RNDZ); 297 if (mpfr_cmp (z, y) != 0) 298 { 299 printf ("Error in mpfr_zeta (2,RNDZ)\n"); 300 exit (1); 301 } 302 mpfr_zeta (z, s, MPFR_RNDU); 303 if (mpfr_cmp (z, y) != 0) 304 { 305 printf ("Error in mpfr_zeta (2,RNDU)\n"); 306 exit (1); 307 } 308 mpfr_zeta (z, s, MPFR_RNDD); 309 mpfr_nexttoinf (y); 310 if (mpfr_cmp (z, y) != 0) 311 { 312 printf ("Error in mpfr_zeta (2,RNDD)\n"); 313 exit (1); 314 } 315 316 mpfr_set_str_binary (s, "0.1100111110100001111110111000110101111001011101000101"); 317 mpfr_set_str_binary (y, "-0.10010111010110000111011111001101100001111011000001010E3"); 318 mpfr_zeta (z, s, MPFR_RNDN); 319 if (mpfr_cmp (z, y) != 0) 320 { 321 printf ("Error in mpfr_zeta (3,RNDN)\n"); 322 exit (1); 323 } 324 mpfr_zeta (z, s, MPFR_RNDD); 325 if (mpfr_cmp (z, y) != 0) 326 { 327 printf ("Error in mpfr_zeta (3,RNDD)\n"); 328 exit (1); 329 } 330 mpfr_nexttozero (y); 331 mpfr_zeta (z, s, MPFR_RNDZ); 332 if (mpfr_cmp (z, y) != 0) 333 { 334 printf ("Error in mpfr_zeta (3,RNDZ)\n"); 335 exit (1); 336 } 337 mpfr_zeta (z, s, MPFR_RNDU); 338 if (mpfr_cmp (z, y) != 0) 339 { 340 printf ("Error in mpfr_zeta (3,RNDU)\n"); 341 exit (1); 342 } 343 344 mpfr_set_str (s, "-400000001", 10, MPFR_RNDZ); 345 mpfr_zeta (z, s, MPFR_RNDN); 346 if (!(mpfr_inf_p (z) && MPFR_SIGN(z) < 0)) 347 { 348 printf ("Error in mpfr_zeta (-400000001)\n"); 349 exit (1); 350 } 351 mpfr_set_str (s, "-400000003", 10, MPFR_RNDZ); 352 mpfr_zeta (z, s, MPFR_RNDN); 353 if (!(mpfr_inf_p (z) && MPFR_SIGN(z) > 0)) 354 { 355 printf ("Error in mpfr_zeta (-400000003)\n"); 356 exit (1); 357 } 358 359 mpfr_set_prec (s, 34); 360 mpfr_set_prec (z, 34); 361 mpfr_set_str_binary (s, "-1.111111100001011110000010001010000e-35"); 362 mpfr_zeta (z, s, MPFR_RNDD); 363 mpfr_set_str_binary (s, "-1.111111111111111111111111111111111e-2"); 364 if (mpfr_cmp (s, z)) 365 { 366 printf ("Error in mpfr_zeta, prec=34, MPFR_RNDD\n"); 367 mpfr_dump (z); 368 exit (1); 369 } 370 371 /* bug found by nightly tests on June 7, 2007 */ 372 mpfr_set_prec (s, 23); 373 mpfr_set_prec (z, 25); 374 mpfr_set_str_binary (s, "-1.0110110110001000000000e-27"); 375 mpfr_zeta (z, s, MPFR_RNDN); 376 mpfr_set_prec (s, 25); 377 mpfr_set_str_binary (s, "-1.111111111111111111111111e-2"); 378 if (mpfr_cmp (s, z)) 379 { 380 printf ("Error in mpfr_zeta, prec=25, MPFR_RNDN\n"); 381 printf ("expected "); mpfr_dump (s); 382 printf ("got "); mpfr_dump (z); 383 exit (1); 384 } 385 386 /* bug reported by Kevin Rauch on 26 Oct 2007 */ 387 mpfr_set_prec (s, 128); 388 mpfr_set_prec (z, 128); 389 mpfr_set_str_binary (s, "-0.1000000000000000000000000000000000000000000000000000000000000001E64"); 390 inex = mpfr_zeta (z, s, MPFR_RNDN); 391 MPFR_ASSERTN (mpfr_inf_p (z) && MPFR_SIGN (z) < 0 && inex < 0); 392 inex = mpfr_zeta (z, s, MPFR_RNDU); 393 mpfr_set_inf (s, -1); 394 mpfr_nextabove (s); 395 MPFR_ASSERTN (mpfr_equal_p (z, s) && inex > 0); 396 397 mpfr_clear (s); 398 mpfr_clear (y); 399 mpfr_clear (z); 400 401 test_generic (2, 70, 5); 402 test2 (); 403 404 tests_end_mpfr (); 405 return 0; 406 } 407