1 /* Test file for gamma function
2
3 Copyright 2001-2023 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba 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 https://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 "mpfr-test.h"
24
25 /* Note: there could be an incorrect test about suspicious overflow
26 (MPFR_SUSPICIOUS_OVERFLOW) for x = 2^(-emax) = 0.5 * 2^(emin+1) in
27 RNDZ or RNDD, but this case is never tested in the generic tests. */
28 #define TEST_FUNCTION mpfr_gamma
29 #define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), 1, RANDS)
30 #include "tgeneric.c"
31
32 static void
33 special (void)
34 {
35 mpfr_t x, y;
36 int inex;
37
38 mpfr_init (x);
39 mpfr_init (y);
40
41 mpfr_set_nan (x);
42 mpfr_gamma (y, x, MPFR_RNDN);
43 if (!mpfr_nan_p (y))
44 {
45 printf ("Error for gamma(NaN)\n");
46 exit (1);
47 }
48
49 mpfr_set_inf (x, -1);
50 mpfr_gamma (y, x, MPFR_RNDN);
51 if (!mpfr_nan_p (y))
52 {
53 printf ("Error for gamma(-Inf)\n");
54 exit (1);
55 }
56
57 mpfr_set_inf (x, 1);
58 mpfr_gamma (y, x, MPFR_RNDN);
59 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0)
60 {
61 printf ("Error for gamma(+Inf)\n");
62 exit (1);
63 }
64
65 mpfr_set_ui (x, 0, MPFR_RNDN);
66 mpfr_gamma (y, x, MPFR_RNDN);
67 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0)
68 {
69 printf ("Error for gamma(+0)\n");
70 exit (1);
71 }
72
73 mpfr_set_ui (x, 0, MPFR_RNDN);
74 mpfr_neg (x, x, MPFR_RNDN);
75 mpfr_gamma (y, x, MPFR_RNDN);
76 if (!mpfr_inf_p (y) || mpfr_sgn (y) > 0)
77 {
78 printf ("Error for gamma(-0)\n");
79 exit (1);
80 }
81
82 mpfr_set_ui (x, 1, MPFR_RNDN);
83 mpfr_gamma (y, x, MPFR_RNDN);
84 if (mpfr_cmp_ui (y, 1))
85 {
86 printf ("Error for gamma(1)\n");
87 exit (1);
88 }
89
90 mpfr_set_si (x, -1, MPFR_RNDN);
91 mpfr_gamma (y, x, MPFR_RNDN);
92 if (!mpfr_nan_p (y))
93 {
94 printf ("Error for gamma(-1)\n");
95 exit (1);
96 }
97
98 mpfr_set_prec (x, 53);
99 mpfr_set_prec (y, 53);
100
101 #define CHECK_X1 "1.0762904832837976166"
102 #define CHECK_Y1 "0.96134843256452096050"
103
104 mpfr_set_str (x, CHECK_X1, 10, MPFR_RNDN);
105 mpfr_gamma (y, x, MPFR_RNDN);
106 mpfr_set_str (x, CHECK_Y1, 10, MPFR_RNDN);
107 if (mpfr_cmp (y, x))
108 {
109 printf ("mpfr_lngamma(" CHECK_X1 ") is wrong:\n"
110 "expected ");
111 mpfr_dump (x);
112 printf ("got ");
113 mpfr_dump (y);
114 exit (1);
115 }
116
117 #define CHECK_X2 "9.23709516716202383435e-01"
118 #define CHECK_Y2 "1.0502315560291053398"
119 mpfr_set_str (x, CHECK_X2, 10, MPFR_RNDN);
120 mpfr_gamma (y, x, MPFR_RNDN);
121 mpfr_set_str (x, CHECK_Y2, 10, MPFR_RNDN);
122 if (mpfr_cmp (y, x))
123 {
124 printf ("mpfr_lngamma(" CHECK_X2 ") is wrong:\n"
125 "expected ");
126 mpfr_dump (x);
127 printf ("got ");
128 mpfr_dump (y);
129 exit (1);
130 }
131
132 mpfr_set_prec (x, 8);
133 mpfr_set_prec (y, 175);
134 mpfr_set_ui (x, 33, MPFR_RNDN);
135 mpfr_gamma (y, x, MPFR_RNDU);
136 mpfr_set_prec (x, 175);
137 mpfr_set_str_binary (x, "0.110010101011010101101000010101010111000110011101001000101011000001100010111001101001011E118");
138 if (mpfr_cmp (x, y))
139 {
140 printf ("Error in mpfr_gamma (1)\n");
141 exit (1);
142 }
143
144 mpfr_set_prec (x, 21);
145 mpfr_set_prec (y, 8);
146 mpfr_set_ui (y, 120, MPFR_RNDN);
147 mpfr_gamma (x, y, MPFR_RNDZ);
148 mpfr_set_prec (y, 21);
149 mpfr_set_str_binary (y, "0.101111101110100110110E654");
150 if (mpfr_cmp (x, y))
151 {
152 printf ("Error in mpfr_gamma (120)\n");
153 printf ("Expected "); mpfr_dump (y);
154 printf ("Got "); mpfr_dump (x);
155 exit (1);
156 }
157
158 mpfr_set_prec (x, 3);
159 mpfr_set_prec (y, 206);
160 mpfr_set_str_binary (x, "0.110e10");
161 inex = mpfr_gamma (y, x, MPFR_RNDN);
162 mpfr_set_prec (x, 206);
163 mpfr_set_str_binary (x, "0.110111100001000001101010010001000111000100000100111000010011100011011111001100011110101000111101101100110001001100110100001001111110000101010000100100011100010011101110000001000010001100010000101001111E6250");
164 if (mpfr_cmp (x, y))
165 {
166 printf ("Error in mpfr_gamma (768)\n");
167 exit (1);
168 }
169 if (inex <= 0)
170 {
171 printf ("Wrong flag for mpfr_gamma (768)\n");
172 exit (1);
173 }
174
175 /* worst case to exercise retry */
176 mpfr_set_prec (x, 1000);
177 mpfr_set_prec (y, 869);
178 mpfr_const_pi (x, MPFR_RNDN);
179 mpfr_gamma (y, x, MPFR_RNDN);
180
181 mpfr_set_prec (x, 4);
182 mpfr_set_prec (y, 4);
183 mpfr_set_str_binary (x, "-0.1100E-66");
184 mpfr_gamma (y, x, MPFR_RNDN);
185 mpfr_set_str_binary (x, "-0.1011E67");
186 if (mpfr_cmp (x, y))
187 {
188 printf ("Error for gamma(-0.1100E-66)\n");
189 exit (1);
190 }
191
192 mpfr_set_prec (x, 2);
193 mpfr_set_prec (y, 2);
194 mpfr_set_ui (x, 2, MPFR_RNDN);
195 mpfr_clear_inexflag ();
196 mpfr_gamma (y, x, MPFR_RNDN);
197 if (mpfr_inexflag_p ())
198 {
199 printf ("Wrong inexact flag for gamma(2)\n");
200 printf ("expected 0, got 1\n");
201 exit (1);
202 }
203
204 mpfr_clear (x);
205 mpfr_clear (y);
206 }
207
208 static void
209 special_overflow (void)
210 {
211 mpfr_t x, y;
212 mpfr_exp_t emin, emax;
213 int inex;
214
215 emin = mpfr_get_emin ();
216 emax = mpfr_get_emax ();
217
218 set_emin (-125);
219 set_emax (128);
220
221 mpfr_init2 (x, 24);
222 mpfr_init2 (y, 24);
223 mpfr_set_str_binary (x, "0.101100100000000000110100E7");
224 mpfr_gamma (y, x, MPFR_RNDN);
225 if (!mpfr_inf_p (y))
226 {
227 printf ("Overflow error.\n");
228 mpfr_dump (y);
229 exit (1);
230 }
231
232 /* problem mentioned by Kenneth Wilder, 18 Aug 2005 */
233 mpfr_set_prec (x, 29);
234 mpfr_set_prec (y, 29);
235 mpfr_set_str (x, "-200000000.5", 10, MPFR_RNDN); /* exact */
236 mpfr_gamma (y, x, MPFR_RNDN);
237 if (!(mpfr_zero_p (y) && MPFR_IS_NEG (y)))
238 {
239 printf ("Error for gamma(-200000000.5)\n");
240 printf ("expected -0");
241 printf ("got ");
242 mpfr_dump (y);
243 exit (1);
244 }
245
246 mpfr_set_prec (x, 53);
247 mpfr_set_prec (y, 53);
248 mpfr_set_str (x, "-200000000.1", 10, MPFR_RNDN);
249 mpfr_gamma (y, x, MPFR_RNDN);
250 if (!(mpfr_zero_p (y) && MPFR_IS_NEG (y)))
251 {
252 printf ("Error for gamma(-200000000.1), prec=53\n");
253 printf ("expected -0");
254 printf ("got ");
255 mpfr_dump (y);
256 exit (1);
257 }
258
259 /* another problem mentioned by Kenneth Wilder, 29 Aug 2005 */
260 mpfr_set_prec (x, 333);
261 mpfr_set_prec (y, 14);
262 mpfr_set_str (x, "-2.0000000000000000000000005", 10, MPFR_RNDN);
263 mpfr_gamma (y, x, MPFR_RNDN);
264 mpfr_set_prec (x, 14);
265 mpfr_set_str_binary (x, "-11010011110001E66");
266 if (mpfr_cmp (x, y))
267 {
268 printf ("Error for gamma(-2.0000000000000000000000005)\n");
269 printf ("expected "); mpfr_dump (x);
270 printf ("got "); mpfr_dump (y);
271 exit (1);
272 }
273
274 /* another tests from Kenneth Wilder, 31 Aug 2005 */
275 set_emax (200);
276 set_emin (-200);
277 mpfr_set_prec (x, 38);
278 mpfr_set_prec (y, 54);
279 mpfr_set_str_binary (x, "0.11101111011100111101001001010110101001E-166");
280 mpfr_gamma (y, x, MPFR_RNDN);
281 mpfr_set_prec (x, 54);
282 mpfr_set_str_binary (x, "0.100010001101100001110110001010111111010000100101011E167");
283 if (mpfr_cmp (x, y))
284 {
285 printf ("Error for gamma (test 1)\n");
286 printf ("expected "); mpfr_dump (x);
287 printf ("got "); mpfr_dump (y);
288 exit (1);
289 }
290
291 set_emax (1000);
292 set_emin (-2000);
293 mpfr_set_prec (x, 38);
294 mpfr_set_prec (y, 71);
295 mpfr_set_str_binary (x, "10101011011100001111111000010111110010E-1034");
296 /* 184083777010*2^(-1034) */
297 mpfr_gamma (y, x, MPFR_RNDN);
298 mpfr_set_prec (x, 71);
299 mpfr_set_str_binary (x, "10111111001000011110010001000000000000110011110000000011101011111111100E926");
300 /* 1762885132679550982140*2^926 */
301 if (mpfr_cmp (x, y))
302 {
303 printf ("Error for gamma (test 2)\n");
304 printf ("expected "); mpfr_dump (x);
305 printf ("got "); mpfr_dump (y);
306 exit (1);
307 }
308
309 mpfr_set_prec (x, 38);
310 mpfr_set_prec (y, 88);
311 mpfr_set_str_binary (x, "10111100111001010000100001100100100101E-104");
312 /* 202824096037*2^(-104) */
313 mpfr_gamma (y, x, MPFR_RNDN);
314 mpfr_set_prec (x, 88);
315 mpfr_set_str_binary (x, "1010110101111000111010111100010110101010100110111000001011000111000011101100001101110010E-21");
316 /* 209715199999500283894743922*2^(-21) */
317 if (mpfr_cmp (x, y))
318 {
319 printf ("Error for gamma (test 3)\n");
320 printf ("expected "); mpfr_dump (x);
321 printf ("got "); mpfr_dump (y);
322 exit (1);
323 }
324
325 mpfr_set_prec (x, 171);
326 mpfr_set_prec (y, 38);
327 mpfr_set_str (x, "-2993155353253689176481146537402947624254601559176535", 10,
328 MPFR_RNDN);
329 mpfr_div_2ui (x, x, 170, MPFR_RNDN);
330 mpfr_gamma (y, x, MPFR_RNDN);
331 mpfr_set_prec (x, 38);
332 mpfr_set_str (x, "201948391737", 10, MPFR_RNDN);
333 mpfr_mul_2ui (x, x, 92, MPFR_RNDN);
334 if (mpfr_cmp (x, y))
335 {
336 printf ("Error for gamma (test 5)\n");
337 printf ("expected "); mpfr_dump (x);
338 printf ("got "); mpfr_dump (y);
339 exit (1);
340 }
341
342 set_emin (-500000);
343 mpfr_set_prec (x, 337);
344 mpfr_set_prec (y, 38);
345 mpfr_set_str (x, "-30000.000000000000000000000000000000000000000000001", 10,
346 MPFR_RNDN);
347 mpfr_gamma (y, x, MPFR_RNDN);
348 mpfr_set_prec (x, 38);
349 mpfr_set_str (x, "-3.623795987425E-121243", 10, MPFR_RNDN);
350 if (mpfr_cmp (x, y))
351 {
352 printf ("Error for gamma (test 7)\n");
353 printf ("expected "); mpfr_dump (x);
354 printf ("got "); mpfr_dump (y);
355 exit (1);
356 }
357
358 /* was producing infinite loop */
359 set_emin (emin);
360 mpfr_set_prec (x, 71);
361 mpfr_set_prec (y, 71);
362 mpfr_set_str (x, "-200000000.1", 10, MPFR_RNDN);
363 mpfr_gamma (y, x, MPFR_RNDN);
364 if (!(mpfr_zero_p (y) && MPFR_IS_NEG (y)))
365 {
366 printf ("Error for gamma (test 8)\n");
367 printf ("expected "); mpfr_dump (x);
368 printf ("got "); mpfr_dump (y);
369 exit (1);
370 }
371
372 set_emax (1073741821);
373 mpfr_set_prec (x, 29);
374 mpfr_set_prec (y, 29);
375 mpfr_set_str (x, "423786866", 10, MPFR_RNDN);
376 mpfr_gamma (y, x, MPFR_RNDN);
377 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0)
378 {
379 printf ("Error for gamma(423786866)\n");
380 exit (1);
381 }
382
383 /* check exact result */
384 mpfr_set_prec (x, 2);
385 mpfr_set_ui (x, 3, MPFR_RNDN);
386 inex = mpfr_gamma (x, x, MPFR_RNDN);
387 if (inex != 0 || mpfr_cmp_ui (x, 2) != 0)
388 {
389 printf ("Error for gamma(3)\n");
390 exit (1);
391 }
392
393 set_emax (1024);
394 mpfr_set_prec (x, 53);
395 mpfr_set_prec (y, 53);
396 mpfr_set_str_binary (x, "101010110100110011111010000110001000111100000110101E-43");
397 mpfr_gamma (x, x, MPFR_RNDU);
398 mpfr_set_str_binary (y, "110000011110001000111110110101011110000100001111111E971");
399 if (mpfr_cmp (x, y) != 0)
400 {
401 printf ("Error for gamma(4)\n");
402 printf ("expected "); mpfr_dump (y);
403 printf ("got "); mpfr_dump (x);
404 exit (1);
405 }
406
407 set_emin (emin);
408 set_emax (emax);
409
410 /* bug found by Kevin Rauch, 26 Oct 2007 */
411 mpfr_set_str (x, "1e19", 10, MPFR_RNDN);
412 inex = mpfr_gamma (x, x, MPFR_RNDN);
413 MPFR_ASSERTN(mpfr_inf_p (x) && inex > 0);
414
415 mpfr_clear (y);
416 mpfr_clear (x);
417 }
418
419 /* test gamma on some integral values (from Christopher Creutzig). */
420 static void
421 gamma_integer (void)
422 {
423 mpz_t n;
424 mpfr_t x, y;
425 unsigned int i;
426
427 mpz_init (n);
428 mpfr_init2 (x, 149);
429 mpfr_init2 (y, 149);
430
431 for (i = 0; i < 100; i++)
432 {
433 mpz_fac_ui (n, i);
434 mpfr_set_ui (x, i+1, MPFR_RNDN);
435 mpfr_gamma (y, x, MPFR_RNDN);
436 mpfr_set_z (x, n, MPFR_RNDN);
437 if (!mpfr_equal_p (x, y))
438 {
439 printf ("Error for gamma(%u)\n", i+1);
440 printf ("expected "); mpfr_dump (x);
441 printf ("got "); mpfr_dump (y);
442 exit (1);
443 }
444 }
445 mpfr_clear (y);
446 mpfr_clear (x);
447 mpz_clear (n);
448 }
449
450 /* bug found by Kevin Rauch */
451 static void
452 test20071231 (void)
453 {
454 mpfr_t x;
455 int inex;
456 mpfr_exp_t emin;
457
458 emin = mpfr_get_emin ();
459 set_emin (-1000000);
460
461 mpfr_init2 (x, 21);
462 mpfr_set_str (x, "-1000001.5", 10, MPFR_RNDN);
463 inex = mpfr_gamma (x, x, MPFR_RNDN);
464 MPFR_ASSERTN(MPFR_IS_ZERO(x) && MPFR_IS_POS(x) && inex < 0);
465 mpfr_clear (x);
466
467 set_emin (emin);
468
469 mpfr_init2 (x, 53);
470 mpfr_set_str (x, "-1000000001.5", 10, MPFR_RNDN);
471 inex = mpfr_gamma (x, x, MPFR_RNDN);
472 MPFR_ASSERTN(MPFR_IS_ZERO(x) && MPFR_IS_POS(x) && inex < 0);
473 mpfr_clear (x);
474 }
475
476 /* bug found by Stathis in mpfr_gamma, only occurs on 32-bit machines;
477 the second test is for 64-bit machines. This bug reappeared due to
478 r8159. */
479 static void
480 test20100709 (void)
481 {
482 mpfr_t x, y, z;
483 int sign;
484 int inex;
485 mpfr_exp_t emin;
486
487 mpfr_init2 (x, 100);
488 mpfr_init2 (y, 32);
489 mpfr_init2 (z, 32);
490 mpfr_set_str (x, "-4.6308260837372266e+07", 10, MPFR_RNDN);
491 mpfr_set_ui (y, 0, MPFR_RNDN);
492 mpfr_nextabove (y);
493 mpfr_log (y, y, MPFR_RNDD);
494 mpfr_const_log2 (z, MPFR_RNDU);
495 mpfr_sub (y, y, z, MPFR_RNDD); /* log(MIN/2) = log(MIN) - log(2) */
496 mpfr_lgamma (z, &sign, x, MPFR_RNDU);
497 MPFR_ASSERTN (sign == -1);
498 MPFR_ASSERTN (mpfr_less_p (z, y)); /* thus underflow */
499 inex = mpfr_gamma (x, x, MPFR_RNDN);
500 MPFR_ASSERTN (MPFR_IS_ZERO(x) && MPFR_IS_NEG(x) && inex > 0);
501 mpfr_clear (x);
502 mpfr_clear (y);
503 mpfr_clear (z);
504
505 /* Similar test for 64-bit machines (also valid with a 32-bit exponent,
506 but will not trigger the bug). */
507 emin = mpfr_get_emin ();
508 set_emin (MPFR_EMIN_MIN);
509 mpfr_init2 (x, 100);
510 mpfr_init2 (y, 32);
511 mpfr_init2 (z, 32);
512 mpfr_set_str (x, "-90.6308260837372266e+15", 10, MPFR_RNDN);
513 mpfr_set_ui (y, 0, MPFR_RNDN);
514 mpfr_nextabove (y);
515 mpfr_log (y, y, MPFR_RNDD);
516 mpfr_const_log2 (z, MPFR_RNDU);
517 mpfr_sub (y, y, z, MPFR_RNDD); /* log(MIN/2) = log(MIN) - log(2) */
518 mpfr_lgamma (z, &sign, x, MPFR_RNDU);
519 MPFR_ASSERTN (sign == -1);
520 MPFR_ASSERTN (mpfr_less_p (z, y)); /* thus underflow */
521 inex = mpfr_gamma (x, x, MPFR_RNDN);
522 MPFR_ASSERTN (MPFR_IS_ZERO(x) && MPFR_IS_NEG(x) && inex > 0);
523 mpfr_clear (x);
524 mpfr_clear (y);
525 mpfr_clear (z);
526 set_emin (emin);
527 }
528
529 /* bug found by Giridhar Tammana */
530 static void
531 test20120426 (void)
532 {
533 mpfr_t xa, xb;
534 int i;
535 mpfr_exp_t emin;
536
537 mpfr_init2 (xa, 53);
538 mpfr_init2 (xb, 53);
539 mpfr_set_si_2exp (xb, -337, -1, MPFR_RNDN);
540 emin = mpfr_get_emin ();
541 set_emin (-1073);
542 i = mpfr_gamma (xa, xb, MPFR_RNDN);
543 i = mpfr_subnormalize (xa, i, MPFR_RNDN); /* new ternary value */
544 mpfr_set_str (xb, "-9.5737343987585366746184749943e-304", 10, MPFR_RNDN);
545 if (!((i > 0) && (mpfr_cmp (xa, xb) == 0)))
546 {
547 printf ("Error in test20120426, i=%d\n", i);
548 printf ("expected ");
549 mpfr_dump (xb);
550 printf ("got ");
551 mpfr_dump (xa);
552 exit (1);
553 }
554 set_emin (emin);
555 mpfr_clear (xa);
556 mpfr_clear (xb);
557 }
558
559 static void
560 exprange (void)
561 {
562 mpfr_exp_t emin, emax;
563 mpfr_t x, y, z;
564 int inex1, inex2;
565 unsigned int flags1, flags2;
566
567 emin = mpfr_get_emin ();
568 emax = mpfr_get_emax ();
569
570 mpfr_init2 (x, 16);
571 mpfr_inits2 (8, y, z, (mpfr_ptr) 0);
572
573 mpfr_set_ui_2exp (x, 5, -1, MPFR_RNDN);
574 mpfr_clear_flags ();
575 inex1 = mpfr_gamma (y, x, MPFR_RNDN);
576 flags1 = __gmpfr_flags;
577 MPFR_ASSERTN (mpfr_inexflag_p ());
578 set_emin (0);
579 mpfr_clear_flags ();
580 inex2 = mpfr_gamma (z, x, MPFR_RNDN);
581 flags2 = __gmpfr_flags;
582 MPFR_ASSERTN (mpfr_inexflag_p ());
583 set_emin (emin);
584 if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z))
585 {
586 printf ("Error in exprange (test1)\n");
587 printf ("x = ");
588 mpfr_dump (x);
589 printf ("Expected inex1 = %d, flags1 = %u, ", VSIGN (inex1), flags1);
590 mpfr_dump (y);
591 printf ("Got inex2 = %d, flags2 = %u, ", VSIGN (inex2), flags2);
592 mpfr_dump (z);
593 exit (1);
594 }
595
596 mpfr_set_ui_2exp (x, 32769, -60, MPFR_RNDN);
597 mpfr_clear_flags ();
598 inex1 = mpfr_gamma (y, x, MPFR_RNDD);
599 flags1 = __gmpfr_flags;
600 MPFR_ASSERTN (mpfr_inexflag_p ());
601 set_emax (45);
602 mpfr_clear_flags ();
603 inex2 = mpfr_gamma (z, x, MPFR_RNDD);
604 flags2 = __gmpfr_flags;
605 MPFR_ASSERTN (mpfr_inexflag_p ());
606 set_emax (emax);
607 if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z))
608 {
609 printf ("Error in exprange (test2)\n");
610 printf ("x = ");
611 mpfr_dump (x);
612 printf ("Expected inex1 = %d, flags1 = %u, ", VSIGN (inex1), flags1);
613 mpfr_dump (y);
614 printf ("Got inex2 = %d, flags2 = %u, ", VSIGN (inex2), flags2);
615 mpfr_dump (z);
616 exit (1);
617 }
618
619 set_emax (44);
620 mpfr_clear_flags ();
621 inex1 = mpfr_check_range (y, inex1, MPFR_RNDD);
622 flags1 = __gmpfr_flags;
623 MPFR_ASSERTN (mpfr_inexflag_p ());
624 mpfr_clear_flags ();
625 inex2 = mpfr_gamma (z, x, MPFR_RNDD);
626 flags2 = __gmpfr_flags;
627 MPFR_ASSERTN (mpfr_inexflag_p ());
628 set_emax (emax);
629 if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z))
630 {
631 printf ("Error in exprange (test3)\n");
632 printf ("x = ");
633 mpfr_dump (x);
634 printf ("Expected inex1 = %d, flags1 = %u, ", VSIGN (inex1), flags1);
635 mpfr_dump (y);
636 printf ("Got inex2 = %d, flags2 = %u, ", VSIGN (inex2), flags2);
637 mpfr_dump (z);
638 exit (1);
639 }
640
641 mpfr_set_ui_2exp (x, 1, -60, MPFR_RNDN);
642 mpfr_clear_flags ();
643 inex1 = mpfr_gamma (y, x, MPFR_RNDD);
644 flags1 = __gmpfr_flags;
645 MPFR_ASSERTN (mpfr_inexflag_p ());
646 set_emax (60);
647 mpfr_clear_flags ();
648 inex2 = mpfr_gamma (z, x, MPFR_RNDD);
649 flags2 = __gmpfr_flags;
650 MPFR_ASSERTN (mpfr_inexflag_p ());
651 set_emax (emax);
652 if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z))
653 {
654 printf ("Error in exprange (test4)\n");
655 printf ("x = ");
656 mpfr_dump (x);
657 printf ("Expected inex1 = %d, flags1 = %u, ", VSIGN (inex1), flags1);
658 mpfr_dump (y);
659 printf ("Got inex2 = %d, flags2 = %u, ", VSIGN (inex2), flags2);
660 mpfr_dump (z);
661 exit (1);
662 }
663
664 MPFR_ASSERTN (MPFR_EMIN_MIN == - MPFR_EMAX_MAX);
665 set_emin (MPFR_EMIN_MIN);
666 set_emax (MPFR_EMAX_MAX);
667 mpfr_set_ui (x, 0, MPFR_RNDN);
668 mpfr_nextabove (x); /* x = 2^(emin - 1) */
669 mpfr_set_inf (y, 1);
670 inex1 = 1;
671 flags1 = MPFR_FLAGS_INEXACT | MPFR_FLAGS_OVERFLOW;
672 mpfr_clear_flags ();
673 /* MPFR_RNDU: overflow, infinity since 1/x = 2^(emax + 1) */
674 inex2 = mpfr_gamma (z, x, MPFR_RNDU);
675 flags2 = __gmpfr_flags;
676 MPFR_ASSERTN (mpfr_inexflag_p ());
677 if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z))
678 {
679 printf ("Error in exprange (test5)\n");
680 printf ("x = ");
681 mpfr_dump (x);
682 printf ("Expected inex1 = %d, flags1 = %u, ", VSIGN (inex1), flags1);
683 mpfr_dump (y);
684 printf ("Got inex2 = %d, flags2 = %u, ", VSIGN (inex2), flags2);
685 mpfr_dump (z);
686 exit (1);
687 }
688 mpfr_clear_flags ();
689 /* MPFR_RNDN: overflow, infinity since 1/x = 2^(emax + 1) */
690 inex2 = mpfr_gamma (z, x, MPFR_RNDN);
691 flags2 = __gmpfr_flags;
692 MPFR_ASSERTN (mpfr_inexflag_p ());
693 if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z))
694 {
695 printf ("Error in exprange (test6)\n");
696 printf ("x = ");
697 mpfr_dump (x);
698 printf ("Expected inex1 = %d, flags1 = %u, ", VSIGN (inex1), flags1);
699 mpfr_dump (y);
700 printf ("Got inex2 = %d, flags2 = %u, ", VSIGN (inex2), flags2);
701 mpfr_dump (z);
702 exit (1);
703 }
704 mpfr_nextbelow (y);
705 inex1 = -1;
706 mpfr_clear_flags ();
707 /* MPFR_RNDD: overflow, maxnum since 1/x = 2^(emax + 1) */
708 inex2 = mpfr_gamma (z, x, MPFR_RNDD);
709 flags2 = __gmpfr_flags;
710 MPFR_ASSERTN (mpfr_inexflag_p ());
711 if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z))
712 {
713 printf ("Error in exprange (test7)\n");
714 printf ("x = ");
715 mpfr_dump (x);
716 printf ("Expected inex1 = %d, flags1 = %u, ", VSIGN (inex1), flags1);
717 mpfr_dump (y);
718 printf ("Got inex2 = %d, flags2 = %u, ", VSIGN (inex2), flags2);
719 mpfr_dump (z);
720 exit (1);
721 }
722 mpfr_mul_2ui (x, x, 1, MPFR_RNDN); /* x = 2^emin */
723 mpfr_set_inf (y, 1);
724 inex1 = 1;
725 flags1 = MPFR_FLAGS_INEXACT | MPFR_FLAGS_OVERFLOW;
726 mpfr_clear_flags ();
727 /* MPFR_RNDU: overflow, infinity since 1/x = 2^emax */
728 inex2 = mpfr_gamma (z, x, MPFR_RNDU);
729 flags2 = __gmpfr_flags;
730 MPFR_ASSERTN (mpfr_inexflag_p ());
731 if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z))
732 {
733 printf ("Error in exprange (test8)\n");
734 printf ("x = ");
735 mpfr_dump (x);
736 printf ("Expected inex1 = %d, flags1 = %u, ", VSIGN (inex1), flags1);
737 mpfr_dump (y);
738 printf ("Got inex2 = %d, flags2 = %u, ", VSIGN (inex2), flags2);
739 mpfr_dump (z);
740 exit (1);
741 }
742 mpfr_clear_flags ();
743 /* MPFR_RNDN: overflow, infinity since 1/x = 2^emax */
744 inex2 = mpfr_gamma (z, x, MPFR_RNDN);
745 flags2 = __gmpfr_flags;
746 MPFR_ASSERTN (mpfr_inexflag_p ());
747 if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z))
748 {
749 printf ("Error in exprange (test9)\n");
750 printf ("x = ");
751 mpfr_dump (x);
752 printf ("Expected inex1 = %d, flags1 = %u, ", VSIGN (inex1), flags1);
753 mpfr_dump (y);
754 printf ("Got inex2 = %d, flags2 = %u, ", VSIGN (inex2), flags2);
755 mpfr_dump (z);
756 exit (1);
757 }
758 mpfr_nextbelow (y);
759 inex1 = -1;
760 flags1 = MPFR_FLAGS_INEXACT;
761 mpfr_clear_flags ();
762 /* MPFR_RNDD: no overflow, maxnum since 1/x = 2^emax and euler > 0 */
763 inex2 = mpfr_gamma (z, x, MPFR_RNDD);
764 flags2 = __gmpfr_flags;
765 MPFR_ASSERTN (mpfr_inexflag_p ());
766 if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z))
767 {
768 printf ("Error in exprange (test10)\n");
769 printf ("x = ");
770 mpfr_dump (x);
771 printf ("Expected inex1 = %d, flags1 = %u, ", VSIGN (inex1), flags1);
772 mpfr_dump (y);
773 printf ("Got inex2 = %d, flags2 = %u, ", VSIGN (inex2), flags2);
774 mpfr_dump (z);
775 exit (1);
776 }
777 set_emin (emin);
778 set_emax (emax);
779
780 mpfr_clears (x, y, z, (mpfr_ptr) 0);
781 }
782
783 static int
784 tiny_aux (int stop, mpfr_exp_t e)
785 {
786 mpfr_t x, y, z;
787 int r, s, spm, inex, err = 0;
788 int expected_dir[2][5] = { { 1, -1, 1, -1, 1 }, { 1, 1, 1, -1, -1 } };
789 mpfr_exp_t saved_emax;
790
791 saved_emax = mpfr_get_emax ();
792
793 mpfr_init2 (x, 32);
794 mpfr_inits2 (8, y, z, (mpfr_ptr) 0);
795
796 mpfr_set_ui_2exp (x, 1, e, MPFR_RNDN);
797 spm = 1;
798 for (s = 0; s < 2; s++)
799 {
800 RND_LOOP_NO_RNDF (r)
801 {
802 mpfr_rnd_t rr = (mpfr_rnd_t) r;
803 mpfr_exp_t exponent, emax;
804
805 /* Exponent of the rounded value in unbounded exponent range. */
806 exponent = expected_dir[s][r] < 0 && s == 0 ? - e : 1 - e;
807
808 for (emax = exponent - 1; emax <= exponent; emax++)
809 {
810 unsigned int flags, expected_flags = MPFR_FLAGS_INEXACT;
811 int overflow, expected_inex = expected_dir[s][r];
812
813 if (emax > MPFR_EMAX_MAX)
814 break;
815 set_emax (emax);
816
817 mpfr_clear_flags ();
818 inex = mpfr_gamma (y, x, rr);
819 flags = __gmpfr_flags;
820 mpfr_clear_flags ();
821 mpfr_set_si_2exp (z, spm, - e, MPFR_RNDU);
822 overflow = mpfr_overflow_p ();
823 /* z is 1/x - euler rounded toward +inf */
824
825 if (overflow && rr == MPFR_RNDN && s == 1)
826 expected_inex = -1;
827
828 if (expected_inex < 0)
829 mpfr_nextbelow (z); /* 1/x - euler rounded toward -inf */
830
831 if (exponent > emax)
832 expected_flags |= MPFR_FLAGS_OVERFLOW;
833
834 if (!(mpfr_equal_p (y, z) && flags == expected_flags
835 && SAME_SIGN (inex, expected_inex)))
836 {
837 printf ("Error in tiny for s = %d, r = %s, emax = %"
838 MPFR_EXP_FSPEC "d%s\n on ",
839 s, mpfr_print_rnd_mode (rr), (mpfr_eexp_t) emax,
840 exponent > emax ? " (overflow)" : "");
841 mpfr_dump (x);
842 printf (" expected inex = %2d, ", expected_inex);
843 mpfr_dump (z);
844 printf (" got inex = %2d, ", VSIGN (inex));
845 mpfr_dump (y);
846 printf (" expected flags = %u, got %u\n",
847 expected_flags, flags);
848 if (stop)
849 exit (1);
850 err = 1;
851 }
852 }
853 }
854 mpfr_neg (x, x, MPFR_RNDN);
855 spm = - spm;
856 }
857
858 mpfr_clears (x, y, z, (mpfr_ptr) 0);
859 set_emax (saved_emax);
860 return err;
861 }
862
863 static void
864 tiny (int stop)
865 {
866 mpfr_exp_t emin;
867 int err = 0;
868
869 emin = mpfr_get_emin ();
870
871 /* Note: in r7499, exponent -17 will select the generic code (in
872 tiny_aux, x has precision 32), while the other exponent values
873 will select special code for tiny values. */
874 err |= tiny_aux (stop, -17);
875 err |= tiny_aux (stop, -999);
876 err |= tiny_aux (stop, mpfr_get_emin ());
877
878 if (emin != MPFR_EMIN_MIN)
879 {
880 set_emin (MPFR_EMIN_MIN);
881 err |= tiny_aux (stop, MPFR_EMIN_MIN);
882 set_emin (emin);
883 }
884
885 if (err)
886 exit (1);
887 }
888
889 /* Test mpfr_gamma in precision p1 by comparing it with exp(lgamma(x))
890 computing with a working precision p2. Assume that x is not an
891 integer <= 2. */
892 static void
893 exp_lgamma (mpfr_ptr x, mpfr_prec_t p1, mpfr_prec_t p2)
894 {
895 mpfr_t yd, yu, zd, zu;
896 int inexd, inexu, sign;
897 int underflow = -1, overflow = -1; /* -1: we don't know */
898 int got_underflow, got_overflow;
899
900 if (mpfr_integer_p (x) && mpfr_cmp_si (x, 2) <= 0)
901 {
902 printf ("Warning! x is an integer <= 2 in exp_lgamma: ");
903 mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN); putchar ('\n');
904 return;
905 }
906 mpfr_inits2 (p2, yd, yu, (mpfr_ptr) 0);
907 inexd = mpfr_lgamma (yd, &sign, x, MPFR_RNDD);
908 mpfr_set (yu, yd, MPFR_RNDN); /* exact */
909 if (inexd)
910 mpfr_nextabove (yu);
911 mpfr_clear_flags ();
912 mpfr_exp (yd, yd, MPFR_RNDD);
913 if (! mpfr_underflow_p ())
914 underflow = 0;
915 if (mpfr_overflow_p ())
916 overflow = 1;
917 mpfr_clear_flags ();
918 mpfr_exp (yu, yu, MPFR_RNDU);
919 if (mpfr_underflow_p ())
920 underflow = 1;
921 if (! mpfr_overflow_p ())
922 overflow = 0;
923 if (sign < 0)
924 {
925 mpfr_neg (yd, yd, MPFR_RNDN); /* exact */
926 mpfr_neg (yu, yu, MPFR_RNDN); /* exact */
927 mpfr_swap (yd, yu);
928 }
929 /* yd < Gamma(x) < yu (strict inequalities since x != 1 and x != 2) */
930 mpfr_inits2 (p1, zd, zu, (mpfr_ptr) 0);
931 mpfr_clear_flags ();
932 inexd = mpfr_gamma (zd, x, MPFR_RNDD); /* zd <= Gamma(x) < yu */
933 got_underflow = underflow == -1 ? -1 : !! mpfr_underflow_p ();
934 got_overflow = overflow == -1 ? -1 : !! mpfr_overflow_p ();
935 if (! mpfr_less_p (zd, yu) || inexd > 0 ||
936 got_underflow != underflow ||
937 got_overflow != overflow)
938 {
939 printf ("Error in exp_lgamma on x = ");
940 mpfr_out_str (stdout, 16, 0, x, MPFR_RNDN); putchar ('\n');
941 printf ("yu = ");
942 mpfr_dump (yu);
943 printf ("zd = ");
944 mpfr_dump (zd);
945 printf ("got inexd = %d, expected <= 0\n", inexd);
946 printf ("got underflow = %d, expected %d\n", got_underflow, underflow);
947 printf ("got overflow = %d, expected %d\n", got_overflow, overflow);
948 exit (1);
949 }
950 mpfr_clear_flags ();
951 inexu = mpfr_gamma (zu, x, MPFR_RNDU); /* zu >= Gamma(x) > yd */
952 got_underflow = underflow == -1 ? -1 : !! mpfr_underflow_p ();
953 got_overflow = overflow == -1 ? -1 : !! mpfr_overflow_p ();
954 if (! mpfr_greater_p (zu, yd) || inexu < 0 ||
955 got_underflow != underflow ||
956 got_overflow != overflow)
957 {
958 printf ("Error in exp_lgamma on x = ");
959 mpfr_out_str (stdout, 16, 0, x, MPFR_RNDN); putchar ('\n');
960 printf ("yd = ");
961 mpfr_dump (yd);
962 printf ("zu = ");
963 mpfr_dump (zu);
964 printf ("got inexu = %d, expected >= 0\n", inexu);
965 printf ("got underflow = %d, expected %d\n", got_underflow, underflow);
966 printf ("got overflow = %d, expected %d\n", got_overflow, overflow);
967 exit (1);
968 }
969 if (mpfr_equal_p (zd, zu))
970 {
971 if (inexd != 0 || inexu != 0)
972 {
973 printf ("Error in exp_lgamma on x = ");
974 mpfr_out_str (stdout, 16, 0, x, MPFR_RNDN); putchar ('\n');
975 printf ("zd = zu, thus exact, but inexd = %d and inexu = %d\n",
976 inexd, inexu);
977 exit (1);
978 }
979 MPFR_ASSERTN (got_underflow == 0);
980 MPFR_ASSERTN (got_overflow == 0);
981 }
982 else if (inexd == 0 || inexu == 0)
983 {
984 printf ("Error in exp_lgamma on x = ");
985 mpfr_out_str (stdout, 16, 0, x, MPFR_RNDN); putchar ('\n');
986 printf ("zd != zu, thus inexact, but inexd = %d and inexu = %d\n",
987 inexd, inexu);
988 exit (1);
989 }
990 mpfr_clears (yd, yu, zd, zu, (mpfr_ptr) 0);
991 }
992
993 static void
994 exp_lgamma_tests (void)
995 {
996 mpfr_t x;
997 mpfr_exp_t emin, emax;
998 int i;
999
1000 emin = mpfr_get_emin ();
1001 emax = mpfr_get_emax ();
1002 set_emin (MPFR_EMIN_MIN);
1003 set_emax (MPFR_EMAX_MAX);
1004
1005 mpfr_init2 (x, 96);
1006 for (i = 3; i <= 8; i++)
1007 {
1008 mpfr_set_ui (x, i, MPFR_RNDN);
1009 exp_lgamma (x, 53, 64);
1010 mpfr_nextbelow (x);
1011 exp_lgamma (x, 53, 64);
1012 mpfr_nextabove (x);
1013 mpfr_nextabove (x);
1014 exp_lgamma (x, 53, 64);
1015 }
1016 mpfr_set_str (x, "1.7", 10, MPFR_RNDN);
1017 exp_lgamma (x, 53, 64);
1018 mpfr_set_str (x, "-4.6308260837372266e+07", 10, MPFR_RNDN);
1019 exp_lgamma (x, 53, 64);
1020 mpfr_set_str (x, "-90.6308260837372266e+15", 10, MPFR_RNDN);
1021 exp_lgamma (x, 53, 64);
1022 /* The following test gives a large positive result < +Inf */
1023 mpfr_set_str (x, "1.2b13fc45a92dea1@14", 16, MPFR_RNDN);
1024 exp_lgamma (x, 53, 64);
1025 /* Idem for a large negative result > -Inf */
1026 mpfr_set_str (x, "-1.2b13fc45a92de81@14", 16, MPFR_RNDN);
1027 exp_lgamma (x, 53, 64);
1028 /* The following two tests trigger an endless loop in r8186
1029 on 64-bit machines (64-bit exponent). The second one (due
1030 to undetected overflow) is a direct consequence of the
1031 first one, due to the call of Gamma(2-x) if x < 1. */
1032 mpfr_set_str (x, "1.2b13fc45a92dec8@14", 16, MPFR_RNDN);
1033 exp_lgamma (x, 53, 64);
1034 mpfr_set_str (x, "-1.2b13fc45a92dea8@14", 16, MPFR_RNDN);
1035 exp_lgamma (x, 53, 64);
1036 /* Similar tests (overflow threshold) for 32-bit machines. */
1037 mpfr_set_str (x, "2ab68d8.657542f855111c61", 16, MPFR_RNDN);
1038 exp_lgamma (x, 12, 64);
1039 mpfr_set_str (x, "-2ab68d6.657542f855111c61", 16, MPFR_RNDN);
1040 exp_lgamma (x, 12, 64);
1041 /* The following test is an overflow on 32-bit and 64-bit machines.
1042 Revision r8189 fails on 64-bit machines as the flag is unset. */
1043 mpfr_set_str (x, "1.2b13fc45a92ded8@14", 16, MPFR_RNDN);
1044 exp_lgamma (x, 53, 64);
1045 /* On the following tests, with r8196, one gets an underflow on
1046 32-bit machines, while a normal result is expected (see FIXME
1047 in gamma.c:382). */
1048 mpfr_set_str (x, "-2ab68d6.657542f855111c6104", 16, MPFR_RNDN);
1049 exp_lgamma (x, 12, 64); /* failure on 32-bit machines */
1050 mpfr_set_str (x, "-12b13fc45a92deb.1c6c5bc964", 16, MPFR_RNDN);
1051 exp_lgamma (x, 12, 64); /* failure on 64-bit machines */
1052 mpfr_clear (x);
1053
1054 set_emin (emin);
1055 set_emax (emax);
1056 }
1057
1058 /* Bug reported by Frithjof Blomquist on May 19, 2020.
1059 For the record, this bug was present since r8981
1060 (in mpfr_bernoulli_internal, den was wrongly reset to 1 in case
1061 of failure in Ziv's loop). The bug only occurred up from r8986
1062 where the initial precision was reduced, but was potentially
1063 present in any case of failure of Ziv's loop. */
1064 static void
1065 bug20200519 (void)
1066 {
1067 mpfr_prec_t prec = 25093;
1068 mpfr_t x, y, z, d;
1069 double dd;
1070 size_t min_memory_limit, old_memory_limit;
1071
1072 old_memory_limit = tests_memory_limit;
1073 min_memory_limit = 24000000;
1074 if (tests_memory_limit > 0 && tests_memory_limit < min_memory_limit)
1075 tests_memory_limit = min_memory_limit;
1076
1077 mpfr_init2 (x, prec);
1078 mpfr_init2 (y, prec);
1079 mpfr_init2 (z, prec + 100);
1080 mpfr_init2 (d, 24);
1081 mpfr_set_d (x, 2.5, MPFR_RNDN);
1082 mpfr_gamma (y, x, MPFR_RNDN);
1083 mpfr_gamma (z, x, MPFR_RNDN);
1084 mpfr_sub (d, y, z, MPFR_RNDN);
1085 mpfr_mul_2si (d, d, prec - mpfr_get_exp (y), MPFR_RNDN);
1086 dd = mpfr_get_d (d, MPFR_RNDN);
1087 if (dd < -0.5 || 0.5 < dd)
1088 {
1089 printf ("Error in bug20200519: dd=%f\n", dd);
1090 exit (1);
1091 }
1092 mpfr_clear (x);
1093 mpfr_clear (y);
1094 mpfr_clear (z);
1095 mpfr_clear (d);
1096
1097 tests_memory_limit = old_memory_limit;
1098 }
1099
1100 int
1101 main (int argc, char *argv[])
1102 {
1103 tests_start_mpfr ();
1104
1105 if (argc == 3) /* tgamma x prec: print gamma(x) to prec bits */
1106 {
1107 mpfr_prec_t p = atoi (argv[2]);
1108 mpfr_t x;
1109 mpfr_init2 (x, p);
1110 mpfr_set_str (x, argv[1], 10, MPFR_RNDN);
1111 mpfr_gamma (x, x, MPFR_RNDN);
1112 mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN);
1113 printf ("\n");
1114 mpfr_clear (x);
1115 return 0;
1116 }
1117
1118 special ();
1119 special_overflow ();
1120 exprange ();
1121 tiny (argc == 1);
1122 test_generic (MPFR_PREC_MIN, 100, 2);
1123 gamma_integer ();
1124 test20071231 ();
1125 test20100709 ();
1126 test20120426 ();
1127 exp_lgamma_tests ();
1128
1129 data_check ("data/gamma", mpfr_gamma, "mpfr_gamma");
1130
1131 /* this test takes about one minute */
1132 if (getenv ("MPFR_CHECK_EXPENSIVE") != NULL &&
1133 getenv ("MPFR_CHECK_LARGEMEM") != NULL)
1134 bug20200519 ();
1135
1136 tests_end_mpfr ();
1137 return 0;
1138 }