diff -Naurd mpfr-2.2.0-p14/add1.c mpfr-2.2.0-p15/add1.c --- mpfr-2.2.0-p14/add1.c 2005-06-08 10:18:24.000000000 +0000 +++ mpfr-2.2.0-p15/add1.c 2006-08-23 13:17:45.000000000 +0000 @@ -537,5 +537,5 @@ end_of_add: MPFR_TMP_FREE(marker); - return inex; + MPFR_RET (inex); } diff -Naurd mpfr-2.2.0-p14/add1sp.c mpfr-2.2.0-p15/add1sp.c --- mpfr-2.2.0-p14/add1sp.c 2005-09-02 14:22:40.000000000 +0000 +++ mpfr-2.2.0-p15/add1sp.c 2006-08-23 13:17:45.000000000 +0000 @@ -375,5 +375,5 @@ MPFR_SET_SAME_SIGN(a,b); MPFR_TMP_FREE(marker); - return inexact*MPFR_INT_SIGN(a); + MPFR_RET (inexact * MPFR_INT_SIGN (a)); } diff -Naurd mpfr-2.2.0-p14/exceptions.c mpfr-2.2.0-p15/exceptions.c --- mpfr-2.2.0-p14/exceptions.c 2005-08-18 17:03:06.000000000 +0000 +++ mpfr-2.2.0-p15/exceptions.c 2006-08-23 13:17:45.000000000 +0000 @@ -214,7 +214,7 @@ if (MPFR_UNLIKELY( exp > __gmpfr_emax) ) return mpfr_overflow(x, rnd_mode, MPFR_SIGN(x)); } - return t; /* propagate inexact ternary value, unlike most functions */ + MPFR_RET (t); /* propagate inexact ternary value, unlike most functions */ } #undef mpfr_underflow_p diff -Naurd mpfr-2.2.0-p14/exp.c mpfr-2.2.0-p15/exp.c --- mpfr-2.2.0-p14/exp.c 2005-09-13 14:39:20.000000000 +0000 +++ mpfr-2.2.0-p15/exp.c 2006-08-23 13:17:45.000000000 +0000 @@ -128,5 +128,5 @@ inexact = mpfr_check_range (y, inexact, rnd_mode); } - return inexact; + MPFR_RET (inexact); } diff -Naurd mpfr-2.2.0-p14/exp2.c mpfr-2.2.0-p15/exp2.c --- mpfr-2.2.0-p14/exp2.c 2005-08-18 17:03:06.000000000 +0000 +++ mpfr-2.2.0-p15/exp2.c 2006-08-23 13:17:45.000000000 +0000 @@ -31,6 +31,8 @@ mpfr_exp2 (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) { int inexact; + long xint; + mpfr_t xfrac; MPFR_SAVE_EXPO_DECL (expo); if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) @@ -58,9 +60,8 @@ /* since the smallest representable non-zero float is 1/2*2^__gmpfr_emin, if x < __gmpfr_emin - 1, the result is either 1/2*2^__gmpfr_emin or 0 */ - MPFR_ASSERTD (MPFR_EMIN_MIN - 2 >= LONG_MIN); - - if (mpfr_cmp_si_2exp (x, __gmpfr_emin - 1, 0) < 0) + MPFR_ASSERTN (MPFR_EMIN_MIN >= LONG_MIN + 2); + if (MPFR_UNLIKELY (mpfr_cmp_si (x, __gmpfr_emin - 1) < 0)) { mp_rnd_t rnd2 = rnd_mode; /* in round to nearest mode, round to zero when x <= __gmpfr_emin-2 */ @@ -70,64 +71,71 @@ return mpfr_underflow (y, rnd2, 1); } - if (mpfr_integer_p (x)) /* we know that x >= 2^(emin-1) */ - { - long xd; + MPFR_ASSERTN (MPFR_EMAX_MAX <= LONG_MAX); + if (MPFR_UNLIKELY (mpfr_cmp_si (x, __gmpfr_emax) >= 0)) + return mpfr_overflow (y, rnd_mode, 1); - MPFR_ASSERTD (MPFR_EMAX_MAX <= LONG_MAX); - if (mpfr_cmp_si_2exp (x, __gmpfr_emax, 0) > 0) - return mpfr_overflow (y, rnd_mode, 1); + /* We now know that emin - 1 <= x < emax. */ - xd = mpfr_get_si (x, GMP_RNDN); + MPFR_SAVE_EXPO_MARK (expo); - mpfr_set_ui (y, 1, GMP_RNDZ); - return mpfr_mul_2si (y, y, xd, rnd_mode); - } + xint = mpfr_get_si (x, GMP_RNDZ); + mpfr_init2 (xfrac, MPFR_PREC (x)); + mpfr_sub_si (xfrac, x, xint, GMP_RNDN); /* exact */ - MPFR_SAVE_EXPO_MARK (expo); + if (MPFR_IS_ZERO (xfrac)) + { + mpfr_set_ui (y, 1, GMP_RNDN); + inexact = 0; + } + else + { + /* Declaration of the intermediary variable */ + mpfr_t t; - /* General case */ - { - /* Declaration of the intermediary variable */ - mpfr_t t; + /* Declaration of the size variable */ + mp_prec_t Ny = MPFR_PREC(y); /* target precision */ + mp_prec_t Nt; /* working precision */ + mp_exp_t err; /* error */ + MPFR_ZIV_DECL (loop); - /* Declaration of the size variable */ - mp_prec_t Ny = MPFR_PREC(y); /* target precision */ - mp_prec_t Nt; /* working precision */ - mp_exp_t err; /* error */ - MPFR_ZIV_DECL (loop); + /* compute the precision of intermediary variable */ + /* the optimal number of bits : see algorithms.tex */ + Nt = Ny + 5 + MPFR_INT_CEIL_LOG2 (Ny); - /* compute the precision of intermediary variable */ - /* the optimal number of bits : see algorithms.tex */ - Nt = Ny + 5 + MPFR_INT_CEIL_LOG2 (Ny); + /* initialise of intermediary variable */ + mpfr_init2 (t, Nt); - /* initialise of intermediary variable */ - mpfr_init2 (t, Nt); + /* First computation */ + MPFR_ZIV_INIT (loop, Nt); + for (;;) + { + /* compute exp(x*ln(2))*/ + mpfr_const_log2 (t, GMP_RNDU); /* ln(2) */ + mpfr_mul (t, xfrac, t, GMP_RNDU); /* xfrac * ln(2) */ + err = Nt - (MPFR_GET_EXP (t) + 2); /* Estimate of the error */ + mpfr_exp (t, t, GMP_RNDN); /* exp(xfrac * ln(2)) */ - /* First computation */ - MPFR_ZIV_INIT (loop, Nt); - for (;;) - { - /* compute exp(x*ln(2))*/ - mpfr_const_log2 (t, GMP_RNDU); /* ln(2) */ - mpfr_mul (t, x, t, GMP_RNDU); /* x*ln(2) */ - err = Nt - (MPFR_GET_EXP (t) + 2); /* Estimate of the error */ - mpfr_exp (t, t, GMP_RNDN); /* exp(x*ln(2))*/ + if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, Ny, rnd_mode))) + break; - if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, Ny, rnd_mode))) - break; + /* Actualisation of the precision */ + MPFR_ZIV_NEXT (loop, Nt); + mpfr_set_prec (t, Nt); + } + MPFR_ZIV_FREE (loop); - /* Actualisation of the precision */ - MPFR_ZIV_NEXT (loop, Nt); - mpfr_set_prec (t, Nt); - } - MPFR_ZIV_FREE (loop); + inexact = mpfr_set (y, t, rnd_mode); - inexact = mpfr_set (y, t, rnd_mode); + mpfr_clear (t); + } - mpfr_clear (t); - } + mpfr_clear (xfrac); + mpfr_clear_flags (); + mpfr_mul_2si (y, y, xint, GMP_RNDN); /* exact or overflow */ + /* Note: We can have an overflow only when t was rounded up to 2. */ + MPFR_ASSERTD (!mpfr_overflow_p () || inexact > 0); + MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags); MPFR_SAVE_EXPO_FREE (expo); - return mpfr_check_range (y, inexact, rnd_mode); } diff -Naurd mpfr-2.2.0-p14/expm1.c mpfr-2.2.0-p15/expm1.c --- mpfr-2.2.0-p14/expm1.c 2006-03-16 17:47:51.000000000 +0000 +++ mpfr-2.2.0-p15/expm1.c 2006-08-23 13:17:45.000000000 +0000 @@ -19,6 +19,8 @@ the Free Software Foundation, Inc., 51 Franklin Place, Fifth Floor, Boston, MA 02110-1301, USA. */ +#include + #define MPFR_NEED_LONGLONG_H #include "mpfr-impl.h" @@ -30,6 +32,7 @@ mpfr_expm1 (mpfr_ptr y, mpfr_srcptr x , mp_rnd_t rnd_mode) { int inexact; + mp_exp_t ex; MPFR_SAVE_EXPO_DECL (expo); if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) @@ -60,8 +63,39 @@ } } - /* exp(x)-1 = x +x^2/2 + ... so the error is < 2^(2*EXP(x)-1) */ - MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -MPFR_GET_EXP (x)+1,1,rnd_mode,{}); + ex = MPFR_GET_EXP (x); + if (ex < 0) + { + /* For -1 < x < 0, abs(expm1(x)-x) < x^2/2. + For 0 < x < 1, abs(expm1(x)-x) < x^2. */ + if (MPFR_IS_POS (x)) + MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, - ex, 1, rnd_mode, {}); + else + MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, 1 - ex, 0, rnd_mode, {}); + } + + if (MPFR_IS_NEG (x) && ex > 5) /* x <= -32 */ + { + mpfr_t minus_one, t; + mp_exp_t err; + + mpfr_init2 (minus_one, 2); + mpfr_init2 (t, 64); + mpfr_set_si (minus_one, -1, GMP_RNDN); + mpfr_const_log2 (t, GMP_RNDU); /* round upward since x is negative */ + mpfr_div (t, x, t, GMP_RNDU); /* > x / ln(2) */ + err = mpfr_cmp_si (t, MPFR_EMIN_MIN >= -LONG_MAX ? + MPFR_EMIN_MIN : -LONG_MAX) <= 0 ? + - (MPFR_EMIN_MIN >= -LONG_MAX ? MPFR_EMIN_MIN : -LONG_MAX) : + - mpfr_get_si (t, GMP_RNDU); + /* exp(x) = 2^(x/ln(2)) + <= 2^max(MPFR_EMIN_MIN,-LONG_MAX,ceil(x/ln(2)+epsilon)) + with epsilon > 0 */ + mpfr_clear (t); + MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, minus_one, err, 0, rnd_mode, + { mpfr_clear (minus_one); }); + mpfr_clear (minus_one); + } MPFR_SAVE_EXPO_MARK (expo); /* General case */ @@ -80,8 +114,8 @@ /* if |x| is smaller than 2^(-e), we will loose about e bits in the subtraction exp(x) - 1 */ - if (MPFR_EXP(x) < 0) - Nt += -MPFR_EXP(x); + if (ex < 0) + Nt += - ex; /* initialize auxiliary variable */ mpfr_init2 (t, Nt); diff -Naurd mpfr-2.2.0-p14/lngamma.c mpfr-2.2.0-p15/lngamma.c --- mpfr-2.2.0-p14/lngamma.c 2005-09-29 11:27:04.000000000 +0000 +++ mpfr-2.2.0-p15/lngamma.c 2006-08-23 13:17:45.000000000 +0000 @@ -154,11 +154,16 @@ } /* if x < 0 and -2k-1 <= x <= -2k, then lngamma(x) = NaN */ - if (MPFR_IS_NEG(z0) && ((mpfr_get_si (z0, GMP_RNDZ) % 2) == 0 - || mpfr_integer_p (z0))) + if (MPFR_IS_NEG (z0)) { - MPFR_SET_NAN (y); - MPFR_RET_NAN; + MPFR_SAVE_EXPO_MARK (expo); + if (mpfr_get_si (z0, GMP_RNDZ) % 2 == 0 || mpfr_integer_p (z0)) + { + MPFR_SAVE_EXPO_FREE (expo); + MPFR_SET_NAN (y); + MPFR_RET_NAN; + } + MPFR_SAVE_EXPO_FREE (expo); } #endif diff -Naurd mpfr-2.2.0-p14/log1p.c mpfr-2.2.0-p15/log1p.c --- mpfr-2.2.0-p14/log1p.c 2006-03-16 17:47:51.000000000 +0000 +++ mpfr-2.2.0-p15/log1p.c 2006-08-23 13:17:45.000000000 +0000 @@ -29,6 +29,7 @@ mpfr_log1p (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) { int comp, inexact; + mp_exp_t ex; MPFR_SAVE_EXPO_DECL (expo); if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) @@ -62,8 +63,16 @@ } } - /* log(1+x) = x-x^2/2 + ... so the error is < 2^(2*EXP(x)-1) */ - MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -MPFR_GET_EXP (x)+1,0,rnd_mode,{}); + ex = MPFR_GET_EXP (x); + if (ex < 0) /* -0.5 < x < 0.5 */ + { + /* For x > 0, abs(log(1+x)-x) < x^2/2. + For x > -0.5, abs(log(1+x)-x) < x^2. */ + if (MPFR_IS_POS (x)) + MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, - ex - 1, 0, rnd_mode, {}); + else + MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, - ex, 1, rnd_mode, {}); + } comp = mpfr_cmp_si (x, -1); /* log1p(x) is undefined for x < -1 */ diff -Naurd mpfr-2.2.0-p14/mpfr-impl.h mpfr-2.2.0-p15/mpfr-impl.h --- mpfr-2.2.0-p14/mpfr-impl.h 2006-05-26 21:00:45.000000000 +0000 +++ mpfr-2.2.0-p15/mpfr-impl.h 2006-08-23 13:17:45.000000000 +0000 @@ -915,7 +915,8 @@ /* Speed up final checking */ #define mpfr_check_range(x,t,r) \ (MPFR_LIKELY (MPFR_EXP (x) >= __gmpfr_emin && MPFR_EXP (x) <= __gmpfr_emax) \ - ? (t) : mpfr_check_range(x,t,r)) + ? ((t) ? (__gmpfr_flags |= MPFR_FLAGS_INEXACT, (t)) : 0) \ + : mpfr_check_range(x,t,r)) /****************************************************** diff -Naurd mpfr-2.2.0-p14/mul.c mpfr-2.2.0-p15/mul.c --- mpfr-2.2.0-p14/mul.c 2005-09-11 22:57:03.000000000 +0000 +++ mpfr-2.2.0-p15/mul.c 2006-08-23 13:17:45.000000000 +0000 @@ -155,7 +155,7 @@ MPFR_SET_EXP (a, ax2); MPFR_SET_SIGN(a, sign_product); } - return inexact; + MPFR_RET (inexact); } int @@ -508,5 +508,5 @@ rnd_mode = GMP_RNDZ; return mpfr_underflow (a, rnd_mode, sign); } - return inexact; + MPFR_RET (inexact); } diff -Naurd mpfr-2.2.0-p14/sub1.c mpfr-2.2.0-p15/sub1.c --- mpfr-2.2.0-p14/sub1.c 2005-08-18 16:35:14.000000000 +0000 +++ mpfr-2.2.0-p15/sub1.c 2006-08-23 13:17:45.000000000 +0000 @@ -107,9 +107,9 @@ if (MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (a))) { mpfr_nexttozero (a); - return -MPFR_INT_SIGN (a); + MPFR_RET (- MPFR_INT_SIGN (a)); } - return MPFR_INT_SIGN (a); + MPFR_RET (MPFR_INT_SIGN (a)); } else { @@ -137,7 +137,7 @@ mpfr_nexttozero (a); inexact = -MPFR_INT_SIGN (a); } - return inexact; + MPFR_RET (inexact); } } @@ -531,5 +531,5 @@ #endif /* check that result is msb-normalized */ MPFR_ASSERTD(ap[an-1] > ~ap[an-1]); - return inexact * MPFR_INT_SIGN(a); + MPFR_RET (inexact * MPFR_INT_SIGN (a)); } diff -Naurd mpfr-2.2.0-p14/sub1sp.c mpfr-2.2.0-p15/sub1sp.c --- mpfr-2.2.0-p14/sub1sp.c 2005-09-02 14:22:41.000000000 +0000 +++ mpfr-2.2.0-p15/sub1sp.c 2006-08-23 13:17:45.000000000 +0000 @@ -791,7 +791,5 @@ MPFR_SET_EXP (a, bx); MPFR_TMP_FREE(marker); - - return inexact*MPFR_INT_SIGN(a); + MPFR_RET (inexact * MPFR_INT_SIGN (a)); } - diff -Naurd mpfr-2.2.0-p14/tests/tests.c mpfr-2.2.0-p15/tests/tests.c --- mpfr-2.2.0-p14/tests/tests.c 2005-08-18 17:03:14.000000000 +0000 +++ mpfr-2.2.0-p15/tests/tests.c 2006-08-23 13:17:45.000000000 +0000 @@ -289,7 +289,8 @@ } /* Open a file in the src directory - can't use fopen directly */ -FILE *src_fopen (const char *filename, const char *mode) +FILE * +src_fopen (const char *filename, const char *mode) { const char *srcdir = getenv ("srcdir"); char *buffer; @@ -309,7 +310,8 @@ return f; } -void set_emin (mp_exp_t exponent) +void +set_emin (mp_exp_t exponent) { if (mpfr_set_emin (exponent)) { @@ -318,7 +320,8 @@ } } -void set_emax (mp_exp_t exponent) +void +set_emax (mp_exp_t exponent) { if (mpfr_set_emax (exponent)) { diff -Naurd mpfr-2.2.0-p14/tests/texp2.c mpfr-2.2.0-p15/tests/texp2.c --- mpfr-2.2.0-p14/tests/texp2.c 2005-08-18 16:35:17.000000000 +0000 +++ mpfr-2.2.0-p15/tests/texp2.c 2006-08-23 13:17:45.000000000 +0000 @@ -22,6 +22,7 @@ #include #include +#include #include "mpfr-test.h" @@ -32,6 +33,7 @@ special_overflow (void) { mpfr_t x, y; + int inex; set_emin (-125); set_emax (128); @@ -40,11 +42,12 @@ mpfr_init2 (y, 24); mpfr_set_str_binary (x, "0.101100100000000000110100E15"); - mpfr_exp2 (y, x, GMP_RNDN); - if (!mpfr_inf_p(y)) + inex = mpfr_exp2 (y, x, GMP_RNDN); + if (!mpfr_inf_p (y) || inex <= 0) { - printf("Overflow error.\n"); + printf ("Overflow error.\n"); mpfr_dump (y); + printf ("inex = %d\n", inex); exit (1); } @@ -54,6 +57,80 @@ set_emax (MPFR_EMAX_MAX); } +static void +emax_m_eps (void) +{ + if (mpfr_get_emax () <= LONG_MAX) + { + mpfr_t x, y; + int inex, ov; + + mpfr_init2 (x, 64); + mpfr_init2 (y, 8); + mpfr_set_si (x, mpfr_get_emax (), GMP_RNDN); + + mpfr_clear_flags (); + inex = mpfr_exp2 (y, x, GMP_RNDN); + ov = mpfr_overflow_p (); + if (!ov || !mpfr_inf_p (y) || inex <= 0) + { + printf ("Overflow error for x = emax, GMP_RNDN.\n"); + mpfr_dump (y); + printf ("inex = %d, %soverflow\n", inex, ov ? "" : "no "); + exit (1); + } + + mpfr_nextbelow (x); + + mpfr_clear_flags (); + inex = mpfr_exp2 (y, x, GMP_RNDN); + ov = mpfr_overflow_p (); + if (!ov || !mpfr_inf_p (y) || inex <= 0) + { + printf ("Overflow error for x = emax - eps, GMP_RNDN.\n"); + mpfr_dump (y); + printf ("inex = %d, %soverflow\n", inex, ov ? "" : "no "); + exit (1); + } + + mpfr_clear_flags (); + inex = mpfr_exp2 (y, x, GMP_RNDD); + ov = mpfr_overflow_p (); + if (ov || mpfr_inf_p (y) || inex >= 0 || + (mpfr_nextabove (y), !mpfr_inf_p (y))) + { + printf ("Overflow error for x = emax - eps, GMP_RNDD.\n"); + mpfr_dump (y); + printf ("inex = %d, %soverflow\n", inex, ov ? "" : "no "); + exit (1); + } + + mpfr_clear (x); + mpfr_clear (y); + } +} + +static void +exp_range (void) +{ + mpfr_t x; + + set_emin (3); + mpfr_init2 (x, 8); + mpfr_set_ui (x, 5, GMP_RNDN); + mpfr_exp2 (x, x, GMP_RNDN); + set_emin (MPFR_EMIN_MIN); + if (mpfr_nan_p (x) || mpfr_cmp_ui (x, 32) != 0) + { + printf ("Error in mpfr_exp2 for x = 5, with emin = 3\n"); + printf ("Expected 32, got "); + mpfr_out_str (stdout, 2, 0, x, GMP_RNDN); + printf ("\n"); + exit (1); + } + mpfr_clear (x); +} + int main (int argc, char *argv[]) { @@ -63,6 +140,8 @@ tests_start_mpfr (); special_overflow (); + emax_m_eps (); + exp_range (); mpfr_init (x); mpfr_init (y); @@ -154,6 +233,17 @@ exit (1); } + mpfr_set_prec (x, 40); + mpfr_set_prec (y, 40); + mpfr_set_str (x, "10000000000.5", 10, GMP_RNDN); + mpfr_clear_flags (); + mpfr_exp2 (y, x, GMP_RNDN); + if (!(MPFR_IS_INF (y) && MPFR_IS_POS (y) && mpfr_overflow_p ())) + { + printf ("exp2(10000000000.5) should overflow.\n"); + exit (1); + } + test_generic (2, 100, 100); mpfr_clear (x); diff -Naurd mpfr-2.2.0-p14/tests/tlog1p.c mpfr-2.2.0-p15/tests/tlog1p.c --- mpfr-2.2.0-p14/tests/tlog1p.c 2005-08-18 17:03:16.000000000 +0000 +++ mpfr-2.2.0-p15/tests/tlog1p.c 2006-08-23 13:17:45.000000000 +0000 @@ -88,12 +88,40 @@ mpfr_clear (x); } +static void +other (void) +{ + mpfr_t x, y; + + /* Bug reported by Guillaume Melquiond on 2006-08-14. */ + mpfr_init2 (x, 53); + mpfr_set_str (x, "-1.5e4f72873ed9a@-100", 16, GMP_RNDN); + mpfr_init2 (y, 57); + mpfr_log1p (y, x, GMP_RNDU); + if (mpfr_cmp (x, y) != 0) + { + printf ("Error in tlog1p for x = "); + mpfr_out_str (stdout, 16, 0, x, GMP_RNDN); + printf (", rnd = GMP_RNDU\nExpected "); + mpfr_out_str (stdout, 16, 15, x, GMP_RNDN); + printf ("\nGot "); + mpfr_out_str (stdout, 16, 15, y, GMP_RNDN); + printf ("\n"); + exit (1); + } + + mpfr_clear (y); + mpfr_clear (x); + return; +} + int main (int argc, char *argv[]) { tests_start_mpfr (); special (); + other (); test_generic (2, 100, 50);