1/* Test file for mpfr_factorial.
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.
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#define TEST_FUNCTION mpfr_fac_ui
29
30static void
31special (void)
32{
33  mpfr_t x, y;
34  int inex;
35
36  mpfr_init (x);
37  mpfr_init (y);
38
39  mpfr_set_prec (x, 21);
40  mpfr_set_prec (y, 21);
41  mpfr_fac_ui (x, 119, MPFR_RNDZ);
42  mpfr_set_str_binary (y, "0.101111101110100110110E654");
43  if (mpfr_cmp (x, y))
44    {
45      printf ("Error in mpfr_fac_ui (119)\n");
46      exit (1);
47    }
48
49  mpfr_set_prec (y, 206);
50  inex = mpfr_fac_ui (y, 767, MPFR_RNDN);
51  mpfr_set_prec (x, 206);
52  mpfr_set_str_binary (x, "0.110111100001000001101010010001000111000100000100111000010011100011011111001100011110101000111101101100110001001100110100001001111110000101010000100100011100010011101110000001000010001100010000101001111E6250");
53  if (mpfr_cmp (x, y))
54    {
55      printf ("Error in mpfr_fac_ui (767)\n");
56      exit (1);
57    }
58  if (inex <= 0)
59    {
60      printf ("Wrong flag for mpfr_fac_ui (767)\n");
61      exit (1);
62    }
63
64  mpfr_set_prec (y, 202);
65  mpfr_fac_ui (y, 69, MPFR_RNDU);
66
67  mpfr_clear (x);
68  mpfr_clear (y);
69}
70
71static void
72test_int (void)
73{
74  unsigned long n0 = 1, n1 = 80, n;
75  mpz_t f;
76  mpfr_t x, y;
77  mpfr_prec_t prec_f, p;
78  int r;
79  int inex1, inex2;
80
81  mpz_init (f);
82  mpfr_init (x);
83  mpfr_init (y);
84
85  mpz_fac_ui (f, n0 - 1);
86  for (n = n0; n <= n1; n++)
87    {
88      mpz_mul_ui (f, f, n); /* f = n! */
89      prec_f = mpz_sizeinbase (f, 2) - mpz_scan1 (f, 0);
90      for (p = MPFR_PREC_MIN; p <= prec_f; p++)
91        {
92          mpfr_set_prec (x, p);
93          mpfr_set_prec (y, p);
94          for (r = 0; r < MPFR_RND_MAX; r++)
95            {
96              inex1 = mpfr_fac_ui (x, n, (mpfr_rnd_t) r);
97              inex2 = mpfr_set_z (y, f, (mpfr_rnd_t) r);
98              if (mpfr_cmp (x, y))
99                {
100                  printf ("Error for n=%lu prec=%lu rnd=%s\n",
101                          n, (unsigned long) p, mpfr_print_rnd_mode ((mpfr_rnd_t) r));
102                  exit (1);
103                }
104              if ((inex1 < 0 && inex2 >= 0) || (inex1 == 0 && inex2 != 0)
105                  || (inex1 > 0 && inex2 <= 0))
106                {
107                  printf ("Wrong inexact flag for n=%lu prec=%lu rnd=%s\n",
108                          n, (unsigned long) p, mpfr_print_rnd_mode ((mpfr_rnd_t) r));
109                  exit (1);
110                }
111            }
112        }
113    }
114
115  mpz_clear (f);
116  mpfr_clear (x);
117  mpfr_clear (y);
118}
119
120static void
121overflowed_fac0 (void)
122{
123  mpfr_t x, y;
124  int inex, rnd, err = 0;
125  mpfr_exp_t old_emax;
126
127  old_emax = mpfr_get_emax ();
128
129  mpfr_init2 (x, 8);
130  mpfr_init2 (y, 8);
131
132  mpfr_set_ui (y, 1, MPFR_RNDN);
133  mpfr_nextbelow (y);
134  set_emax (0);  /* 1 is not representable. */
135  RND_LOOP (rnd)
136    {
137      mpfr_clear_flags ();
138      inex = mpfr_fac_ui (x, 0, (mpfr_rnd_t) rnd);
139      if (! mpfr_overflow_p ())
140        {
141          printf ("Error in overflowed_fac0 (rnd = %s):\n"
142                  "  The overflow flag is not set.\n",
143                  mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
144          err = 1;
145        }
146      if (rnd == MPFR_RNDZ || rnd == MPFR_RNDD)
147        {
148          if (inex >= 0)
149            {
150              printf ("Error in overflowed_fac0 (rnd = %s):\n"
151                      "  The inexact value must be negative.\n",
152                      mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
153              err = 1;
154            }
155          if (! mpfr_equal_p (x, y))
156            {
157              printf ("Error in overflowed_fac0 (rnd = %s):\n"
158                      "  Got ", mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
159              mpfr_print_binary (x);
160              printf (" instead of 0.11111111E0.\n");
161              err = 1;
162            }
163        }
164      else
165        {
166          if (inex <= 0)
167            {
168              printf ("Error in overflowed_fac0 (rnd = %s):\n"
169                      "  The inexact value must be positive.\n",
170                      mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
171              err = 1;
172            }
173          if (! (mpfr_inf_p (x) && MPFR_SIGN (x) > 0))
174            {
175              printf ("Error in overflowed_fac0 (rnd = %s):\n"
176                      "  Got ", mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
177              mpfr_print_binary (x);
178              printf (" instead of +Inf.\n");
179              err = 1;
180            }
181        }
182    }
183  set_emax (old_emax);
184
185  if (err)
186    exit (1);
187  mpfr_clear (x);
188  mpfr_clear (y);
189}
190
191int
192main (int argc, char *argv[])
193{
194  unsigned int prec, err, yprec, n, k, zeros;
195  int rnd;
196  mpfr_t x, y, z, t;
197  int inexact;
198
199  tests_start_mpfr ();
200
201  special ();
202
203  test_int ();
204
205  mpfr_init (x);
206  mpfr_init (y);
207  mpfr_init (z);
208  mpfr_init (t);
209
210  mpfr_fac_ui (y, 0, MPFR_RNDN);
211
212  if (mpfr_cmp_ui (y, 1))
213    {
214      printf ("mpfr_fac_ui(0) does not give 1\n");
215      exit (1);
216    }
217
218  for (prec = 2; prec <= 100; prec++)
219    {
220      mpfr_set_prec (x, prec);
221      mpfr_set_prec (z, prec);
222      mpfr_set_prec (t, prec);
223      yprec = prec + 10;
224      mpfr_set_prec (y, yprec);
225
226      for (n = 0; n < 50; n++)
227        for (rnd = 0; rnd < MPFR_RND_MAX; rnd++)
228          {
229            inexact = mpfr_fac_ui (y, n, (mpfr_rnd_t) rnd);
230            err = (rnd == MPFR_RNDN) ? yprec + 1 : yprec;
231            if (mpfr_can_round (y, err, (mpfr_rnd_t) rnd, (mpfr_rnd_t) rnd, prec))
232              {
233                mpfr_set (t, y, (mpfr_rnd_t) rnd);
234                inexact = mpfr_fac_ui (z, n, (mpfr_rnd_t) rnd);
235                /* fact(n) ends with floor(n/2)+floor(n/4)+... zeros */
236                for (k=n/2, zeros=0; k; k >>= 1)
237                  zeros += k;
238                if (MPFR_EXP(y) <= (mpfr_exp_t) (prec + zeros))
239                  /* result should be exact */
240                  {
241                    if (inexact)
242                      {
243                        printf ("Wrong inexact flag: expected exact\n");
244                        exit (1);
245                      }
246                  }
247                else /* result is inexact */
248                  {
249                    if (!inexact)
250                      {
251                        printf ("Wrong inexact flag: expected inexact\n");
252                        printf ("n=%u prec=%u\n", n, prec);
253                        mpfr_print_binary(z); puts ("");
254                        exit (1);
255                      }
256                  }
257                if (mpfr_cmp (t, z))
258                  {
259                    printf ("results differ for x=");
260                    mpfr_out_str (stdout, 2, prec, x, MPFR_RNDN);
261                    printf (" prec=%u rnd_mode=%s\n", prec,
262                            mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
263                    printf ("   got ");
264                    mpfr_out_str (stdout, 2, prec, z, MPFR_RNDN);
265                    puts ("");
266                    printf ("   expected ");
267                    mpfr_out_str (stdout, 2, prec, t, MPFR_RNDN);
268                    puts ("");
269                    printf ("   approximation was ");
270                    mpfr_print_binary (y);
271                    puts ("");
272                    exit (1);
273                  }
274              }
275          }
276    }
277
278  mpfr_clear (x);
279  mpfr_clear (y);
280  mpfr_clear (z);
281  mpfr_clear (t);
282
283  overflowed_fac0 ();
284
285  tests_end_mpfr ();
286  return 0;
287}
288