1 /* Test file for the various power functions
2
3 Copyright 2008-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 /* Note: some tests of the other tpow* test files could be moved there.
24 The main goal of this test file is to test _all_ the power functions
25 on special values, to make sure that they are consistent and give the
26 expected result, in particular because such special cases are handled
27 in different ways in each function. */
28
29 /* Execute with at least an argument to report all the errors found by
30 comparisons. */
31
32 #include "mpfr-test.h"
33
34 /* Behavior of cmpres (called by test_others):
35 * 0: stop as soon as an error is found.
36 * 1: report all errors found by test_others.
37 * -1: the 1 is changed to this value as soon as an error has been found.
38 */
39 static int all_cmpres_errors;
40
41 /* Non-zero when extended exponent range */
42 static int ext = 0;
43
44 static const char *val[] =
45 { "min", "min+", "max", "@NaN@", "-@Inf@", "-4", "-3", "-2", "-1.5",
46 "-1", "-0.5", "-0", "0", "0.5", "1", "1.5", "2", "3", "4", "@Inf@" };
47
48 static void
49 err (const char *s, int i, int j, int rnd, mpfr_srcptr z, int inex)
50 {
51 puts (s);
52 if (ext)
53 puts ("extended exponent range");
54 printf ("x = %s, y = %s, %s\n", val[i], val[j],
55 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
56 printf ("z = ");
57 mpfr_out_str (stdout, 10, 0, z, MPFR_RNDN);
58 printf ("\ninex = %d\n", inex);
59 exit (1);
60 }
61
62 /* Arguments:
63 * spx: non-zero if px is a stringm zero if px is a MPFR number.
64 * px: value of x (string or MPFR number).
65 * sy: value of y (string).
66 * rnd: rounding mode.
67 * z1: expected result (null pointer if unknown pure FP value).
68 * inex1: expected ternary value (if z1 is not a null pointer).
69 * z2: computed result.
70 * inex2: computed ternary value.
71 * flags1: expected flags (computed flags in __gmpfr_flags).
72 * s1, s2: strings about the context.
73 */
74 static void
75 cmpres (int spx, const void *px, const char *sy, mpfr_rnd_t rnd,
76 mpfr_srcptr z1, int inex1, mpfr_srcptr z2, int inex2,
77 unsigned int flags1, const char *s1, const char *s2)
78 {
79 unsigned int flags2 = __gmpfr_flags;
80
81 if (flags1 == flags2)
82 {
83 /* Note: the test on the sign of z1 and z2 is needed
84 in case they are both zeros. */
85 if (z1 == NULL)
86 {
87 if (MPFR_IS_PURE_FP (z2))
88 return;
89 }
90 else if (SAME_SIGN (inex1, inex2) && SAME_VAL (z1, z2))
91 return;
92 }
93
94 printf ("Error in %s\nwith %s%s\nx = ", s1, s2,
95 ext ? ", extended exponent range" : "");
96 if (spx)
97 printf ("%s, ", (char *) px);
98 else
99 {
100 mpfr_out_str (stdout, 16, 0, (mpfr_ptr) px, MPFR_RNDN);
101 puts (",");
102 }
103 printf ("y = %s, %s\n", sy, mpfr_print_rnd_mode (rnd));
104 printf ("Expected ");
105 if (z1 == NULL)
106 {
107 printf ("pure FP value, flags =");
108 flags_out (flags1);
109 }
110 else
111 {
112 mpfr_out_str (stdout, 16, 0, z1, MPFR_RNDN);
113 printf (", inex = %d,\n flags =", VSIGN (inex1));
114 flags_out (flags1);
115 }
116 printf ("Got ");
117 mpfr_out_str (stdout, 16, 0, z2, MPFR_RNDN);
118 printf (", inex = %d,\n flags =", VSIGN (inex2));
119 flags_out (flags2);
120 if (all_cmpres_errors != 0)
121 all_cmpres_errors = -1;
122 else
123 exit (1);
124 }
125
126 static int
127 is_odd (mpfr_srcptr x)
128 {
129 /* works only with the values from val[] */
130 return mpfr_integer_p (x) && mpfr_fits_slong_p (x, MPFR_RNDN) &&
131 (mpfr_get_si (x, MPFR_RNDN) & 1);
132 }
133
134 /* Compare the result (z1,inex1) of mpfr_pow with all flags cleared
135 with those of mpfr_pow with all flags set and of the other power
136 functions. Arguments x and y are the input values; sx and sy are
137 their string representations (sx may be null); rnd contains the
138 rounding mode; s is a string containing the function that called
139 test_others. */
140 static void
141 test_others (const void *sx, const char *sy, mpfr_rnd_t rnd,
142 mpfr_srcptr x, mpfr_srcptr y, mpfr_srcptr z1,
143 int inex1, unsigned int flags, const char *s)
144 {
145 mpfr_t z2;
146 int inex2;
147 int spx = sx != NULL;
148
149 if (!spx)
150 sx = x;
151
152 mpfr_init2 (z2, mpfr_get_prec (z1));
153
154 __gmpfr_flags = MPFR_FLAGS_ALL;
155 inex2 = mpfr_pow (z2, x, y, rnd);
156 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
157 s, "mpfr_pow, flags set");
158
159 /* If y is an integer that fits in an unsigned long and is not -0,
160 we can test mpfr_pow_ui. */
161 if (MPFR_IS_POS (y) && mpfr_integer_p (y) &&
162 mpfr_fits_ulong_p (y, MPFR_RNDN))
163 {
164 unsigned long yy = mpfr_get_ui (y, MPFR_RNDN);
165
166 mpfr_clear_flags ();
167 inex2 = mpfr_pow_ui (z2, x, yy, rnd);
168 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
169 s, "mpfr_pow_ui, flags cleared");
170 __gmpfr_flags = MPFR_FLAGS_ALL;
171 inex2 = mpfr_pow_ui (z2, x, yy, rnd);
172 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
173 s, "mpfr_pow_ui, flags set");
174
175 /* If x is an integer that fits in an unsigned long and is not -0,
176 we can also test mpfr_ui_pow_ui. */
177 if (MPFR_IS_POS (x) && mpfr_integer_p (x) &&
178 mpfr_fits_ulong_p (x, MPFR_RNDN))
179 {
180 unsigned long xx = mpfr_get_ui (x, MPFR_RNDN);
181
182 mpfr_clear_flags ();
183 inex2 = mpfr_ui_pow_ui (z2, xx, yy, rnd);
184 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
185 s, "mpfr_ui_pow_ui, flags cleared");
186 __gmpfr_flags = MPFR_FLAGS_ALL;
187 inex2 = mpfr_ui_pow_ui (z2, xx, yy, rnd);
188 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
189 s, "mpfr_ui_pow_ui, flags set");
190 }
191 }
192
193 /* If y is an integer but not -0 and not huge, we can test mpfr_pow_z,
194 and possibly mpfr_pow_si (and possibly mpfr_ui_div). */
195 if (MPFR_IS_ZERO (y) ? MPFR_IS_POS (y) :
196 (mpfr_integer_p (y) && MPFR_GET_EXP (y) < 256))
197 {
198 mpz_t yyy;
199
200 /* If y fits in a long, we can test mpfr_pow_si. */
201 if (mpfr_fits_slong_p (y, MPFR_RNDN))
202 {
203 long yy = mpfr_get_si (y, MPFR_RNDN);
204
205 mpfr_clear_flags ();
206 inex2 = mpfr_pow_si (z2, x, yy, rnd);
207 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
208 s, "mpfr_pow_si, flags cleared");
209 __gmpfr_flags = MPFR_FLAGS_ALL;
210 inex2 = mpfr_pow_si (z2, x, yy, rnd);
211 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
212 s, "mpfr_pow_si, flags set");
213
214 /* If y = -1, we can test mpfr_ui_div. */
215 if (yy == -1)
216 {
217 mpfr_clear_flags ();
218 inex2 = mpfr_ui_div (z2, 1, x, rnd);
219 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
220 s, "mpfr_ui_div, flags cleared");
221 __gmpfr_flags = MPFR_FLAGS_ALL;
222 inex2 = mpfr_ui_div (z2, 1, x, rnd);
223 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
224 s, "mpfr_ui_div, flags set");
225 }
226
227 /* If y = 2, we can test mpfr_sqr. */
228 if (yy == 2)
229 {
230 mpfr_clear_flags ();
231 inex2 = mpfr_sqr (z2, x, rnd);
232 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
233 s, "mpfr_sqr, flags cleared");
234 __gmpfr_flags = MPFR_FLAGS_ALL;
235 inex2 = mpfr_sqr (z2, x, rnd);
236 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
237 s, "mpfr_sqr, flags set");
238 }
239 }
240
241 /* Test mpfr_pow_z. */
242 mpz_init (yyy);
243 mpfr_get_z (yyy, y, MPFR_RNDN);
244 mpfr_clear_flags ();
245 inex2 = mpfr_pow_z (z2, x, yyy, rnd);
246 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
247 s, "mpfr_pow_z, flags cleared");
248 __gmpfr_flags = MPFR_FLAGS_ALL;
249 inex2 = mpfr_pow_z (z2, x, yyy, rnd);
250 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
251 s, "mpfr_pow_z, flags set");
252 mpz_clear (yyy);
253 }
254
255 /* If y = 0.5, we can test mpfr_sqrt, except if x is -0 or -Inf (because
256 the rule for mpfr_pow on these special values is different). */
257 if (MPFR_IS_PURE_FP (y) && mpfr_cmp_str1 (y, "0.5") == 0 &&
258 ! ((MPFR_IS_ZERO (x) || MPFR_IS_INF (x)) && MPFR_IS_NEG (x)))
259 {
260 mpfr_clear_flags ();
261 inex2 = mpfr_sqrt (z2, x, rnd);
262 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
263 s, "mpfr_sqrt, flags cleared");
264 __gmpfr_flags = MPFR_FLAGS_ALL;
265 inex2 = mpfr_sqrt (z2, x, rnd);
266 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
267 s, "mpfr_sqrt, flags set");
268 }
269
270 /* If y = -0.5, we can test mpfr_rec_sqrt, except if x = -Inf
271 (because the rule for mpfr_pow on -Inf is different). */
272 if (MPFR_IS_PURE_FP (y) && mpfr_cmp_str1 (y, "-0.5") == 0 &&
273 ! (MPFR_IS_INF (x) && MPFR_IS_NEG (x)))
274 {
275 mpfr_clear_flags ();
276 inex2 = mpfr_rec_sqrt (z2, x, rnd);
277 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
278 s, "mpfr_rec_sqrt, flags cleared");
279 __gmpfr_flags = MPFR_FLAGS_ALL;
280 inex2 = mpfr_rec_sqrt (z2, x, rnd);
281 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
282 s, "mpfr_rec_sqrt, flags set");
283 }
284
285 /* If x is an integer that fits in an unsigned long and is not -0,
286 we can test mpfr_ui_pow. */
287 if (MPFR_IS_POS (x) && mpfr_integer_p (x) &&
288 mpfr_fits_ulong_p (x, MPFR_RNDN))
289 {
290 unsigned long xx = mpfr_get_ui (x, MPFR_RNDN);
291
292 mpfr_clear_flags ();
293 inex2 = mpfr_ui_pow (z2, xx, y, rnd);
294 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
295 s, "mpfr_ui_pow, flags cleared");
296 __gmpfr_flags = MPFR_FLAGS_ALL;
297 inex2 = mpfr_ui_pow (z2, xx, y, rnd);
298 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
299 s, "mpfr_ui_pow, flags set");
300
301 /* If x = 2, we can test mpfr_exp2. */
302 if (xx == 2)
303 {
304 mpfr_clear_flags ();
305 inex2 = mpfr_exp2 (z2, y, rnd);
306 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
307 s, "mpfr_exp2, flags cleared");
308 __gmpfr_flags = MPFR_FLAGS_ALL;
309 inex2 = mpfr_exp2 (z2, y, rnd);
310 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
311 s, "mpfr_exp2, flags set");
312 }
313
314 /* If x = 10, we can test mpfr_exp10. */
315 if (xx == 10)
316 {
317 mpfr_clear_flags ();
318 inex2 = mpfr_exp10 (z2, y, rnd);
319 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
320 s, "mpfr_exp10, flags cleared");
321 __gmpfr_flags = MPFR_FLAGS_ALL;
322 inex2 = mpfr_exp10 (z2, y, rnd);
323 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
324 s, "mpfr_exp10, flags set");
325 }
326 }
327
328 mpfr_clear (z2);
329 }
330
331 static int
332 my_setstr (mpfr_ptr t, const char *s)
333 {
334 if (strcmp (s, "min") == 0)
335 {
336 mpfr_setmin (t, mpfr_get_emin ());
337 MPFR_SET_POS (t);
338 return 0;
339 }
340 if (strcmp (s, "min+") == 0)
341 {
342 mpfr_setmin (t, mpfr_get_emin ());
343 MPFR_SET_POS (t);
344 mpfr_nextabove (t);
345 return 0;
346 }
347 if (strcmp (s, "max") == 0)
348 {
349 mpfr_setmax (t, mpfr_get_emax ());
350 MPFR_SET_POS (t);
351 return 0;
352 }
353 return mpfr_set_str (t, s, 10, MPFR_RNDN);
354 }
355
356 static void
357 tst (void)
358 {
359 int sv = numberof (val);
360 int i, j;
361 int rnd;
362 mpfr_t x, y, z, tmp;
363
364 mpfr_inits2 (53, x, y, z, tmp, (mpfr_ptr) 0);
365
366 for (i = 0; i < sv; i++)
367 for (j = 0; j < sv; j++)
368 RND_LOOP_NO_RNDF (rnd)
369 {
370 int exact, inex;
371 unsigned int flags;
372
373 if (my_setstr (x, val[i]) || my_setstr (y, val[j]))
374 {
375 printf ("internal error for (%d,%d,%d)\n", i, j, rnd);
376 exit (1);
377 }
378 mpfr_clear_flags ();
379 inex = mpfr_pow (z, x, y, (mpfr_rnd_t) rnd);
380 flags = __gmpfr_flags;
381 if (! MPFR_IS_NAN (z) && mpfr_nanflag_p ())
382 err ("got NaN flag without NaN value", i, j, rnd, z, inex);
383 if (MPFR_IS_NAN (z) && ! mpfr_nanflag_p ())
384 err ("got NaN value without NaN flag", i, j, rnd, z, inex);
385 if (inex != 0 && ! mpfr_inexflag_p ())
386 err ("got non-zero ternary value without inexact flag",
387 i, j, rnd, z, inex);
388 if (inex == 0 && mpfr_inexflag_p ())
389 err ("got null ternary value with inexact flag",
390 i, j, rnd, z, inex);
391 if (i >= 3 && j >= 3)
392 {
393 if (mpfr_underflow_p ())
394 err ("got underflow", i, j, rnd, z, inex);
395 if (mpfr_overflow_p ())
396 err ("got overflow", i, j, rnd, z, inex);
397 exact = MPFR_IS_SINGULAR (z) ||
398 (mpfr_mul_2ui (tmp, z, 16, MPFR_RNDN), mpfr_integer_p (tmp));
399 if (exact && inex != 0)
400 err ("got exact value with ternary flag different from 0",
401 i, j, rnd, z, inex);
402 if (! exact && inex == 0)
403 err ("got inexact value with ternary flag equal to 0",
404 i, j, rnd, z, inex);
405 }
406 if (MPFR_IS_ZERO (x) && ! MPFR_IS_NAN (y) && MPFR_NOTZERO (y))
407 {
408 if (MPFR_IS_NEG (y) && ! MPFR_IS_INF (z))
409 err ("expected an infinity", i, j, rnd, z, inex);
410 if (MPFR_IS_POS (y) && ! MPFR_IS_ZERO (z))
411 err ("expected a zero", i, j, rnd, z, inex);
412 if ((MPFR_IS_NEG (x) && is_odd (y)) ^ MPFR_IS_NEG (z))
413 err ("wrong sign", i, j, rnd, z, inex);
414 }
415 if (! MPFR_IS_NAN (x) && mpfr_cmp_si (x, -1) == 0)
416 {
417 /* x = -1 */
418 if (! (MPFR_IS_INF (y) || mpfr_integer_p (y)) &&
419 ! MPFR_IS_NAN (z))
420 err ("expected NaN", i, j, rnd, z, inex);
421 if ((MPFR_IS_INF (y) || (mpfr_integer_p (y) && ! is_odd (y)))
422 && ! mpfr_equal_p (z, __gmpfr_one))
423 err ("expected 1", i, j, rnd, z, inex);
424 if (is_odd (y) &&
425 (MPFR_IS_NAN (z) || mpfr_cmp_si (z, -1) != 0))
426 err ("expected -1", i, j, rnd, z, inex);
427 }
428 if ((mpfr_equal_p (x, __gmpfr_one) || MPFR_IS_ZERO (y)) &&
429 ! mpfr_equal_p (z, __gmpfr_one))
430 err ("expected 1", i, j, rnd, z, inex);
431 if (MPFR_IS_PURE_FP (x) && MPFR_IS_NEG (x) &&
432 MPFR_IS_FP (y) && ! mpfr_integer_p (y) &&
433 ! MPFR_IS_NAN (z))
434 err ("expected NaN", i, j, rnd, z, inex);
435 if (MPFR_IS_INF (y) && MPFR_NOTZERO (x))
436 {
437 int cmpabs1 = mpfr_cmpabs (x, __gmpfr_one);
438
439 if ((MPFR_IS_NEG (y) ? (cmpabs1 < 0) : (cmpabs1 > 0)) &&
440 ! (MPFR_IS_POS (z) && MPFR_IS_INF (z)))
441 err ("expected +Inf", i, j, rnd, z, inex);
442 if ((MPFR_IS_NEG (y) ? (cmpabs1 > 0) : (cmpabs1 < 0)) &&
443 ! (MPFR_IS_POS (z) && MPFR_IS_ZERO (z)))
444 err ("expected +0", i, j, rnd, z, inex);
445 }
446 if (MPFR_IS_INF (x) && ! MPFR_IS_NAN (y) && MPFR_NOTZERO (y))
447 {
448 if (MPFR_IS_POS (y) && ! MPFR_IS_INF (z))
449 err ("expected an infinity", i, j, rnd, z, inex);
450 if (MPFR_IS_NEG (y) && ! MPFR_IS_ZERO (z))
451 err ("expected a zero", i, j, rnd, z, inex);
452 if ((MPFR_IS_NEG (x) && is_odd (y)) ^ MPFR_IS_NEG (z))
453 err ("wrong sign", i, j, rnd, z, inex);
454 }
455 test_others (val[i], val[j], (mpfr_rnd_t) rnd, x, y, z, inex, flags,
456 "tst");
457 }
458 mpfr_clears (x, y, z, tmp, (mpfr_ptr) 0);
459 }
460
461 static void
462 underflow_up1 (void)
463 {
464 mpfr_t delta, x, y, z, z0;
465 mpfr_exp_t n;
466 int inex;
467 int rnd;
468 int i;
469
470 n = mpfr_get_emin ();
471 if (n < LONG_MIN)
472 return;
473
474 mpfr_init2 (delta, 2);
475 inex = mpfr_set_ui_2exp (delta, 1, -2, MPFR_RNDN);
476 MPFR_ASSERTN (inex == 0);
477
478 mpfr_init2 (x, 8);
479 inex = mpfr_set_ui (x, 2, MPFR_RNDN);
480 MPFR_ASSERTN (inex == 0);
481
482 mpfr_init2 (y, sizeof (long) * CHAR_BIT + 4);
483 inex = mpfr_set_si (y, n, MPFR_RNDN);
484 MPFR_ASSERTN (inex == 0);
485
486 mpfr_init2 (z0, 2);
487 mpfr_set_ui (z0, 0, MPFR_RNDN);
488
489 mpfr_init2 (z, 32);
490
491 for (i = 0; i <= 12; i++)
492 {
493 unsigned int flags = 0;
494 char sy[256]; /* larger than needed, for maintainability */
495
496 /* Test 2^(emin - i/4).
497 * --> Underflow iff i > 4.
498 * --> Zero in MPFR_RNDN iff i >= 8.
499 */
500
501 if (i != 0 && i != 4)
502 flags |= MPFR_FLAGS_INEXACT;
503 if (i > 4)
504 flags |= MPFR_FLAGS_UNDERFLOW;
505
506 sprintf (sy, "emin - %d/4", i);
507
508 RND_LOOP_NO_RNDF (rnd)
509 {
510 int zero;
511
512 zero = (i > 4 && (rnd == MPFR_RNDZ || rnd == MPFR_RNDD)) ||
513 (i >= 8 && rnd == MPFR_RNDN);
514
515 mpfr_clear_flags ();
516 inex = mpfr_pow (z, x, y, (mpfr_rnd_t) rnd);
517 cmpres (1, "2", sy, (mpfr_rnd_t) rnd, zero ? z0 : (mpfr_ptr) NULL,
518 -1, z, inex, flags, "underflow_up1", "mpfr_pow");
519 test_others ("2", sy, (mpfr_rnd_t) rnd, x, y, z, inex, flags,
520 "underflow_up1");
521 }
522
523 inex = mpfr_sub (y, y, delta, MPFR_RNDN);
524 MPFR_ASSERTN (inex == 0);
525 }
526
527 mpfr_clears (delta, x, y, z, z0, (mpfr_ptr) 0);
528 }
529
530 /* With pow.c r5497, the following test fails on a 64-bit Linux machine
531 * due to a double-rounding problem when rescaling the result:
532 * Error with underflow_up2 and extended exponent range
533 * x = 7.fffffffffffffff0@-1,
534 * y = 4611686018427387904, MPFR_RNDN
535 * Expected 1.0000000000000000@-1152921504606846976, inex = 1, flags = 9
536 * Got 0, inex = -1, flags = 9
537 * With pow_ui.c r5423, the following test fails on a 64-bit Linux machine
538 * as underflows and overflows are not handled correctly (the approximation
539 * error is ignored):
540 * Error with mpfr_pow_ui, flags cleared
541 * x = 7.fffffffffffffff0@-1,
542 * y = 4611686018427387904, MPFR_RNDN
543 * Expected 1.0000000000000000@-1152921504606846976, inex = 1, flags = 9
544 * Got 0, inex = -1, flags = 9
545 */
546 static void
547 underflow_up2 (void)
548 {
549 mpfr_t x, y, z, z0, eps;
550 mpfr_exp_t n;
551 int inex;
552 int rnd;
553
554 n = 1 - mpfr_get_emin ();
555 MPFR_ASSERTN (n > 1);
556 if (n > ULONG_MAX)
557 return;
558
559 mpfr_init2 (eps, 2);
560 mpfr_set_ui_2exp (eps, 1, -1, MPFR_RNDN); /* 1/2 */
561 mpfr_div_ui (eps, eps, n, MPFR_RNDZ); /* 1/(2n) rounded toward zero */
562
563 mpfr_init2 (x, sizeof (unsigned long) * CHAR_BIT + 1);
564 inex = mpfr_ui_sub (x, 1, eps, MPFR_RNDN);
565 MPFR_ASSERTN (inex == 0); /* since n < 2^(size_of_long_in_bits) */
566 inex = mpfr_div_2ui (x, x, 1, MPFR_RNDN); /* 1/2 - eps/2 exactly */
567 MPFR_ASSERTN (inex == 0);
568
569 mpfr_init2 (y, sizeof (unsigned long) * CHAR_BIT);
570 inex = mpfr_set_ui (y, n, MPFR_RNDN);
571 MPFR_ASSERTN (inex == 0);
572
573 /* 0 < eps < 1 / (2n), thus (1 - eps)^n > 1/2,
574 and 1/2 (1/2)^n < (1/2 - eps/2)^n < (1/2)^n. */
575 mpfr_inits2 (64, z, z0, (mpfr_ptr) 0);
576 RND_LOOP_NO_RNDF (rnd)
577 {
578 unsigned int ufinex = MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT;
579 int expected_inex;
580 char sy[256];
581
582 mpfr_set_ui (z0, 0, MPFR_RNDN);
583 expected_inex = rnd == MPFR_RNDN || rnd == MPFR_RNDU || rnd == MPFR_RNDA ?
584 (mpfr_nextabove (z0), 1) : -1;
585 sprintf (sy, "%lu", (unsigned long) n);
586
587 mpfr_clear_flags ();
588 inex = mpfr_pow (z, x, y, (mpfr_rnd_t) rnd);
589 cmpres (0, x, sy, (mpfr_rnd_t) rnd, z0, expected_inex, z, inex, ufinex,
590 "underflow_up2", "mpfr_pow");
591 test_others (NULL, sy, (mpfr_rnd_t) rnd, x, y, z, inex, ufinex,
592 "underflow_up2");
593 }
594
595 mpfr_clears (x, y, z, z0, eps, (mpfr_ptr) 0);
596 }
597
598 static void
599 underflow_up3 (void)
600 {
601 mpfr_t x, y, z, z0;
602 int inex;
603 int rnd;
604 int i;
605
606 mpfr_init2 (x, 64);
607 mpfr_init2 (y, sizeof (mpfr_exp_t) * CHAR_BIT);
608 mpfr_init2 (z, 32);
609 mpfr_init2 (z0, 2);
610
611 inex = mpfr_set_exp_t (y, mpfr_get_emin () - 2, MPFR_RNDN);
612 MPFR_ASSERTN (inex == 0);
613 for (i = -1; i <= 1; i++)
614 RND_LOOP_NO_RNDF (rnd)
615 {
616 unsigned int ufinex = MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT;
617 int expected_inex;
618
619 mpfr_set_ui (x, 2, MPFR_RNDN);
620 if (i < 0)
621 mpfr_nextbelow (x);
622 if (i > 0)
623 mpfr_nextabove (x);
624 /* x = 2 + i * eps, y = emin - 2, x^y ~= 2^(emin - 2) */
625
626 expected_inex = rnd == MPFR_RNDU || rnd == MPFR_RNDA
627 || (rnd == MPFR_RNDN && i < 0) ? 1 : -1;
628
629 mpfr_set_ui (z0, 0, MPFR_RNDN);
630 if (expected_inex > 0)
631 mpfr_nextabove (z0);
632
633 mpfr_clear_flags ();
634 inex = mpfr_pow (z, x, y, (mpfr_rnd_t) rnd);
635 cmpres (0, x, "emin - 2", (mpfr_rnd_t) rnd, z0, expected_inex, z, inex,
636 ufinex, "underflow_up3", "mpfr_pow");
637 test_others (NULL, "emin - 2", (mpfr_rnd_t) rnd, x, y, z, inex, ufinex,
638 "underflow_up3");
639 }
640
641 mpfr_clears (x, y, z, z0, (mpfr_ptr) 0);
642 }
643
644 static void
645 underflow_up (void)
646 {
647 underflow_up1 ();
648 underflow_up2 ();
649 underflow_up3 ();
650 }
651
652 static void
653 overflow_inv (void)
654 {
655 mpfr_t x, y, z;
656 int precx;
657 int s, t;
658 int inex;
659 int rnd;
660
661 mpfr_init2 (y, 2);
662 mpfr_init2 (z, 8);
663
664 mpfr_set_si (y, -1, MPFR_RNDN);
665 for (precx = 10; precx <= 100; precx += 90)
666 {
667 const char *sp = precx == 10 ?
668 "overflow_inv (precx = 10)" : "overflow_inv (precx = 100)";
669
670 mpfr_init2 (x, precx);
671 for (s = -1; s <= 1; s += 2)
672 {
673 inex = mpfr_set_si_2exp (x, s, - mpfr_get_emax (), MPFR_RNDN);
674 MPFR_ASSERTN (inex == 0);
675 for (t = 0; t <= 5; t++)
676 {
677 /* If precx = 10:
678 * x = s * 2^(-emax) * (1 + t * 2^(-9)), so that
679 * 1/x = s * 2^emax * (1 - t * 2^(-9) + eps) with eps > 0.
680 * Values of (1/x) / 2^emax and overflow condition for x > 0:
681 * t = 0: 1 o: always
682 * t = 1: 0.11111111 100000000011... o: MPFR_RNDN and MPFR_RNDU
683 * t = 2: 0.11111111 000000001111... o: MPFR_RNDU
684 * t = 3: 0.11111110 100000100011... o: never
685 *
686 * If precx = 100:
687 * t = 0: always overflow
688 * t > 0: overflow for MPFR_RNDN and MPFR_RNDU.
689 */
690 RND_LOOP_NO_RNDF (rnd)
691 {
692 int inf, overflow;
693 mpfr_rnd_t rnd2;
694
695 if (rnd == MPFR_RNDA)
696 rnd2 = s < 0 ? MPFR_RNDD : MPFR_RNDU;
697 else
698 rnd2 = (mpfr_rnd_t) rnd;
699
700 overflow = t == 0 ||
701 ((mpfr_rnd_t) rnd == MPFR_RNDN &&
702 (precx > 10 || t == 1)) ||
703 (rnd2 == (s < 0 ? MPFR_RNDD : MPFR_RNDU) &&
704 (precx > 10 || t <= 2));
705 inf = overflow &&
706 ((mpfr_rnd_t) rnd == MPFR_RNDN ||
707 rnd2 == (s < 0 ? MPFR_RNDD : MPFR_RNDU));
708 mpfr_clear_flags ();
709 inex = mpfr_pow (z, x, y, (mpfr_rnd_t) rnd);
710 if (overflow ^ !! mpfr_overflow_p ())
711 {
712 printf ("Bad overflow flag in %s\nfor mpfr_pow%s\n"
713 "s = %d, t = %d, %s\n", sp,
714 ext ? ", extended exponent range" : "",
715 s, t, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
716 exit (1);
717 }
718 if (overflow && (inf ^ !! MPFR_IS_INF (z)))
719 {
720 printf ("Bad value in %s\nfor mpfr_pow%s\n"
721 "s = %d, t = %d, %s\nGot ", sp,
722 ext ? ", extended exponent range" : "",
723 s, t, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
724 mpfr_out_str (stdout, 16, 0, z, MPFR_RNDN);
725 printf (" instead of %s value.\n",
726 inf ? "infinite" : "finite");
727 exit (1);
728 }
729 test_others (NULL, "-1", (mpfr_rnd_t) rnd, x, y, z,
730 inex, __gmpfr_flags, sp);
731 } /* RND_LOOP */
732 mpfr_nexttoinf (x);
733 } /* for (t = ...) */
734 } /* for (s = ...) */
735 mpfr_clear (x);
736 } /* for (precx = ...) */
737
738 mpfr_clears (y, z, (mpfr_ptr) 0);
739 }
740
741 static void
742 alltst (void)
743 {
744 mpfr_exp_t emin, emax;
745
746 ext = 0;
747 tst ();
748 underflow_up ();
749 overflow_inv ();
750
751 emin = mpfr_get_emin ();
752 emax = mpfr_get_emax ();
753 set_emin (MPFR_EMIN_MIN);
754 set_emax (MPFR_EMAX_MAX);
755 if (mpfr_get_emin () != emin || mpfr_get_emax () != emax)
756 {
757 ext = 1;
758 tst ();
759 underflow_up ();
760 overflow_inv ();
761 set_emin (emin);
762 set_emax (emax);
763 }
764 }
765
766 int
767 main (int argc, char *argv[])
768 {
769 tests_start_mpfr ();
770 all_cmpres_errors = argc > 1;
771 alltst ();
772 tests_end_mpfr ();
773 return all_cmpres_errors < 0;
774 }