1/* mpfr_tgamma -- test file for gamma function 2 3Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc. 4Contributed by the Arenaire and Cacao projects, INRIA. 5 6This file is part of the GNU MPFR Library, and was contributed by Mathieu Dutour. 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 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 35static void 36special (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 199static void 200special_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). */ 411static void 412gamma_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 */ 442static void 443test20071231 (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 */ 468static void 469test20100709 (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 481static void 482exprange (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 705static int 706tiny_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 785static void 786tiny (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 811int 812main (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