1/* mpfr_tlgamma -- test file for lgamma function 2 3Copyright 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. 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 int 29mpfr_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 55static void 56special (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 383static int 384mpfr_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 391int 392main (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