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