1/* mpfr_tlgamma -- test file for lgamma function 2 3Copyright 2005-2023 Free Software Foundation, Inc. 4Contributed by the AriC and Caramba projects, INRIA. 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 20https://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 "mpfr-test.h" 24 25static int 26mpfr_lgamma_nosign (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd) 27{ 28 int inex, sign; 29 30 inex = mpfr_lgamma (y, &sign, x, rnd); 31 if (!MPFR_IS_SINGULAR (y)) 32 { 33 MPFR_ASSERTN (sign == 1 || sign == -1); 34 if (sign == -1 && (rnd == MPFR_RNDN || rnd == MPFR_RNDZ)) 35 { 36 mpfr_neg (y, y, MPFR_RNDN); 37 inex = -inex; 38 /* This is a way to check with the generic tests, that the value 39 returned in the sign variable is consistent, but warning! The 40 tested function depends on the rounding mode: it is 41 * lgamma(x) = log(|Gamma(x)|) in MPFR_RNDD and MPFR_RNDU; 42 * lgamma(x) * sign(Gamma(x)) in MPFR_RNDN and MPFR_RNDZ. */ 43 } 44 } 45 46 return inex; 47} 48 49#define TEST_FUNCTION mpfr_lgamma_nosign 50#include "tgeneric.c" 51 52static void 53special (void) 54{ 55 mpfr_t x, y; 56 int inex; 57 int sign; 58 mpfr_exp_t emin, emax; 59 60 mpfr_init (x); 61 mpfr_init (y); 62 63 mpfr_set_nan (x); 64 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 65 if (!mpfr_nan_p (y)) 66 { 67 printf ("Error for lgamma(NaN)\n"); 68 exit (1); 69 } 70 71 mpfr_set_inf (x, -1); 72 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 73 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0) 74 { 75 printf ("Error for lgamma(-Inf)\n"); 76 exit (1); 77 } 78 79 mpfr_set_inf (x, 1); 80 sign = -17; 81 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 82 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0 || sign != 1) 83 { 84 printf ("Error for lgamma(+Inf)\n"); 85 exit (1); 86 } 87 88 mpfr_set_ui (x, 0, MPFR_RNDN); 89 sign = -17; 90 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 91 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0 || sign != 1) 92 { 93 printf ("Error for lgamma(+0)\n"); 94 exit (1); 95 } 96 97 mpfr_set_ui (x, 0, MPFR_RNDN); 98 mpfr_neg (x, x, MPFR_RNDN); 99 sign = -17; 100 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 101 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0 || sign != -1) 102 { 103 printf ("Error for lgamma(-0)\n"); 104 exit (1); 105 } 106 107 mpfr_set_ui (x, 1, MPFR_RNDN); 108 sign = -17; 109 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 110 if (MPFR_IS_NAN (y) || mpfr_cmp_ui (y, 0) || MPFR_IS_NEG (y) || sign != 1) 111 { 112 printf ("Error for lgamma(1)\n"); 113 exit (1); 114 } 115 116 mpfr_set_si (x, -1, MPFR_RNDN); 117 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 118 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0) 119 { 120 printf ("Error for lgamma(-1)\n"); 121 exit (1); 122 } 123 124 mpfr_set_ui (x, 2, MPFR_RNDN); 125 sign = -17; 126 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 127 if (MPFR_IS_NAN (y) || mpfr_cmp_ui (y, 0) || MPFR_IS_NEG (y) || sign != 1) 128 { 129 printf ("Error for lgamma(2)\n"); 130 exit (1); 131 } 132 133 mpfr_set_prec (x, 53); 134 mpfr_set_prec (y, 53); 135 136#define CHECK_X1 "1.0762904832837976166" 137#define CHECK_Y1 "-0.039418362817587634939" 138 139 mpfr_set_str (x, CHECK_X1, 10, MPFR_RNDN); 140 sign = -17; 141 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 142 mpfr_set_str (x, CHECK_Y1, 10, MPFR_RNDN); 143 if (mpfr_equal_p (y, x) == 0 || sign != 1) 144 { 145 printf ("mpfr_lgamma(" CHECK_X1 ") is wrong:\n" 146 "expected "); 147 mpfr_dump (x); 148 printf ("got "); 149 mpfr_dump (y); 150 exit (1); 151 } 152 153#define CHECK_X2 "9.23709516716202383435e-01" 154#define CHECK_Y2 "0.049010669407893718563" 155 mpfr_set_str (x, CHECK_X2, 10, MPFR_RNDN); 156 sign = -17; 157 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 158 mpfr_set_str (x, CHECK_Y2, 10, MPFR_RNDN); 159 if (mpfr_equal_p (y, x) == 0 || sign != 1) 160 { 161 printf ("mpfr_lgamma(" CHECK_X2 ") is wrong:\n" 162 "expected "); 163 mpfr_dump (x); 164 printf ("got "); 165 mpfr_dump (y); 166 exit (1); 167 } 168 169 mpfr_set_prec (x, 8); 170 mpfr_set_prec (y, 175); 171 mpfr_set_ui (x, 33, MPFR_RNDN); 172 sign = -17; 173 mpfr_lgamma (y, &sign, x, MPFR_RNDU); 174 mpfr_set_prec (x, 175); 175 mpfr_set_str_binary (x, "0.1010001100011101101011001101110010100001000001000001110011000001101100001111001001000101011011100100010101011110100111110101010100010011010010000101010111001100011000101111E7"); 176 if (mpfr_equal_p (x, y) == 0 || sign != 1) 177 { 178 printf ("Error in mpfr_lgamma (1)\n"); 179 exit (1); 180 } 181 182 mpfr_set_prec (x, 21); 183 mpfr_set_prec (y, 8); 184 mpfr_set_ui (y, 120, MPFR_RNDN); 185 sign = -17; 186 mpfr_lgamma (x, &sign, y, MPFR_RNDZ); 187 mpfr_set_prec (y, 21); 188 mpfr_set_str_binary (y, "0.111000101000001100101E9"); 189 if (mpfr_equal_p (x, y) == 0 || sign != 1) 190 { 191 printf ("Error in mpfr_lgamma (120)\n"); 192 printf ("Expected "); mpfr_dump (y); 193 printf ("Got "); mpfr_dump (x); 194 exit (1); 195 } 196 197 mpfr_set_prec (x, 3); 198 mpfr_set_prec (y, 206); 199 mpfr_set_str_binary (x, "0.110e10"); 200 sign = -17; 201 inex = mpfr_lgamma (y, &sign, x, MPFR_RNDN); 202 mpfr_set_prec (x, 206); 203 mpfr_set_str_binary (x, "0.10000111011000000011100010101001100110001110000111100011000100100110110010001011011110101001111011110110000001010100111011010000000011100110110101100111000111010011110010000100010111101010001101000110101001E13"); 204 if (mpfr_equal_p (x, y) == 0 || sign != 1) 205 { 206 printf ("Error in mpfr_lgamma (768)\n"); 207 exit (1); 208 } 209 if (inex >= 0) 210 { 211 printf ("Wrong flag for mpfr_lgamma (768)\n"); 212 exit (1); 213 } 214 215 mpfr_set_prec (x, 4); 216 mpfr_set_prec (y, 4); 217 mpfr_set_str_binary (x, "0.1100E-66"); 218 sign = -17; 219 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 220 mpfr_set_str_binary (x, "0.1100E6"); 221 if (mpfr_equal_p (x, y) == 0 || sign != 1) 222 { 223 printf ("Error for lgamma(0.1100E-66)\n"); 224 printf ("Expected "); 225 mpfr_dump (x); 226 printf ("Got "); 227 mpfr_dump (y); 228 exit (1); 229 } 230 231 mpfr_set_prec (x, 256); 232 mpfr_set_prec (y, 32); 233 mpfr_set_si_2exp (x, -1, 200, MPFR_RNDN); 234 mpfr_add_ui (x, x, 1, MPFR_RNDN); 235 mpfr_div_2ui (x, x, 1, MPFR_RNDN); 236 sign = -17; 237 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 238 mpfr_set_prec (x, 32); 239 mpfr_set_str_binary (x, "-0.10001000111011111011000010100010E207"); 240 if (mpfr_equal_p (x, y) == 0 || sign != 1) 241 { 242 printf ("Error for lgamma(-2^199+0.5)\n"); 243 printf ("Got "); 244 mpfr_dump (y); 245 printf ("instead of "); 246 mpfr_dump (x); 247 exit (1); 248 } 249 250 mpfr_set_prec (x, 256); 251 mpfr_set_prec (y, 32); 252 mpfr_set_si_2exp (x, -1, 200, MPFR_RNDN); 253 mpfr_sub_ui (x, x, 1, MPFR_RNDN); 254 mpfr_div_2ui (x, x, 1, MPFR_RNDN); 255 sign = -17; 256 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 257 mpfr_set_prec (x, 32); 258 mpfr_set_str_binary (x, "-0.10001000111011111011000010100010E207"); 259 if (mpfr_equal_p (x, y) == 0 || sign != -1) 260 { 261 printf ("Error for lgamma(-2^199-0.5)\n"); 262 printf ("Got "); 263 mpfr_dump (y); 264 printf ("with sign %d instead of ", sign); 265 mpfr_dump (x); 266 printf ("with sign -1.\n"); 267 exit (1); 268 } 269 270 mpfr_set_prec (x, 10); 271 mpfr_set_prec (y, 10); 272 mpfr_set_str_binary (x, "-0.1101111000E-3"); 273 inex = mpfr_lgamma (y, &sign, x, MPFR_RNDN); 274 mpfr_set_str_binary (x, "10.01001011"); 275 if (mpfr_equal_p (x, y) == 0 || sign != -1 || inex >= 0) 276 { 277 printf ("Error for lgamma(-0.1101111000E-3)\n"); 278 printf ("Got "); 279 mpfr_dump (y); 280 printf ("instead of "); 281 mpfr_dump (x); 282 printf ("with sign %d instead of -1 (inex=%d).\n", sign, inex); 283 exit (1); 284 } 285 286 mpfr_set_prec (x, 18); 287 mpfr_set_prec (y, 28); 288 mpfr_set_str_binary (x, "-1.10001101010001101e-196"); 289 inex = mpfr_lgamma (y, &sign, x, MPFR_RNDN); 290 mpfr_set_prec (x, 28); 291 mpfr_set_str_binary (x, "0.100001110110101011011010011E8"); 292 MPFR_ASSERTN (mpfr_equal_p (x, y) && inex < 0); 293 294 /* values reported by Kaveh Ghazi on 14 Jul 2007, where mpfr_lgamma() 295 takes forever */ 296#define VAL1 "-0.11100001001010110111001010001001001011110100110000110E-55" 297#define OUT1 "100110.01000000010111001110110101110101001001100110111" 298#define VAL2 "-0.11100001001010110111001010001001001011110011111111100E-55" 299#define OUT2 "100110.0100000001011100111011010111010100100110011111" 300#define VAL3 "-0.11100001001010110111001010001001001001110101101010100E-55" 301#define OUT3 "100110.01000000010111001110110101110101001011110111011" 302#define VAL4 "-0.10001111110110110100100100000000001111110001001001011E-57" 303#define OUT4 "101000.0001010111110011101101000101111111010001100011" 304#define VAL5 "-0.10001111110110110100100100000000001111011111100001000E-57" 305#define OUT5 "101000.00010101111100111011010001011111110100111000001" 306#define VAL6 "-0.10001111110110110100100100000000001111011101100011001E-57" 307#define OUT6 "101000.0001010111110011101101000101111111010011101111" 308 309 mpfr_set_prec (x, 53); 310 mpfr_set_prec (y, 53); 311 312 mpfr_set_str_binary (x, VAL1); 313 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 314 mpfr_set_str_binary (x, OUT1); 315 MPFR_ASSERTN(sign == -1 && mpfr_equal_p(x, y)); 316 317 mpfr_set_str_binary (x, VAL2); 318 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 319 mpfr_set_str_binary (x, OUT2); 320 MPFR_ASSERTN(sign == -1 && mpfr_equal_p (x, y)); 321 322 mpfr_set_str_binary (x, VAL3); 323 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 324 mpfr_set_str_binary (x, OUT3); 325 MPFR_ASSERTN(sign == -1 && mpfr_equal_p (x, y)); 326 327 mpfr_set_str_binary (x, VAL4); 328 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 329 mpfr_set_str_binary (x, OUT4); 330 MPFR_ASSERTN(sign == -1 && mpfr_equal_p (x, y)); 331 332 mpfr_set_str_binary (x, VAL5); 333 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 334 mpfr_set_str_binary (x, OUT5); 335 MPFR_ASSERTN(sign == -1 && mpfr_equal_p (x, y)); 336 337 mpfr_set_str_binary (x, VAL6); 338 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 339 mpfr_set_str_binary (x, OUT6); 340 MPFR_ASSERTN(sign == -1 && mpfr_equal_p (x, y)); 341 342 /* further test from Kaveh Ghazi */ 343 mpfr_set_str_binary (x, "-0.10011010101001010010001110010111010111011101010111001E-53"); 344 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 345 mpfr_set_str_binary (x, "100101.00111101101010000000101010111010001111001101111"); 346 MPFR_ASSERTN(sign == -1 && mpfr_equal_p (x, y)); 347 348 /* bug found by Kevin Rauch on 26 Oct 2007 */ 349 emin = mpfr_get_emin (); 350 emax = mpfr_get_emax (); 351 set_emin (-1000000000); 352 set_emax (1000000000); 353 mpfr_set_ui (x, 1, MPFR_RNDN); 354 mpfr_lgamma (x, &sign, x, MPFR_RNDN); 355 MPFR_ASSERTN(mpfr_get_emin () == -1000000000); 356 MPFR_ASSERTN(mpfr_get_emax () == 1000000000); 357 set_emin (emin); 358 set_emax (emax); 359 360 /* two other bugs reported by Kevin Rauch on 27 Oct 2007 */ 361 mpfr_set_prec (x, 128); 362 mpfr_set_prec (y, 128); 363 mpfr_set_str_binary (x, "0.11000110011110111111110010100110000000000000000000000000000000000000000000000000000000000000000001000011000110100100110111101010E-765689"); 364 inex = mpfr_lgamma (y, &sign, x, MPFR_RNDN); 365 mpfr_set_str_binary (x, "10000001100100101111011011010000111010001001110000111010011000101001011111011111110011011010110100101111110111001001010100011101E-108"); 366 MPFR_ASSERTN(inex < 0 && mpfr_cmp (y, x) == 0 && sign > 0); 367 368 mpfr_set_prec (x, 128); 369 mpfr_set_prec (y, 256); 370 mpfr_set_str_binary (x, "0.1011111111111111100000111011111E-31871"); 371 inex = mpfr_lgamma (y, &sign, x, MPFR_RNDN); 372 mpfr_set_prec (x, 256); 373 mpfr_set_str (x, "AC9729B83707E6797612D0D76DAF42B1240A677FF1B6E3783FD4E53037143B1P-237", 16, MPFR_RNDN); 374 MPFR_ASSERTN(inex < 0 && mpfr_cmp (y, x) == 0 && sign > 0); 375 376 mpfr_clear (x); 377 mpfr_clear (y); 378} 379 380static int 381mpfr_lgamma1 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t r) 382{ 383 int sign; 384 385 return mpfr_lgamma (y, &sign, x, r); 386} 387 388/* Since r10377, the following test causes a "too much memory" error 389 when MPFR is built with Debian's tcc 0.9.27~git20151227.933c223-1 390 on x86_64. The problem came from __gmpfr_ceil_log2, now fixed in 391 r10443 (according to the integer promotion rules, this appeared to 392 be a bug in tcc, not in MPFR; however, relying on such an obscure 393 rule was not a good idea). */ 394static void 395tcc_bug20160606 (void) 396{ 397 mpfr_t x, y; 398 int sign; 399 400 mpfr_init2 (x, 53); 401 mpfr_init2 (y, 53); 402 mpfr_set_ui_2exp (x, 1, -1, MPFR_RNDN); 403 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 404 mpfr_clear (x); 405 mpfr_clear (y); 406} 407 408/* With r12088, mpfr_lgamma is much too slow with a reduced emax that 409 yields an overflow, even though this case is easier. In practice, 410 this test will hang. */ 411static void 412bug20180110 (void) 413{ 414 mpfr_exp_t emax, e; 415 mpfr_t x, y, z; 416 mpfr_flags_t flags, eflags; 417 int i, inex, sign; 418 419 emax = mpfr_get_emax (); 420 421 mpfr_init2 (x, 2); 422 mpfr_inits2 (64, y, z, (mpfr_ptr) 0); 423 eflags = MPFR_FLAGS_INEXACT | MPFR_FLAGS_OVERFLOW; 424 425 for (i = 10; i <= 30; i++) 426 { 427 mpfr_set_si_2exp (x, -1, -(1L << i), MPFR_RNDN); /* -2^(-2^i) */ 428 mpfr_lgamma (y, &sign, x, MPFR_RNDZ); 429 e = mpfr_get_exp (y); 430 set_emax (e - 1); 431 mpfr_clear_flags (); 432 inex = mpfr_lgamma (y, &sign, x, MPFR_RNDZ); 433 flags = __gmpfr_flags; 434 mpfr_set_inf (z, 1); 435 mpfr_nextbelow (z); 436 set_emax (emax); 437 if (! (mpfr_equal_p (y, z) && SAME_SIGN (inex, -1) && flags == eflags)) 438 { 439 printf ("Error in bug20180110 for i = %d:\n", i); 440 printf ("Expected "); 441 mpfr_dump (z); 442 printf ("with inex = %d and flags =", -1); 443 flags_out (eflags); 444 printf ("Got "); 445 mpfr_dump (y); 446 printf ("with inex = %d and flags =", inex); 447 flags_out (flags); 448 exit (1); 449 } 450 } 451 452 mpfr_clears (x, y, z, (mpfr_ptr) 0); 453} 454 455int 456main (void) 457{ 458 tests_start_mpfr (); 459 460 tcc_bug20160606 (); 461 bug20180110 (); 462 463 special (); 464 test_generic (MPFR_PREC_MIN, 100, 2); 465 466 data_check ("data/lgamma", mpfr_lgamma1, "mpfr_lgamma"); 467 468 tests_end_mpfr (); 469 return 0; 470} 471