--- /dev/null
+/* Test file for the various power functions
+
+Copyright 2008, 2009 Free Software Foundation, Inc.
+Contributed by the Arenaire and Cacao projects, INRIA.
+
+This file is part of the GNU MPFR Library.
+
+The GNU MPFR Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
+option) any later version.
+
+The GNU MPFR Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MPFR Library; see the file COPYING.LIB. If not, write to
+the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
+MA 02110-1301, USA. */
+
+/* Note: some tests of the other tpow* test files could be moved there.
+ The main goal of this test file is to test _all_ the power functions
+ on special values, to make sure that they are consistent and give the
+ expected result, in particular because such special cases are handled
+ in different ways in each function. */
+
+/* Execute with at least an argument to report all the errors found by
+ comparisons. */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <limits.h>
+
+#include "mpfr-test.h"
+
+/* Behavior of cmpres (called by test_others):
+ * 0: stop as soon as an error is found.
+ * 1: report all errors found by test_others.
+ * -1: the 1 is changed to this value as soon as an error has been found.
+ */
+static int all_cmpres_errors;
+
+/* Non-zero when extended exponent range */
+static int ext = 0;
+
+static char *val[] =
+ { "min", "min+", "max", "@NaN@", "-@Inf@", "-4", "-3", "-2", "-1.5",
+ "-1", "-0.5", "-0", "0", "0.5", "1", "1.5", "2", "3", "4", "@Inf@" };
+
+static void
+err (const char *s, int i, int j, int rnd, mpfr_srcptr z, int inex)
+{
+ puts (s);
+ if (ext)
+ puts ("extended exponent range");
+ printf ("x = %s, y = %s, %s\n", val[i], val[j],
+ mpfr_print_rnd_mode ((mp_rnd_t) rnd));
+ printf ("z = ");
+ mpfr_out_str (stdout, 10, 0, z, GMP_RNDN);
+ printf ("\ninex = %d\n", inex);
+ exit (1);
+}
+
+/* Arguments:
+ * spx: non-zero if px is a stringm zero if px is a MPFR number.
+ * px: value of x (string or MPFR number).
+ * sy: value of y (string).
+ * rnd: rounding mode.
+ * z1: expected result (null pointer if unknown pure FP value).
+ * inex1: expected ternary value (if z1 is not a null pointer).
+ * z2: computed result.
+ * inex2: computed ternary value.
+ * flags1: expected flags (computed flags in __gmpfr_flags).
+ * s1, s2: strings about the context.
+ */
+static void
+cmpres (int spx, const void *px, const char *sy, mp_rnd_t rnd,
+ mpfr_srcptr z1, int inex1, mpfr_srcptr z2, int inex2,
+ unsigned int flags1, const char *s1, const char *s2)
+{
+ unsigned int flags2 = __gmpfr_flags;
+
+ if (flags1 == flags2)
+ {
+ /* Note: the test on the sign of z1 and z2 is needed
+ in case they are both zeros. */
+ if (z1 == NULL)
+ {
+ if (MPFR_IS_PURE_FP (z2))
+ return;
+ }
+ else if (SAME_SIGN (inex1, inex2) &&
+ ((MPFR_IS_NAN (z1) && MPFR_IS_NAN (z2)) ||
+ ((MPFR_IS_NEG (z1) ^ MPFR_IS_NEG (z2)) == 0 &&
+ mpfr_equal_p (z1, z2))))
+ return;
+ }
+
+ printf ("Error in %s\nwith %s%s\nx = ", s1, s2,
+ ext ? ", extended exponent range" : "");
+ if (spx)
+ printf ("%s, ", (char *) px);
+ else
+ {
+ mpfr_out_str (stdout, 16, 0, (mpfr_ptr) px, GMP_RNDN);
+ puts (",");
+ }
+ printf ("y = %s, %s\n", sy, mpfr_print_rnd_mode (rnd));
+ printf ("Expected ");
+ if (z1 == NULL)
+ {
+ printf ("pure FP value, flags = %u\n", flags1);
+ }
+ else
+ {
+ mpfr_out_str (stdout, 16, 0, z1, GMP_RNDN);
+ printf (", inex = %d, flags = %u\n", SIGN (inex1), flags1);
+ }
+ printf ("Got ");
+ mpfr_out_str (stdout, 16, 0, z2, GMP_RNDN);
+ printf (", inex = %d, flags = %u\n", SIGN (inex2), flags2);
+ if (all_cmpres_errors != 0)
+ all_cmpres_errors = -1;
+ else
+ exit (1);
+}
+
+static int
+is_odd (mpfr_srcptr x)
+{
+ /* works only with the values from val[] */
+ return mpfr_integer_p (x) && mpfr_fits_slong_p (x, GMP_RNDN) &&
+ (mpfr_get_si (x, GMP_RNDN) & 1);
+}
+
+/* Compare the result (z1,inex1) of mpfr_pow with all flags cleared
+ with those of mpfr_pow with all flags set and of the other power
+ functions. Arguments x and y are the input values; sx and sy are
+ their string representations (sx may be null); rnd contains the
+ rounding mode; s is a string containing the function that called
+ test_others. */
+static void
+test_others (const void *sx, const char *sy, mp_rnd_t rnd,
+ mpfr_srcptr x, mpfr_srcptr y, mpfr_srcptr z1,
+ int inex1, unsigned int flags, const char *s)
+{
+ mpfr_t z2;
+ int inex2;
+ int spx = sx != NULL;
+
+ if (!spx)
+ sx = x;
+
+ mpfr_init2 (z2, mpfr_get_prec (z1));
+
+ __gmpfr_flags = MPFR_FLAGS_ALL;
+ inex2 = mpfr_pow (z2, x, y, rnd);
+ cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
+ s, "mpfr_pow, flags set");
+
+ /* If y is an integer that fits in an unsigned long and is not -0,
+ we can test mpfr_pow_ui. */
+ if (MPFR_IS_POS (y) && mpfr_integer_p (y) &&
+ mpfr_fits_ulong_p (y, GMP_RNDN))
+ {
+ unsigned long yy = mpfr_get_ui (y, GMP_RNDN);
+
+ mpfr_clear_flags ();
+ inex2 = mpfr_pow_ui (z2, x, yy, rnd);
+ cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
+ s, "mpfr_pow_ui, flags cleared");
+ __gmpfr_flags = MPFR_FLAGS_ALL;
+ inex2 = mpfr_pow_ui (z2, x, yy, rnd);
+ cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
+ s, "mpfr_pow_ui, flags set");
+
+ /* If x is an integer that fits in an unsigned long and is not -0,
+ we can also test mpfr_ui_pow_ui. */
+ if (MPFR_IS_POS (x) && mpfr_integer_p (x) &&
+ mpfr_fits_ulong_p (x, GMP_RNDN))
+ {
+ unsigned long xx = mpfr_get_ui (x, GMP_RNDN);
+
+ mpfr_clear_flags ();
+ inex2 = mpfr_ui_pow_ui (z2, xx, yy, rnd);
+ cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
+ s, "mpfr_ui_pow_ui, flags cleared");
+ __gmpfr_flags = MPFR_FLAGS_ALL;
+ inex2 = mpfr_ui_pow_ui (z2, xx, yy, rnd);
+ cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
+ s, "mpfr_ui_pow_ui, flags set");
+ }
+ }
+
+ /* If y is an integer but not -0 and not huge, we can test mpfr_pow_z,
+ and possibly mpfr_pow_si (and possibly mpfr_ui_div). */
+ if (MPFR_IS_ZERO (y) ? MPFR_IS_POS (y) :
+ (mpfr_integer_p (y) && MPFR_GET_EXP (y) < 256))
+ {
+ mpz_t yyy;
+
+ /* If y fits in a long, we can test mpfr_pow_si. */
+ if (mpfr_fits_slong_p (y, GMP_RNDN))
+ {
+ long yy = mpfr_get_si (y, GMP_RNDN);
+
+ mpfr_clear_flags ();
+ inex2 = mpfr_pow_si (z2, x, yy, rnd);
+ cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
+ s, "mpfr_pow_si, flags cleared");
+ __gmpfr_flags = MPFR_FLAGS_ALL;
+ inex2 = mpfr_pow_si (z2, x, yy, rnd);
+ cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
+ s, "mpfr_pow_si, flags set");
+
+ /* If y = -1, we can test mpfr_ui_div. */
+ if (yy == -1)
+ {
+ mpfr_clear_flags ();
+ inex2 = mpfr_ui_div (z2, 1, x, rnd);
+ cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
+ s, "mpfr_ui_div, flags cleared");
+ __gmpfr_flags = MPFR_FLAGS_ALL;
+ inex2 = mpfr_ui_div (z2, 1, x, rnd);
+ cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
+ s, "mpfr_ui_div, flags set");
+ }
+
+ /* If y = 2, we can test mpfr_sqr. */
+ if (yy == 2)
+ {
+ mpfr_clear_flags ();
+ inex2 = mpfr_sqr (z2, x, rnd);
+ cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
+ s, "mpfr_sqr, flags cleared");
+ __gmpfr_flags = MPFR_FLAGS_ALL;
+ inex2 = mpfr_sqr (z2, x, rnd);
+ cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
+ s, "mpfr_sqr, flags set");
+ }
+ }
+
+ /* Test mpfr_pow_z. */
+ mpz_init (yyy);
+ mpfr_get_z (yyy, y, GMP_RNDN);
+ mpfr_clear_flags ();
+ inex2 = mpfr_pow_z (z2, x, yyy, rnd);
+ cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
+ s, "mpfr_pow_z, flags cleared");
+ __gmpfr_flags = MPFR_FLAGS_ALL;
+ inex2 = mpfr_pow_z (z2, x, yyy, rnd);
+ cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
+ s, "mpfr_pow_z, flags set");
+ mpz_clear (yyy);
+ }
+
+ /* If y = 0.5, we can test mpfr_sqrt, except if x is -0 or -Inf (because
+ the rule for mpfr_pow on these special values is different). */
+ if (MPFR_IS_PURE_FP (y) && mpfr_cmp_str1 (y, "0.5") == 0 &&
+ ! ((MPFR_IS_ZERO (x) || MPFR_IS_INF (x)) && MPFR_IS_NEG (x)))
+ {
+ mpfr_clear_flags ();
+ inex2 = mpfr_sqrt (z2, x, rnd);
+ cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
+ s, "mpfr_sqrt, flags cleared");
+ __gmpfr_flags = MPFR_FLAGS_ALL;
+ inex2 = mpfr_sqrt (z2, x, rnd);
+ cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
+ s, "mpfr_sqrt, flags set");
+ }
+
+#if MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)
+ /* If y = -0.5, we can test mpfr_rec_sqrt, except if x = -Inf
+ (because the rule for mpfr_pow on -Inf is different). */
+ if (MPFR_IS_PURE_FP (y) && mpfr_cmp_str1 (y, "-0.5") == 0 &&
+ ! (MPFR_IS_INF (x) && MPFR_IS_NEG (x)))
+ {
+ mpfr_clear_flags ();
+ inex2 = mpfr_rec_sqrt (z2, x, rnd);
+ cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
+ s, "mpfr_rec_sqrt, flags cleared");
+ __gmpfr_flags = MPFR_FLAGS_ALL;
+ inex2 = mpfr_rec_sqrt (z2, x, rnd);
+ cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
+ s, "mpfr_rec_sqrt, flags set");
+ }
+#endif
+
+ /* If x is an integer that fits in an unsigned long and is not -0,
+ we can test mpfr_ui_pow. */
+ if (MPFR_IS_POS (x) && mpfr_integer_p (x) &&
+ mpfr_fits_ulong_p (x, GMP_RNDN))
+ {
+ unsigned long xx = mpfr_get_ui (x, GMP_RNDN);
+
+ mpfr_clear_flags ();
+ inex2 = mpfr_ui_pow (z2, xx, y, rnd);
+ cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
+ s, "mpfr_ui_pow, flags cleared");
+ __gmpfr_flags = MPFR_FLAGS_ALL;
+ inex2 = mpfr_ui_pow (z2, xx, y, rnd);
+ cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
+ s, "mpfr_ui_pow, flags set");
+
+ /* If x = 2, we can test mpfr_exp2. */
+ if (xx == 2)
+ {
+ mpfr_clear_flags ();
+ inex2 = mpfr_exp2 (z2, y, rnd);
+ cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
+ s, "mpfr_exp2, flags cleared");
+ __gmpfr_flags = MPFR_FLAGS_ALL;
+ inex2 = mpfr_exp2 (z2, y, rnd);
+ cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
+ s, "mpfr_exp2, flags set");
+ }
+
+ /* If x = 10, we can test mpfr_exp10. */
+ if (xx == 10)
+ {
+ mpfr_clear_flags ();
+ inex2 = mpfr_exp10 (z2, y, rnd);
+ cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
+ s, "mpfr_exp10, flags cleared");
+ __gmpfr_flags = MPFR_FLAGS_ALL;
+ inex2 = mpfr_exp10 (z2, y, rnd);
+ cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
+ s, "mpfr_exp10, flags set");
+ }
+ }
+
+ mpfr_clear (z2);
+}
+
+static int
+my_setstr (mpfr_ptr t, const char *s)
+{
+ if (strcmp (s, "min") == 0)
+ {
+ mpfr_setmin (t, mpfr_get_emin ());
+ MPFR_SET_POS (t);
+ return 0;
+ }
+ if (strcmp (s, "min+") == 0)
+ {
+ mpfr_setmin (t, mpfr_get_emin ());
+ MPFR_SET_POS (t);
+ mpfr_nextabove (t);
+ return 0;
+ }
+ if (strcmp (s, "max") == 0)
+ {
+ mpfr_setmax (t, mpfr_get_emax ());
+ MPFR_SET_POS (t);
+ return 0;
+ }
+ return mpfr_set_str (t, s, 10, GMP_RNDN);
+}
+
+static void
+tst (void)
+{
+ int sv = sizeof (val) / sizeof (*val);
+ int i, j;
+ int rnd;
+ mpfr_t x, y, z, tmp;
+
+ mpfr_inits2 (53, x, y, z, tmp, (mpfr_ptr) 0);
+
+ for (i = 0; i < sv; i++)
+ for (j = 0; j < sv; j++)
+ RND_LOOP (rnd)
+ {
+ int exact, inex;
+ unsigned int flags;
+
+ if (my_setstr (x, val[i]) || my_setstr (y, val[j]))
+ {
+ printf ("internal error for (%d,%d,%d)\n", i, j, rnd);
+ exit (1);
+ }
+ mpfr_clear_flags ();
+ inex = mpfr_pow (z, x, y, (mp_rnd_t) rnd);
+ flags = __gmpfr_flags;
+ if (! MPFR_IS_NAN (z) && mpfr_nanflag_p ())
+ err ("got NaN flag without NaN value", i, j, rnd, z, inex);
+ if (MPFR_IS_NAN (z) && ! mpfr_nanflag_p ())
+ err ("got NaN value without NaN flag", i, j, rnd, z, inex);
+ if (inex != 0 && ! mpfr_inexflag_p ())
+ err ("got non-zero ternary value without inexact flag",
+ i, j, rnd, z, inex);
+ if (inex == 0 && mpfr_inexflag_p ())
+ err ("got null ternary value with inexact flag",
+ i, j, rnd, z, inex);
+ if (i >= 3 && j >= 3)
+ {
+ if (mpfr_underflow_p ())
+ err ("got underflow", i, j, rnd, z, inex);
+ if (mpfr_overflow_p ())
+ err ("got overflow", i, j, rnd, z, inex);
+ exact = MPFR_IS_SINGULAR (z) ||
+ (mpfr_mul_2ui (tmp, z, 16, GMP_RNDN), mpfr_integer_p (tmp));
+ if (exact && inex != 0)
+ err ("got exact value with ternary flag different from 0",
+ i, j, rnd, z, inex);
+ if (! exact && inex == 0)
+ err ("got inexact value with ternary flag equal to 0",
+ i, j, rnd, z, inex);
+ }
+ if (MPFR_IS_ZERO (x) && ! MPFR_IS_NAN (y) && MPFR_NOTZERO (y))
+ {
+ if (MPFR_IS_NEG (y) && ! MPFR_IS_INF (z))
+ err ("expected an infinity", i, j, rnd, z, inex);
+ if (MPFR_IS_POS (y) && ! MPFR_IS_ZERO (z))
+ err ("expected a zero", i, j, rnd, z, inex);
+ if ((MPFR_IS_NEG (x) && is_odd (y)) ^ MPFR_IS_NEG (z))
+ err ("wrong sign", i, j, rnd, z, inex);
+ }
+ if (! MPFR_IS_NAN (x) && mpfr_cmp_si (x, -1) == 0)
+ {
+ /* x = -1 */
+ if (! (MPFR_IS_INF (y) || mpfr_integer_p (y)) &&
+ ! MPFR_IS_NAN (z))
+ err ("expected NaN", i, j, rnd, z, inex);
+ if ((MPFR_IS_INF (y) || (mpfr_integer_p (y) && ! is_odd (y)))
+ && ! mpfr_equal_p (z, __gmpfr_one))
+ err ("expected 1", i, j, rnd, z, inex);
+ if (is_odd (y) &&
+ (MPFR_IS_NAN (z) || mpfr_cmp_si (z, -1) != 0))
+ err ("expected -1", i, j, rnd, z, inex);
+ }
+ if ((mpfr_equal_p (x, __gmpfr_one) || MPFR_IS_ZERO (y)) &&
+ ! mpfr_equal_p (z, __gmpfr_one))
+ err ("expected 1", i, j, rnd, z, inex);
+ if (MPFR_IS_PURE_FP (x) && MPFR_IS_NEG (x) &&
+ MPFR_IS_FP (y) && ! mpfr_integer_p (y) &&
+ ! MPFR_IS_NAN (z))
+ err ("expected NaN", i, j, rnd, z, inex);
+ if (MPFR_IS_INF (y) && MPFR_NOTZERO (x))
+ {
+ int cmpabs1 = mpfr_cmpabs (x, __gmpfr_one);
+
+ if ((MPFR_IS_NEG (y) ? (cmpabs1 < 0) : (cmpabs1 > 0)) &&
+ ! (MPFR_IS_POS (z) && MPFR_IS_INF (z)))
+ err ("expected +Inf", i, j, rnd, z, inex);
+ if ((MPFR_IS_NEG (y) ? (cmpabs1 > 0) : (cmpabs1 < 0)) &&
+ ! (MPFR_IS_POS (z) && MPFR_IS_ZERO (z)))
+ err ("expected +0", i, j, rnd, z, inex);
+ }
+ if (MPFR_IS_INF (x) && ! MPFR_IS_NAN (y) && MPFR_NOTZERO (y))
+ {
+ if (MPFR_IS_POS (y) && ! MPFR_IS_INF (z))
+ err ("expected an infinity", i, j, rnd, z, inex);
+ if (MPFR_IS_NEG (y) && ! MPFR_IS_ZERO (z))
+ err ("expected a zero", i, j, rnd, z, inex);
+ if ((MPFR_IS_NEG (x) && is_odd (y)) ^ MPFR_IS_NEG (z))
+ err ("wrong sign", i, j, rnd, z, inex);
+ }
+ test_others (val[i], val[j], (mp_rnd_t) rnd, x, y, z, inex, flags,
+ "tst");
+ }
+ mpfr_clears (x, y, z, tmp, (mpfr_ptr) 0);
+}
+
+static void
+underflow_up1 (void)
+{
+ mpfr_t delta, x, y, z, z0;
+ mp_exp_t n;
+ int inex;
+ int rnd;
+ int i;
+
+ n = mpfr_get_emin ();
+ if (n < LONG_MIN)
+ return;
+
+ mpfr_init2 (delta, 2);
+ inex = mpfr_set_ui_2exp (delta, 1, -2, GMP_RNDN);
+ MPFR_ASSERTN (inex == 0);
+
+ mpfr_init2 (x, 8);
+ inex = mpfr_set_ui (x, 2, GMP_RNDN);
+ MPFR_ASSERTN (inex == 0);
+
+ mpfr_init2 (y, sizeof (long) * CHAR_BIT + 4);
+ inex = mpfr_set_si (y, n, GMP_RNDN);
+ MPFR_ASSERTN (inex == 0);
+
+ mpfr_init2 (z0, 2);
+ mpfr_set_ui (z0, 0, GMP_RNDN);
+
+ mpfr_init2 (z, 32);
+
+ for (i = 0; i <= 12; i++)
+ {
+ unsigned int flags = 0;
+ char sy[16];
+
+ /* Test 2^(emin - i/4).
+ * --> Underflow iff i > 4.
+ * --> Zero in GMP_RNDN iff i >= 8.
+ */
+
+ if (i != 0 && i != 4)
+ flags |= MPFR_FLAGS_INEXACT;
+ if (i > 4)
+ flags |= MPFR_FLAGS_UNDERFLOW;
+
+ sprintf (sy, "emin - %d/4", i);
+
+ RND_LOOP (rnd)
+ {
+ int zero;
+
+ zero = (i > 4 && (rnd == GMP_RNDZ || rnd == GMP_RNDD)) ||
+ (i >= 8 && rnd == GMP_RNDN);
+
+ mpfr_clear_flags ();
+ inex = mpfr_pow (z, x, y, (mp_rnd_t) rnd);
+ cmpres (1, "2", sy, (mp_rnd_t) rnd, zero ? z0 : (mpfr_ptr) NULL,
+ -1, z, inex, flags, "underflow_up1", "mpfr_pow");
+ test_others ("2", sy, (mp_rnd_t) rnd, x, y, z, inex, flags,
+ "underflow_up1");
+ }
+
+ inex = mpfr_sub (y, y, delta, GMP_RNDN);
+ MPFR_ASSERTN (inex == 0);
+ }
+
+ mpfr_clears (delta, x, y, z, z0, (mpfr_ptr) 0);
+}
+
+/* With pow.c r5497, the following test fails on a 64-bit Linux machine
+ * due to a double-rounding problem when rescaling the result:
+ * Error with underflow_up2 and extended exponent range
+ * x = 7.fffffffffffffff0@-1,
+ * y = 4611686018427387904, GMP_RNDN
+ * Expected 1.0000000000000000@-1152921504606846976, inex = 1, flags = 9
+ * Got 0, inex = -1, flags = 9
+ * With pow_ui.c r5423, the following test fails on a 64-bit Linux machine
+ * as underflows and overflows are not handled correctly (the approximation
+ * error is ignored):
+ * Error with mpfr_pow_ui, flags cleared
+ * x = 7.fffffffffffffff0@-1,
+ * y = 4611686018427387904, GMP_RNDN
+ * Expected 1.0000000000000000@-1152921504606846976, inex = 1, flags = 9
+ * Got 0, inex = -1, flags = 9
+ */
+static void
+underflow_up2 (void)
+{
+ mpfr_t x, y, z, z0, eps;
+ mp_exp_t n;
+ int inex;
+ int rnd;
+
+ n = 1 - mpfr_get_emin ();
+ MPFR_ASSERTN (n > 1);
+ if (n > ULONG_MAX)
+ return;
+
+ mpfr_init2 (eps, 2);
+ mpfr_set_ui_2exp (eps, 1, -1, GMP_RNDN); /* 1/2 */
+ mpfr_div_ui (eps, eps, n, GMP_RNDZ); /* 1/(2n) rounded toward zero */
+
+ mpfr_init2 (x, sizeof (unsigned long) * CHAR_BIT + 1);
+ inex = mpfr_ui_sub (x, 1, eps, GMP_RNDN);
+ MPFR_ASSERTN (inex == 0); /* since n < 2^(size_of_long_in_bits) */
+ inex = mpfr_div_2ui (x, x, 1, GMP_RNDN); /* 1/2 - eps/2 exactly */
+ MPFR_ASSERTN (inex == 0);
+
+ mpfr_init2 (y, sizeof (unsigned long) * CHAR_BIT);
+ inex = mpfr_set_ui (y, n, GMP_RNDN);
+ MPFR_ASSERTN (inex == 0);
+
+ /* 0 < eps < 1 / (2n), thus (1 - eps)^n > 1/2,
+ and 1/2 (1/2)^n < (1/2 - eps/2)^n < (1/2)^n. */
+ mpfr_inits2 (64, z, z0, (mpfr_ptr) 0);
+ RND_LOOP (rnd)
+ {
+ unsigned int ufinex = MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT;
+ int expected_inex;
+ char sy[256];
+
+ mpfr_set_ui (z0, 0, GMP_RNDN);
+ expected_inex = rnd == GMP_RNDN || rnd == GMP_RNDU ?
+ (mpfr_nextabove (z0), 1) : -1;
+ sprintf (sy, "%lu", (unsigned long) n);
+
+ mpfr_clear_flags ();
+ inex = mpfr_pow (z, x, y, (mp_rnd_t) rnd);
+ cmpres (0, x, sy, (mp_rnd_t) rnd, z0, expected_inex, z, inex, ufinex,
+ "underflow_up2", "mpfr_pow");
+ test_others (NULL, sy, (mp_rnd_t) rnd, x, y, z, inex, ufinex,
+ "underflow_up2");
+ }
+
+ mpfr_clears (x, y, z, z0, eps, (mpfr_ptr) 0);
+}
+
+static void
+underflow_up3 (void)
+{
+ mpfr_t x, y, z, z0;
+ int inex;
+ int rnd;
+ int i;
+
+ mpfr_init2 (x, 64);
+ mpfr_init2 (y, sizeof (mp_exp_t) * CHAR_BIT);
+ mpfr_init2 (z, 32);
+ mpfr_init2 (z0, 2);
+
+ inex = mpfr_set_exp_t (y, mpfr_get_emin () - 2, GMP_RNDN);
+ MPFR_ASSERTN (inex == 0);
+ for (i = -1; i <= 1; i++)
+ RND_LOOP (rnd)
+ {
+ unsigned int ufinex = MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT;
+ int expected_inex;
+
+ mpfr_set_ui (x, 2, GMP_RNDN);
+ if (i < 0)
+ mpfr_nextbelow (x);
+ if (i > 0)
+ mpfr_nextabove (x);
+ /* x = 2 + i * eps, y = emin - 2, x^y ~= 2^(emin - 2) */
+
+ expected_inex = rnd == GMP_RNDU || (rnd == GMP_RNDN && i < 0) ?
+ 1 : -1;
+
+ mpfr_set_ui (z0, 0, GMP_RNDN);
+ if (expected_inex > 0)
+ mpfr_nextabove (z0);
+
+ mpfr_clear_flags ();
+ inex = mpfr_pow (z, x, y, (mp_rnd_t) rnd);
+ cmpres (0, x, "emin - 2", (mp_rnd_t) rnd, z0, expected_inex, z, inex,
+ ufinex, "underflow_up3", "mpfr_pow");
+ test_others (NULL, "emin - 2", (mp_rnd_t) rnd, x, y, z, inex, ufinex,
+ "underflow_up3");
+ }
+
+ mpfr_clears (x, y, z, z0, (mpfr_ptr) 0);
+}
+
+static void
+underflow_up (void)
+{
+ underflow_up1 ();
+ underflow_up2 ();
+ underflow_up3 ();
+}
+
+static void
+overflow_inv (void)
+{
+ mpfr_t x, y, z;
+ int precx;
+ int s, t;
+ int inex;
+ int rnd;
+
+ mpfr_init2 (y, 2);
+ mpfr_init2 (z, 8);
+
+ mpfr_set_si (y, -1, GMP_RNDN);
+ for (precx = 10; precx <= 100; precx += 90)
+ {
+ const char *sp = precx == 10 ?
+ "overflow_inv (precx = 10)" : "overflow_inv (precx = 100)";
+
+ mpfr_init2 (x, precx);
+ for (s = -1; s <= 1; s += 2)
+ {
+ inex = mpfr_set_si_2exp (x, s, - mpfr_get_emax (), GMP_RNDN);
+ MPFR_ASSERTN (inex == 0);
+ for (t = 0; t <= 5; t++)
+ {
+ /* If precx = 10:
+ * x = s * 2^(-emax) * (1 + t * 2^(-9)), so that
+ * 1/x = s * 2^emax * (1 - t * 2^(-9) + eps) with eps > 0.
+ * Values of (1/x) / 2^emax and overflow condition for x > 0:
+ * t = 0: 1 o: always
+ * t = 1: 0.11111111 100000000011... o: GMP_RNDN and GMP_RNDU
+ * t = 2: 0.11111111 000000001111... o: GMP_RNDU
+ * t = 3: 0.11111110 100000100011... o: never
+ *
+ * If precx = 100:
+ * t = 0: always overflow
+ * t > 0: overflow for GMP_RNDN and GMP_RNDU.
+ */
+ RND_LOOP (rnd)
+ {
+ int inf, overflow;
+
+ overflow = t == 0 ||
+ ((mp_rnd_t) rnd == GMP_RNDN && (precx > 10 || t == 1)) ||
+ ((mp_rnd_t) rnd == (s < 0 ? GMP_RNDD : GMP_RNDU) &&
+ (precx > 10 || t <= 2));
+ inf = overflow &&
+ ((mp_rnd_t) rnd == GMP_RNDN ||
+ (mp_rnd_t) rnd == (s < 0 ? GMP_RNDD : GMP_RNDU));
+ mpfr_clear_flags ();
+ inex = mpfr_pow (z, x, y, (mp_rnd_t) rnd);
+ if (overflow ^ !! mpfr_overflow_p ())
+ {
+ printf ("Bad overflow flag in %s\nfor mpfr_pow%s\n"
+ "s = %d, t = %d, %s\n", sp,
+ ext ? ", extended exponent range" : "",
+ s, t, mpfr_print_rnd_mode ((mp_rnd_t) rnd));
+ exit (1);
+ }
+ if (overflow && (inf ^ !! MPFR_IS_INF (z)))
+ {
+ printf ("Bad value in %s\nfor mpfr_pow%s\n"
+ "s = %d, t = %d, %s\nGot ", sp,
+ ext ? ", extended exponent range" : "",
+ s, t, mpfr_print_rnd_mode ((mp_rnd_t) rnd));
+ mpfr_out_str (stdout, 16, 0, z, GMP_RNDN);
+ printf (" instead of %s value.\n",
+ inf ? "infinite" : "finite");
+ exit (1);
+ }
+ test_others (NULL, "-1", (mp_rnd_t) rnd, x, y, z,
+ inex, __gmpfr_flags, sp);
+ } /* RND_LOOP */
+ mpfr_nexttoinf (x);
+ } /* for (t = ...) */
+ } /* for (s = ...) */
+ mpfr_clear (x);
+ } /* for (precx = ...) */
+
+ mpfr_clears (y, z, (mpfr_ptr) 0);
+}
+
+static void
+alltst (void)
+{
+ mp_exp_t emin, emax;
+
+ ext = 0;
+ tst ();
+ underflow_up ();
+ overflow_inv ();
+
+ emin = mpfr_get_emin ();
+ emax = mpfr_get_emax ();
+ set_emin (MPFR_EMIN_MIN);
+ set_emax (MPFR_EMAX_MAX);
+ if (mpfr_get_emin () != emin || mpfr_get_emax () != emax)
+ {
+ ext = 1;
+ tst ();
+ underflow_up ();
+ overflow_inv ();
+ set_emin (emin);
+ set_emax (emax);
+ }
+}
+
+int
+main (int argc, char *argv[])
+{
+ tests_start_mpfr ();
+ all_cmpres_errors = argc > 1;
+ alltst ();
+ tests_end_mpfr ();
+ return all_cmpres_errors < 0;
+}
projects / msp430-gcc.git / blobdiff