1 /* mpfr_tlgamma -- test file for lgamma function 2 3 Copyright 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. 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 int 29 mpfr_lgamma_nosign (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd) 30 { 31 int inex, sign; 32 33 inex = mpfr_lgamma (y, &sign, x, rnd); 34 if (!MPFR_IS_SINGULAR (y)) 35 { 36 MPFR_ASSERTN (sign == 1 || sign == -1); 37 if (sign == -1 && (rnd == MPFR_RNDN || rnd == MPFR_RNDZ)) 38 { 39 mpfr_neg (y, y, MPFR_RNDN); 40 inex = -inex; 41 /* This is a way to check with the generic tests, that the value 42 returned in the sign variable is consistent, but warning! The 43 tested function depends on the rounding mode: it is 44 * lgamma(x) = log(|Gamma(x)|) in MPFR_RNDD and MPFR_RNDU; 45 * lgamma(x) * sign(Gamma(x)) in MPFR_RNDN and MPFR_RNDZ. */ 46 } 47 } 48 49 return inex; 50 } 51 52 #define TEST_FUNCTION mpfr_lgamma_nosign 53 #include "tgeneric.c" 54 55 static void 56 special (void) 57 { 58 mpfr_t x, y; 59 int inex; 60 int sign; 61 mpfr_exp_t emin, emax; 62 63 mpfr_init (x); 64 mpfr_init (y); 65 66 mpfr_set_nan (x); 67 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 68 if (!mpfr_nan_p (y)) 69 { 70 printf ("Error for lgamma(NaN)\n"); 71 exit (1); 72 } 73 74 mpfr_set_inf (x, -1); 75 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 76 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0) 77 { 78 printf ("Error for lgamma(-Inf)\n"); 79 exit (1); 80 } 81 82 mpfr_set_inf (x, 1); 83 sign = -17; 84 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 85 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0 || sign != 1) 86 { 87 printf ("Error for lgamma(+Inf)\n"); 88 exit (1); 89 } 90 91 mpfr_set_ui (x, 0, MPFR_RNDN); 92 sign = -17; 93 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 94 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0 || sign != 1) 95 { 96 printf ("Error for lgamma(+0)\n"); 97 exit (1); 98 } 99 100 mpfr_set_ui (x, 0, MPFR_RNDN); 101 mpfr_neg (x, x, MPFR_RNDN); 102 sign = -17; 103 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 104 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0 || sign != -1) 105 { 106 printf ("Error for lgamma(-0)\n"); 107 exit (1); 108 } 109 110 mpfr_set_ui (x, 1, MPFR_RNDN); 111 sign = -17; 112 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 113 if (MPFR_IS_NAN (y) || mpfr_cmp_ui (y, 0) || MPFR_IS_NEG (y) || sign != 1) 114 { 115 printf ("Error for lgamma(1)\n"); 116 exit (1); 117 } 118 119 mpfr_set_si (x, -1, MPFR_RNDN); 120 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 121 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0) 122 { 123 printf ("Error for lgamma(-1)\n"); 124 exit (1); 125 } 126 127 mpfr_set_ui (x, 2, MPFR_RNDN); 128 sign = -17; 129 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 130 if (MPFR_IS_NAN (y) || mpfr_cmp_ui (y, 0) || MPFR_IS_NEG (y) || sign != 1) 131 { 132 printf ("Error for lgamma(2)\n"); 133 exit (1); 134 } 135 136 mpfr_set_prec (x, 53); 137 mpfr_set_prec (y, 53); 138 139 #define CHECK_X1 "1.0762904832837976166" 140 #define CHECK_Y1 "-0.039418362817587634939" 141 142 mpfr_set_str (x, CHECK_X1, 10, MPFR_RNDN); 143 sign = -17; 144 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 145 mpfr_set_str (x, CHECK_Y1, 10, MPFR_RNDN); 146 if (mpfr_equal_p (y, x) == 0 || sign != 1) 147 { 148 printf ("mpfr_lgamma("CHECK_X1") is wrong:\n" 149 "expected "); 150 mpfr_print_binary (x); putchar ('\n'); 151 printf ("got "); 152 mpfr_print_binary (y); putchar ('\n'); 153 exit (1); 154 } 155 156 #define CHECK_X2 "9.23709516716202383435e-01" 157 #define CHECK_Y2 "0.049010669407893718563" 158 mpfr_set_str (x, CHECK_X2, 10, MPFR_RNDN); 159 sign = -17; 160 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 161 mpfr_set_str (x, CHECK_Y2, 10, MPFR_RNDN); 162 if (mpfr_equal_p (y, x) == 0 || sign != 1) 163 { 164 printf ("mpfr_lgamma("CHECK_X2") is wrong:\n" 165 "expected "); 166 mpfr_print_binary (x); putchar ('\n'); 167 printf ("got "); 168 mpfr_print_binary (y); putchar ('\n'); 169 exit (1); 170 } 171 172 mpfr_set_prec (x, 8); 173 mpfr_set_prec (y, 175); 174 mpfr_set_ui (x, 33, MPFR_RNDN); 175 sign = -17; 176 mpfr_lgamma (y, &sign, x, MPFR_RNDU); 177 mpfr_set_prec (x, 175); 178 mpfr_set_str_binary (x, "0.1010001100011101101011001101110010100001000001000001110011000001101100001111001001000101011011100100010101011110100111110101010100010011010010000101010111001100011000101111E7"); 179 if (mpfr_equal_p (x, y) == 0 || sign != 1) 180 { 181 printf ("Error in mpfr_lgamma (1)\n"); 182 exit (1); 183 } 184 185 mpfr_set_prec (x, 21); 186 mpfr_set_prec (y, 8); 187 mpfr_set_ui (y, 120, MPFR_RNDN); 188 sign = -17; 189 mpfr_lgamma (x, &sign, y, MPFR_RNDZ); 190 mpfr_set_prec (y, 21); 191 mpfr_set_str_binary (y, "0.111000101000001100101E9"); 192 if (mpfr_equal_p (x, y) == 0 || sign != 1) 193 { 194 printf ("Error in mpfr_lgamma (120)\n"); 195 printf ("Expected "); mpfr_print_binary (y); puts (""); 196 printf ("Got "); mpfr_print_binary (x); puts (""); 197 exit (1); 198 } 199 200 mpfr_set_prec (x, 3); 201 mpfr_set_prec (y, 206); 202 mpfr_set_str_binary (x, "0.110e10"); 203 sign = -17; 204 inex = mpfr_lgamma (y, &sign, x, MPFR_RNDN); 205 mpfr_set_prec (x, 206); 206 mpfr_set_str_binary (x, "0.10000111011000000011100010101001100110001110000111100011000100100110110010001011011110101001111011110110000001010100111011010000000011100110110101100111000111010011110010000100010111101010001101000110101001E13"); 207 if (mpfr_equal_p (x, y) == 0 || sign != 1) 208 { 209 printf ("Error in mpfr_lgamma (768)\n"); 210 exit (1); 211 } 212 if (inex >= 0) 213 { 214 printf ("Wrong flag for mpfr_lgamma (768)\n"); 215 exit (1); 216 } 217 218 mpfr_set_prec (x, 4); 219 mpfr_set_prec (y, 4); 220 mpfr_set_str_binary (x, "0.1100E-66"); 221 sign = -17; 222 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 223 mpfr_set_str_binary (x, "0.1100E6"); 224 if (mpfr_equal_p (x, y) == 0 || sign != 1) 225 { 226 printf ("Error for lgamma(0.1100E-66)\n"); 227 printf ("Expected "); 228 mpfr_dump (x); 229 printf ("Got "); 230 mpfr_dump (y); 231 exit (1); 232 } 233 234 mpfr_set_prec (x, 256); 235 mpfr_set_prec (y, 32); 236 mpfr_set_si_2exp (x, -1, 200, MPFR_RNDN); 237 mpfr_add_ui (x, x, 1, MPFR_RNDN); 238 mpfr_div_2ui (x, x, 1, MPFR_RNDN); 239 sign = -17; 240 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 241 mpfr_set_prec (x, 32); 242 mpfr_set_str_binary (x, "-0.10001000111011111011000010100010E207"); 243 if (mpfr_equal_p (x, y) == 0 || sign != 1) 244 { 245 printf ("Error for lgamma(-2^199+0.5)\n"); 246 printf ("Got "); 247 mpfr_dump (y); 248 printf ("instead of "); 249 mpfr_dump (x); 250 exit (1); 251 } 252 253 mpfr_set_prec (x, 256); 254 mpfr_set_prec (y, 32); 255 mpfr_set_si_2exp (x, -1, 200, MPFR_RNDN); 256 mpfr_sub_ui (x, x, 1, MPFR_RNDN); 257 mpfr_div_2ui (x, x, 1, MPFR_RNDN); 258 sign = -17; 259 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 260 mpfr_set_prec (x, 32); 261 mpfr_set_str_binary (x, "-0.10001000111011111011000010100010E207"); 262 if (mpfr_equal_p (x, y) == 0 || sign != -1) 263 { 264 printf ("Error for lgamma(-2^199-0.5)\n"); 265 printf ("Got "); 266 mpfr_dump (y); 267 printf ("with sign %d instead of ", sign); 268 mpfr_dump (x); 269 printf ("with sign -1.\n"); 270 exit (1); 271 } 272 273 mpfr_set_prec (x, 10); 274 mpfr_set_prec (y, 10); 275 mpfr_set_str_binary (x, "-0.1101111000E-3"); 276 inex = mpfr_lgamma (y, &sign, x, MPFR_RNDN); 277 mpfr_set_str_binary (x, "10.01001011"); 278 if (mpfr_equal_p (x, y) == 0 || sign != -1 || inex >= 0) 279 { 280 printf ("Error for lgamma(-0.1101111000E-3)\n"); 281 printf ("Got "); 282 mpfr_dump (y); 283 printf ("instead of "); 284 mpfr_dump (x); 285 printf ("with sign %d instead of -1 (inex=%d).\n", sign, inex); 286 exit (1); 287 } 288 289 mpfr_set_prec (x, 18); 290 mpfr_set_prec (y, 28); 291 mpfr_set_str_binary (x, "-1.10001101010001101e-196"); 292 inex = mpfr_lgamma (y, &sign, x, MPFR_RNDN); 293 mpfr_set_prec (x, 28); 294 mpfr_set_str_binary (x, "0.100001110110101011011010011E8"); 295 MPFR_ASSERTN (mpfr_equal_p (x, y) && inex < 0); 296 297 /* values reported by Kaveh Ghazi on 14 Jul 2007, where mpfr_lgamma() 298 takes forever */ 299 #define VAL1 "-0.11100001001010110111001010001001001011110100110000110E-55" 300 #define OUT1 "100110.01000000010111001110110101110101001001100110111" 301 #define VAL2 "-0.11100001001010110111001010001001001011110011111111100E-55" 302 #define OUT2 "100110.0100000001011100111011010111010100100110011111" 303 #define VAL3 "-0.11100001001010110111001010001001001001110101101010100E-55" 304 #define OUT3 "100110.01000000010111001110110101110101001011110111011" 305 #define VAL4 "-0.10001111110110110100100100000000001111110001001001011E-57" 306 #define OUT4 "101000.0001010111110011101101000101111111010001100011" 307 #define VAL5 "-0.10001111110110110100100100000000001111011111100001000E-57" 308 #define OUT5 "101000.00010101111100111011010001011111110100111000001" 309 #define VAL6 "-0.10001111110110110100100100000000001111011101100011001E-57" 310 #define OUT6 "101000.0001010111110011101101000101111111010011101111" 311 312 mpfr_set_prec (x, 53); 313 mpfr_set_prec (y, 53); 314 315 mpfr_set_str_binary (x, VAL1); 316 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 317 mpfr_set_str_binary (x, OUT1); 318 MPFR_ASSERTN(sign == -1 && mpfr_equal_p(x, y)); 319 320 mpfr_set_str_binary (x, VAL2); 321 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 322 mpfr_set_str_binary (x, OUT2); 323 MPFR_ASSERTN(sign == -1 && mpfr_equal_p (x, y)); 324 325 mpfr_set_str_binary (x, VAL3); 326 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 327 mpfr_set_str_binary (x, OUT3); 328 MPFR_ASSERTN(sign == -1 && mpfr_equal_p (x, y)); 329 330 mpfr_set_str_binary (x, VAL4); 331 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 332 mpfr_set_str_binary (x, OUT4); 333 MPFR_ASSERTN(sign == -1 && mpfr_equal_p (x, y)); 334 335 mpfr_set_str_binary (x, VAL5); 336 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 337 mpfr_set_str_binary (x, OUT5); 338 MPFR_ASSERTN(sign == -1 && mpfr_equal_p (x, y)); 339 340 mpfr_set_str_binary (x, VAL6); 341 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 342 mpfr_set_str_binary (x, OUT6); 343 MPFR_ASSERTN(sign == -1 && mpfr_equal_p (x, y)); 344 345 /* further test from Kaveh Ghazi */ 346 mpfr_set_str_binary (x, "-0.10011010101001010010001110010111010111011101010111001E-53"); 347 mpfr_lgamma (y, &sign, x, MPFR_RNDN); 348 mpfr_set_str_binary (x, "100101.00111101101010000000101010111010001111001101111"); 349 MPFR_ASSERTN(sign == -1 && mpfr_equal_p (x, y)); 350 351 /* bug found by Kevin Rauch on 26 Oct 2007 */ 352 emin = mpfr_get_emin (); 353 emax = mpfr_get_emax (); 354 mpfr_set_emin (-1000000000); 355 mpfr_set_emax (1000000000); 356 mpfr_set_ui (x, 1, MPFR_RNDN); 357 mpfr_lgamma (x, &sign, x, MPFR_RNDN); 358 MPFR_ASSERTN(mpfr_get_emin () == -1000000000); 359 MPFR_ASSERTN(mpfr_get_emax () == 1000000000); 360 mpfr_set_emin (emin); 361 mpfr_set_emax (emax); 362 363 /* two other bugs reported by Kevin Rauch on 27 Oct 2007 */ 364 mpfr_set_prec (x, 128); 365 mpfr_set_prec (y, 128); 366 mpfr_set_str_binary (x, "0.11000110011110111111110010100110000000000000000000000000000000000000000000000000000000000000000001000011000110100100110111101010E-765689"); 367 inex = mpfr_lgamma (y, &sign, x, MPFR_RNDN); 368 mpfr_set_str_binary (x, "10000001100100101111011011010000111010001001110000111010011000101001011111011111110011011010110100101111110111001001010100011101E-108"); 369 MPFR_ASSERTN(inex < 0 && mpfr_cmp (y, x) == 0 && sign > 0); 370 371 mpfr_set_prec (x, 128); 372 mpfr_set_prec (y, 256); 373 mpfr_set_str_binary (x, "0.1011111111111111100000111011111E-31871"); 374 inex = mpfr_lgamma (y, &sign, x, MPFR_RNDN); 375 mpfr_set_prec (x, 256); 376 mpfr_set_str (x, "AC9729B83707E6797612D0D76DAF42B1240A677FF1B6E3783FD4E53037143B1P-237", 16, MPFR_RNDN); 377 MPFR_ASSERTN(inex < 0 && mpfr_cmp (y, x) == 0 && sign > 0); 378 379 mpfr_clear (x); 380 mpfr_clear (y); 381 } 382 383 static int 384 mpfr_lgamma1 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t r) 385 { 386 int sign; 387 388 return mpfr_lgamma (y, &sign, x, r); 389 } 390 391 int 392 main (void) 393 { 394 tests_start_mpfr (); 395 396 special (); 397 test_generic (2, 100, 2); 398 399 data_check ("data/lgamma", mpfr_lgamma1, "mpfr_lgamma"); 400 401 tests_end_mpfr (); 402 return 0; 403 } 404