OLD | NEW |
| (Empty) |
1 /* Test file for mpfr_factorial. | |
2 | |
3 Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Fou
ndation, Inc. | |
4 Contributed by the Arenaire and Cacao 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 2.1 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.LIB. If not, write to | |
20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, | |
21 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 | |
30 static void | |
31 special (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, GMP_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, GMP_RNDN); | |
51 mpfr_set_prec (x, 206); | |
52 mpfr_set_str_binary (x, "0.110111100001000001101010010001000111000100000100111
00001001110001101111100110001111010100011110110110011000100110011010000100111111
0000101010000100100011100010011101110000001000010001100010000101001111E6250"); | |
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, GMP_RNDU); | |
66 | |
67 mpfr_clear (x); | |
68 mpfr_clear (y); | |
69 } | |
70 | |
71 static void | |
72 test_int (void) | |
73 { | |
74 unsigned long n0 = 1, n1 = 80, n; | |
75 mpz_t f; | |
76 mpfr_t x, y; | |
77 mp_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 < GMP_RND_MAX; r++) | |
95 { | |
96 inex1 = mpfr_fac_ui (x, n, (mp_rnd_t) r); | |
97 inex2 = mpfr_set_z (y, f, (mp_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 ((mp_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 ((mp_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 | |
120 static void | |
121 overflowed_fac0 (void) | |
122 { | |
123 mpfr_t x, y; | |
124 int inex, rnd, err = 0; | |
125 mp_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, GMP_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, (mp_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 ((mp_rnd_t) rnd)); | |
144 err = 1; | |
145 } | |
146 if (rnd == GMP_RNDZ || rnd == GMP_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 ((mp_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 ((mp_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 ((mp_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 ((mp_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 | |
191 int | |
192 main (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, GMP_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 < GMP_RND_MAX; rnd++) | |
228 { | |
229 inexact = mpfr_fac_ui (y, n, (mp_rnd_t) rnd); | |
230 err = (rnd == GMP_RNDN) ? yprec + 1 : yprec; | |
231 if (mpfr_can_round (y, err, (mp_rnd_t) rnd, (mp_rnd_t) rnd, prec)) | |
232 { | |
233 mpfr_set (t, y, (mp_rnd_t) rnd); | |
234 inexact = mpfr_fac_ui (z, n, (mp_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) <= (mp_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, GMP_RNDN); | |
261 printf (" prec=%u rnd_mode=%s\n", prec, | |
262 mpfr_print_rnd_mode ((mp_rnd_t) rnd)); | |
263 printf (" got "); | |
264 mpfr_out_str (stdout, 2, prec, z, GMP_RNDN); | |
265 puts (""); | |
266 printf (" expected "); | |
267 mpfr_out_str (stdout, 2, prec, t, GMP_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 } | |
OLD | NEW |