1 /* Test file for mpfr_pow, mpfr_pow_ui and mpfr_pow_si.
2
3 Copyright 2000-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 #define _MPFR_NO_DEPRECATED_ROOT
24 #define MPFR_NEED_INTMAX_H
25 #include "mpfr-test.h"
26
27 #ifdef CHECK_EXTERNAL
28 static int
test_pow(mpfr_ptr a,mpfr_srcptr b,mpfr_srcptr c,mpfr_rnd_t rnd_mode)29 test_pow (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
30 {
31 int res;
32 int ok = rnd_mode == MPFR_RNDN && mpfr_number_p (b) && mpfr_number_p (c)
33 && mpfr_get_prec (a) >= 53;
34 if (ok)
35 {
36 mpfr_print_raw (b);
37 printf (" ");
38 mpfr_print_raw (c);
39 }
40 res = mpfr_pow (a, b, c, rnd_mode);
41 if (ok)
42 {
43 printf (" ");
44 mpfr_print_raw (a);
45 printf ("\n");
46 }
47 return res;
48 }
49 #else
50 #define test_pow mpfr_pow
51 #endif
52
53 #define TEST_FUNCTION test_pow
54 #define TWO_ARGS
55 #define TEST_RANDOM_POS 16 /* the 2nd argument is negative with prob. 16/512 */
56 #define TGENERIC_NOWARNING 1
57 #include "tgeneric.c"
58
59 #define TEST_FUNCTION mpfr_pow_ui
60 #define INTEGER_TYPE unsigned long
61 #define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), 1, RANDS)
62 #define INT_RAND_FUNCTION() \
63 (randlimb () % 16 == 0 ? randulong () : (unsigned long) (randlimb () % 32))
64 #include "tgeneric_ui.c"
65
66 #define TEST_FUNCTION mpfr_pow_si
67 #define INTEGER_TYPE long
68 #define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), 1, RANDS)
69 #define INT_RAND_FUNCTION() \
70 (randlimb () % 16 == 0 ? randlong () : (long) (randlimb () % 31) - 15)
71 #define test_generic_ui test_generic_si
72 #include "tgeneric_ui.c"
73
74 #define DEFN(N) \
75 static int powu##N (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd) \
76 { return mpfr_pow_ui (y, x, N, rnd); } \
77 static int pows##N (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd) \
78 { return mpfr_pow_si (y, x, N, rnd); } \
79 static int powm##N (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd) \
80 { return mpfr_pow_si (y, x, -(N), rnd); } \
81 static int root##N (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd) \
82 { return RAND_BOOL () ? \
83 mpfr_root (y, x, N, rnd) : mpfr_rootn_ui (y, x, N, rnd); } \
84 static int rootm##N (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd) \
85 { return mpfr_rootn_si (y, x, -(N), rnd); }
86
87
88 DEFN(2)
89 DEFN(3)
90 DEFN(4)
91 DEFN(5)
92 DEFN(17)
93 DEFN(120)
94
95 static void
check_pow_ui(void)96 check_pow_ui (void)
97 {
98 mpfr_t a, b;
99 unsigned long n;
100 int res;
101
102 mpfr_init2 (a, 53);
103 mpfr_init2 (b, 53);
104
105 /* check in-place operations */
106 mpfr_set_str (b, "0.6926773", 10, MPFR_RNDN);
107 mpfr_pow_ui (a, b, 10, MPFR_RNDN);
108 mpfr_pow_ui (b, b, 10, MPFR_RNDN);
109 if (mpfr_cmp (a, b))
110 {
111 printf ("Error for mpfr_pow_ui (b, b, ...)\n");
112 exit (1);
113 }
114
115 /* check large exponents */
116 mpfr_set_ui (b, 1, MPFR_RNDN);
117 mpfr_pow_ui (a, b, 4294967295UL, MPFR_RNDN);
118
119 mpfr_set_inf (a, -1);
120 mpfr_pow_ui (a, a, 4049053855UL, MPFR_RNDN);
121 if (!mpfr_inf_p (a) || (mpfr_sgn (a) >= 0))
122 {
123 printf ("Error for (-Inf)^4049053855\n");
124 exit (1);
125 }
126
127 mpfr_set_inf (a, -1);
128 mpfr_pow_ui (a, a, (unsigned long) 30002752, MPFR_RNDN);
129 if (!mpfr_inf_p (a) || (mpfr_sgn (a) <= 0))
130 {
131 printf ("Error for (-Inf)^30002752\n");
132 exit (1);
133 }
134
135 /* Check underflow */
136 mpfr_set_str_binary (a, "1E-1");
137 res = mpfr_pow_ui (a, a, -mpfr_get_emin (), MPFR_RNDN);
138 if (MPFR_GET_EXP (a) != mpfr_get_emin () + 1)
139 {
140 printf ("Error for (1e-1)^MPFR_EMAX_MAX\n");
141 mpfr_dump (a);
142 exit (1);
143 }
144
145 mpfr_set_str_binary (a, "1E-10");
146 res = mpfr_pow_ui (a, a, -mpfr_get_emin (), MPFR_RNDZ);
147 if (MPFR_NOTZERO (a))
148 {
149 printf ("Error for (1e-10)^MPFR_EMAX_MAX\n");
150 mpfr_dump (a);
151 exit (1);
152 }
153
154 /* Check overflow */
155 mpfr_set_str_binary (a, "1E10");
156 res = mpfr_pow_ui (a, a, ULONG_MAX, MPFR_RNDN);
157 if (!MPFR_IS_INF (a) || MPFR_IS_NEG (a))
158 {
159 printf ("Error for (1e10)^ULONG_MAX\n");
160 exit (1);
161 }
162
163 /* Bug in pow_ui.c from r3214 to r5107: if x = y (same mpfr_t argument),
164 the input argument is negative, n is odd, an overflow or underflow
165 occurs, and the temporary result res is positive, then the result
166 gets a wrong sign (positive instead of negative). */
167 mpfr_set_str_binary (a, "-1E10");
168 n = (ULONG_MAX ^ (ULONG_MAX >> 1)) + 1;
169 res = mpfr_pow_ui (a, a, n, MPFR_RNDN);
170 if (!MPFR_IS_INF (a) || MPFR_IS_POS (a))
171 {
172 printf ("Error for (-1e10)^%lu, expected -Inf,\ngot ", n);
173 mpfr_dump (a);
174 exit (1);
175 }
176
177 /* Check 0 */
178 MPFR_SET_ZERO (a);
179 MPFR_SET_POS (a);
180 mpfr_set_si (b, -1, MPFR_RNDN);
181 res = mpfr_pow_ui (b, a, 1, MPFR_RNDN);
182 if (res != 0 || MPFR_IS_NEG (b))
183 {
184 printf ("Error for (0+)^1\n");
185 exit (1);
186 }
187 MPFR_SET_ZERO (a);
188 MPFR_SET_NEG (a);
189 mpfr_set_ui (b, 1, MPFR_RNDN);
190 res = mpfr_pow_ui (b, a, 5, MPFR_RNDN);
191 if (res != 0 || MPFR_IS_POS (b))
192 {
193 printf ("Error for (0-)^5\n");
194 exit (1);
195 }
196 MPFR_SET_ZERO (a);
197 MPFR_SET_NEG (a);
198 mpfr_set_si (b, -1, MPFR_RNDN);
199 res = mpfr_pow_ui (b, a, 6, MPFR_RNDN);
200 if (res != 0 || MPFR_IS_NEG (b))
201 {
202 printf ("Error for (0-)^6\n");
203 exit (1);
204 }
205
206 mpfr_set_prec (a, 122);
207 mpfr_set_prec (b, 122);
208 mpfr_set_str_binary (a, "0.10000010010000111101001110100101101010011110011100001111000001001101000110011001001001001011001011010110110110101000111011E1");
209 mpfr_set_str_binary (b, "0.11111111100101001001000001000001100011100000001110111111100011111000111011100111111111110100011000111011000100100011001011E51290375");
210 mpfr_pow_ui (a, a, 2026876995UL, MPFR_RNDU);
211 if (mpfr_cmp (a, b) != 0)
212 {
213 printf ("Error for x^2026876995\n");
214 exit (1);
215 }
216
217 mpfr_set_prec (a, 29);
218 mpfr_set_prec (b, 29);
219 mpfr_set_str_binary (a, "1.0000000000000000000000001111");
220 mpfr_set_str_binary (b, "1.1001101111001100111001010111e165");
221 mpfr_pow_ui (a, a, 2055225053, MPFR_RNDZ);
222 if (mpfr_cmp (a, b) != 0)
223 {
224 printf ("Error for x^2055225053\n");
225 printf ("Expected ");
226 mpfr_out_str (stdout, 2, 0, b, MPFR_RNDN);
227 printf ("\nGot ");
228 mpfr_out_str (stdout, 2, 0, a, MPFR_RNDN);
229 printf ("\n");
230 exit (1);
231 }
232
233 /* worst case found by Vincent Lefevre, 25 Nov 2006 */
234 mpfr_set_prec (a, 53);
235 mpfr_set_prec (b, 53);
236 mpfr_set_str_binary (a, "1.0000010110000100001000101101101001011101101011010111");
237 mpfr_set_str_binary (b, "1.0000110111101111011010110100001100010000001010110100E1");
238 mpfr_pow_ui (a, a, 35, MPFR_RNDN);
239 if (mpfr_cmp (a, b) != 0)
240 {
241 printf ("Error in mpfr_pow_ui for worst case (1)\n");
242 printf ("Expected ");
243 mpfr_out_str (stdout, 2, 0, b, MPFR_RNDN);
244 printf ("\nGot ");
245 mpfr_out_str (stdout, 2, 0, a, MPFR_RNDN);
246 printf ("\n");
247 exit (1);
248 }
249 /* worst cases found on 2006-11-26 */
250 mpfr_set_str_binary (a, "1.0110100111010001101001010111001110010100111111000011");
251 mpfr_set_str_binary (b, "1.1111010011101110001111010110000101110000110110101100E17");
252 mpfr_pow_ui (a, a, 36, MPFR_RNDD);
253 if (mpfr_cmp (a, b) != 0)
254 {
255 printf ("Error in mpfr_pow_ui for worst case (2)\n");
256 printf ("Expected ");
257 mpfr_out_str (stdout, 2, 0, b, MPFR_RNDN);
258 printf ("\nGot ");
259 mpfr_out_str (stdout, 2, 0, a, MPFR_RNDN);
260 printf ("\n");
261 exit (1);
262 }
263 mpfr_set_str_binary (a, "1.1001010100001110000110111111100011011101110011000100");
264 mpfr_set_str_binary (b, "1.1100011101101101100010110001000001110001111110010001E23");
265 mpfr_pow_ui (a, a, 36, MPFR_RNDU);
266 if (mpfr_cmp (a, b) != 0)
267 {
268 printf ("Error in mpfr_pow_ui for worst case (3)\n");
269 printf ("Expected ");
270 mpfr_out_str (stdout, 2, 0, b, MPFR_RNDN);
271 printf ("\nGot ");
272 mpfr_out_str (stdout, 2, 0, a, MPFR_RNDN);
273 printf ("\n");
274 exit (1);
275 }
276
277 mpfr_clear (a);
278 mpfr_clear (b);
279 }
280
281 static void
check_pow_si(void)282 check_pow_si (void)
283 {
284 mpfr_t x;
285
286 mpfr_init (x);
287
288 mpfr_set_nan (x);
289 mpfr_pow_si (x, x, -1, MPFR_RNDN);
290 MPFR_ASSERTN(mpfr_nan_p (x));
291
292 mpfr_set_inf (x, 1);
293 mpfr_pow_si (x, x, -1, MPFR_RNDN);
294 MPFR_ASSERTN(MPFR_IS_ZERO (x) && MPFR_IS_POS (x));
295
296 mpfr_set_inf (x, -1);
297 mpfr_pow_si (x, x, -1, MPFR_RNDN);
298 MPFR_ASSERTN(MPFR_IS_ZERO (x) && MPFR_IS_NEG (x));
299
300 mpfr_set_inf (x, -1);
301 mpfr_pow_si (x, x, -2, MPFR_RNDN);
302 MPFR_ASSERTN(MPFR_IS_ZERO (x) && MPFR_IS_POS (x));
303
304 mpfr_set_ui (x, 0, MPFR_RNDN);
305 mpfr_pow_si (x, x, -1, MPFR_RNDN);
306 MPFR_ASSERTN(mpfr_inf_p (x) && mpfr_sgn (x) > 0);
307
308 mpfr_set_ui (x, 0, MPFR_RNDN);
309 mpfr_neg (x, x, MPFR_RNDN);
310 mpfr_pow_si (x, x, -1, MPFR_RNDN);
311 MPFR_ASSERTN(mpfr_inf_p (x) && mpfr_sgn (x) < 0);
312
313 mpfr_set_ui (x, 0, MPFR_RNDN);
314 mpfr_neg (x, x, MPFR_RNDN);
315 mpfr_pow_si (x, x, -2, MPFR_RNDN);
316 MPFR_ASSERTN(mpfr_inf_p (x) && mpfr_sgn (x) > 0);
317
318 mpfr_set_si (x, 2, MPFR_RNDN);
319 mpfr_pow_si (x, x, LONG_MAX, MPFR_RNDN); /* 2^LONG_MAX */
320 if (LONG_MAX > mpfr_get_emax () - 1) /* LONG_MAX + 1 > emax */
321 {
322 MPFR_ASSERTN (mpfr_inf_p (x));
323 }
324 else
325 {
326 MPFR_ASSERTN (mpfr_cmp_si_2exp (x, 1, (mpfr_exp_t) LONG_MAX));
327 }
328
329 mpfr_set_si (x, 2, MPFR_RNDN);
330 mpfr_pow_si (x, x, LONG_MIN, MPFR_RNDN); /* 2^LONG_MIN */
331 if (LONG_MIN + 1 < mpfr_get_emin ())
332 {
333 MPFR_ASSERTN (mpfr_zero_p (x));
334 }
335 else
336 {
337 MPFR_ASSERTN (mpfr_cmp_si_2exp (x, 1, (mpfr_exp_t) LONG_MIN));
338 }
339
340 mpfr_set_si (x, 2, MPFR_RNDN);
341 mpfr_pow_si (x, x, LONG_MIN + 1, MPFR_RNDN); /* 2^(LONG_MIN+1) */
342 if (mpfr_nan_p (x))
343 {
344 printf ("Error in pow_si(2, LONG_MIN+1): got NaN\n");
345 exit (1);
346 }
347 if (LONG_MIN + 2 < mpfr_get_emin ())
348 {
349 MPFR_ASSERTN (mpfr_zero_p (x));
350 }
351 else
352 {
353 MPFR_ASSERTN (mpfr_cmp_si_2exp (x, 1, (mpfr_exp_t) (LONG_MIN + 1)));
354 }
355
356 mpfr_set_si_2exp (x, 1, -1, MPFR_RNDN); /* 0.5 */
357 mpfr_pow_si (x, x, LONG_MIN, MPFR_RNDN); /* 2^(-LONG_MIN) */
358 if (LONG_MIN < 1 - mpfr_get_emax ()) /* 1 - LONG_MIN > emax */
359 {
360 MPFR_ASSERTN (mpfr_inf_p (x));
361 }
362 else
363 {
364 MPFR_ASSERTN (mpfr_cmp_si_2exp (x, 2, (mpfr_exp_t) - (LONG_MIN + 1)));
365 }
366
367 mpfr_clear (x);
368 }
369
370 /* check the IEEE 754-2019 special rules for pown */
371 static void
check_pown_ieee754_2019(void)372 check_pown_ieee754_2019 (void)
373 {
374 #ifdef _MPFR_H_HAVE_INTMAX_T
375 mpfr_t x;
376
377 mpfr_init2 (x, 5); /* ensures 17 is exact */
378
379 /* pown (x, 0) is 1 if x is not a signaling NaN: in MPFR we decide to
380 return 1 for a NaN */
381 mpfr_set_nan (x);
382 mpfr_pown (x, x, 0, MPFR_RNDN);
383 MPFR_ASSERTN(mpfr_cmp_ui0 (x, 1) == 0);
384 mpfr_set_inf (x, 1);
385 mpfr_pown (x, x, 0, MPFR_RNDN);
386 MPFR_ASSERTN(mpfr_cmp_ui0 (x, 1) == 0);
387 mpfr_set_inf (x, -1);
388 mpfr_pown (x, x, 0, MPFR_RNDN);
389 MPFR_ASSERTN(mpfr_cmp_ui0 (x, 1) == 0);
390 mpfr_set_zero (x, 1);
391 mpfr_pown (x, x, 0, MPFR_RNDN);
392 MPFR_ASSERTN(mpfr_cmp_ui0 (x, 1) == 0);
393 mpfr_set_zero (x, -1);
394 mpfr_pown (x, x, 0, MPFR_RNDN);
395 MPFR_ASSERTN(mpfr_cmp_ui0 (x, 1) == 0);
396 mpfr_set_si (x, 17, MPFR_RNDN);
397 mpfr_pown (x, x, 0, MPFR_RNDN);
398 MPFR_ASSERTN(mpfr_cmp_ui0 (x, 1) == 0);
399 mpfr_set_si (x, -17, MPFR_RNDN);
400 mpfr_pown (x, x, 0, MPFR_RNDN);
401 MPFR_ASSERTN(mpfr_cmp_ui0 (x, 1) == 0);
402
403 /* pown (±0, n) is ±∞ and signals the divideByZero exception for odd n < 0 */
404 mpfr_set_zero (x, 1);
405 mpfr_clear_divby0 ();
406 mpfr_pown (x, x, -17, MPFR_RNDN);
407 MPFR_ASSERTN(mpfr_inf_p (x) && mpfr_sgn (x) > 0 && mpfr_divby0_p ());
408 mpfr_set_zero (x, -1);
409 mpfr_clear_divby0 ();
410 mpfr_pown (x, x, -17, MPFR_RNDN);
411 MPFR_ASSERTN(mpfr_inf_p (x) && mpfr_sgn (x) < 0 && mpfr_divby0_p ());
412
413 /* pown (±0, n) is +∞ and signals the divideByZero exception for even n < 0*/
414 mpfr_set_zero (x, 1);
415 mpfr_clear_divby0 ();
416 mpfr_pown (x, x, -42, MPFR_RNDN);
417 MPFR_ASSERTN(mpfr_inf_p (x) && mpfr_sgn (x) > 0 && mpfr_divby0_p ());
418 mpfr_set_zero (x, -1);
419 mpfr_clear_divby0 ();
420 mpfr_pown (x, x, -42, MPFR_RNDN);
421 MPFR_ASSERTN(mpfr_inf_p (x) && mpfr_sgn (x) > 0 && mpfr_divby0_p ());
422
423 /* pown (±0, n) is +0 for even n > 0 */
424 mpfr_set_zero (x, 1);
425 mpfr_pown (x, x, 42, MPFR_RNDN);
426 MPFR_ASSERTN(mpfr_zero_p (x) && mpfr_signbit (x) == 0);
427 mpfr_set_zero (x, -1);
428 mpfr_pown (x, x, 42, MPFR_RNDN);
429 MPFR_ASSERTN(mpfr_zero_p (x) && mpfr_signbit (x) == 0);
430
431 /* pown (±0, n) is ±0 for odd n > 0 */
432 mpfr_set_zero (x, 1);
433 mpfr_pown (x, x, 17, MPFR_RNDN);
434 MPFR_ASSERTN(mpfr_zero_p (x) && mpfr_signbit (x) == 0);
435 mpfr_set_zero (x, -1);
436 mpfr_pown (x, x, 17, MPFR_RNDN);
437 MPFR_ASSERTN(mpfr_zero_p (x) && mpfr_signbit (x) == 1);
438
439 /* pown (+∞, n) is +∞ for n > 0 */
440 mpfr_set_inf (x, 1);
441 mpfr_pown (x, x, 17, MPFR_RNDN);
442 MPFR_ASSERTN(mpfr_inf_p (x) && mpfr_sgn (x) > 0);
443
444 /* pown (−∞, n) is −∞ for odd n > 0 */
445 mpfr_set_inf (x, -1);
446 mpfr_pown (x, x, 17, MPFR_RNDN);
447 MPFR_ASSERTN(mpfr_inf_p (x) && mpfr_sgn (x) < 0);
448
449 /* pown (−∞, n) is +∞ for even n > 0 */
450 mpfr_set_inf (x, -1);
451 mpfr_pown (x, x, 42, MPFR_RNDN);
452 MPFR_ASSERTN(mpfr_inf_p (x) && mpfr_sgn (x) > 0);
453
454 /* pown (+∞, n) is +0 for n < 0 */
455 mpfr_set_inf (x, 1);
456 mpfr_pown (x, x, -17, MPFR_RNDN);
457 MPFR_ASSERTN(mpfr_zero_p (x) && mpfr_signbit (x) == 0);
458 mpfr_set_inf (x, 1);
459 mpfr_pown (x, x, -42, MPFR_RNDN);
460 MPFR_ASSERTN(mpfr_zero_p (x) && mpfr_signbit (x) == 0);
461
462 /* pown (−∞, n) is −0 for odd n < 0 */
463 mpfr_set_inf (x, -1);
464 mpfr_pown (x, x, -17, MPFR_RNDN);
465 MPFR_ASSERTN(mpfr_zero_p (x) && mpfr_signbit (x) != 0);
466
467 /* pown (−∞, n) is +0 for even n < 0 */
468 mpfr_set_inf (x, -1);
469 mpfr_pown (x, x, -42, MPFR_RNDN);
470 MPFR_ASSERTN(mpfr_zero_p (x) && mpfr_signbit (x) == 0);
471
472 mpfr_clear (x);
473 #endif
474 }
475
476 static void
check_special_pow_si(void)477 check_special_pow_si (void)
478 {
479 mpfr_t a, b;
480 mpfr_exp_t emin;
481
482 mpfr_init (a);
483 mpfr_init (b);
484 mpfr_set_str (a, "2E100000000", 10, MPFR_RNDN);
485 mpfr_set_si (b, -10, MPFR_RNDN);
486 test_pow (b, a, b, MPFR_RNDN);
487 if (MPFR_NOTZERO(b))
488 {
489 printf("Pow(2E10000000, -10) failed\n");
490 mpfr_dump (a);
491 mpfr_dump (b);
492 exit(1);
493 }
494
495 emin = mpfr_get_emin ();
496 set_emin (-10);
497 mpfr_set_si (a, -2, MPFR_RNDN);
498 mpfr_pow_si (b, a, -10000, MPFR_RNDN);
499 if (MPFR_NOTZERO (b))
500 {
501 printf ("Pow_so (1, -10000) doesn't underflow if emin=-10.\n");
502 mpfr_dump (a);
503 mpfr_dump (b);
504 exit (1);
505 }
506 set_emin (emin);
507 mpfr_clear (a);
508 mpfr_clear (b);
509 }
510
511 static void
pow_si_long_min(void)512 pow_si_long_min (void)
513 {
514 mpfr_t x, y, z;
515 int inex;
516
517 mpfr_inits2 (sizeof(long) * CHAR_BIT + 32, x, y, z, (mpfr_ptr) 0);
518 mpfr_set_si_2exp (x, 3, -1, MPFR_RNDN); /* 1.5 */
519
520 inex = mpfr_set_si (y, LONG_MIN, MPFR_RNDN);
521 MPFR_ASSERTN (inex == 0);
522 mpfr_nextbelow (y);
523 mpfr_pow (y, x, y, MPFR_RNDD);
524
525 inex = mpfr_set_si (z, LONG_MIN, MPFR_RNDN);
526 MPFR_ASSERTN (inex == 0);
527 mpfr_nextabove (z);
528 mpfr_pow (z, x, z, MPFR_RNDU);
529
530 mpfr_pow_si (x, x, LONG_MIN, MPFR_RNDN); /* 1.5^LONG_MIN */
531 if (mpfr_cmp (x, y) < 0 || mpfr_cmp (x, z) > 0)
532 {
533 printf ("Error in pow_si_long_min\n");
534 exit (1);
535 }
536
537 mpfr_clears (x, y, z, (mpfr_ptr) 0);
538 }
539
540 static void
check_inexact(mpfr_prec_t p)541 check_inexact (mpfr_prec_t p)
542 {
543 mpfr_t x, y, z, t;
544 unsigned long u;
545 mpfr_prec_t q;
546 int inexact, cmp;
547 int rnd;
548
549 mpfr_init2 (x, p);
550 mpfr_init (y);
551 mpfr_init (z);
552 mpfr_init (t);
553 mpfr_urandomb (x, RANDS);
554 u = RAND_BOOL ();
555 for (q = MPFR_PREC_MIN; q <= p; q++)
556 RND_LOOP_NO_RNDF(rnd)
557 {
558 mpfr_set_prec (y, q);
559 mpfr_set_prec (z, q + 10);
560 mpfr_set_prec (t, q);
561 inexact = mpfr_pow_ui (y, x, u, (mpfr_rnd_t) rnd);
562 cmp = mpfr_pow_ui (z, x, u, (mpfr_rnd_t) rnd);
563 /* Note: that test makes no sense for RNDF, since according to the
564 reference manual, if the mpfr_can_round() call succeeds, one would
565 have to use RNDN in the mpfr_set() call below, which might give a
566 different value for t that the value y obtained above. */
567 if (mpfr_can_round (z, q + 10, (mpfr_rnd_t) rnd, (mpfr_rnd_t) rnd, q))
568 {
569 cmp = mpfr_set (t, z, (mpfr_rnd_t) rnd) || cmp;
570 if (mpfr_cmp (y, t))
571 {
572 printf ("results differ for u=%lu rnd=%s\n",
573 u, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
574 printf ("x="); mpfr_dump (x);
575 printf ("y="); mpfr_dump (y);
576 printf ("t="); mpfr_dump (t);
577 printf ("z="); mpfr_dump (z);
578 exit (1);
579 }
580 if (((inexact == 0) && (cmp != 0)) ||
581 ((inexact != 0) && (cmp == 0)))
582 {
583 printf ("Wrong inexact flag for p=%u, q=%u, rnd=%s\n",
584 (unsigned int) p, (unsigned int) q,
585 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
586 printf ("expected %d, got %d\n", cmp, inexact);
587 printf ("u=%lu x=", u); mpfr_dump (x);
588 printf ("y="); mpfr_dump (y);
589 exit (1);
590 }
591 }
592 }
593
594 /* check exact power */
595 mpfr_set_prec (x, p);
596 mpfr_set_prec (y, p);
597 mpfr_set_prec (z, p);
598 mpfr_set_ui (x, 4, MPFR_RNDN);
599 mpfr_set_str (y, "0.5", 10, MPFR_RNDN);
600 test_pow (z, x, y, MPFR_RNDZ);
601
602 mpfr_clear (x);
603 mpfr_clear (y);
604 mpfr_clear (z);
605 mpfr_clear (t);
606 }
607
608 static void
special(void)609 special (void)
610 {
611 mpfr_t x, y, z, t;
612 mpfr_exp_t emin, emax;
613 int inex;
614
615 mpfr_init2 (x, 53);
616 mpfr_init2 (y, 53);
617 mpfr_init2 (z, 53);
618 mpfr_init2 (t, 2);
619
620 mpfr_set_ui (x, 2, MPFR_RNDN);
621 mpfr_pow_si (x, x, -2, MPFR_RNDN);
622 if (mpfr_cmp_ui_2exp (x, 1, -2))
623 {
624 printf ("Error in pow_si(x,x,-2) for x=2\n");
625 exit (1);
626 }
627 mpfr_set_ui (x, 2, MPFR_RNDN);
628 mpfr_set_si (y, -2, MPFR_RNDN);
629 test_pow (x, x, y, MPFR_RNDN);
630 if (mpfr_cmp_ui_2exp (x, 1, -2))
631 {
632 printf ("Error in pow(x,x,y) for x=2, y=-2\n");
633 exit (1);
634 }
635
636 mpfr_set_prec (x, 2);
637 mpfr_set_str_binary (x, "1.0e-1");
638 mpfr_set_prec (y, 53);
639 mpfr_set_str_binary (y, "0.11010110011100101010110011001010100111000001000101110E-1");
640 mpfr_set_prec (z, 2);
641 test_pow (z, x, y, MPFR_RNDZ);
642 mpfr_set_str_binary (x, "1.0e-1");
643 if (mpfr_cmp (x, z))
644 {
645 printf ("Error in mpfr_pow (1)\n");
646 exit (1);
647 }
648
649 mpfr_set_prec (x, 64);
650 mpfr_set_prec (y, 64);
651 mpfr_set_prec (z, 64);
652 mpfr_set_prec (t, 64);
653 mpfr_set_str_binary (x, "0.111011000111100000111010000101010100110011010000011");
654 mpfr_set_str_binary (y, "0.111110010100110000011101100011010111000010000100101");
655 mpfr_set_str_binary (t, "0.1110110011110110001000110100100001001111010011111000010000011001");
656
657 test_pow (z, x, y, MPFR_RNDN);
658 if (mpfr_cmp (z, t))
659 {
660 printf ("Error in mpfr_pow for prec=64, rnd=MPFR_RNDN\n");
661 exit (1);
662 }
663
664 mpfr_set_prec (x, 53);
665 mpfr_set_prec (y, 53);
666 mpfr_set_prec (z, 53);
667 mpfr_set_str (x, "5.68824667828621954868e-01", 10, MPFR_RNDN);
668 mpfr_set_str (y, "9.03327850535952658895e-01", 10, MPFR_RNDN);
669 test_pow (z, x, y, MPFR_RNDZ);
670 if (mpfr_cmp_str1 (z, "0.60071044650456473235"))
671 {
672 printf ("Error in mpfr_pow for prec=53, rnd=MPFR_RNDZ\n");
673 exit (1);
674 }
675
676 mpfr_set_prec (t, 2);
677 mpfr_set_prec (x, 30);
678 mpfr_set_prec (y, 30);
679 mpfr_set_prec (z, 30);
680 mpfr_set_str (x, "1.00000000001010111110001111011e1", 2, MPFR_RNDN);
681 mpfr_set_str (t, "-0.5", 10, MPFR_RNDN);
682 test_pow (z, x, t, MPFR_RNDN);
683 mpfr_set_str (y, "1.01101001111010101110000101111e-1", 2, MPFR_RNDN);
684 if (mpfr_cmp (z, y))
685 {
686 printf ("Error in mpfr_pow for prec=30, rnd=MPFR_RNDN\n");
687 exit (1);
688 }
689
690 mpfr_set_prec (x, 21);
691 mpfr_set_prec (y, 21);
692 mpfr_set_prec (z, 21);
693 mpfr_set_str (x, "1.11111100100001100101", 2, MPFR_RNDN);
694 test_pow (z, x, t, MPFR_RNDZ);
695 mpfr_set_str (y, "1.01101011010001100000e-1", 2, MPFR_RNDN);
696 if (mpfr_cmp (z, y))
697 {
698 printf ("Error in mpfr_pow for prec=21, rnd=MPFR_RNDZ\n");
699 printf ("Expected ");
700 mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN);
701 printf ("\nGot ");
702 mpfr_out_str (stdout, 2, 0, z, MPFR_RNDN);
703 printf ("\n");
704 exit (1);
705 }
706
707 /* From https://web.archive.org/web/20050824044408/http://www.terra.es/personal9/ismaeljc/hall.htm */
708 mpfr_set_prec (x, 113);
709 mpfr_set_prec (y, 2);
710 mpfr_set_prec (z, 169);
711 mpfr_set_str1 (x, "6078673043126084065007902175846955");
712 mpfr_set_ui_2exp (y, 3, -1, MPFR_RNDN);
713 test_pow (z, x, y, MPFR_RNDZ);
714 if (mpfr_cmp_str1 (z, "473928882491000966028828671876527456070714790264144"))
715 {
716 printf ("Error in mpfr_pow for 6078673043126084065007902175846955");
717 printf ("^(3/2), MPFR_RNDZ\nExpected ");
718 printf ("4.73928882491000966028828671876527456070714790264144e50");
719 printf ("\nGot ");
720 mpfr_out_str (stdout, 10, 0, z, MPFR_RNDN);
721 printf ("\n");
722 exit (1);
723 }
724 test_pow (z, x, y, MPFR_RNDU);
725 if (mpfr_cmp_str1 (z, "473928882491000966028828671876527456070714790264145"))
726 {
727 printf ("Error in mpfr_pow for 6078673043126084065007902175846955");
728 printf ("^(3/2), MPFR_RNDU\nExpected ");
729 printf ("4.73928882491000966028828671876527456070714790264145e50");
730 printf ("\nGot ");
731 mpfr_out_str (stdout, 10, 0, z, MPFR_RNDN);
732 printf ("\n");
733 exit (1);
734 }
735
736 mpfr_set_inf (x, 1);
737 mpfr_set_prec (y, 2);
738 mpfr_set_str_binary (y, "1E10");
739 test_pow (z, x, y, MPFR_RNDN);
740 MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_POS(z));
741 mpfr_set_inf (x, -1);
742 test_pow (z, x, y, MPFR_RNDN);
743 MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_POS(z));
744 mpfr_set_prec (y, 10);
745 mpfr_set_str_binary (y, "1.000000001E9");
746 test_pow (z, x, y, MPFR_RNDN);
747 MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_NEG(z));
748 mpfr_set_str_binary (y, "1.000000001E8");
749 test_pow (z, x, y, MPFR_RNDN);
750 MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_POS(z));
751
752 mpfr_set_inf (x, -1);
753 mpfr_set_prec (y, 2 * mp_bits_per_limb);
754 mpfr_set_ui (y, 1, MPFR_RNDN);
755 mpfr_mul_2ui (y, y, mp_bits_per_limb - 1, MPFR_RNDN);
756 /* y = 2^(mp_bits_per_limb - 1) */
757 test_pow (z, x, y, MPFR_RNDN);
758 MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_POS(z));
759 mpfr_nextabove (y);
760 test_pow (z, x, y, MPFR_RNDN);
761 /* y = 2^(mp_bits_per_limb - 1) + epsilon */
762 MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_POS(z));
763 mpfr_nextbelow (y);
764 mpfr_div_2ui (y, y, 1, MPFR_RNDN);
765 mpfr_nextabove (y);
766 test_pow (z, x, y, MPFR_RNDN);
767 /* y = 2^(mp_bits_per_limb - 2) + epsilon */
768 MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_POS(z));
769
770 mpfr_set_si (x, -1, MPFR_RNDN);
771 mpfr_set_prec (y, 2);
772 mpfr_set_str_binary (y, "1E10");
773 test_pow (z, x, y, MPFR_RNDN);
774 MPFR_ASSERTN(mpfr_cmp_ui0 (z, 1) == 0);
775
776 /* Check (-0)^(17.0001) */
777 mpfr_set_prec (x, 6);
778 mpfr_set_prec (y, 640);
779 MPFR_SET_ZERO (x); MPFR_SET_NEG (x);
780 mpfr_set_ui (y, 17, MPFR_RNDN); mpfr_nextabove (y);
781 test_pow (z, x, y, MPFR_RNDN);
782 MPFR_ASSERTN (MPFR_IS_ZERO (z) && MPFR_IS_POS (z));
783
784 /* Bugs reported by Kevin Rauch on 29 Oct 2007 */
785 emin = mpfr_get_emin ();
786 emax = mpfr_get_emax ();
787 set_emin (-1000000);
788 set_emax ( 1000000);
789 mpfr_set_prec (x, 64);
790 mpfr_set_prec (y, 64);
791 mpfr_set_prec (z, 64);
792 mpfr_set_str (x, "-0.5", 10, MPFR_RNDN);
793 mpfr_set_str (y, "-0.ffffffffffffffff", 16, MPFR_RNDN);
794 mpfr_set_exp (y, mpfr_get_emax ());
795 inex = mpfr_pow (z, x, y, MPFR_RNDN);
796 /* (-0.5)^(-n) = 1/2^n for n even */
797 MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_POS (z) && inex > 0);
798
799 /* (-1)^(-n) = 1 for n even */
800 mpfr_set_str (x, "-1", 10, MPFR_RNDN);
801 inex = mpfr_pow (z, x, y, MPFR_RNDN);
802 MPFR_ASSERTN(mpfr_cmp_ui0 (z, 1) == 0 && inex == 0);
803
804 /* (-1)^n = 1 for n even */
805 mpfr_set_str (x, "-1", 10, MPFR_RNDN);
806 mpfr_neg (y, y, MPFR_RNDN);
807 inex = mpfr_pow (z, x, y, MPFR_RNDN);
808 MPFR_ASSERTN(mpfr_cmp_ui0 (z, 1) == 0 && inex == 0);
809
810 /* (-1.5)^n = +Inf for n even */
811 mpfr_set_str (x, "-1.5", 10, MPFR_RNDN);
812 inex = mpfr_pow (z, x, y, MPFR_RNDN);
813 MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_POS (z) && inex > 0);
814
815 /* (-n)^1.5 = NaN for n even */
816 mpfr_neg (y, y, MPFR_RNDN);
817 mpfr_set_str (x, "1.5", 10, MPFR_RNDN);
818 inex = mpfr_pow (z, y, x, MPFR_RNDN);
819 MPFR_ASSERTN(mpfr_nan_p (z));
820
821 /* x^(-1.5) = NaN for x small < 0 */
822 mpfr_set_str (x, "-0.8", 16, MPFR_RNDN);
823 mpfr_set_exp (x, mpfr_get_emin ());
824 mpfr_set_str (y, "-1.5", 10, MPFR_RNDN);
825 inex = mpfr_pow (z, x, y, MPFR_RNDN);
826 MPFR_ASSERTN(mpfr_nan_p (z));
827
828 set_emin (emin);
829 set_emax (emax);
830 mpfr_clear (x);
831 mpfr_clear (y);
832 mpfr_clear (z);
833 mpfr_clear (t);
834 }
835
836 static void
particular_cases(void)837 particular_cases (void)
838 {
839 mpfr_t t[11], r, r2;
840 mpz_t z;
841 long si;
842
843 static const char *name[11] = {
844 "NaN", "+inf", "-inf", "+0", "-0", "+1", "-1", "+2", "-2", "+0.5", "-0.5"};
845 int i, j;
846 int error = 0;
847
848 mpz_init (z);
849
850 for (i = 0; i < 11; i++)
851 mpfr_init2 (t[i], 2);
852 mpfr_init2 (r, 6);
853 mpfr_init2 (r2, 6);
854
855 mpfr_set_nan (t[0]);
856 mpfr_set_inf (t[1], 1);
857 mpfr_set_ui (t[3], 0, MPFR_RNDN);
858 mpfr_set_ui (t[5], 1, MPFR_RNDN);
859 mpfr_set_ui (t[7], 2, MPFR_RNDN);
860 mpfr_div_2ui (t[9], t[5], 1, MPFR_RNDN);
861 for (i = 1; i < 11; i += 2)
862 mpfr_neg (t[i+1], t[i], MPFR_RNDN);
863
864 for (i = 0; i < 11; i++)
865 for (j = 0; j < 11; j++)
866 {
867 double d;
868 int p;
869 static const int q[11][11] = {
870 /* NaN +inf -inf +0 -0 +1 -1 +2 -2 +0.5 -0.5 */
871 /* NaN */ { 0, 0, 0, 128, 128, 0, 0, 0, 0, 0, 0 },
872 /* +inf */ { 0, 1, 2, 128, 128, 1, 2, 1, 2, 1, 2 },
873 /* -inf */ { 0, 1, 2, 128, 128, -1, -2, 1, 2, 1, 2 },
874 /* +0 */ { 0, 2, 1, 128, 128, 2, 1, 2, 1, 2, 1 },
875 /* -0 */ { 0, 2, 1, 128, 128, -2, -1, 2, 1, 2, 1 },
876 /* +1 */ {128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128 },
877 /* -1 */ { 0, 128, 128, 128, 128,-128,-128, 128, 128, 0, 0 },
878 /* +2 */ { 0, 1, 2, 128, 128, 256, 64, 512, 32, 180, 90 },
879 /* -2 */ { 0, 1, 2, 128, 128,-256, -64, 512, 32, 0, 0 },
880 /* +0.5 */ { 0, 2, 1, 128, 128, 64, 256, 32, 512, 90, 180 },
881 /* -0.5 */ { 0, 2, 1, 128, 128, -64,-256, 32, 512, 0, 0 }
882 };
883 /* This define is used to make the following table readable */
884 #define N MPFR_FLAGS_NAN
885 #define I MPFR_FLAGS_INEXACT
886 #define D MPFR_FLAGS_DIVBY0
887 static const unsigned int f[11][11] = {
888 /* NaN +inf -inf +0 -0 +1 -1 +2 -2 +0.5 -0.5 */
889 /* NaN */ { N, N, N, 0, 0, N, N, N, N, N, N },
890 /* +inf */ { N, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
891 /* -inf */ { N, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
892 /* +0 */ { N, 0, 0, 0, 0, 0, D, 0, D, 0, D },
893 /* -0 */ { N, 0, 0, 0, 0, 0, D, 0, D, 0, D },
894 /* +1 */ { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
895 /* -1 */ { N, 0, 0, 0, 0, 0, 0, 0, 0, N, N },
896 /* +2 */ { N, 0, 0, 0, 0, 0, 0, 0, 0, I, I },
897 /* -2 */ { N, 0, 0, 0, 0, 0, 0, 0, 0, N, N },
898 /* +0.5 */ { N, 0, 0, 0, 0, 0, 0, 0, 0, I, I },
899 /* -0.5 */ { N, 0, 0, 0, 0, 0, 0, 0, 0, N, N }
900 };
901 #undef N
902 #undef I
903 #undef D
904 mpfr_clear_flags ();
905 test_pow (r, t[i], t[j], MPFR_RNDN);
906 p = mpfr_nan_p (r) ? 0 : mpfr_inf_p (r) ? 1 :
907 mpfr_zero_p (r) ? 2 :
908 (d = mpfr_get_d (r, MPFR_RNDN), (int) (ABS(d) * 128.0));
909 if (p != 0 && MPFR_IS_NEG (r))
910 p = -p;
911 if (p != q[i][j])
912 {
913 printf ("Error in mpfr_pow for (%s)^(%s) (%d,%d):\n"
914 "got %d instead of %d\n",
915 name[i], name[j], i, j, p, q[i][j]);
916 mpfr_dump (r);
917 error = 1;
918 }
919 if (__gmpfr_flags != f[i][j])
920 {
921 printf ("Error in mpfr_pow for (%s)^(%s) (%d,%d):\n"
922 "Flags = %u instead of expected %u\n",
923 name[i], name[j], i, j,
924 (unsigned int) __gmpfr_flags, f[i][j]);
925 mpfr_dump (r);
926 error = 1;
927 }
928 /* Perform the same tests with pow_z & pow_si & pow_ui
929 if t[j] is an integer */
930 if (mpfr_integer_p (t[j]))
931 {
932 /* mpfr_pow_z */
933 mpfr_clear_flags ();
934 mpfr_get_z (z, t[j], MPFR_RNDN);
935 mpfr_pow_z (r, t[i], z, MPFR_RNDN);
936 p = mpfr_nan_p (r) ? 0 : mpfr_inf_p (r) ? 1 :
937 mpfr_zero_p (r) ? 2 :
938 (d = mpfr_get_d (r, MPFR_RNDN), (int) (ABS(d) * 128.0));
939 if (p != 0 && MPFR_IS_NEG (r))
940 p = -p;
941 if (p != q[i][j])
942 {
943 printf ("Error in mpfr_pow_z for (%s)^(%s) (%d,%d):\n"
944 "got %d instead of %d\n",
945 name[i], name[j], i, j, p, q[i][j]);
946 mpfr_dump (r);
947 error = 1;
948 }
949 if (__gmpfr_flags != f[i][j])
950 {
951 printf ("Error in mpfr_pow_z for (%s)^(%s) (%d,%d):\n"
952 "Flags = %u instead of expected %u\n",
953 name[i], name[j], i, j,
954 (unsigned int) __gmpfr_flags, f[i][j]);
955 mpfr_dump (r);
956 error = 1;
957 }
958 /* mpfr_pow_si */
959 mpfr_clear_flags ();
960 si = mpfr_get_si (t[j], MPFR_RNDN);
961 mpfr_pow_si (r, t[i], si, MPFR_RNDN);
962 p = mpfr_nan_p (r) ? 0 : mpfr_inf_p (r) ? 1 :
963 mpfr_zero_p (r) ? 2 :
964 (d = mpfr_get_d (r, MPFR_RNDN), (int) (ABS(d) * 128.0));
965 if (p != 0 && MPFR_IS_NEG (r))
966 p = -p;
967 if (p != q[i][j])
968 {
969 printf ("Error in mpfr_pow_si for (%s)^(%s) (%d,%d):\n"
970 "got %d instead of %d\n",
971 name[i], name[j], i, j, p, q[i][j]);
972 mpfr_dump (r);
973 error = 1;
974 }
975 if (__gmpfr_flags != f[i][j])
976 {
977 printf ("Error in mpfr_pow_si for (%s)^(%s) (%d,%d):\n"
978 "Flags = %u instead of expected %u\n",
979 name[i], name[j], i, j,
980 (unsigned int) __gmpfr_flags, f[i][j]);
981 mpfr_dump (r);
982 error = 1;
983 }
984 /* if si >= 0, test mpfr_pow_ui */
985 if (si >= 0)
986 {
987 mpfr_clear_flags ();
988 mpfr_pow_ui (r, t[i], si, MPFR_RNDN);
989 p = mpfr_nan_p (r) ? 0 : mpfr_inf_p (r) ? 1 :
990 mpfr_zero_p (r) ? 2 :
991 (d = mpfr_get_d (r, MPFR_RNDN), (int) (ABS(d) * 128.0));
992 if (p != 0 && MPFR_IS_NEG (r))
993 p = -p;
994 if (p != q[i][j])
995 {
996 printf ("Error in mpfr_pow_ui for (%s)^(%s) (%d,%d):\n"
997 "got %d instead of %d\n",
998 name[i], name[j], i, j, p, q[i][j]);
999 mpfr_dump (r);
1000 error = 1;
1001 }
1002 if (__gmpfr_flags != f[i][j])
1003 {
1004 printf ("Error in mpfr_pow_ui for (%s)^(%s) (%d,%d):\n"
1005 "Flags = %u instead of expected %u\n",
1006 name[i], name[j], i, j,
1007 (unsigned int) __gmpfr_flags, f[i][j]);
1008 mpfr_dump (r);
1009 error = 1;
1010 }
1011 }
1012 } /* integer_p */
1013 /* Perform the same tests with mpfr_ui_pow */
1014 if (mpfr_integer_p (t[i]) && MPFR_IS_POS (t[i]))
1015 {
1016 /* mpfr_ui_pow */
1017 mpfr_clear_flags ();
1018 si = mpfr_get_si (t[i], MPFR_RNDN);
1019 mpfr_ui_pow (r, si, t[j], MPFR_RNDN);
1020 p = mpfr_nan_p (r) ? 0 : mpfr_inf_p (r) ? 1 :
1021 mpfr_zero_p (r) ? 2 :
1022 (d = mpfr_get_d (r, MPFR_RNDN), (int) (ABS(d) * 128.0));
1023 if (p != 0 && MPFR_IS_NEG (r))
1024 p = -p;
1025 if (p != q[i][j])
1026 {
1027 printf ("Error in mpfr_ui_pow for (%s)^(%s) (%d,%d):\n"
1028 "got %d instead of %d\n",
1029 name[i], name[j], i, j, p, q[i][j]);
1030 mpfr_dump (r);
1031 error = 1;
1032 }
1033 if (__gmpfr_flags != f[i][j])
1034 {
1035 printf ("Error in mpfr_ui_pow for (%s)^(%s) (%d,%d):\n"
1036 "Flags = %u instead of expected %u\n",
1037 name[i], name[j], i, j,
1038 (unsigned int) __gmpfr_flags, f[i][j]);
1039 mpfr_dump (r);
1040 error = 1;
1041 }
1042 }
1043 }
1044
1045 for (i = 0; i < 11; i++)
1046 mpfr_clear (t[i]);
1047 mpfr_clear (r);
1048 mpfr_clear (r2);
1049 mpz_clear (z);
1050
1051 if (error)
1052 exit (1);
1053 }
1054
1055 static void
underflows(void)1056 underflows (void)
1057 {
1058 mpfr_t x, y, z;
1059 int err = 0;
1060 int inexact;
1061 int i;
1062 mpfr_exp_t emin;
1063
1064 mpfr_init2 (x, 64);
1065 mpfr_init2 (y, 64);
1066
1067 mpfr_set_ui (x, 1, MPFR_RNDN);
1068 mpfr_set_exp (x, mpfr_get_emin());
1069
1070 for (i = 3; i < 10; i++)
1071 {
1072 mpfr_set_ui (y, i, MPFR_RNDN);
1073 mpfr_div_2ui (y, y, 1, MPFR_RNDN);
1074 test_pow (y, x, y, MPFR_RNDN);
1075 if (MPFR_NOTZERO (y))
1076 {
1077 printf ("Error in mpfr_pow for ");
1078 mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN);
1079 printf (" ^ (%d/2)\nGot ", i);
1080 mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN);
1081 printf (" instead of 0.\n");
1082 exit (1);
1083 }
1084 }
1085
1086 mpfr_init2 (z, 55);
1087 mpfr_set_str (x, "0.110011010011101001110001110100010000110111101E0",
1088 2, MPFR_RNDN);
1089 mpfr_set_str (y, "0.101110010011111001011010100011011100111110011E40",
1090 2, MPFR_RNDN);
1091 mpfr_clear_flags ();
1092 inexact = mpfr_pow (z, x, y, MPFR_RNDU);
1093 if (!mpfr_underflow_p ())
1094 {
1095 printf ("Underflow flag is not set for special underflow test.\n");
1096 err = 1;
1097 }
1098 if (inexact <= 0)
1099 {
1100 printf ("Ternary value is wrong for special underflow test.\n");
1101 err = 1;
1102 }
1103 mpfr_set_ui (x, 0, MPFR_RNDN);
1104 mpfr_nextabove (x);
1105 if (mpfr_cmp (x, z) != 0)
1106 {
1107 printf ("Wrong value for special underflow test.\nGot ");
1108 mpfr_out_str (stdout, 2, 0, z, MPFR_RNDN);
1109 printf ("\ninstead of ");
1110 mpfr_out_str (stdout, 2, 2, x, MPFR_RNDN);
1111 printf ("\n");
1112 err = 1;
1113 }
1114 if (err)
1115 exit (1);
1116
1117 /* MPFR currently (2006-08-19) segfaults on the following code (and
1118 possibly makes other programs crash due to the lack of memory),
1119 because y is converted into an mpz_t, and the required precision
1120 is too high. */
1121 mpfr_set_prec (x, 2);
1122 mpfr_set_prec (y, 2);
1123 mpfr_set_prec (z, 12);
1124 mpfr_set_ui_2exp (x, 3, -2, MPFR_RNDN);
1125 mpfr_set_ui_2exp (y, 1, mpfr_get_emax () - 1, MPFR_RNDN);
1126 mpfr_clear_flags ();
1127 mpfr_pow (z, x, y, MPFR_RNDN);
1128 if (!mpfr_underflow_p () || MPFR_NOTZERO (z))
1129 {
1130 printf ("Underflow test with large y fails.\n");
1131 exit (1);
1132 }
1133
1134 emin = mpfr_get_emin ();
1135 set_emin (-256);
1136 mpfr_set_prec (x, 2);
1137 mpfr_set_prec (y, 2);
1138 mpfr_set_prec (z, 12);
1139 mpfr_set_ui_2exp (x, 3, -2, MPFR_RNDN);
1140 mpfr_set_ui_2exp (y, 1, 38, MPFR_RNDN);
1141 mpfr_clear_flags ();
1142 inexact = mpfr_pow (z, x, y, MPFR_RNDN);
1143 if (!mpfr_underflow_p () || MPFR_NOTZERO (z) || inexact >= 0)
1144 {
1145 printf ("Bad underflow detection for 0.75^(2^38). Obtained:\n"
1146 "Underflow flag... %-3s (should be 'yes')\n"
1147 "Zero result...... %-3s (should be 'yes')\n"
1148 "Inexact value.... %-3d (should be negative)\n",
1149 mpfr_underflow_p () ? "yes" : "no",
1150 MPFR_IS_ZERO (z) ? "yes" : "no", inexact);
1151 exit (1);
1152 }
1153 set_emin (emin);
1154
1155 emin = mpfr_get_emin ();
1156 set_emin (-256);
1157 mpfr_set_prec (x, 2);
1158 mpfr_set_prec (y, 40);
1159 mpfr_set_prec (z, 12);
1160 mpfr_set_ui_2exp (x, 3, -1, MPFR_RNDN);
1161 mpfr_set_si_2exp (y, -1, 38, MPFR_RNDN);
1162 for (i = 0; i < 4; i++)
1163 {
1164 if (i == 2)
1165 mpfr_neg (x, x, MPFR_RNDN);
1166 mpfr_clear_flags ();
1167 inexact = mpfr_pow (z, x, y, MPFR_RNDN);
1168 if (!mpfr_underflow_p () || MPFR_NOTZERO (z) ||
1169 (i == 3 ? (inexact <= 0) : (inexact >= 0)))
1170 {
1171 printf ("Bad underflow detection for (");
1172 mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN);
1173 printf (")^(-2^38-%d). Obtained:\n"
1174 "Overflow flag.... %-3s (should be 'no')\n"
1175 "Underflow flag... %-3s (should be 'yes')\n"
1176 "Zero result...... %-3s (should be 'yes')\n"
1177 "Inexact value.... %-3d (should be %s)\n", i,
1178 mpfr_overflow_p () ? "yes" : "no",
1179 mpfr_underflow_p () ? "yes" : "no",
1180 MPFR_IS_ZERO (z) ? "yes" : "no", inexact,
1181 i == 3 ? "positive" : "negative");
1182 exit (1);
1183 }
1184 inexact = mpfr_sub_ui (y, y, 1, MPFR_RNDN);
1185 MPFR_ASSERTN (inexact == 0);
1186 }
1187 set_emin (emin);
1188
1189 mpfr_clears (x, y, z, (mpfr_ptr) 0);
1190 }
1191
1192 static void
overflows(void)1193 overflows (void)
1194 {
1195 mpfr_t a, b;
1196
1197 /* bug found by Ming J. Tsai <mingjt@delvron.us>, 4 Oct 2003 */
1198
1199 mpfr_init_set_str (a, "5.1e32", 10, MPFR_RNDN);
1200 mpfr_init (b);
1201
1202 test_pow (b, a, a, MPFR_RNDN);
1203 if (!(mpfr_inf_p (b) && mpfr_sgn (b) > 0))
1204 {
1205 printf ("Error for a^a for a=5.1e32\n");
1206 printf ("Expected +Inf, got ");
1207 mpfr_out_str (stdout, 10, 0, b, MPFR_RNDN);
1208 printf ("\n");
1209 exit (1);
1210 }
1211
1212 mpfr_clear(a);
1213 mpfr_clear(b);
1214 }
1215
1216 static void
overflows2(void)1217 overflows2 (void)
1218 {
1219 mpfr_t x, y, z;
1220 mpfr_exp_t emin, emax;
1221 int e;
1222
1223 /* x^y in reduced exponent range, where x = 2^b and y is not an integer
1224 (so that mpfr_pow_z is not used). */
1225
1226 emin = mpfr_get_emin ();
1227 emax = mpfr_get_emax ();
1228 set_emin (-128);
1229
1230 mpfr_inits2 (16, x, y, z, (mpfr_ptr) 0);
1231
1232 mpfr_set_si_2exp (x, 1, -64, MPFR_RNDN); /* 2^(-64) */
1233 mpfr_set_si_2exp (y, -1, -1, MPFR_RNDN); /* -0.5 */
1234 for (e = 2; e <= 32; e += 17)
1235 {
1236 set_emax (e);
1237 mpfr_clear_flags ();
1238 mpfr_pow (z, x, y, MPFR_RNDN);
1239 if (MPFR_IS_NEG (z) || ! mpfr_inf_p (z))
1240 {
1241 printf ("Error in overflows2 (e = %d): expected +Inf, got ", e);
1242 mpfr_dump (z);
1243 exit (1);
1244 }
1245 if (__gmpfr_flags != (MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT))
1246 {
1247 printf ("Error in overflows2 (e = %d): bad flags (%u)\n",
1248 e, (unsigned int) __gmpfr_flags);
1249 exit (1);
1250 }
1251 }
1252
1253 mpfr_clears (x, y, z, (mpfr_ptr) 0);
1254
1255 set_emin (emin);
1256 set_emax (emax);
1257 }
1258
1259 static void
overflows3(void)1260 overflows3 (void)
1261 {
1262 /* x^y where x = 2^b, y is not an integer (so that mpfr_pow_z is not used)
1263 and b * y = emax in the extended exponent range. If emax is divisible
1264 by 3, we choose x = 2^(-2*emax/3) and y = -3/2.
1265 Test also with nextbelow(x). */
1266
1267 if (MPFR_EMAX_MAX % 3 == 0)
1268 {
1269 mpfr_t x, y, z, t;
1270 mpfr_exp_t emin, emax;
1271 unsigned int flags;
1272 int i;
1273
1274 emin = mpfr_get_emin ();
1275 emax = mpfr_get_emax ();
1276 set_emin (MPFR_EMIN_MIN);
1277 set_emax (MPFR_EMAX_MAX);
1278
1279 mpfr_inits2 (16, x, y, z, t, (mpfr_ptr) 0);
1280
1281 mpfr_set_si_2exp (x, 1, -2 * (MPFR_EMAX_MAX / 3), MPFR_RNDN);
1282 for (i = 0; i <= 1; i++)
1283 {
1284 mpfr_set_si_2exp (y, -3, -1, MPFR_RNDN);
1285 mpfr_clear_flags ();
1286 mpfr_pow (z, x, y, MPFR_RNDN);
1287 if (MPFR_IS_NEG (z) || ! mpfr_inf_p (z))
1288 {
1289 printf ("Error in overflows3 (RNDN, i = %d): expected +Inf,"
1290 " got ", i);
1291 mpfr_dump (z);
1292 exit (1);
1293 }
1294 if (__gmpfr_flags != (MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT))
1295 {
1296 printf ("Error in overflows3 (RNDN, i = %d): bad flags (%u)\n",
1297 i, (unsigned int) __gmpfr_flags);
1298 exit (1);
1299 }
1300
1301 mpfr_clear_flags ();
1302 mpfr_pow (z, x, y, MPFR_RNDZ);
1303 flags = __gmpfr_flags;
1304 mpfr_set (t, z, MPFR_RNDN);
1305 mpfr_nextabove (t);
1306 if (MPFR_IS_NEG (z) || mpfr_inf_p (z) || ! mpfr_inf_p (t))
1307 {
1308 printf ("Error in overflows3 (RNDZ, i = %d):\nexpected ", i);
1309 mpfr_set_inf (t, 1);
1310 mpfr_nextbelow (t);
1311 mpfr_dump (t);
1312 printf ("got ");
1313 mpfr_dump (z);
1314 exit (1);
1315 }
1316 if (flags != (MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT))
1317 {
1318 printf ("Error in overflows3 (RNDZ, i = %d): bad flags (%u)\n",
1319 i, flags);
1320 exit (1);
1321 }
1322 mpfr_nextbelow (x);
1323 }
1324
1325 mpfr_clears (x, y, z, t, (mpfr_ptr) 0);
1326
1327 set_emin (emin);
1328 set_emax (emax);
1329 }
1330 }
1331
1332 static void
x_near_one(void)1333 x_near_one (void)
1334 {
1335 mpfr_t x, y, z;
1336 int inex;
1337
1338 mpfr_init2 (x, 32);
1339 mpfr_init2 (y, 4);
1340 mpfr_init2 (z, 33);
1341
1342 mpfr_set_ui (x, 1, MPFR_RNDN);
1343 mpfr_nextbelow (x);
1344 mpfr_set_ui_2exp (y, 11, -2, MPFR_RNDN);
1345 inex = mpfr_pow (z, x, y, MPFR_RNDN);
1346 if (mpfr_cmp_str (z, "0.111111111111111111111111111111011E0", 2, MPFR_RNDN)
1347 || inex <= 0)
1348 {
1349 printf ("Failure in x_near_one, got inex = %d and\nz = ", inex);
1350 mpfr_dump (z);
1351 }
1352
1353 mpfr_clears (x, y, z, (mpfr_ptr) 0);
1354 }
1355
1356 static int
mpfr_pow275(mpfr_ptr y,mpfr_srcptr x,mpfr_rnd_t r)1357 mpfr_pow275 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t r)
1358 {
1359 mpfr_t z;
1360 int inex;
1361
1362 mpfr_init2 (z, 4);
1363 mpfr_set_ui_2exp (z, 11, -2, MPFR_RNDN);
1364 inex = mpfr_pow (y, x, z, MPFR_RNDN);
1365 mpfr_clear (z);
1366 return inex;
1367 }
1368
1369 /* Bug found by Kevin P. Rauch */
1370 static void
bug20071103(void)1371 bug20071103 (void)
1372 {
1373 mpfr_t x, y, z;
1374 mpfr_exp_t emin, emax;
1375
1376 emin = mpfr_get_emin ();
1377 emax = mpfr_get_emax ();
1378 set_emin (-1000000);
1379 set_emax ( 1000000);
1380
1381 mpfr_inits2 (64, x, y, z, (mpfr_ptr) 0);
1382 mpfr_set_si_2exp (x, -3, -1, MPFR_RNDN); /* x = -1.5 */
1383 mpfr_set_str (y, "-0.ffffffffffffffff", 16, MPFR_RNDN);
1384 mpfr_set_exp (y, mpfr_get_emax ());
1385 mpfr_clear_flags ();
1386 mpfr_pow (z, x, y, MPFR_RNDN);
1387 MPFR_ASSERTN (mpfr_zero_p (z) && MPFR_IS_POS (z) &&
1388 __gmpfr_flags == (MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT));
1389 mpfr_clears (x, y, z, (mpfr_ptr) 0);
1390
1391 set_emin (emin);
1392 set_emax (emax);
1393 }
1394
1395 /* Bug found by Kevin P. Rauch */
1396 static void
bug20071104(void)1397 bug20071104 (void)
1398 {
1399 mpfr_t x, y, z;
1400 mpfr_exp_t emin, emax;
1401 int inex;
1402
1403 emin = mpfr_get_emin ();
1404 emax = mpfr_get_emax ();
1405 set_emin (-1000000);
1406 set_emax ( 1000000);
1407
1408 mpfr_inits2 (20, x, y, z, (mpfr_ptr) 0);
1409 mpfr_set_ui (x, 0, MPFR_RNDN);
1410 mpfr_nextbelow (x); /* x = -2^(emin-1) */
1411 mpfr_set_si (y, -2, MPFR_RNDN); /* y = -2 */
1412 mpfr_clear_flags ();
1413 inex = mpfr_pow (z, x, y, MPFR_RNDN);
1414 if (! mpfr_inf_p (z) || MPFR_IS_NEG (z))
1415 {
1416 printf ("Error in bug20071104: expected +Inf, got ");
1417 mpfr_dump (z);
1418 exit (1);
1419 }
1420 if (inex <= 0)
1421 {
1422 printf ("Error in bug20071104: bad ternary value (%d)\n", inex);
1423 exit (1);
1424 }
1425 if (__gmpfr_flags != (MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT))
1426 {
1427 printf ("Error in bug20071104: bad flags (%u)\n",
1428 (unsigned int) __gmpfr_flags);
1429 exit (1);
1430 }
1431 mpfr_clears (x, y, z, (mpfr_ptr) 0);
1432
1433 set_emin (emin);
1434 set_emax (emax);
1435 }
1436
1437 /* Bug found by Kevin P. Rauch */
1438 static void
bug20071127(void)1439 bug20071127 (void)
1440 {
1441 mpfr_t x, y, z;
1442 int i, tern;
1443 mpfr_exp_t emin, emax;
1444
1445 emin = mpfr_get_emin ();
1446 emax = mpfr_get_emax ();
1447 set_emin (-1000000);
1448 set_emax ( 1000000);
1449
1450 mpfr_init2 (x, 128);
1451 mpfr_init2 (y, 128);
1452 mpfr_init2 (z, 128);
1453
1454 mpfr_set_str (x, "0.80000000000000000000000000000001", 16, MPFR_RNDN);
1455
1456 for (i = 1; i < 9; i *= 2)
1457 {
1458 mpfr_set_str (y, "8000000000000000", 16, MPFR_RNDN);
1459 mpfr_add_si (y, y, i, MPFR_RNDN);
1460 tern = mpfr_pow (z, x, y, MPFR_RNDN);
1461 MPFR_ASSERTN (mpfr_zero_p (z) && MPFR_IS_POS (z) && tern < 0);
1462 }
1463
1464 mpfr_clear (x);
1465 mpfr_clear (y);
1466 mpfr_clear (z);
1467
1468 set_emin (emin);
1469 set_emax (emax);
1470 }
1471
1472 /* Bug found by Kevin P. Rauch */
1473 static void
bug20071128(void)1474 bug20071128 (void)
1475 {
1476 mpfr_t max_val, x, y, z;
1477 int i, tern;
1478 mpfr_exp_t emin, emax;
1479
1480 emin = mpfr_get_emin ();
1481 emax = mpfr_get_emax ();
1482 set_emin (-1000000);
1483 set_emax ( 1000000);
1484
1485 mpfr_init2 (max_val, 64);
1486 mpfr_init2 (x, 64);
1487 mpfr_init2 (y, 64);
1488 mpfr_init2 (z, 64);
1489
1490 mpfr_set_str (max_val, "0.ffffffffffffffff", 16, MPFR_RNDN);
1491 mpfr_set_exp (max_val, mpfr_get_emax ());
1492
1493 mpfr_neg (x, max_val, MPFR_RNDN);
1494
1495 /* on 64-bit machines */
1496 for (i = 41; i < 45; i++)
1497 {
1498 mpfr_set_si_2exp (y, -1, i, MPFR_RNDN);
1499 mpfr_add_si (y, y, 1, MPFR_RNDN);
1500 tern = mpfr_pow (z, x, y, MPFR_RNDN);
1501 MPFR_ASSERTN (mpfr_zero_p (z) && MPFR_IS_NEG (z) && tern > 0);
1502 }
1503
1504 /* on 32-bit machines */
1505 for (i = 9; i < 13; i++)
1506 {
1507 mpfr_set_si_2exp (y, -1, i, MPFR_RNDN);
1508 mpfr_add_si (y, y, 1, MPFR_RNDN);
1509 tern = mpfr_pow (z, x, y, MPFR_RNDN);
1510 MPFR_ASSERTN(mpfr_zero_p (z) && MPFR_IS_NEG (z));
1511 }
1512
1513 mpfr_clear (x);
1514 mpfr_clear (y);
1515 mpfr_clear (z);
1516 mpfr_clear (max_val);
1517
1518 set_emin (emin);
1519 set_emax (emax);
1520 }
1521
1522 /* Bug found by Kevin P. Rauch */
1523 static void
bug20071218(void)1524 bug20071218 (void)
1525 {
1526 mpfr_t x, y, z, t;
1527 int tern;
1528
1529 mpfr_inits2 (64, x, y, z, t, (mpfr_ptr) 0);
1530 mpfr_set_str (x, "0x.80000000000002P-1023", 0, MPFR_RNDN);
1531 mpfr_set_str (y, "100000.000000002", 16, MPFR_RNDN);
1532 mpfr_set_ui (t, 0, MPFR_RNDN);
1533 mpfr_nextabove (t);
1534 tern = mpfr_pow (z, x, y, MPFR_RNDN);
1535 if (mpfr_cmp0 (z, t) != 0)
1536 {
1537 printf ("Error in bug20071218 (1): Expected\n");
1538 mpfr_dump (t);
1539 printf ("Got\n");
1540 mpfr_dump (z);
1541 exit (1);
1542 }
1543 if (tern <= 0)
1544 {
1545 printf ("Error in bug20071218 (1): bad ternary value"
1546 " (%d instead of positive)\n", tern);
1547 exit (1);
1548 }
1549 mpfr_mul_2ui (y, y, 32, MPFR_RNDN);
1550 tern = mpfr_pow (z, x, y, MPFR_RNDN);
1551 if (MPFR_NOTZERO (z) || MPFR_IS_NEG (z))
1552 {
1553 printf ("Error in bug20071218 (2): expected 0, got\n");
1554 mpfr_dump (z);
1555 exit (1);
1556 }
1557 if (tern >= 0)
1558 {
1559 printf ("Error in bug20071218 (2): bad ternary value"
1560 " (%d instead of negative)\n", tern);
1561 exit (1);
1562 }
1563 mpfr_clears (x, y, z, t, (mpfr_ptr) 0);
1564 }
1565
1566 /* With revision 5429, this gives:
1567 * pow.c:43: assertion failed: !mpfr_integer_p (y)
1568 * This is fixed in revision 5432.
1569 */
1570 static void
bug20080721(void)1571 bug20080721 (void)
1572 {
1573 mpfr_t x, y, z, t[2];
1574 int inex;
1575 int rnd;
1576 int err = 0;
1577
1578 /* Note: input values have been chosen in a way to select the
1579 * general case. If mpfr_pow is modified, in particular line
1580 * if (y_is_integer && (MPFR_GET_EXP (y) <= 256))
1581 * make sure that this test still does what we want.
1582 */
1583 mpfr_inits2 (4913, x, y, (mpfr_ptr) 0);
1584 mpfr_inits2 (8, z, t[0], t[1], (mpfr_ptr) 0);
1585 mpfr_set_si (x, -1, MPFR_RNDN);
1586 mpfr_nextbelow (x);
1587 mpfr_set_ui_2exp (y, 1, mpfr_get_prec (y) - 1, MPFR_RNDN);
1588 inex = mpfr_add_ui (y, y, 1, MPFR_RNDN);
1589 MPFR_ASSERTN (inex == 0);
1590 mpfr_set_str_binary (t[0], "-0.10101101e2");
1591 mpfr_set_str_binary (t[1], "-0.10101110e2");
1592 RND_LOOP_NO_RNDF (rnd)
1593 {
1594 int i, inex0;
1595
1596 i = (rnd == MPFR_RNDN || rnd == MPFR_RNDD || rnd == MPFR_RNDA);
1597 inex0 = i ? -1 : 1;
1598 mpfr_clear_flags ();
1599 inex = mpfr_pow (z, x, y, (mpfr_rnd_t) rnd);
1600
1601 if (__gmpfr_flags != MPFR_FLAGS_INEXACT || ! SAME_SIGN (inex, inex0)
1602 || MPFR_IS_NAN (z) || mpfr_cmp (z, t[i]) != 0)
1603 {
1604 unsigned int flags = __gmpfr_flags;
1605
1606 printf ("Error in bug20080721 with %s\n",
1607 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
1608 printf ("expected ");
1609 mpfr_out_str (stdout, 2, 0, t[i], MPFR_RNDN);
1610 printf (", inex = %d, flags = %u\n", inex0,
1611 (unsigned int) MPFR_FLAGS_INEXACT);
1612 printf ("got ");
1613 mpfr_out_str (stdout, 2, 0, z, MPFR_RNDN);
1614 printf (", inex = %d, flags = %u\n", inex, flags);
1615 err = 1;
1616 }
1617 }
1618 mpfr_clears (x, y, z, t[0], t[1], (mpfr_ptr) 0);
1619 if (err)
1620 exit (1);
1621 }
1622
1623 /* The following test fails in r5552 (32-bit and 64-bit). This is due to:
1624 * mpfr_log (t, absx, MPFR_RNDU);
1625 * mpfr_mul (t, y, t, MPFR_RNDU);
1626 * in pow.c, that is supposed to compute an upper bound on exp(y*ln|x|),
1627 * but this is incorrect if y is negative.
1628 */
1629 static void
bug20080820(void)1630 bug20080820 (void)
1631 {
1632 mpfr_exp_t emin;
1633 mpfr_t x, y, z1, z2;
1634
1635 emin = mpfr_get_emin ();
1636 set_emin (MPFR_EMIN_MIN);
1637 mpfr_init2 (x, 80);
1638 mpfr_init2 (y, sizeof (mpfr_exp_t) * CHAR_BIT + 32);
1639 mpfr_init2 (z1, 2);
1640 mpfr_init2 (z2, 80);
1641 mpfr_set_ui (x, 2, MPFR_RNDN);
1642 mpfr_nextbelow (x);
1643 mpfr_set_exp_t (y, mpfr_get_emin () - 2, MPFR_RNDN);
1644 mpfr_nextabove (y);
1645 mpfr_pow (z1, x, y, MPFR_RNDN);
1646 mpfr_pow (z2, x, y, MPFR_RNDN);
1647 /* As x > 0, the rounded value of x^y to nearest in precision p is equal
1648 to 0 iff x^y <= 2^(emin - 2). In particular, this does not depend on
1649 the precision p. Hence the following test. */
1650 if (MPFR_IS_ZERO (z1) && MPFR_NOTZERO (z2))
1651 {
1652 printf ("Error in bug20080820\n");
1653 exit (1);
1654 }
1655 mpfr_clears (x, y, z1, z2, (mpfr_ptr) 0);
1656 set_emin (emin);
1657 }
1658
1659 static void
bug20110320(void)1660 bug20110320 (void)
1661 {
1662 mpfr_exp_t emin;
1663 mpfr_t x, y, z1, z2;
1664 int inex;
1665 unsigned int flags;
1666
1667 emin = mpfr_get_emin ();
1668 set_emin (11);
1669 mpfr_inits2 (2, x, y, z1, z2, (mpfr_ptr) 0);
1670 mpfr_set_ui_2exp (x, 1, 215, MPFR_RNDN);
1671 mpfr_set_ui (y, 1024, MPFR_RNDN);
1672 mpfr_clear_flags ();
1673 inex = mpfr_pow (z1, x, y, MPFR_RNDN);
1674 flags = __gmpfr_flags;
1675 mpfr_set_ui_2exp (z2, 1, 215*1024, MPFR_RNDN);
1676 if (inex != 0 || flags != 0 || ! mpfr_equal_p (z1, z2))
1677 {
1678 printf ("Error in bug20110320\n");
1679 printf ("Expected inex = 0, flags = 0, z = ");
1680 mpfr_dump (z2);
1681 printf ("Got inex = %d, flags = %u, z = ", inex, flags);
1682 mpfr_dump (z1);
1683 exit (1);
1684 }
1685 mpfr_clears (x, y, z1, z2, (mpfr_ptr) 0);
1686 set_emin (emin);
1687 }
1688
1689 static void
tst20140422(void)1690 tst20140422 (void)
1691 {
1692 mpfr_t x, y, z1, z2;
1693 int inex, rnd;
1694 unsigned int flags;
1695
1696 mpfr_inits2 (128, x, y, z1, z2, (mpfr_ptr) 0);
1697 mpfr_set_ui (x, 1296, MPFR_RNDN);
1698 mpfr_set_ui_2exp (y, 3, -2, MPFR_RNDN);
1699 mpfr_set_ui (z2, 216, MPFR_RNDN);
1700 RND_LOOP (rnd)
1701 {
1702 mpfr_clear_flags ();
1703 inex = mpfr_pow (z1, x, y, (mpfr_rnd_t) rnd);
1704 flags = __gmpfr_flags;
1705 if (inex != 0 || flags != 0 || ! mpfr_equal_p (z1, z2))
1706 {
1707 printf ("Error in bug20140422 with %s\n",
1708 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
1709 printf ("Expected inex = 0, flags = 0, z = ");
1710 mpfr_dump (z2);
1711 printf ("Got inex = %d, flags = %u, z = ", inex, flags);
1712 mpfr_dump (z1);
1713 exit (1);
1714 }
1715 }
1716 mpfr_clears (x, y, z1, z2, (mpfr_ptr) 0);
1717 }
1718
1719 static void
coverage(void)1720 coverage (void)
1721 {
1722 mpfr_t x, y, z, t, u;
1723 mpfr_exp_t emin;
1724 int inex;
1725 int i;
1726
1727 /* check a corner case with prec(z)=1 */
1728 mpfr_init2 (x, 2);
1729 mpfr_init2 (y, 8);
1730 mpfr_init2 (z, 1);
1731 mpfr_set_ui_2exp (x, 3, -2, MPFR_RNDN); /* x = 3/4 */
1732 emin = mpfr_get_emin ();
1733 set_emin (-40);
1734 mpfr_set_ui_2exp (y, 199, -1, MPFR_RNDN); /* y = 99.5 */
1735 /* x^y ~ 0.81*2^-41 */
1736 mpfr_clear_underflow ();
1737 inex = mpfr_pow (z, x, y, MPFR_RNDN);
1738 MPFR_ASSERTN(inex > 0);
1739 MPFR_ASSERTN(mpfr_cmp_ui_2exp (z, 1, -41) == 0);
1740 /* there should be no underflow, since with unbounded exponent range,
1741 and a precision of 1 bit, x^y is rounded to 1.0*2^-41 */
1742 MPFR_ASSERTN(!mpfr_underflow_p ());
1743 mpfr_set_ui_2exp (y, 201, -1, MPFR_RNDN); /* y = 100.5 */
1744 /* x^y ~ 0.61*2^-41 */
1745 mpfr_clear_underflow ();
1746 inex = mpfr_pow (z, x, y, MPFR_RNDN);
1747 MPFR_ASSERTN(inex > 0);
1748 MPFR_ASSERTN(mpfr_cmp_ui_2exp (z, 1, -41) == 0);
1749 /* there should be an underflow, since with unbounded exponent range,
1750 and a precision of 1 bit, x^y is rounded to 0.5*2^-41 */
1751 MPFR_ASSERTN(mpfr_underflow_p ());
1752 mpfr_set_ui_2exp (y, 203, -1, MPFR_RNDN); /* y = 101.5 */
1753 /* x^y ~ 0.46*2^-41 */
1754 mpfr_clear_underflow ();
1755 inex = mpfr_pow (z, x, y, MPFR_RNDN);
1756 MPFR_ASSERTN(inex < 0);
1757 MPFR_ASSERTN(mpfr_zero_p (z) && mpfr_signbit (z) == 0);
1758 MPFR_ASSERTN(mpfr_underflow_p ());
1759 mpfr_clears (x, y, z, (mpfr_ptr) 0);
1760 set_emin (emin);
1761
1762 /* test for x = -2, y an odd integer with EXP(y) > i */
1763 mpfr_inits2 (10, t, u, (mpfr_ptr) 0);
1764 mpfr_set_ui_2exp (t, 1, 1UL << 8, MPFR_RNDN);
1765 for (i = 8; i <= 300; i++)
1766 {
1767 mpfr_flags_t flags, flags_u;
1768 int inex_u;
1769
1770 mpfr_mul_si (u, t, -2, MPFR_RNDN); /* u = (-2)^(2^i + 1) */
1771 mpfr_init2 (x, 10);
1772 mpfr_init2 (y, i+1);
1773 mpfr_init2 (z, 10);
1774 mpfr_set_si (x, -2, MPFR_RNDN);
1775 mpfr_set_ui_2exp (y, 1, i, MPFR_RNDN);
1776 mpfr_nextabove (y); /* y = 2^i + 1 */
1777 if (MPFR_IS_INF (u))
1778 {
1779 inex_u = -1;
1780 flags_u = MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT;
1781 }
1782 else
1783 {
1784 inex_u = 0;
1785 flags_u = 0;
1786 }
1787 mpfr_clear_flags ();
1788 inex = mpfr_pow (z, x, y, MPFR_RNDN);
1789 flags = __gmpfr_flags;
1790 if (mpfr_cmp0 (z, u) != 0 ||
1791 ! SAME_SIGN (inex, inex_u) ||
1792 flags != flags_u)
1793 {
1794 printf ("Error in coverage for (-2)^(2^%d + 1):\n", i);
1795 printf ("Expected ");
1796 mpfr_dump (u);
1797 printf (" with inex = %d and flags:", inex_u);
1798 flags_out (flags_u);
1799 printf ("Got ");
1800 mpfr_dump (z);
1801 printf (" with inex = %d and flags:", inex);
1802 flags_out (flags);
1803 exit (1);
1804 }
1805 mpfr_sqr (t, t, MPFR_RNDN);
1806 mpfr_clears (x, y, z, (mpfr_ptr) 0);
1807 }
1808 mpfr_clears (t, u, (mpfr_ptr) 0);
1809
1810 #if _MPFR_EXP_FORMAT >= 3 && _MPFR_PREC_FORMAT == 3 && MPFR_PREC_BITS == 64
1811 /* thus an unsigned long has at least 64 bits and x will be finite */
1812 {
1813 mpfr_exp_t emax;
1814
1815 emax = mpfr_get_emax ();
1816 set_emax ((1UL << 62) - 1);
1817 /* emax = 4611686018427387903 on a 64-bit machine */
1818 mpfr_init2 (x, 65);
1819 mpfr_init2 (y, 65);
1820 mpfr_init2 (z, 64);
1821 mpfr_set_ui_2exp (x, 512, 3074457345618258593UL, MPFR_RNDN);
1822 mpfr_nextbelow (x);
1823 mpfr_set_str_binary (y, "1.1"); /* y = 3/2 */
1824 mpfr_nextabove (y);
1825 inex = mpfr_pow (z, x, y, MPFR_RNDN);
1826 MPFR_ASSERTN(inex > 0);
1827 MPFR_ASSERTN(mpfr_inf_p (z));
1828 mpfr_clears (x, y, z, (mpfr_ptr) 0);
1829 set_emax (emax);
1830 }
1831 #endif
1832 }
1833
1834 static void
check_binary128(void)1835 check_binary128 (void)
1836 {
1837 mpfr_t x, y, z, t;
1838
1839 mpfr_init2 (x, 113);
1840 mpfr_init2 (y, 113);
1841 mpfr_init2 (z, 113);
1842 mpfr_init2 (t, 113);
1843
1844 /* x = 1-2^(-113) */
1845 mpfr_set_ui (x, 1, MPFR_RNDN);
1846 mpfr_nextbelow (x);
1847 /* y = 1.125*2^126 = 9*2^123 */
1848 mpfr_set_ui_2exp (y, 9, 123, MPFR_RNDN);
1849 mpfr_pow (z, x, y, MPFR_RNDN);
1850 /* x^y ~ 3.48e-4003 */
1851 mpfr_set_str (t, "1.16afef53c30899a5c172bb302882p-13296", 16, MPFR_RNDN);
1852 if (! mpfr_equal_p (z, t))
1853 {
1854 printf ("Error in check_binary128\n");
1855 printf ("expected "); mpfr_dump (t);
1856 printf ("got "); mpfr_dump (z);
1857 exit (1);
1858 }
1859
1860 /* x = 5192296858534827628530496329220095/2^112 */
1861 mpfr_set_str (x, "1.fffffffffffffffffffffffffffep-1", 16, MPFR_RNDN);
1862 /* y = -58966440806378323534486035691038613504 */
1863 mpfr_set_str (y, "-1.62e42fefa39ef35793c7673007e5p125", 16, MPFR_RNDN);
1864 mpfr_pow (z, x, y, MPFR_RNDN);
1865 mpfr_set_str (t, "1.fffffffffffffffffffffffff105p16383", 16, MPFR_RNDN);
1866 if (! mpfr_equal_p (z, t))
1867 {
1868 printf ("Error in check_binary128 (2)\n");
1869 printf ("expected "); mpfr_dump (t);
1870 printf ("got "); mpfr_dump (z);
1871 exit (1);
1872 }
1873
1874 mpfr_clear (x);
1875 mpfr_clear (y);
1876 mpfr_clear (z);
1877 mpfr_clear (t);
1878 }
1879
1880 int
main(int argc,char ** argv)1881 main (int argc, char **argv)
1882 {
1883 mpfr_prec_t p;
1884
1885 tests_start_mpfr ();
1886
1887 coverage ();
1888 bug20071127 ();
1889 special ();
1890 particular_cases ();
1891 check_pow_ui ();
1892 check_pow_si ();
1893 check_pown_ieee754_2019 ();
1894 check_special_pow_si ();
1895 pow_si_long_min ();
1896 for (p = MPFR_PREC_MIN; p < 100; p++)
1897 check_inexact (p);
1898 underflows ();
1899 overflows ();
1900 overflows2 ();
1901 overflows3 ();
1902 x_near_one ();
1903 bug20071103 ();
1904 bug20071104 ();
1905 bug20071128 ();
1906 bug20071218 ();
1907 bug20080721 ();
1908 bug20080820 ();
1909 bug20110320 ();
1910 tst20140422 ();
1911 check_binary128 ();
1912
1913 test_generic (MPFR_PREC_MIN, 100, 100);
1914 test_generic_ui (MPFR_PREC_MIN, 100, 100);
1915 test_generic_si (MPFR_PREC_MIN, 100, 100);
1916
1917 data_check ("data/pow275", mpfr_pow275, "mpfr_pow275");
1918
1919 bad_cases (powu2, root2, "mpfr_pow_ui[2]",
1920 8, -256, 255, 4, 128, 800, 40);
1921 bad_cases (pows2, root2, "mpfr_pow_si[2]",
1922 8, -256, 255, 4, 128, 800, 40);
1923 bad_cases (powm2, rootm2, "mpfr_pow_si[-2]",
1924 8, -256, 255, 4, 128, 800, 40);
1925 bad_cases (powu3, root3, "mpfr_pow_ui[3]",
1926 256, -256, 255, 4, 128, 800, 40);
1927 bad_cases (pows3, root3, "mpfr_pow_si[3]",
1928 256, -256, 255, 4, 128, 800, 40);
1929 bad_cases (powm3, rootm3, "mpfr_pow_si[-3]",
1930 256, -256, 255, 4, 128, 800, 40);
1931 bad_cases (powu4, root4, "mpfr_pow_ui[4]",
1932 8, -256, 255, 4, 128, 800, 40);
1933 bad_cases (pows4, root4, "mpfr_pow_si[4]",
1934 8, -256, 255, 4, 128, 800, 40);
1935 bad_cases (powm4, rootm4, "mpfr_pow_si[-4]",
1936 8, -256, 255, 4, 128, 800, 40);
1937 bad_cases (powu5, root5, "mpfr_pow_ui[5]",
1938 256, -256, 255, 4, 128, 800, 40);
1939 bad_cases (pows5, root5, "mpfr_pow_si[5]",
1940 256, -256, 255, 4, 128, 800, 40);
1941 bad_cases (powm5, rootm5, "mpfr_pow_si[-5]",
1942 256, -256, 255, 4, 128, 800, 40);
1943 bad_cases (powu17, root17, "mpfr_pow_ui[17]",
1944 256, -256, 255, 4, 128, 800, 40);
1945 bad_cases (pows17, root17, "mpfr_pow_si[17]",
1946 256, -256, 255, 4, 128, 800, 40);
1947 bad_cases (powm17, rootm17, "mpfr_pow_si[-17]",
1948 256, -256, 255, 4, 128, 800, 40);
1949 bad_cases (powu120, root120, "mpfr_pow_ui[120]",
1950 8, -256, 255, 4, 128, 800, 40);
1951 bad_cases (pows120, root120, "mpfr_pow_si[120]",
1952 8, -256, 255, 4, 128, 800, 40);
1953 bad_cases (powm120, rootm120, "mpfr_pow_si[-120]",
1954 8, -256, 255, 4, 128, 800, 40);
1955
1956 tests_end_mpfr ();
1957 return 0;
1958 }
1959