tests/tlgamma.c

   1 /* mpfr_tlgamma -- test file for lgamma function
   2
   3 Copyright 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012 Free Software Foundation, Inc.
   4 Contributed by the AriC and Caramel 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 }