1 /* mpfr_tgamma -- test file for gamma function 2 3 Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc. 4 Contributed by the Arenaire and Cacao projects, INRIA. 5 6 This file is part of the GNU MPFR Library, and was contributed by Mathieu Dutour. 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 /* Note: there could be an incorrect test about suspicious overflow 29 (MPFR_SUSPICIOUS_OVERFLOW) for x = 2^(-emax) = 0.5 * 2^(emin+1) in 30 RNDZ or RNDD, but this case is never tested in the generic tests. */ 31 #define TEST_FUNCTION mpfr_gamma 32 #define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), 1, RANDS) 33 #include "tgeneric.c" 34 35 static void 36 special (void) 37 { 38 mpfr_t x, y; 39 int inex; 40 41 mpfr_init (x); 42 mpfr_init (y); 43 44 mpfr_set_nan (x); 45 mpfr_gamma (y, x, MPFR_RNDN); 46 if (!mpfr_nan_p (y)) 47 { 48 printf ("Error for gamma(NaN)\n"); 49 exit (1); 50 } 51 52 mpfr_set_inf (x, -1); 53 mpfr_gamma (y, x, MPFR_RNDN); 54 if (!mpfr_nan_p (y)) 55 { 56 printf ("Error for gamma(-Inf)\n"); 57 exit (1); 58 } 59 60 mpfr_set_inf (x, 1); 61 mpfr_gamma (y, x, MPFR_RNDN); 62 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0) 63 { 64 printf ("Error for gamma(+Inf)\n"); 65 exit (1); 66 } 67 68 mpfr_set_ui (x, 0, MPFR_RNDN); 69 mpfr_gamma (y, x, MPFR_RNDN); 70 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0) 71 { 72 printf ("Error for gamma(+0)\n"); 73 exit (1); 74 } 75 76 mpfr_set_ui (x, 0, MPFR_RNDN); 77 mpfr_neg (x, x, MPFR_RNDN); 78 mpfr_gamma (y, x, MPFR_RNDN); 79 if (!mpfr_inf_p (y) || mpfr_sgn (y) > 0) 80 { 81 printf ("Error for gamma(-0)\n"); 82 exit (1); 83 } 84 85 mpfr_set_ui (x, 1, MPFR_RNDN); 86 mpfr_gamma (y, x, MPFR_RNDN); 87 if (mpfr_cmp_ui (y, 1)) 88 { 89 printf ("Error for gamma(1)\n"); 90 exit (1); 91 } 92 93 mpfr_set_si (x, -1, MPFR_RNDN); 94 mpfr_gamma (y, x, MPFR_RNDN); 95 if (!mpfr_nan_p (y)) 96 { 97 printf ("Error for gamma(-1)\n"); 98 exit (1); 99 } 100 101 mpfr_set_prec (x, 53); 102 mpfr_set_prec (y, 53); 103 104 #define CHECK_X1 "1.0762904832837976166" 105 #define CHECK_Y1 "0.96134843256452096050" 106 107 mpfr_set_str (x, CHECK_X1, 10, MPFR_RNDN); 108 mpfr_gamma (y, x, MPFR_RNDN); 109 mpfr_set_str (x, CHECK_Y1, 10, MPFR_RNDN); 110 if (mpfr_cmp (y, x)) 111 { 112 printf ("mpfr_lngamma("CHECK_X1") is wrong:\n" 113 "expected "); 114 mpfr_print_binary (x); putchar ('\n'); 115 printf ("got "); 116 mpfr_print_binary (y); putchar ('\n'); 117 exit (1); 118 } 119 120 #define CHECK_X2 "9.23709516716202383435e-01" 121 #define CHECK_Y2 "1.0502315560291053398" 122 mpfr_set_str (x, CHECK_X2, 10, MPFR_RNDN); 123 mpfr_gamma (y, x, MPFR_RNDN); 124 mpfr_set_str (x, CHECK_Y2, 10, MPFR_RNDN); 125 if (mpfr_cmp (y, x)) 126 { 127 printf ("mpfr_lngamma("CHECK_X2") is wrong:\n" 128 "expected "); 129 mpfr_print_binary (x); putchar ('\n'); 130 printf ("got "); 131 mpfr_print_binary (y); putchar ('\n'); 132 exit (1); 133 } 134 135 mpfr_set_prec (x, 8); 136 mpfr_set_prec (y, 175); 137 mpfr_set_ui (x, 33, MPFR_RNDN); 138 mpfr_gamma (y, x, MPFR_RNDU); 139 mpfr_set_prec (x, 175); 140 mpfr_set_str_binary (x, "0.110010101011010101101000010101010111000110011101001000101011000001100010111001101001011E118"); 141 if (mpfr_cmp (x, y)) 142 { 143 printf ("Error in mpfr_gamma (1)\n"); 144 exit (1); 145 } 146 147 mpfr_set_prec (x, 21); 148 mpfr_set_prec (y, 8); 149 mpfr_set_ui (y, 120, MPFR_RNDN); 150 mpfr_gamma (x, y, MPFR_RNDZ); 151 mpfr_set_prec (y, 21); 152 mpfr_set_str_binary (y, "0.101111101110100110110E654"); 153 if (mpfr_cmp (x, y)) 154 { 155 printf ("Error in mpfr_gamma (120)\n"); 156 printf ("Expected "); mpfr_print_binary (y); puts (""); 157 printf ("Got "); mpfr_print_binary (x); puts (""); 158 exit (1); 159 } 160 161 mpfr_set_prec (x, 3); 162 mpfr_set_prec (y, 206); 163 mpfr_set_str_binary (x, "0.110e10"); 164 inex = mpfr_gamma (y, x, MPFR_RNDN); 165 mpfr_set_prec (x, 206); 166 mpfr_set_str_binary (x, "0.110111100001000001101010010001000111000100000100111000010011100011011111001100011110101000111101101100110001001100110100001001111110000101010000100100011100010011101110000001000010001100010000101001111E6250"); 167 if (mpfr_cmp (x, y)) 168 { 169 printf ("Error in mpfr_gamma (768)\n"); 170 exit (1); 171 } 172 if (inex <= 0) 173 { 174 printf ("Wrong flag for mpfr_gamma (768)\n"); 175 exit (1); 176 } 177 178 /* worst case to exercise retry */ 179 mpfr_set_prec (x, 1000); 180 mpfr_set_prec (y, 869); 181 mpfr_const_pi (x, MPFR_RNDN); 182 mpfr_gamma (y, x, MPFR_RNDN); 183 184 mpfr_set_prec (x, 4); 185 mpfr_set_prec (y, 4); 186 mpfr_set_str_binary (x, "-0.1100E-66"); 187 mpfr_gamma (y, x, MPFR_RNDN); 188 mpfr_set_str_binary (x, "-0.1011E67"); 189 if (mpfr_cmp (x, y)) 190 { 191 printf ("Error for gamma(-0.1100E-66)\n"); 192 exit (1); 193 } 194 195 mpfr_clear (x); 196 mpfr_clear (y); 197 } 198 199 static void 200 special_overflow (void) 201 { 202 mpfr_t x, y; 203 mpfr_exp_t emin, emax; 204 int inex; 205 206 emin = mpfr_get_emin (); 207 emax = mpfr_get_emax (); 208 209 set_emin (-125); 210 set_emax (128); 211 212 mpfr_init2 (x, 24); 213 mpfr_init2 (y, 24); 214 mpfr_set_str_binary (x, "0.101100100000000000110100E7"); 215 mpfr_gamma (y, x, MPFR_RNDN); 216 if (!mpfr_inf_p (y)) 217 { 218 printf ("Overflow error.\n"); 219 mpfr_dump (y); 220 exit (1); 221 } 222 223 /* problem mentioned by Kenneth Wilder, 18 Aug 2005 */ 224 mpfr_set_prec (x, 29); 225 mpfr_set_prec (y, 29); 226 mpfr_set_str (x, "-200000000.5", 10, MPFR_RNDN); /* exact */ 227 mpfr_gamma (y, x, MPFR_RNDN); 228 if (!(mpfr_zero_p (y) && MPFR_SIGN (y) < 0)) 229 { 230 printf ("Error for gamma(-200000000.5)\n"); 231 printf ("expected -0"); 232 printf ("got "); 233 mpfr_dump (y); 234 exit (1); 235 } 236 237 mpfr_set_prec (x, 53); 238 mpfr_set_prec (y, 53); 239 mpfr_set_str (x, "-200000000.1", 10, MPFR_RNDN); 240 mpfr_gamma (y, x, MPFR_RNDN); 241 if (!(mpfr_zero_p (y) && MPFR_SIGN (y) < 0)) 242 { 243 printf ("Error for gamma(-200000000.1), prec=53\n"); 244 printf ("expected -0"); 245 printf ("got "); 246 mpfr_dump (y); 247 exit (1); 248 } 249 250 /* another problem mentioned by Kenneth Wilder, 29 Aug 2005 */ 251 mpfr_set_prec (x, 333); 252 mpfr_set_prec (y, 14); 253 mpfr_set_str (x, "-2.0000000000000000000000005", 10, MPFR_RNDN); 254 mpfr_gamma (y, x, MPFR_RNDN); 255 mpfr_set_prec (x, 14); 256 mpfr_set_str_binary (x, "-11010011110001E66"); 257 if (mpfr_cmp (x, y)) 258 { 259 printf ("Error for gamma(-2.0000000000000000000000005)\n"); 260 printf ("expected "); mpfr_dump (x); 261 printf ("got "); mpfr_dump (y); 262 exit (1); 263 } 264 265 /* another tests from Kenneth Wilder, 31 Aug 2005 */ 266 set_emax (200); 267 set_emin (-200); 268 mpfr_set_prec (x, 38); 269 mpfr_set_prec (y, 54); 270 mpfr_set_str_binary (x, "0.11101111011100111101001001010110101001E-166"); 271 mpfr_gamma (y, x, MPFR_RNDN); 272 mpfr_set_prec (x, 54); 273 mpfr_set_str_binary (x, "0.100010001101100001110110001010111111010000100101011E167"); 274 if (mpfr_cmp (x, y)) 275 { 276 printf ("Error for gamma (test 1)\n"); 277 printf ("expected "); mpfr_dump (x); 278 printf ("got "); mpfr_dump (y); 279 exit (1); 280 } 281 282 set_emax (1000); 283 set_emin (-2000); 284 mpfr_set_prec (x, 38); 285 mpfr_set_prec (y, 71); 286 mpfr_set_str_binary (x, "10101011011100001111111000010111110010E-1034"); 287 /* 184083777010*2^(-1034) */ 288 mpfr_gamma (y, x, MPFR_RNDN); 289 mpfr_set_prec (x, 71); 290 mpfr_set_str_binary (x, "10111111001000011110010001000000000000110011110000000011101011111111100E926"); 291 /* 1762885132679550982140*2^926 */ 292 if (mpfr_cmp (x, y)) 293 { 294 printf ("Error for gamma (test 2)\n"); 295 printf ("expected "); mpfr_dump (x); 296 printf ("got "); mpfr_dump (y); 297 exit (1); 298 } 299 300 mpfr_set_prec (x, 38); 301 mpfr_set_prec (y, 88); 302 mpfr_set_str_binary (x, "10111100111001010000100001100100100101E-104"); 303 /* 202824096037*2^(-104) */ 304 mpfr_gamma (y, x, MPFR_RNDN); 305 mpfr_set_prec (x, 88); 306 mpfr_set_str_binary (x, "1010110101111000111010111100010110101010100110111000001011000111000011101100001101110010E-21"); 307 /* 209715199999500283894743922*2^(-21) */ 308 if (mpfr_cmp (x, y)) 309 { 310 printf ("Error for gamma (test 3)\n"); 311 printf ("expected "); mpfr_dump (x); 312 printf ("got "); mpfr_dump (y); 313 exit (1); 314 } 315 316 mpfr_set_prec (x, 171); 317 mpfr_set_prec (y, 38); 318 mpfr_set_str (x, "-2993155353253689176481146537402947624254601559176535", 10, 319 MPFR_RNDN); 320 mpfr_div_2exp (x, x, 170, MPFR_RNDN); 321 mpfr_gamma (y, x, MPFR_RNDN); 322 mpfr_set_prec (x, 38); 323 mpfr_set_str (x, "201948391737", 10, MPFR_RNDN); 324 mpfr_mul_2exp (x, x, 92, MPFR_RNDN); 325 if (mpfr_cmp (x, y)) 326 { 327 printf ("Error for gamma (test 5)\n"); 328 printf ("expected "); mpfr_dump (x); 329 printf ("got "); mpfr_dump (y); 330 exit (1); 331 } 332 333 set_emin (-500000); 334 mpfr_set_prec (x, 337); 335 mpfr_set_prec (y, 38); 336 mpfr_set_str (x, "-30000.000000000000000000000000000000000000000000001", 10, 337 MPFR_RNDN); 338 mpfr_gamma (y, x, MPFR_RNDN); 339 mpfr_set_prec (x, 38); 340 mpfr_set_str (x, "-3.623795987425E-121243", 10, MPFR_RNDN); 341 if (mpfr_cmp (x, y)) 342 { 343 printf ("Error for gamma (test 7)\n"); 344 printf ("expected "); mpfr_dump (x); 345 printf ("got "); mpfr_dump (y); 346 exit (1); 347 } 348 349 /* was producing infinite loop */ 350 set_emin (emin); 351 mpfr_set_prec (x, 71); 352 mpfr_set_prec (y, 71); 353 mpfr_set_str (x, "-200000000.1", 10, MPFR_RNDN); 354 mpfr_gamma (y, x, MPFR_RNDN); 355 if (!(mpfr_zero_p (y) && MPFR_SIGN (y) < 0)) 356 { 357 printf ("Error for gamma (test 8)\n"); 358 printf ("expected "); mpfr_dump (x); 359 printf ("got "); mpfr_dump (y); 360 exit (1); 361 } 362 363 set_emax (1073741823); 364 mpfr_set_prec (x, 29); 365 mpfr_set_prec (y, 29); 366 mpfr_set_str (x, "423786866", 10, MPFR_RNDN); 367 mpfr_gamma (y, x, MPFR_RNDN); 368 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0) 369 { 370 printf ("Error for gamma(423786866)\n"); 371 exit (1); 372 } 373 374 /* check exact result */ 375 mpfr_set_prec (x, 2); 376 mpfr_set_ui (x, 3, MPFR_RNDN); 377 inex = mpfr_gamma (x, x, MPFR_RNDN); 378 if (inex != 0 || mpfr_cmp_ui (x, 2) != 0) 379 { 380 printf ("Error for gamma(3)\n"); 381 exit (1); 382 } 383 384 mpfr_set_emax (1024); 385 mpfr_set_prec (x, 53); 386 mpfr_set_prec (y, 53); 387 mpfr_set_str_binary (x, "101010110100110011111010000110001000111100000110101E-43"); 388 mpfr_gamma (x, x, MPFR_RNDU); 389 mpfr_set_str_binary (y, "110000011110001000111110110101011110000100001111111E971"); 390 if (mpfr_cmp (x, y) != 0) 391 { 392 printf ("Error for gamma(4)\n"); 393 printf ("expected "); mpfr_dump (y); 394 printf ("got "); mpfr_dump (x); 395 exit (1); 396 } 397 398 set_emin (emin); 399 set_emax (emax); 400 401 /* bug found by Kevin Rauch, 26 Oct 2007 */ 402 mpfr_set_str (x, "1e19", 10, MPFR_RNDN); 403 inex = mpfr_gamma (x, x, MPFR_RNDN); 404 MPFR_ASSERTN(mpfr_inf_p (x) && inex > 0); 405 406 mpfr_clear (y); 407 mpfr_clear (x); 408 } 409 410 /* test gamma on some integral values (from Christopher Creutzig). */ 411 static void 412 gamma_integer (void) 413 { 414 mpz_t n; 415 mpfr_t x, y; 416 unsigned int i; 417 418 mpz_init (n); 419 mpfr_init2 (x, 149); 420 mpfr_init2 (y, 149); 421 422 for (i = 0; i < 100; i++) 423 { 424 mpz_fac_ui (n, i); 425 mpfr_set_ui (x, i+1, MPFR_RNDN); 426 mpfr_gamma (y, x, MPFR_RNDN); 427 mpfr_set_z (x, n, MPFR_RNDN); 428 if (!mpfr_equal_p (x, y)) 429 { 430 printf ("Error for gamma(%u)\n", i+1); 431 printf ("expected "); mpfr_dump (x); 432 printf ("got "); mpfr_dump (y); 433 exit (1); 434 } 435 } 436 mpfr_clear (y); 437 mpfr_clear (x); 438 mpz_clear (n); 439 } 440 441 /* bug found by Kevin Rauch */ 442 static void 443 test20071231 (void) 444 { 445 mpfr_t x; 446 int inex; 447 mpfr_exp_t emin; 448 449 emin = mpfr_get_emin (); 450 mpfr_set_emin (-1000000); 451 452 mpfr_init2 (x, 21); 453 mpfr_set_str (x, "-1000001.5", 10, MPFR_RNDN); 454 inex = mpfr_gamma (x, x, MPFR_RNDN); 455 MPFR_ASSERTN(MPFR_IS_ZERO(x) && MPFR_IS_POS(x) && inex < 0); 456 mpfr_clear (x); 457 458 mpfr_set_emin (emin); 459 460 mpfr_init2 (x, 53); 461 mpfr_set_str (x, "-1000000001.5", 10, MPFR_RNDN); 462 inex = mpfr_gamma (x, x, MPFR_RNDN); 463 MPFR_ASSERTN(MPFR_IS_ZERO(x) && MPFR_IS_POS(x) && inex < 0); 464 mpfr_clear (x); 465 } 466 467 /* bug found by Stathis, only occurs on 32-bit machines */ 468 static void 469 test20100709 (void) 470 { 471 mpfr_t x; 472 int inex; 473 474 mpfr_init2 (x, 100); 475 mpfr_set_str (x, "-4.6308260837372266e+07", 10, MPFR_RNDN); 476 inex = mpfr_gamma (x, x, MPFR_RNDN); 477 MPFR_ASSERTN(MPFR_IS_ZERO(x) && MPFR_IS_NEG(x) && inex > 0); 478 mpfr_clear (x); 479 } 480 481 static void 482 exprange (void) 483 { 484 mpfr_exp_t emin, emax; 485 mpfr_t x, y, z; 486 int inex1, inex2; 487 unsigned int flags1, flags2; 488 489 emin = mpfr_get_emin (); 490 emax = mpfr_get_emax (); 491 492 mpfr_init2 (x, 16); 493 mpfr_inits2 (8, y, z, (mpfr_ptr) 0); 494 495 mpfr_set_ui_2exp (x, 5, -1, MPFR_RNDN); 496 mpfr_clear_flags (); 497 inex1 = mpfr_gamma (y, x, MPFR_RNDN); 498 flags1 = __gmpfr_flags; 499 MPFR_ASSERTN (mpfr_inexflag_p ()); 500 mpfr_set_emin (0); 501 mpfr_clear_flags (); 502 inex2 = mpfr_gamma (z, x, MPFR_RNDN); 503 flags2 = __gmpfr_flags; 504 MPFR_ASSERTN (mpfr_inexflag_p ()); 505 mpfr_set_emin (emin); 506 if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z)) 507 { 508 printf ("Error in exprange (test1)\n"); 509 printf ("x = "); 510 mpfr_dump (x); 511 printf ("Expected inex1 = %d, flags1 = %u, ", SIGN (inex1), flags1); 512 mpfr_dump (y); 513 printf ("Got inex2 = %d, flags2 = %u, ", SIGN (inex2), flags2); 514 mpfr_dump (z); 515 exit (1); 516 } 517 518 mpfr_set_ui_2exp (x, 32769, -60, MPFR_RNDN); 519 mpfr_clear_flags (); 520 inex1 = mpfr_gamma (y, x, MPFR_RNDD); 521 flags1 = __gmpfr_flags; 522 MPFR_ASSERTN (mpfr_inexflag_p ()); 523 mpfr_set_emax (45); 524 mpfr_clear_flags (); 525 inex2 = mpfr_gamma (z, x, MPFR_RNDD); 526 flags2 = __gmpfr_flags; 527 MPFR_ASSERTN (mpfr_inexflag_p ()); 528 mpfr_set_emax (emax); 529 if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z)) 530 { 531 printf ("Error in exprange (test2)\n"); 532 printf ("x = "); 533 mpfr_dump (x); 534 printf ("Expected inex1 = %d, flags1 = %u, ", SIGN (inex1), flags1); 535 mpfr_dump (y); 536 printf ("Got inex2 = %d, flags2 = %u, ", SIGN (inex2), flags2); 537 mpfr_dump (z); 538 exit (1); 539 } 540 541 mpfr_set_emax (44); 542 mpfr_clear_flags (); 543 inex1 = mpfr_check_range (y, inex1, MPFR_RNDD); 544 flags1 = __gmpfr_flags; 545 MPFR_ASSERTN (mpfr_inexflag_p ()); 546 mpfr_clear_flags (); 547 inex2 = mpfr_gamma (z, x, MPFR_RNDD); 548 flags2 = __gmpfr_flags; 549 MPFR_ASSERTN (mpfr_inexflag_p ()); 550 mpfr_set_emax (emax); 551 if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z)) 552 { 553 printf ("Error in exprange (test3)\n"); 554 printf ("x = "); 555 mpfr_dump (x); 556 printf ("Expected inex1 = %d, flags1 = %u, ", SIGN (inex1), flags1); 557 mpfr_dump (y); 558 printf ("Got inex2 = %d, flags2 = %u, ", SIGN (inex2), flags2); 559 mpfr_dump (z); 560 exit (1); 561 } 562 563 mpfr_set_ui_2exp (x, 1, -60, MPFR_RNDN); 564 mpfr_clear_flags (); 565 inex1 = mpfr_gamma (y, x, MPFR_RNDD); 566 flags1 = __gmpfr_flags; 567 MPFR_ASSERTN (mpfr_inexflag_p ()); 568 mpfr_set_emax (60); 569 mpfr_clear_flags (); 570 inex2 = mpfr_gamma (z, x, MPFR_RNDD); 571 flags2 = __gmpfr_flags; 572 MPFR_ASSERTN (mpfr_inexflag_p ()); 573 mpfr_set_emax (emax); 574 if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z)) 575 { 576 printf ("Error in exprange (test4)\n"); 577 printf ("x = "); 578 mpfr_dump (x); 579 printf ("Expected inex1 = %d, flags1 = %u, ", SIGN (inex1), flags1); 580 mpfr_dump (y); 581 printf ("Got inex2 = %d, flags2 = %u, ", SIGN (inex2), flags2); 582 mpfr_dump (z); 583 exit (1); 584 } 585 586 MPFR_ASSERTN (MPFR_EMIN_MIN == - MPFR_EMAX_MAX); 587 mpfr_set_emin (MPFR_EMIN_MIN); 588 mpfr_set_emax (MPFR_EMAX_MAX); 589 mpfr_set_ui (x, 0, MPFR_RNDN); 590 mpfr_nextabove (x); /* x = 2^(emin - 1) */ 591 mpfr_set_inf (y, 1); 592 inex1 = 1; 593 flags1 = MPFR_FLAGS_INEXACT | MPFR_FLAGS_OVERFLOW; 594 mpfr_clear_flags (); 595 /* MPFR_RNDU: overflow, infinity since 1/x = 2^(emax + 1) */ 596 inex2 = mpfr_gamma (z, x, MPFR_RNDU); 597 flags2 = __gmpfr_flags; 598 MPFR_ASSERTN (mpfr_inexflag_p ()); 599 if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z)) 600 { 601 printf ("Error in exprange (test5)\n"); 602 printf ("x = "); 603 mpfr_dump (x); 604 printf ("Expected inex1 = %d, flags1 = %u, ", SIGN (inex1), flags1); 605 mpfr_dump (y); 606 printf ("Got inex2 = %d, flags2 = %u, ", SIGN (inex2), flags2); 607 mpfr_dump (z); 608 exit (1); 609 } 610 mpfr_clear_flags (); 611 /* MPFR_RNDN: overflow, infinity since 1/x = 2^(emax + 1) */ 612 inex2 = mpfr_gamma (z, x, MPFR_RNDN); 613 flags2 = __gmpfr_flags; 614 MPFR_ASSERTN (mpfr_inexflag_p ()); 615 if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z)) 616 { 617 printf ("Error in exprange (test6)\n"); 618 printf ("x = "); 619 mpfr_dump (x); 620 printf ("Expected inex1 = %d, flags1 = %u, ", SIGN (inex1), flags1); 621 mpfr_dump (y); 622 printf ("Got inex2 = %d, flags2 = %u, ", SIGN (inex2), flags2); 623 mpfr_dump (z); 624 exit (1); 625 } 626 mpfr_nextbelow (y); 627 inex1 = -1; 628 mpfr_clear_flags (); 629 /* MPFR_RNDD: overflow, maxnum since 1/x = 2^(emax + 1) */ 630 inex2 = mpfr_gamma (z, x, MPFR_RNDD); 631 flags2 = __gmpfr_flags; 632 MPFR_ASSERTN (mpfr_inexflag_p ()); 633 if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z)) 634 { 635 printf ("Error in exprange (test7)\n"); 636 printf ("x = "); 637 mpfr_dump (x); 638 printf ("Expected inex1 = %d, flags1 = %u, ", SIGN (inex1), flags1); 639 mpfr_dump (y); 640 printf ("Got inex2 = %d, flags2 = %u, ", SIGN (inex2), flags2); 641 mpfr_dump (z); 642 exit (1); 643 } 644 mpfr_mul_2ui (x, x, 1, MPFR_RNDN); /* x = 2^emin */ 645 mpfr_set_inf (y, 1); 646 inex1 = 1; 647 flags1 = MPFR_FLAGS_INEXACT | MPFR_FLAGS_OVERFLOW; 648 mpfr_clear_flags (); 649 /* MPFR_RNDU: overflow, infinity since 1/x = 2^emax */ 650 inex2 = mpfr_gamma (z, x, MPFR_RNDU); 651 flags2 = __gmpfr_flags; 652 MPFR_ASSERTN (mpfr_inexflag_p ()); 653 if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z)) 654 { 655 printf ("Error in exprange (test8)\n"); 656 printf ("x = "); 657 mpfr_dump (x); 658 printf ("Expected inex1 = %d, flags1 = %u, ", SIGN (inex1), flags1); 659 mpfr_dump (y); 660 printf ("Got inex2 = %d, flags2 = %u, ", SIGN (inex2), flags2); 661 mpfr_dump (z); 662 exit (1); 663 } 664 mpfr_clear_flags (); 665 /* MPFR_RNDN: overflow, infinity since 1/x = 2^emax */ 666 inex2 = mpfr_gamma (z, x, MPFR_RNDN); 667 flags2 = __gmpfr_flags; 668 MPFR_ASSERTN (mpfr_inexflag_p ()); 669 if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z)) 670 { 671 printf ("Error in exprange (test9)\n"); 672 printf ("x = "); 673 mpfr_dump (x); 674 printf ("Expected inex1 = %d, flags1 = %u, ", SIGN (inex1), flags1); 675 mpfr_dump (y); 676 printf ("Got inex2 = %d, flags2 = %u, ", SIGN (inex2), flags2); 677 mpfr_dump (z); 678 exit (1); 679 } 680 mpfr_nextbelow (y); 681 inex1 = -1; 682 flags1 = MPFR_FLAGS_INEXACT; 683 mpfr_clear_flags (); 684 /* MPFR_RNDD: no overflow, maxnum since 1/x = 2^emax and euler > 0 */ 685 inex2 = mpfr_gamma (z, x, MPFR_RNDD); 686 flags2 = __gmpfr_flags; 687 MPFR_ASSERTN (mpfr_inexflag_p ()); 688 if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z)) 689 { 690 printf ("Error in exprange (test10)\n"); 691 printf ("x = "); 692 mpfr_dump (x); 693 printf ("Expected inex1 = %d, flags1 = %u, ", SIGN (inex1), flags1); 694 mpfr_dump (y); 695 printf ("Got inex2 = %d, flags2 = %u, ", SIGN (inex2), flags2); 696 mpfr_dump (z); 697 exit (1); 698 } 699 mpfr_set_emin (emin); 700 mpfr_set_emax (emax); 701 702 mpfr_clears (x, y, z, (mpfr_ptr) 0); 703 } 704 705 static int 706 tiny_aux (int stop, mpfr_exp_t e) 707 { 708 mpfr_t x, y, z; 709 int r, s, spm, inex, err = 0; 710 int expected_dir[2][5] = { { 1, -1, 1, -1, 1 }, { 1, 1, 1, -1, -1 } }; 711 mpfr_exp_t saved_emax; 712 713 saved_emax = mpfr_get_emax (); 714 715 mpfr_init2 (x, 32); 716 mpfr_inits2 (8, y, z, (mpfr_ptr) 0); 717 718 mpfr_set_ui_2exp (x, 1, e, MPFR_RNDN); 719 spm = 1; 720 for (s = 0; s < 2; s++) 721 { 722 RND_LOOP(r) 723 { 724 mpfr_rnd_t rr = (mpfr_rnd_t) r; 725 mpfr_exp_t exponent, emax; 726 727 /* Exponent of the rounded value in unbounded exponent range. */ 728 exponent = expected_dir[s][r] < 0 && s == 0 ? - e : 1 - e; 729 730 for (emax = exponent - 1; emax <= exponent; emax++) 731 { 732 unsigned int flags, expected_flags = MPFR_FLAGS_INEXACT; 733 int overflow, expected_inex = expected_dir[s][r]; 734 735 if (emax > MPFR_EMAX_MAX) 736 break; 737 mpfr_set_emax (emax); 738 739 mpfr_clear_flags (); 740 inex = mpfr_gamma (y, x, rr); 741 flags = __gmpfr_flags; 742 mpfr_clear_flags (); 743 mpfr_set_si_2exp (z, spm, - e, MPFR_RNDU); 744 overflow = mpfr_overflow_p (); 745 /* z is 1/x - euler rounded toward +inf */ 746 747 if (overflow && rr == MPFR_RNDN && s == 1) 748 expected_inex = -1; 749 750 if (expected_inex < 0) 751 mpfr_nextbelow (z); /* 1/x - euler rounded toward -inf */ 752 753 if (exponent > emax) 754 expected_flags |= MPFR_FLAGS_OVERFLOW; 755 756 if (!(mpfr_equal_p (y, z) && flags == expected_flags 757 && SAME_SIGN (inex, expected_inex))) 758 { 759 printf ("Error in tiny for s = %d, r = %s, emax = %" 760 MPFR_EXP_FSPEC "d%s\n on ", 761 s, mpfr_print_rnd_mode (rr), emax, 762 exponent > emax ? " (overflow)" : ""); 763 mpfr_dump (x); 764 printf (" expected inex = %2d, ", expected_inex); 765 mpfr_dump (z); 766 printf (" got inex = %2d, ", SIGN (inex)); 767 mpfr_dump (y); 768 printf (" expected flags = %u, got %u\n", 769 expected_flags, flags); 770 if (stop) 771 exit (1); 772 err = 1; 773 } 774 } 775 } 776 mpfr_neg (x, x, MPFR_RNDN); 777 spm = - spm; 778 } 779 780 mpfr_clears (x, y, z, (mpfr_ptr) 0); 781 mpfr_set_emax (saved_emax); 782 return err; 783 } 784 785 static void 786 tiny (int stop) 787 { 788 mpfr_exp_t emin; 789 int err = 0; 790 791 emin = mpfr_get_emin (); 792 793 /* Note: in r7499, exponent -17 will select the generic code (in 794 tiny_aux, x has precision 32), while the other exponent values 795 will select special code for tiny values. */ 796 err |= tiny_aux (stop, -17); 797 err |= tiny_aux (stop, -999); 798 err |= tiny_aux (stop, mpfr_get_emin ()); 799 800 if (emin != MPFR_EMIN_MIN) 801 { 802 mpfr_set_emin (MPFR_EMIN_MIN); 803 err |= tiny_aux (stop, MPFR_EMIN_MIN); 804 mpfr_set_emin (emin); 805 } 806 807 if (err) 808 exit (1); 809 } 810 811 int 812 main (int argc, char *argv[]) 813 { 814 tests_start_mpfr (); 815 816 special (); 817 special_overflow (); 818 exprange (); 819 tiny (argc == 1); 820 test_generic (2, 100, 2); 821 gamma_integer (); 822 test20071231 (); 823 test20100709 (); 824 825 data_check ("data/gamma", mpfr_gamma, "mpfr_gamma"); 826 827 tests_end_mpfr (); 828 return 0; 829 } 830