diff -Naur mpfr-2.1.1/exp_2.c mpfr-2.1.1-p1/exp_2.c --- mpfr-2.1.1/exp_2.c 2004-08-24 14:00:42.000000000 +0000 +++ mpfr-2.1.1-p1/exp_2.c 2005-03-09 07:52:34.226856000 +0000 @@ -1,7 +1,7 @@ /* mpfr_exp_2 -- exponential of a floating-point number using Brent's algorithms in O(n^(1/2)*M(n)) and O(n^(1/3)*M(n)) -Copyright 1999, 2000, 2001, 2002, 2003, 2004 Free Software Foundation, Inc. +Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -147,7 +147,7 @@ mpfr_sub (r, x, r, GMP_RNDU); /* possible cancellation here: the error on r is at most 3*2^(EXP(old_r)-EXP(new_r)) */ - if (MPFR_IS_NEG (r)) + while (MPFR_IS_NEG (r)) { /* initial approximation n was too large */ n--; mpfr_add (r, r, s, GMP_RNDU); diff -Naur mpfr-2.1.1/tests/texp.c mpfr-2.1.1-p1/tests/texp.c --- mpfr-2.1.1/tests/texp.c 2005-01-29 11:36:41.000000000 +0000 +++ mpfr-2.1.1-p1/tests/texp.c 2005-03-05 01:25:57.000000000 +0000 @@ -346,6 +346,21 @@ exit (1); } + /* 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'); + } + mpfr_clear (x); mpfr_clear (y); } diff -Naur mpfr-2.1.1-p1/add1sp.c mpfr-2.1.1-p2/add1sp.c --- mpfr-2.1.1-p1/add1sp.c 2004-05-04 09:05:00.000000000 +0000 +++ mpfr-2.1.1-p2/add1sp.c 2005-03-09 15:11:14.181932000 +0000 @@ -1,7 +1,7 @@ /* mpfr_add1sp -- internal function to perform a "real" addition All the op must have the same precision -Copyright 2004 Free Software Foundation. +Copyright 2004, 2005 Free Software Foundation. Contributed by the Spaces project, INRIA Lorraine. This file is part of the MPFR Library. @@ -101,29 +101,34 @@ else goto add_one_ulp; } - else if (MPFR_UNLIKELY(d >= p)) + else if (MPFR_UNLIKELY (d >= p)) { - if (MPFR_LIKELY(d > p)) + if (MPFR_LIKELY (d > p)) { /* d > p : Copy B in A */ - ap = MPFR_MANT(a); - MPN_COPY (ap, MPFR_MANT(b), n); /* Away: Add 1 Nearest: Trunc Zero: Trunc */ - if (MPFR_LIKELY(rnd_mode==GMP_RNDN)) - { inexact = -1; goto set_exponent; } - MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(b)); - if (rnd_mode==GMP_RNDZ) - { inexact = -1; goto set_exponent; } + if (MPFR_LIKELY (rnd_mode==GMP_RNDN + || MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (b)))) + { + copy_set_exponent: + ap = MPFR_MANT (a); + MPN_COPY (ap, MPFR_MANT(b), n); + inexact = -1; + goto set_exponent; + } else - goto add_one_ulp; + { + copy_add_one_ulp: + ap = MPFR_MANT(a); + MPN_COPY (ap, MPFR_MANT(b), n); + goto add_one_ulp; + } } else { /* d==p : Copy B in A */ - ap = MPFR_MANT(a); - MPN_COPY (ap, MPFR_MANT(b), n); /* Away: Add 1 Nearest: Even Rule if C is a power of 2, else Add 1 Zero: Trunc */ @@ -139,17 +144,16 @@ } while (k>=0 && cp[k]==0); if (MPFR_UNLIKELY(k<0)) /* Power of 2: Even rule */ - if ((ap[0]&(MPFR_LIMB_ONE< 1, and arctanh(+/-1) = +/-Inf */ if (MPFR_EXP(xt) > 0) { - if (MPFR_EXP(xt) == 1) + if (MPFR_EXP(xt) == 1 && mpfr_powerof2_raw (xt)) { - if (mpfr_cmp_ui (xt, 1) || mpfr_cmp_si (xt, -1)) - { - MPFR_SET_INF(y); - MPFR_SET_SAME_SIGN(y, xt); - MPFR_RET(0); - } + MPFR_SET_INF(y); + MPFR_SET_SAME_SIGN(y, xt); + MPFR_RET(0); } MPFR_SET_NAN(y); MPFR_RET_NAN; diff -Naur mpfr-2.1.1-p2/tests/tatanh.c mpfr-2.1.1-p3/tests/tatanh.c --- mpfr-2.1.1-p2/tests/tatanh.c 2004-02-13 09:33:55.000000000 +0000 +++ mpfr-2.1.1-p3/tests/tatanh.c 2005-03-17 00:21:40.347287000 +0000 @@ -1,6 +1,6 @@ /* Test file for mpfr_atanh. -Copyright 2001, 2002, 2003, 2004 Free Software Foundation. +Copyright 2001, 2002, 2003, 2004, 2005 Free Software Foundation. Adapted from tatan.c. This file is part of the MPFR Library. @@ -32,6 +32,7 @@ special (void) { mpfr_t x, y, z; + int i; mpfr_init (x); mpfr_init (y); @@ -53,20 +54,24 @@ exit (1); } - /* atanh(+/-2) = NaN */ - mpfr_set_ui (x, 2, GMP_RNDN); - mpfr_atanh (y, x, GMP_RNDN); - if (!mpfr_nan_p (y)) + /* atanh(+/-x) = NaN if x > 1 */ + for (i = 3; i <= 6; i++) { - printf ("Error: mpfr_atanh(2) <> NaN\n"); - exit (1); - } - mpfr_neg (x, x, GMP_RNDN); - mpfr_atanh (y, x, GMP_RNDN); - if (!mpfr_nan_p (y)) - { - printf ("Error: mpfr_atanh(-2) <> NaN\n"); - exit (1); + mpfr_set_si (x, i, GMP_RNDN); + mpfr_div_2ui (x, x, 1, GMP_RNDN); + mpfr_atanh (y, x, GMP_RNDN); + if (!mpfr_nan_p (y)) + { + printf ("Error: mpfr_atanh(%d/2) <> NaN\n", i); + exit (1); + } + mpfr_neg (x, x, GMP_RNDN); + mpfr_atanh (y, x, GMP_RNDN); + if (!mpfr_nan_p (y)) + { + printf ("Error: mpfr_atanh(-%d/2) <> NaN\n", i); + exit (1); + } } /* atanh(+0) = +0, atanh(-0) = -0 */ diff -Naur mpfr-2.1.1-p3/pow_ui.c mpfr-2.1.1-p4/pow_ui.c --- mpfr-2.1.1-p3/pow_ui.c 2004-02-23 09:43:29.000000000 +0000 +++ mpfr-2.1.1-p4/pow_ui.c 2005-03-30 14:09:16.299081000 +0000 @@ -1,7 +1,8 @@ /* mpfr_pow_ui-- compute the power of a floating-point by a machine integer -Copyright 1999, 2000, 2001, 2002, 2003, 2004 Free Software Foundation, Inc. +Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005 + Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -61,6 +62,10 @@ MPFR_ASSERTD(MPFR_IS_ZERO(y)); /* 0^n = 0 for any n */ MPFR_SET_ZERO(x); + if (MPFR_IS_POS (y) || ((n & 1) == 0)) + MPFR_SET_POS (x); + else + MPFR_SET_NEG (x); MPFR_RET(0); } } @@ -91,7 +96,7 @@ ; mpfr_set_prec (res, prec); inexact = mpfr_set (res, y, rnd1); - err = prec <= (mpfr_prec_t) i ? 0 : prec - (mpfr_prec_t) i; + err = prec <= (mpfr_prec_t) i ? 0 : prec - 1 - (mpfr_prec_t) i; MPFR_ASSERTD (i >= 1); /* now 2^(i-1) <= n < 2^i */ for (i -= 2; i >= 0; i--) diff -Naur mpfr-2.1.1-p3/tests/tpow.c mpfr-2.1.1-p4/tests/tpow.c --- mpfr-2.1.1-p3/tests/tpow.c 2005-01-27 17:34:39.000000000 +0000 +++ mpfr-2.1.1-p4/tests/tpow.c 2005-03-30 14:08:33.308980000 +0000 @@ -84,6 +84,22 @@ exit (1); } + mpfr_set_prec (a, 29); + mpfr_set_prec (b, 29); + mpfr_set_str_binary (a, "1.0000000000000000000000001111"); + mpfr_set_str_binary (b, "1.1001101111001100111001010111e165"); + mpfr_pow_ui (a, a, 2055225053, GMP_RNDZ); + if (mpfr_cmp (a, b) != 0) + { + printf ("Error for x^2055225053\n"); + printf ("Expected "); + mpfr_out_str (stdout, 2, 0, b, GMP_RNDN); + printf ("\nGot "); + mpfr_out_str (stdout, 2, 0, a, GMP_RNDN); + printf ("\n"); + exit (1); + } + mpfr_clear (a); mpfr_clear (b); } diff -Naur mpfr-2.1.1-p4/get_ld.c mpfr-2.1.1-p5/get_ld.c --- mpfr-2.1.1-p4/get_ld.c 2004-02-06 13:27:01.000000000 +0000 +++ mpfr-2.1.1-p5/get_ld.c 2005-04-04 14:47:03.000000000 +0000 @@ -1,7 +1,7 @@ /* mpfr_get_ld -- convert a multiple precision floating-point number to a machine long double -Copyright 2002, 2003, 2004 Free Software Foundation, Inc. +Copyright 2002, 2003, 2004, 2005 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -22,9 +22,14 @@ #include - #include "mpfr-impl.h" +/* The MPFR number rounded to the size MPFR_LDBL_MANT_DIG of a long double + will be scaled so that all its bits are representable in a double. + Therefore the exponent of its LSB must be >= -1074, and the exponent + of the MPFR number must be >= -1074 + MPFR_LDBL_MANT_DIG. */ +#define EMIN (-1074 + MPFR_LDBL_MANT_DIG) + long double mpfr_get_ld (mpfr_srcptr x, mp_rnd_t rnd_mode) { @@ -57,10 +62,10 @@ sh = e - 1023; MPFR_SET_EXP (y, 1023); } - else if (e < -1021) + else if (e < EMIN) { - sh = e + 1021; - MPFR_SET_EXP (y, -1021); + sh = e - EMIN; + MPFR_SET_EXP (y, EMIN); } else { diff -Naur mpfr-2.1.1-p4/tests/tset_ld.c mpfr-2.1.1-p5/tests/tset_ld.c --- mpfr-2.1.1-p4/tests/tset_ld.c 2005-01-27 17:40:41.000000000 +0000 +++ mpfr-2.1.1-p5/tests/tset_ld.c 2005-04-04 13:58:41.000000000 +0000 @@ -87,6 +87,32 @@ } } +static void +test_small (void) +{ + mpfr_t x, y; + long double d; + + mpfr_init2 (x, 64); + mpfr_init2 (y, 64); + + mpfr_set_str (x, "0.1010010100111100110000001110101101000111010110000001111101110011E-16381", 2, GMP_RNDN); + mpfr_get_ld (x, GMP_RNDN); /* infinite loop? */ + mpfr_set_ld (y, d, GMP_RNDN); + mpfr_sub (y, x, y, GMP_RNDN); + mpfr_abs (y, y, GMP_RNDN); + if (mpfr_cmp_str (y, "1E-16434", 2, GMP_RNDN) > 0) + { + printf ("Error with "); + mpfr_out_str (NULL, 16, 0, x, GMP_RNDN); + printf ("\n"); + exit (1); + } + + mpfr_clear (x); + mpfr_clear (y); +} + int main (int argc, char *argv[]) { @@ -183,6 +209,8 @@ mpfr_clear (x); + test_small (); + tests_end_mpfr (); return 0; diff -Naur mpfr-2.1.1-p5/set_ld.c mpfr-2.1.1-p6/set_ld.c --- mpfr-2.1.1-p5/set_ld.c 2004-02-16 14:30:43.000000000 +0000 +++ mpfr-2.1.1-p6/set_ld.c 2005-04-21 14:12:23.000000000 +0000 @@ -1,7 +1,7 @@ /* mpfr_set_ld -- convert a machine long double to a multiple precision floating-point number -Copyright 2002, 2003, 2004 Free Software Foundation, Inc. +Copyright 2002, 2003, 2004, 2005 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -110,21 +110,28 @@ d = d / div9; /* exact */ shift_exp += 512; } - mpfr_set_ui (u, 0, GMP_RNDZ); } else { + long double div9, div10, div11; + div9 = (long double) (double) 7.4583407312002067432909653e-155; + /* div9 = 2^(-2^9) */ + div10 = div9 * div9; /* 2^(-2^10) */ + div11 = div10 * div10; /* 2^(-2^11) if extended precision */ /* since -DBL_MAX <= d <= DBL_MAX, the cast to double should not overflow here */ inexact = mpfr_set_d (u, (double) d, GMP_RNDN); MPFR_ASSERTD(inexact == 0); - if (MPFR_IS_ZERO (u) && (d != (long double) 0.0)) /* underflow */ - { - long double div10, div11, div12, div13; - div10 = (long double) (double) 5.5626846462680034577255e-309; /* 2^(-2^10) */ - div11 = div10 * div10; /* 2^(-2^11) */ - div12 = div11 * div11; /* 2^(-2^12) */ - div13 = div12 * div12; /* 2^(-2^13) */ + if (d != (long double) 0.0 && + ABS(d) < div10 && + div11 != (long double) 0.0 && + div11 / div10 == div10) /* possible underflow */ + { + long double div12, div13; + /* After the divisions, any bit of d must be >= div10, + hence the possible division by div9. */ + div12 = div11 * div11; /* 2^(-2^12) */ + div13 = div12 * div12; /* 2^(-2^13) */ if (ABS(d) <= div13) { d = d / div13; /* exact */ @@ -145,12 +152,20 @@ d = d / div10; /* exact */ shift_exp -= 1024; } + if (ABS(d) <= div9) + { + d = d / div9; /* exact */ + shift_exp -= 512; + } } + else + { + mpfr_add (t, t, u, GMP_RNDN); /* exact */ + if (!mpfr_number_p (t)) + break; + d = d - (long double) mpfr_get_d1 (u); /* exact */ + } } - mpfr_add (t, t, u, GMP_RNDN); /* exact */ - if (!mpfr_number_p (t)) - break; - d = d - (long double) mpfr_get_d1 (u); /* exact */ } /* now t is exactly the input value d */ inexact = mpfr_set (r, t, rnd_mode); diff -Naur mpfr-2.1.1-p5/tests/tset_ld.c mpfr-2.1.1-p6/tests/tset_ld.c --- mpfr-2.1.1-p5/tests/tset_ld.c 2005-04-04 13:58:41.000000000 +0000 +++ mpfr-2.1.1-p6/tests/tset_ld.c 2005-04-21 13:00:38.000000000 +0000 @@ -23,6 +23,9 @@ #include #include #include +#if WITH_FPU_CONTROL +#include +#endif #include "mpfr-test.h" @@ -90,27 +93,42 @@ static void test_small (void) { - mpfr_t x, y; + mpfr_t x, y, z; long double d; mpfr_init2 (x, 64); mpfr_init2 (y, 64); + mpfr_init2 (z, 64); mpfr_set_str (x, "0.1010010100111100110000001110101101000111010110000001111101110011E-16381", 2, GMP_RNDN); - mpfr_get_ld (x, GMP_RNDN); /* infinite loop? */ + d = mpfr_get_ld (x, GMP_RNDN); /* infinite loop? */ mpfr_set_ld (y, d, GMP_RNDN); - mpfr_sub (y, x, y, GMP_RNDN); - mpfr_abs (y, y, GMP_RNDN); - if (mpfr_cmp_str (y, "1E-16434", 2, GMP_RNDN) > 0) + mpfr_sub (z, x, y, GMP_RNDN); + mpfr_abs (z, z, GMP_RNDN); + mpfr_clear_erangeflag (); + /* If long double = double, d should be equal to 0; + in this case, everything is OK. */ + if (d != 0 && (mpfr_cmp_str (z, "1E-16434", 2, GMP_RNDN) > 0 || + mpfr_erangeflag_p ())) { - printf ("Error with "); + printf ("Error with x = "); + mpfr_out_str (NULL, 10, 20, x, GMP_RNDN); + printf (" = "); mpfr_out_str (NULL, 16, 0, x, GMP_RNDN); + printf ("\n -> d = %.20Lg", d); + printf ("\n -> y = "); + mpfr_out_str (NULL, 10, 20, y, GMP_RNDN); + printf (" = "); + mpfr_out_str (NULL, 16, 0, y, GMP_RNDN); + printf ("\n -> |x-y| = "); + mpfr_out_str (NULL, 16, 0, z, GMP_RNDN); printf ("\n"); exit (1); } mpfr_clear (x); mpfr_clear (y); + mpfr_clear (z); } int @@ -120,6 +138,16 @@ mpfr_t x; int i; mp_exp_t emax; +#if WITH_FPU_CONTROL + fpu_control_t cw; + + if (argc > 1) + { + cw = strtol(argv[1], NULL, 0); + printf ("FPU control word: 0x%x\n", (unsigned int) cw); + _FPU_SETCW (cw); + } +#endif check_gcc33_bug (); @@ -170,7 +198,9 @@ check_set_get (-d, x); /* checks the smallest power of two */ - d = 1.0; while ((e = d / 2.0) != (long double) 0.0) d = e; + d = 1.0; + while ((e = d / 2.0) != (long double) 0.0 && e != d) + d = e; check_set_get (d, x); check_set_get (-d, x); diff -Naur mpfr-2.1.1-p6/div_ui.c mpfr-2.1.1-p7/div_ui.c --- mpfr-2.1.1-p6/div_ui.c 2004-03-26 13:27:15.000000000 +0000 +++ mpfr-2.1.1-p7/div_ui.c 2005-05-04 15:55:58.183515000 +0000 @@ -1,6 +1,7 @@ /* mpfr_div_ui -- divide a floating-point number by a machine integer -Copyright 1999, 2000, 2001, 2002, 2003, 2004 Free Software Foundation, Inc. +Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005 + Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -96,40 +97,75 @@ c = mpn_divrem_1 (tmp, 0, xp - dif, yn + 1, c); inexact = (c != 0); + + /* First pass in estimating next bit of the quotient, in case of RNDN * + * In case we just have the right number of bits (postpone this ?), * + * we need to check whether the remainder is more or less than half * + * the divisor. The test must be performed with a subtraction, so as * + * to prevent carries. */ + if (rnd_mode == GMP_RNDN) { - if (2 * c < (mp_limb_t) u) + if (c < (mp_limb_t) u - c) /* We have u > c */ middle = -1; - else if (2 * c > (mp_limb_t) u) + else if (c > (mp_limb_t) u - c) middle = 1; else middle = 0; /* exactly in the middle */ } + + /* If we believe that we are right in the middle or exact, we should check + that we did not neglect any word of x (division large / 1 -> small). */ + for (i=0; ((inexact == 0) || (middle == 0)) && (i < -dif); i++) if (xp[i]) inexact = middle = 1; /* larger than middle */ - if (tmp[yn] == 0) /* high limb is zero */ + /* + If the high limb of the result is 0 (xp[xn-1] < u), remove it. + Otherwise, compute the left shift to be performed to normalize. + In the latter case, we discard some low bits computed. They + contain information useful for the rounding, hence the updating + of middle and inexact. + */ + + if (tmp[yn] == 0) { - tmp--; - sh = 0; + MPN_COPY(yp, tmp, yn); exp -= BITS_PER_MP_LIMB; + sh = 0; } + else + { + count_leading_zeros (sh, tmp[yn]); - /* now we have yn limbs starting from tmp[1], with tmp[yn]<>0 */ + /* shift left to normalize */ + if (sh) + { + mp_limb_t w = tmp[0] << sh; + + mpn_lshift (yp, tmp + 1, yn, sh); + yp[0] += tmp[0] >> (BITS_PER_MP_LIMB - sh); + + if (w > (MPFR_LIMB_ONE << (BITS_PER_MP_LIMB - 1))) + { middle = 1; } + else if (w < (MPFR_LIMB_ONE << (BITS_PER_MP_LIMB - 1))) + { middle = -1; } + else + { middle = (c != 0); } + + inexact = inexact || (w != 0); + exp -= sh; + } + else + { /* this happens only if u == 1 and xp[xn-1] >= + 1<<(BITS_PER_MP_LIMB-1). It might be better to handle the + u == 1 case seperately ? + */ - /* shift left to normalize */ - count_leading_zeros(sh, tmp[yn]); - if (sh) - { - mpn_lshift (yp, tmp + 1, yn, sh); - yp[0] += tmp[0] >> (BITS_PER_MP_LIMB - sh); - middle = middle || ((tmp[0] << sh) != 0); - inexact = inexact || ((tmp[0] << sh) != 0); - exp -= sh; + MPN_COPY (yp, tmp + 1, yn); + } } - else - MPN_COPY(yp, tmp + 1, yn); sh = yn * BITS_PER_MP_LIMB - MPFR_PREC(y); /* it remains sh bits in less significant limb of y */ @@ -160,6 +196,7 @@ default: MPFR_ASSERTD(rnd_mode == GMP_RNDN); + /* we have one more significant bit in yn */ if (sh && d < (MPFR_LIMB_ONE << (sh - 1))) MPFR_RET(-MPFR_INT_SIGN(x)); else if (sh && d > (MPFR_LIMB_ONE << (sh - 1))) @@ -169,12 +206,14 @@ } else /* sh = 0 or d = 1 << (sh-1) */ { - /* we are in the middle if: - (a) sh > 0 and inexact == 0 - (b) sh=0 and middle=1 - */ + /* The first case is "false" even rounding (significant bits + indicate even rounding, but the result is inexact, so up) ; + The second case is the case where middle should be used to + decide the direction of rounding (no further bit computed) ; + The third is the true even rounding. + */ if ((sh && inexact) || (!sh && (middle > 0)) || - (*yp & (MPFR_LIMB_ONE << sh))) + (!inexact && *yp & (MPFR_LIMB_ONE << sh))) { mpfr_add_one_ulp (y, rnd_mode); MPFR_RET(MPFR_INT_SIGN(x)); diff -Naur mpfr-2.1.1-p6/tests/tdiv_ui.c mpfr-2.1.1-p7/tests/tdiv_ui.c --- mpfr-2.1.1-p6/tests/tdiv_ui.c 2005-01-27 17:22:35.000000000 +0000 +++ mpfr-2.1.1-p7/tests/tdiv_ui.c 2005-05-04 16:02:25.803493000 +0000 @@ -219,6 +219,20 @@ } mpfr_clear(x); + mpfr_init2 (x, 32); + mpfr_set_str (x, "1.0000000000000000000000000000111e-2", 2, GMP_RNDD); + mpfr_div_ui (x, x, 275604255, GMP_RNDN); + if (mpfr_cmp_str (x, "1.1111001010101110101010001101111e-31", 2, GMP_RNDD)) + { + printf ("Error for x=1.0000000000000000000000000000111e-2,\n" + "u=275604255, prec=32, rnd_mode=GMP_RNDN\n" + "got "); + mpfr_out_str (stdout, 2, 0, x, GMP_RNDD); + printf ("\nexpected 1.1111001010101110101010001101111e-31\n"); + exit (1); + } + mpfr_clear(x); + tests_end_mpfr (); return 0; } diff -Naur mpfr-2.1.1-p7/hypot.c mpfr-2.1.1-p8/hypot.c --- mpfr-2.1.1-p7/hypot.c 2004-02-14 11:11:31.000000000 +0000 +++ mpfr-2.1.1-p8/hypot.c 2005-05-12 16:06:13.525321000 +0000 @@ -1,6 +1,6 @@ /* mpfr_hypot -- Euclidean distance -Copyright 2001, 2002, 2003, 2004 Free Software Foundation, Inc. +Copyright 2001, 2002, 2003, 2004, 2005 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -84,19 +84,33 @@ or 2^(2*Ey) <= 2^(2*Ex-1-Nz), i.e. 2*diff_exp > Nz. Warning: this is true only for Nx <= Nz, otherwise the trailing bits of x may be already very close to 1/2*ulp(x,Nz)! + If Nx > Nz, then we can notice that it is possible to round on Nx bits + if 2*diff_exp > Nx (see above as if Nz = Nx), therefore on Nz bits. + Hence the condition: 2*diff_exp > MAX(Nz,Nx). */ - if (MPFR_PREC(x) <= Nz && diff_exp > Nz / 2) /* result is |x| or |x|+ulp(|x|,Nz) */ + if (diff_exp > MAX (Nz, MPFR_PREC (x)) / 2) + /* result is |x| or |x|+ulp(|x|,Nz) */ { if (rnd_mode == GMP_RNDU) { - /* if z > abs(x), then it was already rounded up */ - if (mpfr_abs (z, x, rnd_mode) <= 0) + /* If z > abs(x), then it was already rounded up; otherwise + z = abs(x), and we need to add one ulp due to y. */ + if (mpfr_abs (z, x, rnd_mode) == 0) mpfr_add_one_ulp (z, rnd_mode); return 1; } else /* GMP_RNDZ, GMP_RNDD, GMP_RNDN */ { inexact = mpfr_abs (z, x, rnd_mode); + if (inexact == -MPFR_EVEN_INEX) + { + /* In case of tie in GMP_RNDN, one must round away from 0. + Note: One could do better, by using GMP_RNDNA or similar, + but I want to keep the change minimal. Anyway this case + should be quite rare in practice. */ + mpfr_add_one_ulp (z, rnd_mode); + return 1; + } return (inexact) ? inexact : -1; } } diff -Naur mpfr-2.1.1-p7/tests/thypot.c mpfr-2.1.1-p8/tests/thypot.c --- mpfr-2.1.1-p7/tests/thypot.c 2005-01-27 17:32:18.000000000 +0000 +++ mpfr-2.1.1-p8/tests/thypot.c 2005-05-12 13:58:44.034394000 +0000 @@ -106,6 +106,47 @@ mpfr_clear (t); } +static void +test_large_small (void) +{ + mpfr_t x, y, z; + int inexact; + + mpfr_init2 (x, 3); + mpfr_init2 (y, 2); + mpfr_init2 (z, 2); + + mpfr_set_ui_2exp (x, 1, mpfr_get_emax () / 2, GMP_RNDN); + mpfr_set_ui_2exp (y, 1, -1, GMP_RNDN); + inexact = mpfr_hypot (z, x, y, GMP_RNDN); + if (inexact >= 0 || mpfr_cmp (x, z)) + { + printf ("Error 1 in test_large_small\n"); + exit (1); + } + + mpfr_mul_ui (x, x, 5, GMP_RNDN); + inexact = mpfr_hypot (z, x, y, GMP_RNDN); + if (mpfr_cmp (x, z) >= 0) + { + printf ("Error 2 in test_large_small\n"); + printf ("x = "); + mpfr_out_str (stdout, 2, 0, x, GMP_RNDN); + printf ("\n"); + printf ("y = "); + mpfr_out_str (stdout, 2, 0, y, GMP_RNDN); + printf ("\n"); + printf ("z = "); + mpfr_out_str (stdout, 2, 0, z, GMP_RNDN); + printf (" (in precision 2)\n"); + exit (1); + } + + mpfr_clear (x); + mpfr_clear (y); + mpfr_clear (z); +} + int main (int argc, char *argv[]) { @@ -204,6 +245,7 @@ mpfr_clear (t); test_large (); + test_large_small (); tests_end_mpfr (); return 0; diff -Naur mpfr-2.1.1-p8/add1.c mpfr-2.1.1-p9/add1.c --- mpfr-2.1.1-p8/add1.c 2004-03-26 13:27:15.000000000 +0000 +++ mpfr-2.1.1-p9/add1.c 2005-06-08 11:37:27.000000000 +0000 @@ -1,6 +1,6 @@ /* mpfr_add1 -- internal function to perform a "real" addition -Copyright 1999, 2000, 2001, 2002, 2003, 2004 Free Software Foundation. +Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005 Free Software Foundation. Contributed by the Spaces project, INRIA Lorraine. This file is part of the MPFR Library. @@ -523,14 +523,13 @@ add_one_ulp: /* add one unit in last place to a */ if (MPFR_UNLIKELY(mpn_add_1(ap, ap, an, MPFR_LIMB_ONE << sh))) { - /* Case 100000x0 + 1*/ if (MPFR_UNLIKELY(exp == __gmpfr_emax)) - inex = mpfr_set_overflow(a, rnd_mode, MPFR_SIGN(a)); - else { - exp++; - ap[an-1] = MPFR_LIMB_HIGHBIT; + inex = mpfr_set_overflow(a, rnd_mode, MPFR_SIGN(a)); + goto end_of_add; } + exp++; + ap[an-1] = MPFR_LIMB_HIGHBIT; } set_exponent: diff -Naur mpfr-2.1.1-p9/atan.c mpfr-2.1.1-p10/atan.c --- mpfr-2.1.1-p9/atan.c 2004-06-16 12:08:30.000000000 +0000 +++ mpfr-2.1.1-p10/atan.c 2005-07-05 13:00:43.000000000 +0000 @@ -1,6 +1,6 @@ /* mpfr_atan -- arc-tangent of a floating-point number -Copyright 2001, 2002, 2003, 2004 Free Software Foundation. +Copyright 2001, 2002, 2003, 2004, 2005 Free Software Foundation. This file is part of the MPFR Library, and was contributed by Mathieu Dutour. @@ -107,7 +107,8 @@ else /* x is necessarily 0 */ { MPFR_ASSERTD(MPFR_IS_ZERO(x)); - mpfr_set_ui (arctangent, 0, GMP_RNDN); + MPFR_SET_ZERO (arctangent); + MPFR_SET_SAME_SIGN (arctangent, x); return 0; /* exact result */ } } diff -Naur mpfr-2.1.1-p9/tests/tatan.c mpfr-2.1.1-p10/tests/tatan.c --- mpfr-2.1.1-p9/tests/tatan.c 2005-02-04 12:28:53.000000000 +0000 +++ mpfr-2.1.1-p10/tests/tatan.c 2005-07-05 12:57:44.000000000 +0000 @@ -93,19 +93,33 @@ /* atan(+/-0) = +/-0 */ mpfr_set_ui (x, 0, GMP_RNDN); + MPFR_SET_NEG (y); mpfr_atan (y, x, GMP_RNDN); - if (mpfr_cmp_ui (y, 0) || mpfr_sgn (y) < 0) + if (mpfr_cmp_ui (y, 0) || MPFR_IS_NEG (y)) { printf ("Error: mpfr_atan (+0) <> +0\n"); exit (1); } + mpfr_atan (x, x, GMP_RNDN); + if (mpfr_cmp_ui (x, 0) || MPFR_IS_NEG (x)) + { + printf ("Error: mpfr_atan (+0) <> +0 (in place)\n"); + exit (1); + } mpfr_neg (x, x, GMP_RNDN); + MPFR_SET_POS (y); mpfr_atan (y, x, GMP_RNDN); - if (mpfr_cmp_ui (y, 0) || mpfr_sgn (y) > 0) + if (mpfr_cmp_ui (y, 0) || MPFR_IS_POS (y)) { printf ("Error: mpfr_atan (-0) <> -0\n"); exit (1); } + mpfr_atan (x, x, GMP_RNDN); + if (mpfr_cmp_ui (x, 0) || MPFR_IS_POS (x)) + { + printf ("Error: mpfr_atan (-0) <> -0 (in place)\n"); + exit (1); + } mpfr_set_prec (x, 32); mpfr_set_prec (y, 32); diff -Naur mpfr-2.1.1-p10/strtofr.c mpfr-2.1.1-p11/strtofr.c --- mpfr-2.1.1-p10/strtofr.c 2005-01-27 14:47:26.000000000 +0000 +++ mpfr-2.1.1-p11/strtofr.c 2005-07-11 17:38:34.000000000 +0000 @@ -1,6 +1,6 @@ /* mpfr_strtofr -- set a floating-point number from a string -Copyright 2004 Free Software Foundation, Inc. +Copyright 2004, 2005 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -186,6 +186,27 @@ return MPFR_LIKELY (digit < base) ? digit : -1; } +/* Compatible with any locale, but one still assumes that 'a', 'b', 'c', + ..., 'z', and 'A', 'B', 'C', ..., 'Z' are consecutive values (like + in any ASCII-based character set). */ +static int +fast_casecmp (const char *s1, const char *s2) +{ + unsigned char c1, c2; + + do + { + c2 = *(const unsigned char *) s2++; + if (c2 == '\0') + return 0; + c1 = *(const unsigned char *) s1++; + if (c1 >= 'A' && c1 <= 'Z') + c1 = c1 - 'A' + 'a'; + } + while (c1 == c2); + return 1; +} + /* Parse a string and fill pstr. Return the advanced ptr too. It returns: @@ -214,12 +235,12 @@ while (isspace((unsigned char) *str)) str++; /* Can be case-insensitive NAN */ - if (strncasecmp (str, "@nan@", 5) == 0) + if (fast_casecmp (str, "@nan@") == 0) { str += 5; goto set_nan; } - if (base <= 16 && strncasecmp (str, "nan", 3) == 0) + if (base <= 16 && fast_casecmp (str, "nan") == 0) { str += 3; set_nan: @@ -249,17 +270,17 @@ str++; /* Can be case-insensitive INF */ - if (strncasecmp (str, "@inf@", 5) == 0) + if (fast_casecmp (str, "@inf@") == 0) { str += 5; goto set_inf; } - if (base <= 16 && strncasecmp (str, "infinity", 8) == 0) + if (base <= 16 && fast_casecmp (str, "infinity") == 0) { str += 8; goto set_inf; } - if (base <= 16 && strncasecmp (str, "inf", 3) == 0) + if (base <= 16 && fast_casecmp (str, "inf") == 0) { str += 3; set_inf: diff -Naur mpfr-2.1.1-p11/set_ld.c mpfr-2.1.1-p12/set_ld.c --- mpfr-2.1.1-p11/set_ld.c 2005-04-21 14:12:23.000000000 +0000 +++ mpfr-2.1.1-p12/set_ld.c 2005-07-29 14:22:46.000000000 +0000 @@ -46,7 +46,8 @@ mpfr_set_ld (mpfr_ptr r, long double d, mp_rnd_t rnd_mode) { mpfr_t t, u; - int inexact, shift_exp = 0, inexact2 = 0; + int inexact, shift_exp = 0; + long double x; LONGDOUBLE_NAN_ACTION (d, goto nan); @@ -67,11 +68,14 @@ mpfr_init2 (t, MPFR_LDBL_MANT_DIG); mpfr_init2 (u, IEEE_DBL_MANT_DIG); - mpfr_set_ui (t, 0, GMP_RNDN); mpfr_save_emin_emax (); - while (d != (long double) 0.0) + + convert: + x = d; + mpfr_set_ui (t, 0, GMP_RNDN); + while (x != (long double) 0.0) { - if ((d > (long double) DBL_MAX) || ((-d) > (long double) DBL_MAX)) + if ((x > (long double) DBL_MAX) || ((-x) > (long double) DBL_MAX)) { long double div9, div10, div11, div12, div13; @@ -83,31 +87,31 @@ div11 = div10 * div10; /* 2^(2^11) */ div12 = div11 * div11; /* 2^(2^12) */ div13 = div12 * div12; /* 2^(2^13) */ - if (ABS(d) >= div13) + if (ABS(x) >= div13) { - d = d / div13; /* exact */ + x /= div13; /* exact */ shift_exp += 8192; } - if (ABS(d) >= div12) + if (ABS(x) >= div12) { - d = d / div12; /* exact */ + x /= div12; /* exact */ shift_exp += 4096; } - if (ABS(d) >= div11) + if (ABS(x) >= div11) { - d = d / div11; /* exact */ + x /= div11; /* exact */ shift_exp += 2048; } - if (ABS(d) >= div10) + if (ABS(x) >= div10) { - d = d / div10; /* exact */ + x /= div10; /* exact */ shift_exp += 1024; } - /* warning: we may have DBL_MAX=2^1024*(1-2^(-53)) < d < 2^1024, + /* warning: we may have DBL_MAX=2^1024*(1-2^(-53)) < x < 2^1024, therefore we have one extra exponent reduction step */ - if (ABS(d) >= div9) + if (ABS(x) >= div9) { - d = d / div9; /* exact */ + x /= div9; /* exact */ shift_exp += 512; } } @@ -118,63 +122,85 @@ /* div9 = 2^(-2^9) */ div10 = div9 * div9; /* 2^(-2^10) */ div11 = div10 * div10; /* 2^(-2^11) if extended precision */ - /* since -DBL_MAX <= d <= DBL_MAX, the cast to double should not + /* since -DBL_MAX <= x <= DBL_MAX, the cast to double should not overflow here */ - inexact = mpfr_set_d (u, (double) d, GMP_RNDN); + inexact = mpfr_set_d (u, (double) x, GMP_RNDZ); MPFR_ASSERTD(inexact == 0); - if (d != (long double) 0.0 && - ABS(d) < div10 && + if (x != (long double) 0.0 && + ABS(x) < div10 && div11 != (long double) 0.0 && div11 / div10 == div10) /* possible underflow */ { long double div12, div13; - /* After the divisions, any bit of d must be >= div10, + /* After the divisions, any bit of x must be >= div10, hence the possible division by div9. */ div12 = div11 * div11; /* 2^(-2^12) */ div13 = div12 * div12; /* 2^(-2^13) */ - if (ABS(d) <= div13) + if (ABS(x) <= div13) { - d = d / div13; /* exact */ + x /= div13; /* exact */ shift_exp -= 8192; } - if (ABS(d) <= div12) + if (ABS(x) <= div12) { - d = d / div12; /* exact */ + x /= div12; /* exact */ shift_exp -= 4096; } - if (ABS(d) <= div11) + if (ABS(x) <= div11) { - d = d / div11; /* exact */ + x /= div11; /* exact */ shift_exp -= 2048; } - if (ABS(d) <= div10) + if (ABS(x) <= div10) { - d = d / div10; /* exact */ + x /= div10; /* exact */ shift_exp -= 1024; } - if (ABS(d) <= div9) + if (ABS(x) <= div9) { - d = d / div9; /* exact */ + x /= div9; /* exact */ shift_exp -= 512; } } else { - mpfr_add (t, t, u, GMP_RNDN); /* exact */ - if (!mpfr_number_p (t)) - break; - d = d - (long double) mpfr_get_d1 (u); /* exact */ + if (mpfr_add (t, t, u, GMP_RNDZ) != 0) + { + if (!mpfr_number_p (t)) + break; + /* Inexact. This cannot happen unless the C implementation + "lies" on the precision or when long doubles are + implemented with FP expansions like under Mac OS X. */ + if (MPFR_PREC (t) != MPFR_PREC (r) + 1) + { + /* We assume that MPFR_PREC (r) < MPFR_PREC_MAX. + The precision MPFR_PREC (r) + 1 allows us to + deduce the rounding bit and the sticky bit. */ + mpfr_set_prec (t, MPFR_PREC (r) + 1); + goto convert; + } + else + { + mp_limb_t *tp; + int rb_mask; + + /* Since mpfr_add was inexact, the sticky bit is 1. */ + tp = MPFR_MANT (t); + rb_mask = MPFR_LIMB_ONE << + (BITS_PER_MP_LIMB - 1 - + (MPFR_PREC (r) & (BITS_PER_MP_LIMB - 1))); + if (rnd_mode == GMP_RNDN) + rnd_mode = (*tp & rb_mask) ^ MPFR_IS_NEG (t) ? + GMP_RNDU : GMP_RNDD; + *tp |= rb_mask; + break; + } + } + x -= (long double) mpfr_get_d1 (u); /* exact */ } } } - /* now t is exactly the input value d */ - inexact = mpfr_set (r, t, rnd_mode); - if (shift_exp > 0) - inexact2 = mpfr_mul_2exp (r, r, shift_exp, rnd_mode); - else if (shift_exp < 0) - inexact2 = mpfr_div_2exp (r, r, -shift_exp, rnd_mode); - if (inexact2) /* overflow */ - inexact = inexact2; + inexact = mpfr_mul_2si (r, t, shift_exp, rnd_mode); mpfr_clear (t); mpfr_clear (u); mpfr_restore_emin_emax (); diff -Naur mpfr-2.1.1-p11/tests/tset_ld.c mpfr-2.1.1-p12/tests/tset_ld.c --- mpfr-2.1.1-p11/tests/tset_ld.c 2005-04-21 13:00:38.000000000 +0000 +++ mpfr-2.1.1-p12/tests/tset_ld.c 2005-07-29 14:56:38.000000000 +0000 @@ -155,7 +155,7 @@ mpfr_test_init (); tests_machine_prec_long_double (); - mpfr_init2 (x, 113); + mpfr_init2 (x, MPFR_LDBL_MANT_DIG); /* check +0.0 and -0.0 */ d = 0.0; @@ -210,7 +210,7 @@ /* checks that 2^i, 2^i+1 and 2^i-1 are correctly converted */ d = 1.0; - for (i = 1; i <= 113; i++) + for (i = 1; i < MPFR_LDBL_MANT_DIG; i++) { d = 2.0 * d; /* d = 2^i */ check_set_get (d, x);