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