--- /dev/null
+/* Test file for mpfr_exp.
+
+Copyright 1999, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 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. */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <limits.h>
+
+#include "mpfr-test.h"
+
+#ifdef CHECK_EXTERNAL
+static int
+test_exp (mpfr_ptr a, mpfr_srcptr b, mp_rnd_t rnd_mode)
+{
+ int res;
+ int ok = rnd_mode == GMP_RNDN && mpfr_number_p (b) && mpfr_get_prec (a)>=53;
+ if (ok)
+ {
+ mpfr_print_raw (b);
+ }
+ res = mpfr_exp (a, b, rnd_mode);
+ if (ok)
+ {
+ printf (" ");
+ mpfr_print_raw (a);
+ printf ("\n");
+ }
+ return res;
+}
+#else
+#define test_exp mpfr_exp
+#endif
+
+/* returns the number of ulp of error */
+static void
+check3 (const char *op, mp_rnd_t rnd, const char *res)
+{
+ mpfr_t x, y;
+
+ mpfr_inits2 (53, x, y, (mpfr_ptr) 0);
+ /* y negative. If we forget to set the sign in mpfr_exp, we'll see it. */
+ mpfr_set_si (y, -1, GMP_RNDN);
+ mpfr_set_str1 (x, op);
+ test_exp (y, x, rnd);
+ if (mpfr_cmp_str1 (y, res) )
+ {
+ printf ("mpfr_exp failed for x=%s, rnd=%s\n",
+ op, mpfr_print_rnd_mode (rnd));
+ printf ("expected result is %s, got ", res);
+ mpfr_out_str (stdout, 10, 0, y, GMP_RNDN);
+ putchar('\n');
+ exit (1);
+ }
+ mpfr_clears (x, y, (mpfr_ptr) 0);
+}
+
+/* expx is the value of exp(X) rounded towards -infinity */
+static void
+check_worst_case (const char *Xs, const char *expxs)
+{
+ mpfr_t x, y;
+
+ mpfr_inits2 (53, x, y, (mpfr_ptr) 0);
+ mpfr_set_str1(x, Xs);
+ test_exp(y, x, GMP_RNDD);
+ if (mpfr_cmp_str1 (y, expxs))
+ {
+ printf ("exp(x) rounded towards -infinity is wrong\n");
+ exit(1);
+ }
+ mpfr_set_str1(x, Xs);
+ test_exp(x, x, GMP_RNDU);
+ mpfr_nexttoinf (y);
+ if (mpfr_cmp(x,y))
+ {
+ printf ("exp(x) rounded towards +infinity is wrong\n");
+ exit(1);
+ }
+ mpfr_clears (x, y, (mpfr_ptr) 0);
+}
+
+/* worst cases communicated by Jean-Michel Muller and Vincent Lefevre */
+static int
+check_worst_cases (void)
+{
+ mpfr_t x; mpfr_t y;
+
+ mpfr_init(x);
+ mpfr_set_prec (x, 53);
+
+ check_worst_case("4.44089209850062517562e-16", "1.00000000000000022204");
+ check_worst_case("6.39488462184069720009e-14", "1.00000000000006372680");
+ check_worst_case("1.84741111297455401935e-12", "1.00000000000184718907");
+ check_worst_case("1.76177628026265550074e-10", "1.00000000017617751702");
+ check3("1.76177628026265550074e-10", GMP_RNDN, "1.00000000017617773906");
+ check_worst_case("7.54175277499595900852e-10", "1.00000000075417516676");
+ check3("7.54175277499595900852e-10", GMP_RNDN, "1.00000000075417538881");
+ /* bug found by Vincent Lefe`vre on December 8, 1999 */
+ check3("-5.42410311287441459172e+02", GMP_RNDN, "2.7176584868845723e-236");
+ /* further cases communicated by Vincent Lefe`vre on January 27, 2000 */
+ check3("-1.32920285897904911589e-10", GMP_RNDN, "0.999999999867079769622");
+ check3("-1.44037948245738330735e-10", GMP_RNDN, "0.9999999998559621072757");
+ check3("-1.66795910430705305937e-10", GMP_RNDZ, "0.9999999998332040895832");
+ check3("-1.64310953745426656203e-10", GMP_RNDN, "0.9999999998356891017792");
+ check3("-1.38323574826034659172e-10", GMP_RNDZ, "0.9999999998616764251835");
+ check3("-1.23621668465115401498e-10", GMP_RNDZ, "0.9999999998763783315425");
+
+ mpfr_set_prec (x, 601);
+ mpfr_set_str (x, "0.88b6ba510e10450edc258748bc9dfdd466f21b47ed264cdf24aa8f64af1f3fad9ec2301d43c0743f534b5aa20091ff6d352df458ef1ba519811ef6f5b11853534fd8fa32764a0a6d2d0dd20@0", 16, GMP_RNDZ);
+ mpfr_init2 (y, 601);
+ mpfr_exp_2 (y, x, GMP_RNDD);
+ mpfr_exp_3 (x, x, GMP_RNDD);
+ if (mpfr_cmp (x, y))
+ {
+ printf ("mpfr_exp_2 and mpfr_exp_3 differ for prec=601\n");
+ printf ("mpfr_exp_2 gives ");
+ mpfr_out_str (stdout, 2, 0, y, GMP_RNDN);
+ printf ("\nmpfr_exp_3 gives ");
+ mpfr_out_str (stdout, 2, 0, x, GMP_RNDN);
+ printf ("\n");
+ exit (1);
+ }
+
+ mpfr_set_prec (x, 13001);
+ mpfr_set_prec (y, 13001);
+ mpfr_urandomb (x, RANDS);
+ mpfr_exp_3 (y, x, GMP_RNDN);
+ mpfr_exp_2 (x, x, GMP_RNDN);
+ if (mpfr_cmp (x, y))
+ {
+ printf ("mpfr_exp_2 and mpfr_exp_3 differ for prec=13001\n");
+ exit (1);
+ }
+
+ mpfr_clear (x);
+ mpfr_clear (y);
+ return 0;
+}
+
+static void
+compare_exp2_exp3 (mp_prec_t p0, mp_prec_t p1)
+{
+ mpfr_t x, y, z;
+ mp_prec_t prec;
+ mp_rnd_t rnd;
+
+ mpfr_init (x);
+ mpfr_init (y);
+ mpfr_init (z);
+ for (prec = p0; prec <= p1; prec ++)
+ {
+ mpfr_set_prec (x, prec);
+ mpfr_set_prec (y, prec);
+ mpfr_set_prec (z, prec);
+ mpfr_urandomb (x, RANDS);
+ rnd = RND_RAND ();
+ mpfr_exp_2 (y, x, rnd);
+ mpfr_exp_3 (z, x, rnd);
+ if (mpfr_cmp (y,z))
+ {
+ printf ("mpfr_exp_2 and mpfr_exp_3 disagree for rnd=%s and\nx=",
+ mpfr_print_rnd_mode (rnd));
+ mpfr_print_binary (x);
+ puts ("");
+ printf ("mpfr_exp_2 gives ");
+ mpfr_print_binary (y);
+ puts ("");
+ printf ("mpfr_exp_3 gives ");
+ mpfr_print_binary (z);
+ puts ("");
+ exit (1);
+ }
+ }
+
+ mpfr_clear (x);
+ mpfr_clear (y);
+ mpfr_clear (z);
+}
+
+static void
+check_large (void)
+{
+ mpfr_t x, z;
+ mp_prec_t prec;
+
+ /* bug found by Patrick Pe'lissier on 7 Jun 2004 */
+ prec = 203780;
+ mpfr_init2 (x, prec);
+ mpfr_init2 (z, prec);
+ mpfr_set_ui (x, 3, GMP_RNDN);
+ mpfr_sqrt (x, x, GMP_RNDN);
+ mpfr_sub_ui (x, x, 1, GMP_RNDN);
+ mpfr_exp_3 (z, x, GMP_RNDN);
+ mpfr_clear (x);
+ mpfr_clear (z);
+}
+
+#define TEST_FUNCTION test_exp
+#define TEST_RANDOM_EMIN -36
+#define TEST_RANDOM_EMAX 36
+#include "tgeneric.c"
+
+static void
+check_special (void)
+{
+ mpfr_t x, y, z;
+ mp_exp_t emin, emax;
+
+ emin = mpfr_get_emin ();
+ emax = mpfr_get_emax ();
+
+ mpfr_init (x);
+ mpfr_init (y);
+ mpfr_init (z);
+
+ /* check exp(NaN) = NaN */
+ mpfr_set_nan (x);
+ test_exp (y, x, GMP_RNDN);
+ if (!mpfr_nan_p (y))
+ {
+ printf ("Error for exp(NaN)\n");
+ exit (1);
+ }
+
+ /* check exp(+inf) = +inf */
+ mpfr_set_inf (x, 1);
+ test_exp (y, x, GMP_RNDN);
+ if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0)
+ {
+ printf ("Error for exp(+inf)\n");
+ exit (1);
+ }
+
+ /* check exp(-inf) = +0 */
+ mpfr_set_inf (x, -1);
+ test_exp (y, x, GMP_RNDN);
+ if (mpfr_cmp_ui (y, 0) || mpfr_sgn (y) < 0)
+ {
+ printf ("Error for exp(-inf)\n");
+ exit (1);
+ }
+
+ /* Check overflow. Corner case of mpfr_exp_2 */
+ mpfr_set_prec (x, 64);
+ mpfr_set_emax (MPFR_EMAX_DEFAULT);
+ mpfr_set_emin (MPFR_EMIN_DEFAULT);
+ mpfr_set_str (x,
+ "0.1011000101110010000101111111010100001100000001110001100111001101E30",
+ 2, GMP_RNDN);
+ mpfr_exp (x, x, GMP_RNDD);
+ if (mpfr_cmp_str (x,
+".1111111111111111111111111111111111111111111111111111111111111111E1073741823",
+ 2, GMP_RNDN) != 0)
+ {
+ printf ("Wrong overflow detection in mpfr_exp\n");
+ mpfr_dump (x);
+ exit (1);
+ }
+ /* Check underflow. Corner case of mpfr_exp_2 */
+ mpfr_set_str (x,
+"-0.1011000101110010000101111111011111010001110011110111100110101100E30",
+ 2, GMP_RNDN);
+ mpfr_exp (x, x, GMP_RNDN);
+ if (mpfr_cmp_str (x, "0.1E-1073741823", 2, GMP_RNDN) != 0)
+ {
+ printf ("Wrong underflow (1) detection in mpfr_exp\n");
+ mpfr_dump (x);
+ exit (1);
+ }
+ mpfr_set_str (x,
+"-0.1011001101110010000101111111011111010001110011110111100110111101E30",
+ 2, GMP_RNDN);
+ mpfr_exp (x, x, GMP_RNDN);
+ if (mpfr_cmp_ui (x, 0) != 0)
+ {
+ printf ("Wrong underflow (2) detection in mpfr_exp\n");
+ mpfr_dump (x);
+ exit (1);
+ }
+ /* Check overflow. Corner case of mpfr_exp_3 */
+ if (MPFR_PREC_MAX > MPFR_EXP_THRESHOLD + 10)
+ {
+ mpfr_set_prec (x, MPFR_EXP_THRESHOLD + 10);
+ mpfr_set_str (x,
+ "0.1011000101110010000101111111010100001100000001110001100111001101E30",
+ 2, GMP_RNDN);
+ mpfr_clear_overflow ();
+ mpfr_exp (x, x, GMP_RNDD);
+ if (!mpfr_overflow_p ())
+ {
+ printf ("Wrong overflow detection in mpfr_exp_3\n");
+ mpfr_dump (x);
+ exit (1);
+ }
+ /* Check underflow. Corner case of mpfr_exp_3 */
+ mpfr_set_str (x,
+ "-0.1011000101110010000101111111011111010001110011110111100110101100E30",
+ 2, GMP_RNDN);
+ mpfr_clear_underflow ();
+ mpfr_exp (x, x, GMP_RNDN);
+ if (!mpfr_underflow_p ())
+ {
+ printf ("Wrong underflow detection in mpfr_exp_3\n");
+ mpfr_dump (x);
+ exit (1);
+ }
+ mpfr_set_prec (x, 53);
+ }
+
+ /* check overflow */
+ set_emax (10);
+ mpfr_set_ui (x, 7, GMP_RNDN);
+ test_exp (y, x, GMP_RNDN);
+ if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0)
+ {
+ printf ("Error for exp(7) for emax=10\n");
+ exit (1);
+ }
+ set_emax (emax);
+
+ /* check underflow */
+ set_emin (-10);
+ mpfr_set_si (x, -9, GMP_RNDN);
+ test_exp (y, x, GMP_RNDN);
+ if (mpfr_cmp_ui (y, 0) || mpfr_sgn (y) < 0)
+ {
+ printf ("Error for exp(-9) for emin=-10\n");
+ printf ("Expected +0\n");
+ printf ("Got "); mpfr_print_binary (y); puts ("");
+ exit (1);
+ }
+ set_emin (emin);
+
+ /* check case EXP(x) < -precy */
+ mpfr_set_prec (y, 2);
+ mpfr_set_str_binary (x, "-0.1E-3");
+ test_exp (y, x, GMP_RNDD);
+ if (mpfr_cmp_ui_2exp (y, 3, -2))
+ {
+ printf ("Error for exp(-1/16), prec=2, RNDD\n");
+ printf ("expected 0.11, got ");
+ mpfr_dump (y);
+ exit (1);
+ }
+ test_exp (y, x, GMP_RNDZ);
+ if (mpfr_cmp_ui_2exp (y, 3, -2))
+ {
+ printf ("Error for exp(-1/16), prec=2, RNDZ\n");
+ printf ("expected 0.11, got ");
+ mpfr_dump (y);
+ exit (1);
+ }
+ mpfr_set_str_binary (x, "0.1E-3");
+ test_exp (y, x, GMP_RNDN);
+ if (mpfr_cmp_ui (y, 1))
+ {
+ printf ("Error for exp(1/16), prec=2, RNDN\n");
+ exit (1);
+ }
+ test_exp (y, x, GMP_RNDU);
+ if (mpfr_cmp_ui_2exp (y, 3, -1))
+ {
+ printf ("Error for exp(1/16), prec=2, RNDU\n");
+ exit (1);
+ }
+
+ /* bug reported by Franky Backeljauw, 28 Mar 2003 */
+ mpfr_set_prec (x, 53);
+ mpfr_set_prec (y, 53);
+ mpfr_set_str_binary (x, "1.1101011000111101011110000111010010101001101001110111e28");
+ test_exp (y, x, GMP_RNDN);
+
+ mpfr_set_prec (x, 153);
+ mpfr_set_prec (z, 153);
+ mpfr_set_str_binary (x, "1.1101011000111101011110000111010010101001101001110111e28");
+ test_exp (z, x, GMP_RNDN);
+ mpfr_prec_round (z, 53, GMP_RNDN);
+
+ if (mpfr_cmp (y, z))
+ {
+ printf ("Error in mpfr_exp for large argument\n");
+ exit (1);
+ }
+
+ /* corner cases in mpfr_exp_3 */
+ mpfr_set_prec (x, 2);
+ mpfr_set_ui (x, 1, GMP_RNDN);
+ mpfr_set_prec (y, 2);
+ mpfr_exp_3 (y, x, GMP_RNDN);
+
+ /* Check some little things about overflow detection */
+ set_emin (-125);
+ set_emax (128);
+ mpfr_set_prec (x, 107);
+ mpfr_set_prec (y, 107);
+ mpfr_set_str_binary (x, "0.11110000000000000000000000000000000000000000000"
+ "0000000000000000000000000000000000000000000000000000"
+ "00000000E4");
+ test_exp (y, x, GMP_RNDN);
+ if (mpfr_cmp_str (y, "0.11000111100001100110010101111101011010010101010000"
+ "1101110111100010111001011111111000110111001011001101010"
+ "01E22", 2, GMP_RNDN))
+ {
+ printf ("Special overflow error (1)\n");
+ mpfr_dump (y);
+ exit (1);
+ }
+
+ set_emin (emin);
+ set_emax (emax);
+
+ /* Check for overflow producing a segfault with HUGE exponent */
+ mpfr_set_ui (x, 3, GMP_RNDN);
+ mpfr_mul_2ui (x, x, 32, GMP_RNDN);
+ test_exp (y, x, GMP_RNDN); /* Can't test return value: May overflow or not*/
+
+ /* Bug due to wrong approximation of (x)/log2 */
+ mpfr_set_prec (x, 163);
+
+ mpfr_set_str (x, "-4.28ac8fceeadcda06bb56359017b1c81b85b392e7", 16,
+ GMP_RNDN);
+ mpfr_exp (x, x, GMP_RNDN);
+ if (mpfr_cmp_str (x, "3.fffffffffffffffffffffffffffffffffffffffe8@-2",
+ 16, GMP_RNDN))
+ {
+ printf ("Error for x= -4.28ac8fceeadcda06bb56359017b1c81b85b392e7");
+ printf ("expected 3.fffffffffffffffffffffffffffffffffffffffe8@-2");
+ printf ("Got ");
+ mpfr_out_str (stdout, 16, 0, x, GMP_RNDN); putchar ('\n');
+ }
+
+ /* bug found by Guillaume Melquiond, 13 Sep 2005 */
+ mpfr_set_prec (x, 53);
+ mpfr_set_str_binary (x, "-1E-400");
+ mpfr_exp (x, x, GMP_RNDZ);
+ if (mpfr_cmp_ui (x, 1) == 0)
+ {
+ printf ("Error for exp(-2^(-400))\n");
+ exit (1);
+ }
+
+ mpfr_clear (x);
+ mpfr_clear (y);
+ mpfr_clear (z);
+}
+
+/* check sign of inexact flag */
+static void
+check_inexact (void)
+{
+ mpfr_t x, y;
+ int inexact;
+
+ mpfr_init2 (x, 53);
+ mpfr_init2 (y, 53);
+
+ mpfr_set_str_binary (x,
+ "1.0000000000001001000110100100101000001101101011100101e2");
+ inexact = test_exp (y, x, GMP_RNDN);
+ if (inexact <= 0)
+ {
+ printf ("Wrong inexact flag (Got %d instead of 1)\n", inexact);
+ exit (1);
+ }
+
+ mpfr_clear (x);
+ mpfr_clear (y);
+}
+
+static void
+check_exp10(void)
+{
+ mpfr_t x;
+ int inexact;
+
+ mpfr_init2 (x, 200);
+ mpfr_set_ui(x, 4, GMP_RNDN);
+
+ inexact = mpfr_exp10 (x, x, GMP_RNDN);
+ if (mpfr_cmp_ui(x, 10*10*10*10))
+ {
+ printf ("exp10: Wrong returned value\n");
+ exit (1);
+ }
+ if (inexact != 0)
+ {
+ printf ("exp10: Wrong inexact flag\n");
+ exit (1);
+ }
+
+ mpfr_clear (x);
+}
+
+static void
+overflowed_exp0 (void)
+{
+ mpfr_t x, y;
+ int emax, i, inex, rnd, err = 0;
+ mp_exp_t old_emax;
+
+ old_emax = mpfr_get_emax ();
+
+ mpfr_init2 (x, 8);
+ mpfr_init2 (y, 8);
+
+ for (emax = -1; emax <= 0; emax++)
+ {
+ mpfr_set_ui_2exp (y, 1, emax, GMP_RNDN);
+ mpfr_nextbelow (y);
+ set_emax (emax); /* 1 is not representable. */
+ /* and if emax < 0, 1 - eps is not representable either. */
+ for (i = -1; i <= 1; i++)
+ RND_LOOP (rnd)
+ {
+ mpfr_set_si_2exp (x, i, -512 * ABS (i), GMP_RNDN);
+ mpfr_clear_flags ();
+ inex = mpfr_exp (x, x, (mp_rnd_t) rnd);
+ if ((i >= 0 || emax < 0 || rnd == GMP_RNDN || rnd == GMP_RNDU) &&
+ ! mpfr_overflow_p ())
+ {
+ printf ("Error in overflowed_exp0 (i = %d, rnd = %s):\n"
+ " The overflow flag is not set.\n",
+ i, mpfr_print_rnd_mode ((mp_rnd_t) rnd));
+ err = 1;
+ }
+ if (rnd == GMP_RNDZ || rnd == GMP_RNDD)
+ {
+ if (inex >= 0)
+ {
+ printf ("Error in overflowed_exp0 (i = %d, rnd = %s):\n"
+ " The inexact value must be negative.\n",
+ i, mpfr_print_rnd_mode ((mp_rnd_t) rnd));
+ err = 1;
+ }
+ if (! mpfr_equal_p (x, y))
+ {
+ printf ("Error in overflowed_exp0 (i = %d, rnd = %s):\n"
+ " Got ", i, mpfr_print_rnd_mode ((mp_rnd_t) rnd));
+ mpfr_print_binary (x);
+ printf (" instead of 0.11111111E%d.\n", emax);
+ err = 1;
+ }
+ }
+ else
+ {
+ if (inex <= 0)
+ {
+ printf ("Error in overflowed_exp0 (i = %d, rnd = %s):\n"
+ " The inexact value must be positive.\n",
+ i, mpfr_print_rnd_mode ((mp_rnd_t) rnd));
+ err = 1;
+ }
+ if (! (mpfr_inf_p (x) && MPFR_SIGN (x) > 0))
+ {
+ printf ("Error in overflowed_exp0 (i = %d, rnd = %s):\n"
+ " Got ", i, mpfr_print_rnd_mode ((mp_rnd_t) rnd));
+ mpfr_print_binary (x);
+ printf (" instead of +Inf.\n");
+ err = 1;
+ }
+ }
+ }
+ set_emax (old_emax);
+ }
+
+ if (err)
+ exit (1);
+ mpfr_clear (x);
+ mpfr_clear (y);
+}
+
+/* This bug occurs in mpfr_exp_2 on a Linux-64 machine, r5475. */
+static void
+bug20080731 (void)
+{
+ mp_exp_t emin;
+ mpfr_t x, y1, y2;
+ mp_prec_t prec = 64;
+
+ emin = mpfr_get_emin ();
+ set_emin (MPFR_EMIN_MIN);
+
+ mpfr_init2 (x, 200);
+ mpfr_set_str (x, "-2.c5c85fdf473de6af278ece700fcbdabd03cd0cb9ca62d8b62c@7",
+ 16, GMP_RNDN);
+
+ mpfr_init2 (y1, prec);
+ mpfr_exp (y1, x, GMP_RNDU);
+
+ /* Compute the result with a higher internal precision. */
+ mpfr_init2 (y2, 300);
+ mpfr_exp (y2, x, GMP_RNDU);
+ mpfr_prec_round (y2, prec, GMP_RNDU);
+
+ if (mpfr_cmp0 (y1, y2) != 0)
+ {
+ printf ("Error in bug20080731\nExpected ");
+ mpfr_out_str (stdout, 16, 0, y2, GMP_RNDN);
+ printf ("\nGot ");
+ mpfr_out_str (stdout, 16, 0, y1, GMP_RNDN);
+ printf ("\n");
+ exit (1);
+ }
+
+ mpfr_clears (x, y1, y2, (mpfr_ptr) 0);
+ set_emin (emin);
+}
+
+/* Emulate mpfr_exp with mpfr_exp_3 in the general case. */
+static int
+exp_3 (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
+{
+ int inexact;
+
+ inexact = mpfr_exp_3 (y, x, rnd_mode);
+ return mpfr_check_range (y, inexact, rnd_mode);
+}
+
+static void
+underflow_up (int extended_emin)
+{
+ mpfr_t minpos, x, y, t, t2;
+ int precx, precy;
+ int inex;
+ int rnd;
+ int e3;
+ int i, j;
+
+ mpfr_init2 (minpos, 2);
+ mpfr_set_ui (minpos, 0, GMP_RNDN);
+ mpfr_nextabove (minpos);
+
+ /* Let's test values near the underflow boundary.
+ *
+ * Minimum representable positive number: minpos = 2^(emin - 1).
+ * Let's choose an MPFR number x = log(minpos) + eps, with |eps| small
+ * (note: eps cannot be 0, and cannot be a rational number either).
+ * Then exp(x) = minpos * exp(eps) ~= minpos * (1 + eps + eps^2).
+ * We will compute y = rnd(exp(x)) in some rounding mode, precision p.
+ * 1. If eps > 0, then in any rounding mode:
+ * rnd(exp(x)) >= minpos and no underflow.
+ * So, let's take x1 = rndu(log(minpos)) in some precision.
+ * 2. If eps < 0, then exp(x) < minpos and the result will be either 0
+ * or minpos. An underflow always occurs in GMP_RNDZ and GMP_RNDD,
+ * but not necessarily in GMP_RNDN and GMP_RNDU (this is underflow
+ * after rounding in an unbounded exponent range). If -a < eps < -b,
+ * minpos * (1 - a) < exp(x) < minpos * (1 - b + b^2).
+ * - If eps > -2^(-p), no underflow in GMP_RNDU.
+ * - If eps > -2^(-p-1), no underflow in GMP_RNDN.
+ * - If eps < - (2^(-p-1) + 2^(-2p-1)), underflow in GMP_RNDN.
+ * - If eps < - (2^(-p) + 2^(-2p+1)), underflow in GMP_RNDU.
+ * - In GMP_RNDN, result is minpos iff exp(eps) > 1/2, i.e.
+ * - log(2) < eps < ...
+ *
+ * Moreover, since precy < MPFR_EXP_THRESHOLD (to avoid tests that take
+ * too much time), mpfr_exp() always selects mpfr_exp_2(); so, we need
+ * to test mpfr_exp_3() too. This will be done via the e3 variable:
+ * e3 = 0: mpfr_exp(), thus mpfr_exp_2().
+ * e3 = 1: mpfr_exp_3(), via the exp_3() wrapper.
+ * i.e.: inex = e3 ? exp_3 (y, x, rnd) : mpfr_exp (y, x, rnd);
+ */
+
+ /* Case eps > 0. In revision 5461 (trunk) on a 64-bit Linux machine:
+ * Incorrect flags in underflow_up, eps > 0, GMP_RNDN and extended emin
+ * for precx = 96, precy = 16, mpfr_exp_3
+ * Got 9 instead of 8.
+ * Note: testing this case in several precisions for x and y introduces
+ * some useful random. Indeed, the bug is not always triggered.
+ * Fixed in r5469.
+ */
+ for (precx = 16; precx <= 128; precx += 16)
+ {
+ mpfr_init2 (x, precx);
+ mpfr_log (x, minpos, GMP_RNDU);
+ for (precy = 16; precy <= 128; precy += 16)
+ {
+ mpfr_init2 (y, precy);
+
+ for (e3 = 0; e3 <= 1; e3++)
+ {
+ RND_LOOP (rnd)
+ {
+ int err = 0;
+
+ mpfr_clear_flags ();
+ inex = e3 ? exp_3 (y, x, (mpfr_rnd_t) rnd)
+ : mpfr_exp (y, x, (mpfr_rnd_t) rnd);
+ if (__gmpfr_flags != MPFR_FLAGS_INEXACT)
+ {
+ printf ("Incorrect flags in underflow_up, eps > 0, %s",
+ mpfr_print_rnd_mode ((mp_rnd_t) rnd));
+ if (extended_emin)
+ printf (" and extended emin");
+ printf ("\nfor precx = %d, precy = %d, %s\n",
+ precx, precy, e3 ? "mpfr_exp_3" : "mpfr_exp");
+ printf ("Got %u instead of %u.\n", __gmpfr_flags,
+ (unsigned int) MPFR_FLAGS_INEXACT);
+ err = 1;
+ }
+ if (mpfr_cmp0 (y, minpos) < 0)
+ {
+ printf ("Incorrect result in underflow_up, eps > 0, %s",
+ mpfr_print_rnd_mode ((mp_rnd_t) rnd));
+ if (extended_emin)
+ printf (" and extended emin");
+ printf ("\nfor precx = %d, precy = %d, %s\n",
+ precx, precy, e3 ? "mpfr_exp_3" : "mpfr_exp");
+ mpfr_dump (y);
+ err = 1;
+ }
+ MPFR_ASSERTN (inex != 0);
+ if (rnd == GMP_RNDD || rnd == GMP_RNDZ)
+ MPFR_ASSERTN (inex < 0);
+ if (rnd == GMP_RNDU)
+ MPFR_ASSERTN (inex > 0);
+ if (err)
+ exit (1);
+ }
+ }
+
+ mpfr_clear (y);
+ }
+ mpfr_clear (x);
+ }
+
+ /* Case - log(2) < eps < 0 in GMP_RNDN, starting with small-precision x;
+ * only check the result and the ternary value.
+ * Previous to r5453 (trunk), on 32-bit and 64-bit machines, this fails
+ * for precx = 65 and precy = 16, e.g.:
+ * exp_2.c:264: assertion failed: ...
+ * because mpfr_sub (r, x, r, GMP_RNDU); yields a null value. This is
+ * fixed in r5453 by going to next Ziv's iteration.
+ */
+ for (precx = sizeof(mp_exp_t) * CHAR_BIT + 1; precx <= 81; precx += 8)
+ {
+ mpfr_init2 (x, precx);
+ mpfr_log (x, minpos, GMP_RNDD); /* |ulp| <= 1/2 */
+ for (precy = 16; precy <= 128; precy += 16)
+ {
+ mpfr_init2 (y, precy);
+ inex = mpfr_exp (y, x, GMP_RNDN);
+ if (inex <= 0 || mpfr_cmp0 (y, minpos) != 0)
+ {
+ printf ("Error in underflow_up, - log(2) < eps < 0");
+ if (extended_emin)
+ printf (" and extended emin");
+ printf (" for prec = %d\nExpected ", precy);
+ mpfr_out_str (stdout, 16, 0, minpos, GMP_RNDN);
+ printf (" (minimum positive MPFR number) and inex > 0\nGot ");
+ mpfr_out_str (stdout, 16, 0, y, GMP_RNDN);
+ printf ("\nwith inex = %d\n", inex);
+ exit (1);
+ }
+ mpfr_clear (y);
+ }
+ mpfr_clear (x);
+ }
+
+ /* Cases eps ~ -2^(-p) and eps ~ -2^(-p-1). More precisely,
+ * _ for j = 0, eps > -2^(-(p+i)),
+ * _ for j = 1, eps < - (2^(-(p+i)) + 2^(1-2(p+i))),
+ * where i = 0 or 1.
+ */
+ mpfr_inits2 (2, t, t2, (mpfr_ptr) 0);
+ for (precy = 16; precy <= 128; precy += 16)
+ {
+ mpfr_set_ui_2exp (t, 1, - precy, GMP_RNDN); /* 2^(-p) */
+ mpfr_set_ui_2exp (t2, 1, 1 - 2 * precy, GMP_RNDN); /* 2^(-2p+1) */
+ precx = sizeof(mp_exp_t) * CHAR_BIT + 2 * precy + 8;
+ mpfr_init2 (x, precx);
+ mpfr_init2 (y, precy);
+ for (i = 0; i <= 1; i++)
+ {
+ for (j = 0; j <= 1; j++)
+ {
+ if (j == 0)
+ {
+ /* Case eps > -2^(-(p+i)). */
+ mpfr_log (x, minpos, GMP_RNDU);
+ }
+ else /* j == 1 */
+ {
+ /* Case eps < - (2^(-(p+i)) + 2^(1-2(p+i))). */
+ mpfr_log (x, minpos, GMP_RNDD);
+ inex = mpfr_sub (x, x, t2, GMP_RNDN);
+ MPFR_ASSERTN (inex == 0);
+ }
+ inex = mpfr_sub (x, x, t, GMP_RNDN);
+ MPFR_ASSERTN (inex == 0);
+
+ RND_LOOP (rnd)
+ for (e3 = 0; e3 <= 1; e3++)
+ {
+ int err = 0;
+ unsigned int flags;
+
+ flags = MPFR_FLAGS_INEXACT |
+ ((rnd == GMP_RNDU && (i == 1 || j == 0)) ||
+ (rnd == GMP_RNDN && (i == 1 && j == 0)) ?
+ 0 : MPFR_FLAGS_UNDERFLOW);
+ mpfr_clear_flags ();
+ inex = e3 ? exp_3 (y, x, (mpfr_rnd_t) rnd)
+ : mpfr_exp (y, x, (mpfr_rnd_t) rnd);
+ if (__gmpfr_flags != flags)
+ {
+ printf ("Incorrect flags in underflow_up, %s",
+ mpfr_print_rnd_mode ((mp_rnd_t) rnd));
+ if (extended_emin)
+ printf (" and extended emin");
+ printf ("\nfor precx = %d, precy = %d, ",
+ precx, precy);
+ if (j == 0)
+ printf ("eps >~ -2^(-%d)", precy + i);
+ else
+ printf ("eps <~ - (2^(-%d) + 2^(%d))",
+ precy + i, 1 - 2 * (precy + i));
+ printf (", %s\n", e3 ? "mpfr_exp_3" : "mpfr_exp");
+ printf ("Got %u instead of %u.\n",
+ __gmpfr_flags, flags);
+ err = 1;
+ }
+ if (rnd == GMP_RNDU || rnd == GMP_RNDN ?
+ mpfr_cmp0 (y, minpos) != 0 : MPFR_NOTZERO (y))
+ {
+ printf ("Incorrect result in underflow_up, %s",
+ mpfr_print_rnd_mode ((mp_rnd_t) rnd));
+ if (extended_emin)
+ printf (" and extended emin");
+ printf ("\nfor precx = %d, precy = %d, ",
+ precx, precy);
+ if (j == 0)
+ printf ("eps >~ -2^(-%d)", precy + i);
+ else
+ printf ("eps <~ - (2^(-%d) + 2^(%d))",
+ precy + i, 1 - 2 * (precy + i));
+ printf (", %s\n", e3 ? "mpfr_exp_3" : "mpfr_exp");
+ mpfr_dump (y);
+ err = 1;
+ }
+ if (err)
+ exit (1);
+ } /* for (e3 ...) */
+ } /* for (j ...) */
+ mpfr_div_2si (t, t, 1, GMP_RNDN);
+ mpfr_div_2si (t2, t2, 2, GMP_RNDN);
+ } /* for (i ...) */
+ mpfr_clears (x, y, (mpfr_ptr) 0);
+ } /* for (precy ...) */
+ mpfr_clears (t, t2, (mpfr_ptr) 0);
+
+ /* Case exp(eps) ~= 1/2, i.e. eps ~= - log(2).
+ * We choose x0 and x1 with high enough precision such that:
+ * x0 = rndd(rndd(log(minpos)) - rndu(log(2)))
+ * x1 = rndu(rndu(log(minpos)) - rndd(log(2)))
+ * In revision 5507 (trunk) on a 64-bit Linux machine, this fails:
+ * Error in underflow_up, eps >~ - log(2) and extended emin
+ * for precy = 16, mpfr_exp
+ * Expected 1.0@-1152921504606846976 (minimum positive MPFR number),
+ * inex > 0 and flags = 9
+ * Got 0
+ * with inex = -1 and flags = 9
+ * due to a double-rounding problem in mpfr_mul_2si when rescaling
+ * the result.
+ */
+ mpfr_inits2 (sizeof(mp_exp_t) * CHAR_BIT + 64, x, t, (mpfr_ptr) 0);
+ for (i = 0; i <= 1; i++)
+ {
+ mpfr_log (x, minpos, i ? GMP_RNDU : GMP_RNDD);
+ mpfr_const_log2 (t, i ? GMP_RNDD : GMP_RNDU);
+ mpfr_sub (x, x, t, i ? GMP_RNDU : GMP_RNDD);
+ for (precy = 16; precy <= 128; precy += 16)
+ {
+ mpfr_init2 (y, precy);
+ for (e3 = 0; e3 <= 1; e3++)
+ {
+ unsigned int flags, uflags =
+ MPFR_FLAGS_INEXACT | MPFR_FLAGS_UNDERFLOW;
+
+ mpfr_clear_flags ();
+ inex = e3 ? exp_3 (y, x, GMP_RNDN) : mpfr_exp (y, x, GMP_RNDN);
+ flags = __gmpfr_flags;
+ if (flags != uflags ||
+ (i ? (inex <= 0 || mpfr_cmp0 (y, minpos) != 0)
+ : (inex >= 0 || MPFR_NOTZERO (y))))
+ {
+ printf ("Error in underflow_up, eps %c~ - log(2)",
+ i ? '>' : '<');
+ if (extended_emin)
+ printf (" and extended emin");
+ printf ("\nfor precy = %d, %s\nExpected ", precy,
+ e3 ? "mpfr_exp_3" : "mpfr_exp");
+ if (i)
+ {
+ mpfr_out_str (stdout, 16, 0, minpos, GMP_RNDN);
+ printf (" (minimum positive MPFR number),\ninex >");
+ }
+ else
+ {
+ printf ("+0, inex <");
+ }
+ printf (" 0 and flags = %u\nGot ", uflags);
+ mpfr_out_str (stdout, 16, 0, y, GMP_RNDN);
+ printf ("\nwith inex = %d and flags = %u\n", inex, flags);
+ exit (1);
+ }
+ }
+ mpfr_clear (y);
+ }
+ }
+ mpfr_clears (x, t, (mpfr_ptr) 0);
+
+ mpfr_clear (minpos);
+}
+
+static void
+underflow (void)
+{
+ mp_exp_t emin;
+
+ underflow_up (0);
+
+ emin = mpfr_get_emin ();
+ set_emin (MPFR_EMIN_MIN);
+ if (mpfr_get_emin () != emin)
+ {
+ underflow_up (1);
+ set_emin (emin);
+ }
+}
+
+int
+main (int argc, char *argv[])
+{
+ tests_start_mpfr ();
+
+ if (argc > 1)
+ check_large ();
+
+ check_inexact ();
+ check_special ();
+
+ test_generic (2, 100, 100);
+
+ compare_exp2_exp3 (20, 1000);
+ check_worst_cases();
+ check3("0.0", GMP_RNDU, "1.0");
+ check3("-1e-170", GMP_RNDU, "1.0");
+ check3("-1e-170", GMP_RNDN, "1.0");
+ check3("-8.88024741073346941839e-17", GMP_RNDU, "1.0");
+ check3("8.70772839244701057915e-01", GMP_RNDN, "2.38875626491680437269");
+ check3("1.0", GMP_RNDN, "2.71828182845904509080");
+ check3("-3.42135637628104173534e-07", GMP_RNDZ, "0.999999657864420798958");
+ /* worst case for argument reduction, very near from 5*log(2),
+ thanks to Jean-Michel Muller */
+ check3("3.4657359027997265421", GMP_RNDN, "32.0");
+ check3("3.4657359027997265421", GMP_RNDU, "32.0");
+ check3("3.4657359027997265421", GMP_RNDD, "31.999999999999996447");
+ check3("2.26523754332090625496e+01", GMP_RNDD, "6.8833785261699581146e9");
+ check3("1.31478962104089092122e+01", GMP_RNDZ, "5.12930793917860137299e+05");
+ check3("4.25637507920002378103e-01", GMP_RNDU, "1.53056585656161181497e+00");
+ check3("6.26551618962329307459e-16", GMP_RNDU, "1.00000000000000066613e+00");
+ check3("-3.35589513871216568383e-03",GMP_RNDD, "9.96649729583626853291e-01");
+ check3("1.95151388850007272424e+01", GMP_RNDU, "2.98756340674767792225e+08");
+ check3("2.45045953503350730784e+01", GMP_RNDN, "4.38743344916128387451e+10");
+ check3("2.58165606081678085104e+01", GMP_RNDD, "1.62925781879432281494e+11");
+ check3("-2.36539020084338638128e+01",GMP_RNDZ, "5.33630792749924762447e-11");
+ check3("2.39211946135858077866e+01", GMP_RNDU, "2.44817704330214385986e+10");
+ check3("-2.78190533055889162029e+01",GMP_RNDZ, "8.2858803483596879512e-13");
+ check3("2.64028186174889789584e+01", GMP_RNDD, "2.9281844652878973388e11");
+ check3("2.92086338843268329413e+01", GMP_RNDZ, "4.8433797301907177734e12");
+ check3("-2.46355324071459982349e+01",GMP_RNDZ, "1.9995129297760994791e-11");
+ check3("-2.23509444608605427618e+01",GMP_RNDZ, "1.9638492867489702307e-10");
+ check3("-2.41175390197331687148e+01",GMP_RNDD, "3.3564940885530624592e-11");
+ check3("2.46363885231578088053e+01", GMP_RNDU, "5.0055014282693267822e10");
+ check3("111.1263531080090984914932050742208957672119140625", GMP_RNDN, "1.8262572323517295459e48");
+ check3("-3.56196340354684821250e+02",GMP_RNDN, "2.0225297096141478156e-155");
+ check3("6.59678273772710895173e+02", GMP_RNDU, "3.1234469273830195529e286");
+ check3("5.13772529701934331570e+02", GMP_RNDD, "1.3445427121297197752e223");
+ check3("3.57430211008718345056e+02", GMP_RNDD, "1.6981197246857298443e155");
+ check3("3.82001814471465536371e+02", GMP_RNDU, "7.9667300591087367805e165");
+ check3("5.92396038219384422518e+02", GMP_RNDD, "1.880747529554661989e257");
+ check3("-5.02678550462488090034e+02",GMP_RNDU, "4.8919201895446217839e-219");
+ check3("5.30015757134837031117e+02", GMP_RNDD, "1.5237672861171573939e230");
+ check3("5.16239362447650933063e+02", GMP_RNDZ, "1.5845518406744492105e224");
+ check3("6.00812634798592370977e-01", GMP_RNDN, "1.823600119339019443");
+ check_exp10 ();
+
+ bug20080731 ();
+
+ overflowed_exp0 ();
+ underflow ();
+
+ data_check ("data/exp", mpfr_exp, "mpfr_exp");
+ bad_cases (mpfr_exp, mpfr_log, "mpfr_exp", 0, -256, 255, 4, 128, 800, 50);
+
+ tests_end_mpfr ();
+ return 0;
+}
projects / msp430-gcc.git / blobdiff