OLD | NEW |
| (Empty) |
1 /* Test file for mpfr_pow_z -- power function x^z with z a MPZ | |
2 | |
3 Copyright 2005, 2006, 2007, 2008, 2009 Free Software Foundation, 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 <stdlib.h> | |
24 #include <float.h> | |
25 #include <math.h> | |
26 | |
27 #include "mpfr-test.h" | |
28 | |
29 #define ERROR(str) do { printf("Error for "str"\n"); exit (1); } while (0) | |
30 | |
31 static void | |
32 check_special (void) | |
33 { | |
34 mpfr_t x, y; | |
35 mpz_t z; | |
36 int res; | |
37 | |
38 mpfr_init (x); | |
39 mpfr_init (y); | |
40 mpz_init (z); | |
41 | |
42 /* x^0 = 1 except for NAN */ | |
43 mpfr_set_ui (x, 23, GMP_RNDN); | |
44 mpz_set_ui (z, 0); | |
45 res = mpfr_pow_z (y, x, z, GMP_RNDN); | |
46 if (res != 0 || mpfr_cmp_ui (y, 1) != 0) | |
47 ERROR ("23^0"); | |
48 mpfr_set_nan (x); | |
49 res = mpfr_pow_z (y, x, z, GMP_RNDN); | |
50 if (res != 0 || mpfr_nan_p (y) || mpfr_cmp_si (y, 1) != 0) | |
51 ERROR ("NAN^0"); | |
52 mpfr_set_inf (x, 1); | |
53 res = mpfr_pow_z (y, x, z, GMP_RNDN); | |
54 if (res != 0 || mpfr_cmp_ui (y, 1) != 0) | |
55 ERROR ("INF^0"); | |
56 | |
57 /* sINF^N = INF if s==1 or n even if N > 0*/ | |
58 mpz_set_ui (z, 42); | |
59 mpfr_set_inf (x, 1); | |
60 res = mpfr_pow_z (y, x, z, GMP_RNDN); | |
61 if (res != 0 || mpfr_inf_p (y) == 0 || mpfr_sgn (y) <= 0) | |
62 ERROR ("INF^42"); | |
63 mpfr_set_inf (x, -1); | |
64 res = mpfr_pow_z (y, x, z, GMP_RNDN); | |
65 if (res != 0 || mpfr_inf_p (y) == 0 || mpfr_sgn (y) <= 0) | |
66 ERROR ("-INF^42"); | |
67 mpz_set_ui (z, 17); | |
68 mpfr_set_inf (x, 1); | |
69 res = mpfr_pow_z (y, x, z, GMP_RNDN); | |
70 if (res != 0 || mpfr_inf_p (y) == 0 || mpfr_sgn (y) <= 0) | |
71 ERROR ("INF^17"); | |
72 mpfr_set_inf (x, -1); | |
73 res = mpfr_pow_z (y, x, z, GMP_RNDN); | |
74 if (res != 0 || mpfr_inf_p (y) == 0 || mpfr_sgn (y) >= 0) | |
75 ERROR ("-INF^17"); | |
76 | |
77 mpz_set_si (z, -42); | |
78 mpfr_set_inf (x, 1); | |
79 res = mpfr_pow_z (y, x, z, GMP_RNDN); | |
80 if (res != 0 || mpfr_zero_p (y) == 0 || MPFR_SIGN (y) <= 0) | |
81 ERROR ("INF^-42"); | |
82 mpfr_set_inf (x, -1); | |
83 res = mpfr_pow_z (y, x, z, GMP_RNDN); | |
84 if (res != 0 || mpfr_zero_p (y) == 0 || MPFR_SIGN (y) <= 0) | |
85 ERROR ("-INF^-42"); | |
86 mpz_set_si (z, -17); | |
87 mpfr_set_inf (x, 1); | |
88 res = mpfr_pow_z (y, x, z, GMP_RNDN); | |
89 if (res != 0 || mpfr_zero_p (y) == 0 || MPFR_SIGN (y) <= 0) | |
90 ERROR ("INF^-17"); | |
91 mpfr_set_inf (x, -1); | |
92 res = mpfr_pow_z (y, x, z, GMP_RNDN); | |
93 if (res != 0 || mpfr_zero_p (y) == 0 || MPFR_SIGN (y) >= 0) | |
94 ERROR ("-INF^-17"); | |
95 | |
96 /* s0^N = +0 if s==+ or n even if N > 0*/ | |
97 mpz_set_ui (z, 42); | |
98 MPFR_SET_ZERO (x); MPFR_SET_POS (x); | |
99 res = mpfr_pow_z (y, x, z, GMP_RNDN); | |
100 if (res != 0 || mpfr_zero_p (y) == 0 || MPFR_SIGN (y) <= 0) | |
101 ERROR ("+0^42"); | |
102 MPFR_SET_NEG (x); | |
103 res = mpfr_pow_z (y, x, z, GMP_RNDN); | |
104 if (res != 0 || mpfr_zero_p (y) == 0 || MPFR_SIGN (y) <= 0) | |
105 ERROR ("-0^42"); | |
106 mpz_set_ui (z, 17); | |
107 MPFR_SET_POS (x); | |
108 res = mpfr_pow_z (y, x, z, GMP_RNDN); | |
109 if (res != 0 || mpfr_zero_p (y) == 0 || MPFR_SIGN (y) <= 0) | |
110 ERROR ("+0^17"); | |
111 MPFR_SET_NEG (x); | |
112 res = mpfr_pow_z (y, x, z, GMP_RNDN); | |
113 if (res != 0 || mpfr_zero_p (y) == 0 || MPFR_SIGN (y) >= 0) | |
114 ERROR ("-0^17"); | |
115 | |
116 mpz_set_si (z, -42); | |
117 MPFR_SET_ZERO (x); MPFR_SET_POS (x); | |
118 res = mpfr_pow_z (y, x, z, GMP_RNDN); | |
119 if (res != 0 || mpfr_inf_p (y) == 0 || MPFR_SIGN (y) <= 0) | |
120 ERROR ("+0^-42"); | |
121 MPFR_SET_NEG (x); | |
122 res = mpfr_pow_z (y, x, z, GMP_RNDN); | |
123 if (res != 0 || mpfr_inf_p (y) == 0 || MPFR_SIGN (y) <= 0) | |
124 ERROR ("-0^-42"); | |
125 mpz_set_si (z, -17); | |
126 MPFR_SET_POS (x); | |
127 res = mpfr_pow_z (y, x, z, GMP_RNDN); | |
128 if (res != 0 || mpfr_inf_p (y) == 0 || MPFR_SIGN (y) <= 0) | |
129 ERROR ("+0^-17"); | |
130 MPFR_SET_NEG (x); | |
131 res = mpfr_pow_z (y, x, z, GMP_RNDN); | |
132 if (res != 0 || mpfr_inf_p (y) == 0 || MPFR_SIGN (y) >= 0) | |
133 ERROR ("-0^-17"); | |
134 | |
135 mpz_clear (z); | |
136 mpfr_clear (y); | |
137 mpfr_clear (x); | |
138 } | |
139 | |
140 static void | |
141 check_integer (mp_prec_t begin, mp_prec_t end, unsigned long max) | |
142 { | |
143 mpfr_t x, y1, y2; | |
144 mpz_t z; | |
145 unsigned long i, n; | |
146 mp_prec_t p; | |
147 int res1, res2; | |
148 mp_rnd_t rnd; | |
149 | |
150 mpfr_inits2 (begin, x, y1, y2, (mpfr_ptr) 0); | |
151 mpz_init (z); | |
152 for (p = begin ; p < end ; p+=4) | |
153 { | |
154 mpfr_set_prec (x, p); | |
155 mpfr_set_prec (y1, p); | |
156 mpfr_set_prec (y2, p); | |
157 for (i = 0 ; i < max ; i++) | |
158 { | |
159 mpz_urandomb (z, RANDS, GMP_NUMB_BITS); | |
160 if ((i & 1) != 0) | |
161 mpz_neg (z, z); | |
162 mpfr_urandomb (x, RANDS); | |
163 mpfr_mul_2ui (x, x, 1, GMP_RNDN); /* 0 <= x < 2 */ | |
164 rnd = RND_RAND (); | |
165 if (mpz_fits_slong_p (z)) | |
166 { | |
167 n = mpz_get_si (z); | |
168 /* printf ("New test for x=%ld\nCheck Pow_si\n", n); */ | |
169 res1 = mpfr_pow_si (y1, x, n, rnd); | |
170 /* printf ("Check pow_z\n"); */ | |
171 res2 = mpfr_pow_z (y2, x, z, rnd); | |
172 if (mpfr_cmp (y1, y2) != 0) | |
173 { | |
174 printf ("Error for p = %lu, z = %lu, rnd = %s and x = ", | |
175 p, n, mpfr_print_rnd_mode (rnd)); | |
176 mpfr_dump (x); | |
177 printf ("Ypowsi = "); mpfr_dump (y1); | |
178 printf ("Ypowz = "); mpfr_dump (y2); | |
179 exit (1); | |
180 } | |
181 if (res1 != res2) | |
182 { | |
183 printf ("Wrong inexact flags for p = %lu, z = %lu, rnd = %s" | |
184 " and x = ", p, n, mpfr_print_rnd_mode (rnd)); | |
185 mpfr_dump (x); | |
186 printf ("Ypowsi(inex = %2d) = ", res1); mpfr_dump (y1); | |
187 printf ("Ypowz (inex = %2d) = ", res2); mpfr_dump (y2); | |
188 exit (1); | |
189 } | |
190 } | |
191 } /* for i */ | |
192 } /* for p */ | |
193 mpfr_clears (x, y1, y2, (mpfr_ptr) 0); | |
194 mpz_clear (z); | |
195 } | |
196 | |
197 static void | |
198 check_regression (void) | |
199 { | |
200 mpfr_t x, y; | |
201 mpz_t z; | |
202 int res1, res2; | |
203 | |
204 mpz_init_set_ui (z, 2026876995); | |
205 mpfr_init2 (x, 122); | |
206 mpfr_init2 (y, 122); | |
207 | |
208 mpfr_set_str_binary (x, "0.100000100100001111010011101001011010100111100111000
01111000001001101000110011001001001001011001011010110110110101000111011E1"); | |
209 res1 = mpfr_pow_z (y, x, z, GMP_RNDU); | |
210 res2 = mpfr_pow_ui (x, x, 2026876995UL, GMP_RNDU); | |
211 if (mpfr_cmp (x, y) || res1 != res2) | |
212 { | |
213 printf ("Regression (1) tested failed (%d=?%d)\n",res1, res2); | |
214 printf ("pow_ui: "); mpfr_dump (x); | |
215 printf ("pow_z: "); mpfr_dump (y); | |
216 | |
217 exit (1); | |
218 } | |
219 | |
220 mpfr_clear (x); | |
221 mpfr_clear (y); | |
222 mpz_clear (z); | |
223 } | |
224 | |
225 /* Bug found by Kevin P. Rauch */ | |
226 static void | |
227 bug20071104 (void) | |
228 { | |
229 mpfr_t x, y; | |
230 mpz_t z; | |
231 int inex; | |
232 | |
233 mpz_init_set_si (z, -2); | |
234 mpfr_inits2 (20, x, y, (mpfr_ptr) 0); | |
235 mpfr_set_ui (x, 0, GMP_RNDN); | |
236 mpfr_nextbelow (x); /* x = -2^(emin-1) */ | |
237 mpfr_clear_flags (); | |
238 inex = mpfr_pow_z (y, x, z, GMP_RNDN); | |
239 if (! mpfr_inf_p (y) || MPFR_SIGN (y) < 0) | |
240 { | |
241 printf ("Error in bug20071104: expected +Inf, got "); | |
242 mpfr_dump (y); | |
243 exit (1); | |
244 } | |
245 if (inex <= 0) | |
246 { | |
247 printf ("Error in bug20071104: bad ternary value (%d)\n", inex); | |
248 exit (1); | |
249 } | |
250 if (__gmpfr_flags != (MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT)) | |
251 { | |
252 printf ("Error in bug20071104: bad flags (%u)\n", __gmpfr_flags); | |
253 exit (1); | |
254 } | |
255 mpfr_clears (x, y, (mpfr_ptr) 0); | |
256 mpz_clear (z); | |
257 } | |
258 | |
259 static void | |
260 check_overflow (void) | |
261 { | |
262 mpfr_t a; | |
263 mpz_t z; | |
264 unsigned long n; | |
265 int res; | |
266 | |
267 mpfr_init2 (a, 53); | |
268 | |
269 mpfr_set_str_binary (a, "1E10"); | |
270 mpz_init_set_ui (z, ULONG_MAX); | |
271 res = mpfr_pow_z (a, a, z, GMP_RNDN); | |
272 if (!MPFR_IS_INF (a) || MPFR_SIGN (a) < 0) | |
273 { | |
274 printf ("Error for (1e10)^ULONG_MAX\n"); | |
275 exit (1); | |
276 } | |
277 | |
278 /* Bug in pow_z.c up to r5109: if x = y (same mpfr_t argument), the | |
279 input argument is negative and not a power of two, z is positive | |
280 and odd, an overflow or underflow occurs, and the temporary result | |
281 res is positive, then the result gets a wrong sign (positive | |
282 instead of negative). */ | |
283 mpfr_set_str_binary (a, "-1.1E10"); | |
284 n = (ULONG_MAX ^ (ULONG_MAX >> 1)) + 1; | |
285 mpz_set_ui (z, n); | |
286 res = mpfr_pow_z (a, a, z, GMP_RNDN); | |
287 if (!MPFR_IS_INF (a) || MPFR_SIGN (a) > 0) | |
288 { | |
289 printf ("Error for (-1e10)^%lu, expected -Inf,\ngot ", n); | |
290 mpfr_dump (a); | |
291 exit (1); | |
292 } | |
293 | |
294 mpfr_clear (a); | |
295 mpz_clear (z); | |
296 } | |
297 | |
298 /* bug reported by Carl Witty (32-bit architecture) */ | |
299 static void | |
300 bug20080223 (void) | |
301 { | |
302 mpfr_t a, exp, answer; | |
303 | |
304 mpfr_init2 (a, 53); | |
305 mpfr_init2 (exp, 53); | |
306 mpfr_init2 (answer, 53); | |
307 | |
308 mpfr_set_si (exp, -1073741824, GMP_RNDN); | |
309 | |
310 mpfr_set_str (a, "1.999999999", 10, GMP_RNDN); | |
311 /* a = 562949953139837/2^48 */ | |
312 mpfr_pow (answer, a, exp, GMP_RNDN); | |
313 mpfr_set_str_binary (a, "0.110110101111011001110000111111100011101000111011101
E-1073741823"); | |
314 MPFR_ASSERTN(mpfr_cmp0 (answer, a) == 0); | |
315 | |
316 mpfr_clear (a); | |
317 mpfr_clear (exp); | |
318 mpfr_clear (answer); | |
319 } | |
320 | |
321 static void | |
322 bug20080904 (void) | |
323 { | |
324 mpz_t exp; | |
325 mpfr_t a, answer; | |
326 mp_exp_t emin_default; | |
327 | |
328 mpz_init (exp); | |
329 mpfr_init2 (a, 70); | |
330 mpfr_init2 (answer, 70); | |
331 | |
332 emin_default = mpfr_get_emin (); | |
333 mpfr_set_emin (MPFR_EMIN_MIN); | |
334 | |
335 mpz_set_str (exp, "-4eb92f8c7b7bf81e", 16); | |
336 mpfr_set_str_binary (a, "1.110000101110100110100011111000011110111101000011111
001111001010011100"); | |
337 | |
338 mpfr_pow_z (answer, a, exp, GMP_RNDN); | |
339 /* The correct result is near 2^(-2^62), so it underflows when | |
340 MPFR_EMIN_MIN > -2^62 (i.e. with 32 and 64 bits machines). */ | |
341 mpfr_set_str (a, "AA500C0D7A69275DBp-4632850503556296886", 16, GMP_RNDN); | |
342 MPFR_ASSERTN(mpfr_cmp0 (answer, a) == 0); | |
343 | |
344 mpfr_set_emin (emin_default); | |
345 | |
346 mpz_clear (exp); | |
347 mpfr_clear (a); | |
348 mpfr_clear (answer); | |
349 } | |
350 | |
351 int | |
352 main (void) | |
353 { | |
354 tests_start_mpfr (); | |
355 | |
356 check_special (); | |
357 | |
358 check_integer (2, 163, 100); | |
359 check_regression (); | |
360 bug20071104 (); | |
361 bug20080223 (); | |
362 bug20080904 (); | |
363 check_overflow (); | |
364 | |
365 tests_end_mpfr (); | |
366 return 0; | |
367 } | |
OLD | NEW |