1/* mpfr_tgamma -- test file for gamma function
2
3Copyright 2001, 2002, 2003, 2004, 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, and was contributed by Mathieu Dutour.
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
28/* Note: there could be an incorrect test about suspicious overflow
29   (MPFR_SUSPICIOUS_OVERFLOW) for x = 2^(-emax) = 0.5 * 2^(emin+1) in
30   RNDZ or RNDD, but this case is never tested in the generic tests. */
31#define TEST_FUNCTION mpfr_gamma
32#define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), 1, RANDS)
33#include "tgeneric.c"
34
35static void
36special (void)
37{
38  mpfr_t x, y;
39  int inex;
40
41  mpfr_init (x);
42  mpfr_init (y);
43
44  mpfr_set_nan (x);
45  mpfr_gamma (y, x, MPFR_RNDN);
46  if (!mpfr_nan_p (y))
47    {
48      printf ("Error for gamma(NaN)\n");
49      exit (1);
50    }
51
52  mpfr_set_inf (x, -1);
53  mpfr_gamma (y, x, MPFR_RNDN);
54  if (!mpfr_nan_p (y))
55    {
56      printf ("Error for gamma(-Inf)\n");
57      exit (1);
58    }
59
60  mpfr_set_inf (x, 1);
61  mpfr_gamma (y, x, MPFR_RNDN);
62  if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0)
63    {
64      printf ("Error for gamma(+Inf)\n");
65      exit (1);
66    }
67
68  mpfr_set_ui (x, 0, MPFR_RNDN);
69  mpfr_gamma (y, x, MPFR_RNDN);
70  if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0)
71    {
72      printf ("Error for gamma(+0)\n");
73      exit (1);
74    }
75
76  mpfr_set_ui (x, 0, MPFR_RNDN);
77  mpfr_neg (x, x, MPFR_RNDN);
78  mpfr_gamma (y, x, MPFR_RNDN);
79  if (!mpfr_inf_p (y) || mpfr_sgn (y) > 0)
80    {
81      printf ("Error for gamma(-0)\n");
82      exit (1);
83    }
84
85  mpfr_set_ui (x, 1, MPFR_RNDN);
86  mpfr_gamma (y, x, MPFR_RNDN);
87  if (mpfr_cmp_ui (y, 1))
88    {
89      printf ("Error for gamma(1)\n");
90      exit (1);
91    }
92
93  mpfr_set_si (x, -1, MPFR_RNDN);
94  mpfr_gamma (y, x, MPFR_RNDN);
95  if (!mpfr_nan_p (y))
96    {
97      printf ("Error for gamma(-1)\n");
98      exit (1);
99    }
100
101  mpfr_set_prec (x, 53);
102  mpfr_set_prec (y, 53);
103
104#define CHECK_X1 "1.0762904832837976166"
105#define CHECK_Y1 "0.96134843256452096050"
106
107  mpfr_set_str (x, CHECK_X1, 10, MPFR_RNDN);
108  mpfr_gamma (y, x, MPFR_RNDN);
109  mpfr_set_str (x, CHECK_Y1, 10, MPFR_RNDN);
110  if (mpfr_cmp (y, x))
111    {
112      printf ("mpfr_lngamma("CHECK_X1") is wrong:\n"
113              "expected ");
114      mpfr_print_binary (x); putchar ('\n');
115      printf ("got      ");
116      mpfr_print_binary (y); putchar ('\n');
117      exit (1);
118    }
119
120#define CHECK_X2 "9.23709516716202383435e-01"
121#define CHECK_Y2 "1.0502315560291053398"
122  mpfr_set_str (x, CHECK_X2, 10, MPFR_RNDN);
123  mpfr_gamma (y, x, MPFR_RNDN);
124  mpfr_set_str (x, CHECK_Y2, 10, MPFR_RNDN);
125  if (mpfr_cmp (y, x))
126    {
127      printf ("mpfr_lngamma("CHECK_X2") is wrong:\n"
128              "expected ");
129      mpfr_print_binary (x); putchar ('\n');
130      printf ("got      ");
131      mpfr_print_binary (y); putchar ('\n');
132      exit (1);
133    }
134
135  mpfr_set_prec (x, 8);
136  mpfr_set_prec (y, 175);
137  mpfr_set_ui (x, 33, MPFR_RNDN);
138  mpfr_gamma (y, x, MPFR_RNDU);
139  mpfr_set_prec (x, 175);
140  mpfr_set_str_binary (x, "0.110010101011010101101000010101010111000110011101001000101011000001100010111001101001011E118");
141  if (mpfr_cmp (x, y))
142    {
143      printf ("Error in mpfr_gamma (1)\n");
144      exit (1);
145    }
146
147  mpfr_set_prec (x, 21);
148  mpfr_set_prec (y, 8);
149  mpfr_set_ui (y, 120, MPFR_RNDN);
150  mpfr_gamma (x, y, MPFR_RNDZ);
151  mpfr_set_prec (y, 21);
152  mpfr_set_str_binary (y, "0.101111101110100110110E654");
153  if (mpfr_cmp (x, y))
154    {
155      printf ("Error in mpfr_gamma (120)\n");
156      printf ("Expected "); mpfr_print_binary (y); puts ("");
157      printf ("Got      "); mpfr_print_binary (x); puts ("");
158      exit (1);
159    }
160
161  mpfr_set_prec (x, 3);
162  mpfr_set_prec (y, 206);
163  mpfr_set_str_binary (x, "0.110e10");
164  inex = mpfr_gamma (y, x, MPFR_RNDN);
165  mpfr_set_prec (x, 206);
166  mpfr_set_str_binary (x, "0.110111100001000001101010010001000111000100000100111000010011100011011111001100011110101000111101101100110001001100110100001001111110000101010000100100011100010011101110000001000010001100010000101001111E6250");
167  if (mpfr_cmp (x, y))
168    {
169      printf ("Error in mpfr_gamma (768)\n");
170      exit (1);
171    }
172  if (inex <= 0)
173    {
174      printf ("Wrong flag for mpfr_gamma (768)\n");
175      exit (1);
176    }
177
178  /* worst case to exercise retry */
179  mpfr_set_prec (x, 1000);
180  mpfr_set_prec (y, 869);
181  mpfr_const_pi (x, MPFR_RNDN);
182  mpfr_gamma (y, x, MPFR_RNDN);
183
184  mpfr_set_prec (x, 4);
185  mpfr_set_prec (y, 4);
186  mpfr_set_str_binary (x, "-0.1100E-66");
187  mpfr_gamma (y, x, MPFR_RNDN);
188  mpfr_set_str_binary (x, "-0.1011E67");
189  if (mpfr_cmp (x, y))
190    {
191      printf ("Error for gamma(-0.1100E-66)\n");
192      exit (1);
193    }
194
195  mpfr_clear (x);
196  mpfr_clear (y);
197}
198
199static void
200special_overflow (void)
201{
202  mpfr_t x, y;
203  mpfr_exp_t emin, emax;
204  int inex;
205
206  emin = mpfr_get_emin ();
207  emax = mpfr_get_emax ();
208
209  set_emin (-125);
210  set_emax (128);
211
212  mpfr_init2 (x, 24);
213  mpfr_init2 (y, 24);
214  mpfr_set_str_binary (x, "0.101100100000000000110100E7");
215  mpfr_gamma (y, x, MPFR_RNDN);
216  if (!mpfr_inf_p (y))
217    {
218      printf ("Overflow error.\n");
219      mpfr_dump (y);
220      exit (1);
221    }
222
223  /* problem mentioned by Kenneth Wilder, 18 Aug 2005 */
224  mpfr_set_prec (x, 29);
225  mpfr_set_prec (y, 29);
226  mpfr_set_str (x, "-200000000.5", 10, MPFR_RNDN); /* exact */
227  mpfr_gamma (y, x, MPFR_RNDN);
228  if (!(mpfr_zero_p (y) && MPFR_SIGN (y) < 0))
229    {
230      printf ("Error for gamma(-200000000.5)\n");
231      printf ("expected -0");
232      printf ("got      ");
233      mpfr_dump (y);
234      exit (1);
235    }
236
237  mpfr_set_prec (x, 53);
238  mpfr_set_prec (y, 53);
239  mpfr_set_str (x, "-200000000.1", 10, MPFR_RNDN);
240  mpfr_gamma (y, x, MPFR_RNDN);
241  if (!(mpfr_zero_p (y) && MPFR_SIGN (y) < 0))
242    {
243      printf ("Error for gamma(-200000000.1), prec=53\n");
244      printf ("expected -0");
245      printf ("got      ");
246      mpfr_dump (y);
247      exit (1);
248    }
249
250  /* another problem mentioned by Kenneth Wilder, 29 Aug 2005 */
251  mpfr_set_prec (x, 333);
252  mpfr_set_prec (y, 14);
253  mpfr_set_str (x, "-2.0000000000000000000000005", 10, MPFR_RNDN);
254  mpfr_gamma (y, x, MPFR_RNDN);
255  mpfr_set_prec (x, 14);
256  mpfr_set_str_binary (x, "-11010011110001E66");
257  if (mpfr_cmp (x, y))
258    {
259      printf ("Error for gamma(-2.0000000000000000000000005)\n");
260      printf ("expected "); mpfr_dump (x);
261      printf ("got      "); mpfr_dump (y);
262      exit (1);
263    }
264
265  /* another tests from Kenneth Wilder, 31 Aug 2005 */
266  set_emax (200);
267  set_emin (-200);
268  mpfr_set_prec (x, 38);
269  mpfr_set_prec (y, 54);
270  mpfr_set_str_binary (x, "0.11101111011100111101001001010110101001E-166");
271  mpfr_gamma (y, x, MPFR_RNDN);
272  mpfr_set_prec (x, 54);
273  mpfr_set_str_binary (x, "0.100010001101100001110110001010111111010000100101011E167");
274  if (mpfr_cmp (x, y))
275    {
276      printf ("Error for gamma (test 1)\n");
277      printf ("expected "); mpfr_dump (x);
278      printf ("got      "); mpfr_dump (y);
279      exit (1);
280    }
281
282  set_emax (1000);
283  set_emin (-2000);
284  mpfr_set_prec (x, 38);
285  mpfr_set_prec (y, 71);
286  mpfr_set_str_binary (x, "10101011011100001111111000010111110010E-1034");
287  /* 184083777010*2^(-1034) */
288  mpfr_gamma (y, x, MPFR_RNDN);
289  mpfr_set_prec (x, 71);
290  mpfr_set_str_binary (x, "10111111001000011110010001000000000000110011110000000011101011111111100E926");
291  /* 1762885132679550982140*2^926 */
292  if (mpfr_cmp (x, y))
293    {
294      printf ("Error for gamma (test 2)\n");
295      printf ("expected "); mpfr_dump (x);
296      printf ("got      "); mpfr_dump (y);
297      exit (1);
298    }
299
300  mpfr_set_prec (x, 38);
301  mpfr_set_prec (y, 88);
302  mpfr_set_str_binary (x, "10111100111001010000100001100100100101E-104");
303  /* 202824096037*2^(-104) */
304  mpfr_gamma (y, x, MPFR_RNDN);
305  mpfr_set_prec (x, 88);
306  mpfr_set_str_binary (x, "1010110101111000111010111100010110101010100110111000001011000111000011101100001101110010E-21");
307  /* 209715199999500283894743922*2^(-21) */
308  if (mpfr_cmp (x, y))
309    {
310      printf ("Error for gamma (test 3)\n");
311      printf ("expected "); mpfr_dump (x);
312      printf ("got      "); mpfr_dump (y);
313      exit (1);
314    }
315
316  mpfr_set_prec (x, 171);
317  mpfr_set_prec (y, 38);
318  mpfr_set_str (x, "-2993155353253689176481146537402947624254601559176535", 10,
319                MPFR_RNDN);
320  mpfr_div_2exp (x, x, 170, MPFR_RNDN);
321  mpfr_gamma (y, x, MPFR_RNDN);
322  mpfr_set_prec (x, 38);
323  mpfr_set_str (x, "201948391737", 10, MPFR_RNDN);
324  mpfr_mul_2exp (x, x, 92, MPFR_RNDN);
325  if (mpfr_cmp (x, y))
326    {
327      printf ("Error for gamma (test 5)\n");
328      printf ("expected "); mpfr_dump (x);
329      printf ("got      "); mpfr_dump (y);
330      exit (1);
331    }
332
333  set_emin (-500000);
334  mpfr_set_prec (x, 337);
335  mpfr_set_prec (y, 38);
336  mpfr_set_str (x, "-30000.000000000000000000000000000000000000000000001", 10,
337                MPFR_RNDN);
338  mpfr_gamma (y, x, MPFR_RNDN);
339  mpfr_set_prec (x, 38);
340  mpfr_set_str (x, "-3.623795987425E-121243", 10, MPFR_RNDN);
341  if (mpfr_cmp (x, y))
342    {
343      printf ("Error for gamma (test 7)\n");
344      printf ("expected "); mpfr_dump (x);
345      printf ("got      "); mpfr_dump (y);
346      exit (1);
347    }
348
349  /* was producing infinite loop */
350  set_emin (emin);
351  mpfr_set_prec (x, 71);
352  mpfr_set_prec (y, 71);
353  mpfr_set_str (x, "-200000000.1", 10, MPFR_RNDN);
354  mpfr_gamma (y, x, MPFR_RNDN);
355  if (!(mpfr_zero_p (y) && MPFR_SIGN (y) < 0))
356    {
357      printf ("Error for gamma (test 8)\n");
358      printf ("expected "); mpfr_dump (x);
359      printf ("got      "); mpfr_dump (y);
360      exit (1);
361    }
362
363  set_emax (1073741823);
364  mpfr_set_prec (x, 29);
365  mpfr_set_prec (y, 29);
366  mpfr_set_str (x, "423786866", 10, MPFR_RNDN);
367  mpfr_gamma (y, x, MPFR_RNDN);
368  if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0)
369    {
370      printf ("Error for gamma(423786866)\n");
371      exit (1);
372    }
373
374  /* check exact result */
375  mpfr_set_prec (x, 2);
376  mpfr_set_ui (x, 3, MPFR_RNDN);
377  inex = mpfr_gamma (x, x, MPFR_RNDN);
378  if (inex != 0 || mpfr_cmp_ui (x, 2) != 0)
379    {
380      printf ("Error for gamma(3)\n");
381      exit (1);
382    }
383
384  mpfr_set_emax (1024);
385  mpfr_set_prec (x, 53);
386  mpfr_set_prec (y, 53);
387  mpfr_set_str_binary (x, "101010110100110011111010000110001000111100000110101E-43");
388  mpfr_gamma (x, x, MPFR_RNDU);
389  mpfr_set_str_binary (y, "110000011110001000111110110101011110000100001111111E971");
390  if (mpfr_cmp (x, y) != 0)
391    {
392      printf ("Error for gamma(4)\n");
393      printf ("expected "); mpfr_dump (y);
394      printf ("got      "); mpfr_dump (x);
395      exit (1);
396    }
397
398  set_emin (emin);
399  set_emax (emax);
400
401  /* bug found by Kevin Rauch, 26 Oct 2007 */
402  mpfr_set_str (x, "1e19", 10, MPFR_RNDN);
403  inex = mpfr_gamma (x, x, MPFR_RNDN);
404  MPFR_ASSERTN(mpfr_inf_p (x) && inex > 0);
405
406  mpfr_clear (y);
407  mpfr_clear (x);
408}
409
410/* test gamma on some integral values (from Christopher Creutzig). */
411static void
412gamma_integer (void)
413{
414  mpz_t n;
415  mpfr_t x, y;
416  unsigned int i;
417
418  mpz_init (n);
419  mpfr_init2 (x, 149);
420  mpfr_init2 (y, 149);
421
422  for (i = 0; i < 100; i++)
423    {
424      mpz_fac_ui (n, i);
425      mpfr_set_ui (x, i+1, MPFR_RNDN);
426      mpfr_gamma (y, x, MPFR_RNDN);
427      mpfr_set_z (x, n, MPFR_RNDN);
428      if (!mpfr_equal_p (x, y))
429        {
430          printf ("Error for gamma(%u)\n", i+1);
431          printf ("expected "); mpfr_dump (x);
432          printf ("got      "); mpfr_dump (y);
433          exit (1);
434        }
435    }
436  mpfr_clear (y);
437  mpfr_clear (x);
438  mpz_clear (n);
439}
440
441/* bug found by Kevin Rauch */
442static void
443test20071231 (void)
444{
445  mpfr_t x;
446  int inex;
447  mpfr_exp_t emin;
448
449  emin = mpfr_get_emin ();
450  mpfr_set_emin (-1000000);
451
452  mpfr_init2 (x, 21);
453  mpfr_set_str (x, "-1000001.5", 10, MPFR_RNDN);
454  inex = mpfr_gamma (x, x, MPFR_RNDN);
455  MPFR_ASSERTN(MPFR_IS_ZERO(x) && MPFR_IS_POS(x) && inex < 0);
456  mpfr_clear (x);
457
458  mpfr_set_emin (emin);
459
460  mpfr_init2 (x, 53);
461  mpfr_set_str (x, "-1000000001.5", 10, MPFR_RNDN);
462  inex = mpfr_gamma (x, x, MPFR_RNDN);
463  MPFR_ASSERTN(MPFR_IS_ZERO(x) && MPFR_IS_POS(x) && inex < 0);
464  mpfr_clear (x);
465}
466
467/* bug found by Stathis, only occurs on 32-bit machines */
468static void
469test20100709 (void)
470{
471  mpfr_t x;
472  int inex;
473
474  mpfr_init2 (x, 100);
475  mpfr_set_str (x, "-4.6308260837372266e+07", 10, MPFR_RNDN);
476  inex = mpfr_gamma (x, x, MPFR_RNDN);
477  MPFR_ASSERTN(MPFR_IS_ZERO(x) && MPFR_IS_NEG(x) && inex > 0);
478  mpfr_clear (x);
479}
480
481static void
482exprange (void)
483{
484  mpfr_exp_t emin, emax;
485  mpfr_t x, y, z;
486  int inex1, inex2;
487  unsigned int flags1, flags2;
488
489  emin = mpfr_get_emin ();
490  emax = mpfr_get_emax ();
491
492  mpfr_init2 (x, 16);
493  mpfr_inits2 (8, y, z, (mpfr_ptr) 0);
494
495  mpfr_set_ui_2exp (x, 5, -1, MPFR_RNDN);
496  mpfr_clear_flags ();
497  inex1 = mpfr_gamma (y, x, MPFR_RNDN);
498  flags1 = __gmpfr_flags;
499  MPFR_ASSERTN (mpfr_inexflag_p ());
500  mpfr_set_emin (0);
501  mpfr_clear_flags ();
502  inex2 = mpfr_gamma (z, x, MPFR_RNDN);
503  flags2 = __gmpfr_flags;
504  MPFR_ASSERTN (mpfr_inexflag_p ());
505  mpfr_set_emin (emin);
506  if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z))
507    {
508      printf ("Error in exprange (test1)\n");
509      printf ("x = ");
510      mpfr_dump (x);
511      printf ("Expected inex1 = %d, flags1 = %u, ", SIGN (inex1), flags1);
512      mpfr_dump (y);
513      printf ("Got      inex2 = %d, flags2 = %u, ", SIGN (inex2), flags2);
514      mpfr_dump (z);
515      exit (1);
516    }
517
518  mpfr_set_ui_2exp (x, 32769, -60, MPFR_RNDN);
519  mpfr_clear_flags ();
520  inex1 = mpfr_gamma (y, x, MPFR_RNDD);
521  flags1 = __gmpfr_flags;
522  MPFR_ASSERTN (mpfr_inexflag_p ());
523  mpfr_set_emax (45);
524  mpfr_clear_flags ();
525  inex2 = mpfr_gamma (z, x, MPFR_RNDD);
526  flags2 = __gmpfr_flags;
527  MPFR_ASSERTN (mpfr_inexflag_p ());
528  mpfr_set_emax (emax);
529  if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z))
530    {
531      printf ("Error in exprange (test2)\n");
532      printf ("x = ");
533      mpfr_dump (x);
534      printf ("Expected inex1 = %d, flags1 = %u, ", SIGN (inex1), flags1);
535      mpfr_dump (y);
536      printf ("Got      inex2 = %d, flags2 = %u, ", SIGN (inex2), flags2);
537      mpfr_dump (z);
538      exit (1);
539    }
540
541  mpfr_set_emax (44);
542  mpfr_clear_flags ();
543  inex1 = mpfr_check_range (y, inex1, MPFR_RNDD);
544  flags1 = __gmpfr_flags;
545  MPFR_ASSERTN (mpfr_inexflag_p ());
546  mpfr_clear_flags ();
547  inex2 = mpfr_gamma (z, x, MPFR_RNDD);
548  flags2 = __gmpfr_flags;
549  MPFR_ASSERTN (mpfr_inexflag_p ());
550  mpfr_set_emax (emax);
551  if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z))
552    {
553      printf ("Error in exprange (test3)\n");
554      printf ("x = ");
555      mpfr_dump (x);
556      printf ("Expected inex1 = %d, flags1 = %u, ", SIGN (inex1), flags1);
557      mpfr_dump (y);
558      printf ("Got      inex2 = %d, flags2 = %u, ", SIGN (inex2), flags2);
559      mpfr_dump (z);
560      exit (1);
561    }
562
563  mpfr_set_ui_2exp (x, 1, -60, MPFR_RNDN);
564  mpfr_clear_flags ();
565  inex1 = mpfr_gamma (y, x, MPFR_RNDD);
566  flags1 = __gmpfr_flags;
567  MPFR_ASSERTN (mpfr_inexflag_p ());
568  mpfr_set_emax (60);
569  mpfr_clear_flags ();
570  inex2 = mpfr_gamma (z, x, MPFR_RNDD);
571  flags2 = __gmpfr_flags;
572  MPFR_ASSERTN (mpfr_inexflag_p ());
573  mpfr_set_emax (emax);
574  if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z))
575    {
576      printf ("Error in exprange (test4)\n");
577      printf ("x = ");
578      mpfr_dump (x);
579      printf ("Expected inex1 = %d, flags1 = %u, ", SIGN (inex1), flags1);
580      mpfr_dump (y);
581      printf ("Got      inex2 = %d, flags2 = %u, ", SIGN (inex2), flags2);
582      mpfr_dump (z);
583      exit (1);
584    }
585
586  MPFR_ASSERTN (MPFR_EMIN_MIN == - MPFR_EMAX_MAX);
587  mpfr_set_emin (MPFR_EMIN_MIN);
588  mpfr_set_emax (MPFR_EMAX_MAX);
589  mpfr_set_ui (x, 0, MPFR_RNDN);
590  mpfr_nextabove (x);  /* x = 2^(emin - 1) */
591  mpfr_set_inf (y, 1);
592  inex1 = 1;
593  flags1 = MPFR_FLAGS_INEXACT | MPFR_FLAGS_OVERFLOW;
594  mpfr_clear_flags ();
595  /* MPFR_RNDU: overflow, infinity since 1/x = 2^(emax + 1) */
596  inex2 = mpfr_gamma (z, x, MPFR_RNDU);
597  flags2 = __gmpfr_flags;
598  MPFR_ASSERTN (mpfr_inexflag_p ());
599  if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z))
600    {
601      printf ("Error in exprange (test5)\n");
602      printf ("x = ");
603      mpfr_dump (x);
604      printf ("Expected inex1 = %d, flags1 = %u, ", SIGN (inex1), flags1);
605      mpfr_dump (y);
606      printf ("Got      inex2 = %d, flags2 = %u, ", SIGN (inex2), flags2);
607      mpfr_dump (z);
608      exit (1);
609    }
610  mpfr_clear_flags ();
611  /* MPFR_RNDN: overflow, infinity since 1/x = 2^(emax + 1) */
612  inex2 = mpfr_gamma (z, x, MPFR_RNDN);
613  flags2 = __gmpfr_flags;
614  MPFR_ASSERTN (mpfr_inexflag_p ());
615  if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z))
616    {
617      printf ("Error in exprange (test6)\n");
618      printf ("x = ");
619      mpfr_dump (x);
620      printf ("Expected inex1 = %d, flags1 = %u, ", SIGN (inex1), flags1);
621      mpfr_dump (y);
622      printf ("Got      inex2 = %d, flags2 = %u, ", SIGN (inex2), flags2);
623      mpfr_dump (z);
624      exit (1);
625    }
626  mpfr_nextbelow (y);
627  inex1 = -1;
628  mpfr_clear_flags ();
629  /* MPFR_RNDD: overflow, maxnum since 1/x = 2^(emax + 1) */
630  inex2 = mpfr_gamma (z, x, MPFR_RNDD);
631  flags2 = __gmpfr_flags;
632  MPFR_ASSERTN (mpfr_inexflag_p ());
633  if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z))
634    {
635      printf ("Error in exprange (test7)\n");
636      printf ("x = ");
637      mpfr_dump (x);
638      printf ("Expected inex1 = %d, flags1 = %u, ", SIGN (inex1), flags1);
639      mpfr_dump (y);
640      printf ("Got      inex2 = %d, flags2 = %u, ", SIGN (inex2), flags2);
641      mpfr_dump (z);
642      exit (1);
643    }
644  mpfr_mul_2ui (x, x, 1, MPFR_RNDN);  /* x = 2^emin */
645  mpfr_set_inf (y, 1);
646  inex1 = 1;
647  flags1 = MPFR_FLAGS_INEXACT | MPFR_FLAGS_OVERFLOW;
648  mpfr_clear_flags ();
649  /* MPFR_RNDU: overflow, infinity since 1/x = 2^emax */
650  inex2 = mpfr_gamma (z, x, MPFR_RNDU);
651  flags2 = __gmpfr_flags;
652  MPFR_ASSERTN (mpfr_inexflag_p ());
653  if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z))
654    {
655      printf ("Error in exprange (test8)\n");
656      printf ("x = ");
657      mpfr_dump (x);
658      printf ("Expected inex1 = %d, flags1 = %u, ", SIGN (inex1), flags1);
659      mpfr_dump (y);
660      printf ("Got      inex2 = %d, flags2 = %u, ", SIGN (inex2), flags2);
661      mpfr_dump (z);
662      exit (1);
663    }
664  mpfr_clear_flags ();
665  /* MPFR_RNDN: overflow, infinity since 1/x = 2^emax */
666  inex2 = mpfr_gamma (z, x, MPFR_RNDN);
667  flags2 = __gmpfr_flags;
668  MPFR_ASSERTN (mpfr_inexflag_p ());
669  if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z))
670    {
671      printf ("Error in exprange (test9)\n");
672      printf ("x = ");
673      mpfr_dump (x);
674      printf ("Expected inex1 = %d, flags1 = %u, ", SIGN (inex1), flags1);
675      mpfr_dump (y);
676      printf ("Got      inex2 = %d, flags2 = %u, ", SIGN (inex2), flags2);
677      mpfr_dump (z);
678      exit (1);
679    }
680  mpfr_nextbelow (y);
681  inex1 = -1;
682  flags1 = MPFR_FLAGS_INEXACT;
683  mpfr_clear_flags ();
684  /* MPFR_RNDD: no overflow, maxnum since 1/x = 2^emax and euler > 0 */
685  inex2 = mpfr_gamma (z, x, MPFR_RNDD);
686  flags2 = __gmpfr_flags;
687  MPFR_ASSERTN (mpfr_inexflag_p ());
688  if (inex1 != inex2 || flags1 != flags2 || ! mpfr_equal_p (y, z))
689    {
690      printf ("Error in exprange (test10)\n");
691      printf ("x = ");
692      mpfr_dump (x);
693      printf ("Expected inex1 = %d, flags1 = %u, ", SIGN (inex1), flags1);
694      mpfr_dump (y);
695      printf ("Got      inex2 = %d, flags2 = %u, ", SIGN (inex2), flags2);
696      mpfr_dump (z);
697      exit (1);
698    }
699  mpfr_set_emin (emin);
700  mpfr_set_emax (emax);
701
702  mpfr_clears (x, y, z, (mpfr_ptr) 0);
703}
704
705static int
706tiny_aux (int stop, mpfr_exp_t e)
707{
708  mpfr_t x, y, z;
709  int r, s, spm, inex, err = 0;
710  int expected_dir[2][5] = { { 1, -1, 1, -1, 1 }, { 1, 1, 1, -1, -1 } };
711  mpfr_exp_t saved_emax;
712
713  saved_emax = mpfr_get_emax ();
714
715  mpfr_init2 (x, 32);
716  mpfr_inits2 (8, y, z, (mpfr_ptr) 0);
717
718  mpfr_set_ui_2exp (x, 1, e, MPFR_RNDN);
719  spm = 1;
720  for (s = 0; s < 2; s++)
721    {
722      RND_LOOP(r)
723        {
724          mpfr_rnd_t rr = (mpfr_rnd_t) r;
725          mpfr_exp_t exponent, emax;
726
727          /* Exponent of the rounded value in unbounded exponent range. */
728          exponent = expected_dir[s][r] < 0 && s == 0 ? - e : 1 - e;
729
730          for (emax = exponent - 1; emax <= exponent; emax++)
731            {
732              unsigned int flags, expected_flags = MPFR_FLAGS_INEXACT;
733              int overflow, expected_inex = expected_dir[s][r];
734
735              if (emax > MPFR_EMAX_MAX)
736                break;
737              mpfr_set_emax (emax);
738
739              mpfr_clear_flags ();
740              inex = mpfr_gamma (y, x, rr);
741              flags = __gmpfr_flags;
742              mpfr_clear_flags ();
743              mpfr_set_si_2exp (z, spm, - e, MPFR_RNDU);
744              overflow = mpfr_overflow_p ();
745              /* z is 1/x - euler rounded toward +inf */
746
747              if (overflow && rr == MPFR_RNDN && s == 1)
748                expected_inex = -1;
749
750              if (expected_inex < 0)
751                mpfr_nextbelow (z); /* 1/x - euler rounded toward -inf */
752
753              if (exponent > emax)
754                expected_flags |= MPFR_FLAGS_OVERFLOW;
755
756              if (!(mpfr_equal_p (y, z) && flags == expected_flags
757                    && SAME_SIGN (inex, expected_inex)))
758                {
759                  printf ("Error in tiny for s = %d, r = %s, emax = %"
760                          MPFR_EXP_FSPEC "d%s\n  on ",
761                          s, mpfr_print_rnd_mode (rr), emax,
762                          exponent > emax ? " (overflow)" : "");
763                  mpfr_dump (x);
764                  printf ("  expected inex = %2d, ", expected_inex);
765                  mpfr_dump (z);
766                  printf ("  got      inex = %2d, ", SIGN (inex));
767                  mpfr_dump (y);
768                  printf ("  expected flags = %u, got %u\n",
769                          expected_flags, flags);
770                  if (stop)
771                    exit (1);
772                  err = 1;
773                }
774            }
775        }
776      mpfr_neg (x, x, MPFR_RNDN);
777      spm = - spm;
778    }
779
780  mpfr_clears (x, y, z, (mpfr_ptr) 0);
781  mpfr_set_emax (saved_emax);
782  return err;
783}
784
785static void
786tiny (int stop)
787{
788  mpfr_exp_t emin;
789  int err = 0;
790
791  emin = mpfr_get_emin ();
792
793  /* Note: in r7499, exponent -17 will select the generic code (in
794     tiny_aux, x has precision 32), while the other exponent values
795     will select special code for tiny values. */
796  err |= tiny_aux (stop, -17);
797  err |= tiny_aux (stop, -999);
798  err |= tiny_aux (stop, mpfr_get_emin ());
799
800  if (emin != MPFR_EMIN_MIN)
801    {
802      mpfr_set_emin (MPFR_EMIN_MIN);
803      err |= tiny_aux (stop, MPFR_EMIN_MIN);
804      mpfr_set_emin (emin);
805    }
806
807  if (err)
808    exit (1);
809}
810
811int
812main (int argc, char *argv[])
813{
814  tests_start_mpfr ();
815
816  special ();
817  special_overflow ();
818  exprange ();
819  tiny (argc == 1);
820  test_generic (2, 100, 2);
821  gamma_integer ();
822  test20071231 ();
823  test20100709 ();
824
825  data_check ("data/gamma", mpfr_gamma, "mpfr_gamma");
826
827  tests_end_mpfr ();
828  return 0;
829}
830