diff -Naurd mpfr-2.2.0-p4/cache.c mpfr-2.2.0-p5/cache.c --- mpfr-2.2.0-p4/cache.c 2005-08-18 16:35:03.000000000 +0000 +++ mpfr-2.2.0-p5/cache.c 2005-11-23 09:04:29.000000000 +0000 @@ -32,9 +32,9 @@ void mpfr_clear_cache (mpfr_cache_t cache) { - if (MPFR_PREC(cache->x) != 0) + if (MPFR_PREC (cache->x) != 0) mpfr_clear (cache->x); - MPFR_PREC(cache->x) = 0; + MPFR_PREC (cache->x) = 0; } int @@ -47,45 +47,58 @@ MPFR_SAVE_EXPO_MARK (expo); - /* Check if the cache has been already filled */ - if (MPFR_UNLIKELY(pold == 0)) - mpfr_init2 (cache->x, prec); + if (MPFR_UNLIKELY (prec > pold)) + { + /* No previous result in the cache or the precision of the + previous result is not sufficient. */ - /* Check if we can round with the previous result */ - else if (MPFR_LIKELY(prec <= pold)) - goto round; + if (MPFR_UNLIKELY (pold == 0)) /* No previous result. */ + mpfr_init2 (cache->x, prec); - /* Update the cache */ - pold = prec /*MPFR_PREC_MIN + prec + __gmpfr_ceil_exp2 (prec)*/; - MPFR_ASSERTD (pold >= prec); - mpfr_prec_round (cache->x, pold, GMP_RNDN); - cache->inexact = (*cache->func) (cache->x, GMP_RNDN); + /* Update the cache. */ + pold = prec; + mpfr_prec_round (cache->x, pold, GMP_RNDN); + cache->inexact = (*cache->func) (cache->x, GMP_RNDN); + } - round: - /* First check if the cache has the exact value (Unlikely) - Else the exact value is between (assuming x=cache->x > 0) - x and x+ulp(x) if cache->inexact < 0 - x-ulp(x) and x if cache->inexact > 0 - and abs(x-exact) <= ulp(x)/2 */ - MPFR_ASSERTD (MPFR_IS_POS(cache->x)); /* TODO...*/ - /* We must use nextbelow instead of sub_one_ulp, since we know - that the exact value is < 1/2ulp(x) (We want sub_demi_ulp(x)). - Can't use mpfr_set since we need the even flag. */ + /* First, check if the cache has the exact value (unlikely). + Else the exact value is between (assuming x=cache->x > 0): + x and x+ulp(x) if cache->inexact < 0, + x-ulp(x) and x if cache->inexact > 0, + and abs(x-exact) <= ulp(x)/2. */ + MPFR_ASSERTN (MPFR_IS_POS (cache->x)); /* TODO... */ sign = MPFR_SIGN (cache->x); MPFR_SET_EXP (dest, MPFR_GET_EXP (cache->x)); MPFR_SET_SIGN (dest, sign); - MPFR_RNDRAW_EVEN (inexact, dest, - MPFR_MANT (cache->x), MPFR_PREC (cache->x), rnd, sign, - if (MPFR_UNLIKELY ( ++MPFR_EXP (dest) > __gmpfr_emax)) - mpfr_overflow (dest, rnd, sign) ); - /* inexact = mpfr_set (dest, cache->x, rnd); */ - if (MPFR_LIKELY(cache->inexact != 0)) + MPFR_RNDRAW_GEN (inexact, dest, + MPFR_MANT (cache->x), MPFR_PREC (cache->x), rnd, sign, + if (MPFR_UNLIKELY (cache->inexact == 0)) + { + if ((sp[0] & ulp) == 0) + { + inexact = -sign; + goto trunc_doit; + } + else + goto addoneulp; + } + else if (cache->inexact < 0) + goto addoneulp; + else + { + inexact = -sign; + goto trunc_doit; + }, + if (MPFR_UNLIKELY (++MPFR_EXP (dest) > __gmpfr_emax)) + mpfr_overflow (dest, rnd, sign); + ); + if (MPFR_LIKELY (cache->inexact != 0)) { switch (rnd) { case GMP_RNDZ: case GMP_RNDD: - if (MPFR_UNLIKELY(inexact == 0)) + if (MPFR_UNLIKELY (inexact == 0)) { inexact = cache->inexact; if (inexact > 0) @@ -93,7 +106,7 @@ } break; case GMP_RNDU: - if (MPFR_UNLIKELY(inexact == 0)) + if (MPFR_UNLIKELY (inexact == 0)) { inexact = cache->inexact; if (inexact < 0) @@ -101,16 +114,7 @@ } break; default: /* GMP_RNDN */ - if (MPFR_UNLIKELY(inexact == MPFR_EVEN_INEX || - inexact == -MPFR_EVEN_INEX)) - { - if (cache->inexact < 0) - mpfr_nextabove (dest); - else - mpfr_nextbelow (dest); - inexact = -inexact; - } - else if (MPFR_UNLIKELY(inexact == 0)) + if (MPFR_UNLIKELY(inexact == 0)) inexact = cache->inexact; break; } diff -Naurd mpfr-2.2.0-p4/hypot.c mpfr-2.2.0-p5/hypot.c --- mpfr-2.2.0-p4/hypot.c 2005-06-06 13:39:39.000000000 +0000 +++ mpfr-2.2.0-p5/hypot.c 2005-11-23 09:04:29.000000000 +0000 @@ -73,6 +73,7 @@ Ey = MPFR_GET_EXP (y); diff_exp = (mp_exp_unsigned_t) Ex - Ey; + Nx = MPFR_PREC (x); /* Precision of input variable */ Nz = MPFR_PREC (z); /* Precision of output variable */ /* we have x < 2^Ex thus x^2 < 2^(2*Ex), @@ -87,10 +88,10 @@ 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 (diff_exp > MAX (Nz, MPFR_PREC (x)) / 2) + if (diff_exp > MAX (Nz, Nx) / 2) /* result is |x| or |x|+ulp(|x|,Nz) */ { - if (rnd_mode == GMP_RNDU) + if (MPFR_UNLIKELY (rnd_mode == GMP_RNDU)) { /* If z > abs(x), then it was already rounded up; otherwise z = abs(x), and we need to add one ulp due to y. */ @@ -100,14 +101,27 @@ } else /* GMP_RNDZ, GMP_RNDD, GMP_RNDN */ { - inexact = mpfr_abs (z, x, rnd_mode); - return (inexact) ? inexact : -1; + if (MPFR_LIKELY (Nz >= Nx)) + { + mpfr_abs (z, x, rnd_mode); /* exact */ + return -1; + } + else + { + MPFR_SET_EXP (z, Ex); + MPFR_SET_SIGN (z, 1); + MPFR_RNDRAW_GEN (inexact, z, MPFR_MANT (x), Nx, rnd_mode, 1, + goto addoneulp, + if (MPFR_UNLIKELY (++MPFR_EXP (z) > __gmpfr_emax)) + return mpfr_overflow (z, rnd_mode, 1); + ); + return inexact ? inexact : -1; + } } } /* General case */ - Nx = MPFR_PREC(x); /* Precision of input variable */ Ny = MPFR_PREC(y); /* Precision of input variable */ /* compute the working precision -- see algorithms.ps */ diff -Naurd mpfr-2.2.0-p4/mpfr-impl.h mpfr-2.2.0-p5/mpfr-impl.h --- mpfr-2.2.0-p4/mpfr-impl.h 2005-09-11 22:13:24.000000000 +0000 +++ mpfr-2.2.0-p5/mpfr-impl.h 2005-11-23 09:04:29.000000000 +0000 @@ -184,6 +184,14 @@ # define MPFR_THREAD_ATTR #endif +/* Cache struct */ +struct __gmpfr_cache_s { + mpfr_t x; + int inexact; + int (*func)(mpfr_ptr, mpfr_rnd_t); +}; +typedef struct __gmpfr_cache_s mpfr_cache_t[1]; + #if defined (__cplusplus) extern "C" { #endif @@ -924,113 +932,19 @@ ******************************************************/ /* - * Round Mantissa (`srcp`, `sprec`) to mpfr_t `dest` using rounding mode `rnd` - * assuming dest's sign is `sign`. - * Execute OVERFLOW_HANDLER in case of overflow when rounding (Power 2 case) + * Note: due to the labels, one cannot use a macro MPFR_RNDRAW* more than + * once in a function (otherwise these labels would not be unique). */ -#define MPFR_RNDRAW(inexact, dest, srcp, sprec, rnd, sign, OVERFLOW_HANDLER)\ - do { \ - mp_size_t dests, srcs; \ - mp_limb_t *destp; \ - mp_prec_t destprec, srcprec; \ - \ - /* Check Trivial Case when Dest Mantissa has more bits than source */ \ - srcprec = sprec; \ - destprec = MPFR_PREC (dest); \ - destp = MPFR_MANT (dest); \ - if (MPFR_UNLIKELY (destprec >= srcprec)) \ - { \ - srcs = (srcprec + BITS_PER_MP_LIMB-1)/BITS_PER_MP_LIMB; \ - dests = (destprec + BITS_PER_MP_LIMB-1)/BITS_PER_MP_LIMB - srcs; \ - MPN_COPY (destp + dests, srcp, srcs); \ - MPN_ZERO (destp, dests); \ - inexact = 0; \ - } \ - else \ - { \ - /* Non trivial case: rounding needed */ \ - mp_prec_t sh; \ - mp_limb_t *sp; \ - mp_limb_t rb, sb, ulp; \ - \ - /* Compute Position and shift */ \ - srcs = (srcprec + BITS_PER_MP_LIMB-1)/BITS_PER_MP_LIMB; \ - dests = (destprec + BITS_PER_MP_LIMB-1)/BITS_PER_MP_LIMB; \ - MPFR_UNSIGNED_MINUS_MODULO (sh, destprec); \ - sp = srcp + srcs - dests; \ - \ - /* General case when prec % BITS_PER_MP_LIMB != 0 */ \ - if (MPFR_LIKELY (sh != 0)) \ - { \ - mp_limb_t mask; \ - /* Compute Rounding Bit and Sticky Bit */ \ - mask = MPFR_LIMB_ONE << (sh-1); \ - rb = sp[0] & mask; \ - sb = sp[0] & (mask-1); \ - if (MPFR_UNLIKELY (sb == 0)) \ - { /* TODO: Improve it */ \ - mp_limb_t *tmp; \ - mp_size_t n; \ - for (tmp = sp, n = srcs - dests ; n != 0 && sb == 0 ; n--) \ - sb = *--tmp; \ - } \ - ulp = 2*mask; \ - } \ - else /* sh == 0 */ \ - { \ - MPFR_ASSERTD (dests < srcs); \ - /* Compute Rounding Bit and Sticky Bit */ \ - rb = sp[-1] & MPFR_LIMB_HIGHBIT; \ - sb = sp[-1] & (MPFR_LIMB_HIGHBIT-1); \ - if (MPFR_UNLIKELY (sb == 0)) \ - { \ - mp_limb_t *tmp; \ - mp_size_t n; \ - for (tmp = sp-1, n = srcs - dests-1 ; n!=0 && sb==0 ; n--) \ - sb = *--tmp; \ - } \ - ulp = MPFR_LIMB_ONE; \ - } \ - /* Rounding */ \ - if (MPFR_LIKELY (rnd == GMP_RNDN)) \ - { \ - if (rb == 0 || MPFR_UNLIKELY (sb == 0 && (sp[0] & ulp) == 0)) \ - { \ - trunc: \ - inexact = MPFR_LIKELY ((sb | rb) != 0) ? -sign : 0; \ - MPN_COPY (destp, sp, dests); \ - destp[0] &= ~(ulp-1); \ - } \ - else \ - { \ - addoneulp: \ - if (MPFR_UNLIKELY (mpn_add_1 (destp, sp, dests, ulp))) \ - { \ - destp[dests-1] = MPFR_LIMB_HIGHBIT; \ - OVERFLOW_HANDLER; \ - } \ - destp[0] &= ~(ulp-1); \ - inexact = sign; \ - } \ - } \ - else \ - { /* Not Rounding to Nearest */ \ - if (MPFR_LIKELY (MPFR_IS_LIKE_RNDZ (rnd, MPFR_IS_NEG_SIGN (sign)))\ - || MPFR_UNLIKELY ((sb | rb) == 0)) \ - goto trunc; \ - else \ - goto addoneulp; \ - } \ - } \ - } while (0) /* - * Round Mantissa (`srcp`, `sprec`) to mpfr_t `dest` using rounding mode `rnd` - * assuming dest's sign is `sign`. - * Execute OVERFLOW_HANDLER in case of overflow when rounding (Power 2 case) - * Return MPFR_EVEN_INEX in case of EVEN rounding + * Round mantissa (srcp, sprec) to mpfr_t dest using rounding mode rnd + * assuming dest's sign is sign. + * In rounding to nearest mode, execute MIDDLE_HANDLER when the value + * is the middle of two consecutive numbers in dest precision. + * Execute OVERFLOW_HANDLER in case of overflow when rounding. */ -#define MPFR_RNDRAW_EVEN(inexact, dest, srcp, sprec, rnd, sign, OVERFLOW_HANDLER)\ +#define MPFR_RNDRAW_GEN(inexact, dest, srcp, sprec, rnd, sign, \ + MIDDLE_HANDLER, OVERFLOW_HANDLER) \ do { \ mp_size_t dests, srcs; \ mp_limb_t *destp; \ @@ -1105,19 +1019,8 @@ destp[0] &= ~(ulp-1); \ } \ else if (MPFR_UNLIKELY (sb == 0)) \ - { \ - /* EVEN rounding */ \ - if ((sp[0] & ulp) == 0) \ - { \ - MPFR_ASSERTD (rb != 0); \ - inexact = -MPFR_EVEN_INEX*sign; \ - goto trunc_doit; \ - } \ - else \ - { \ - inexact = MPFR_EVEN_INEX*sign; \ - goto addoneulp_doit; \ - } \ + { /* Middle of two consecutive representable numbers */ \ + MIDDLE_HANDLER; \ } \ else \ { \ @@ -1133,16 +1036,58 @@ } \ } \ else \ - { /* Not Rounding to Nearest */ \ - if (MPFR_LIKELY (MPFR_IS_LIKE_RNDZ (rnd, MPFR_IS_NEG_SIGN (sign)))\ - || MPFR_UNLIKELY ((sb | rb) == 0)) \ + { /* Directed rounding mode */ \ + if (MPFR_LIKELY (MPFR_IS_LIKE_RNDZ (rnd, \ + MPFR_IS_NEG_SIGN (sign)))) \ goto trunc; \ + else if (MPFR_UNLIKELY ((sb | rb) == 0)) \ + { \ + inexact = 0; \ + goto trunc_doit; \ + } \ else \ goto addoneulp; \ } \ } \ } while (0) +/* + * Round mantissa (srcp, sprec) to mpfr_t dest using rounding mode rnd + * assuming dest's sign is sign. + * Execute OVERFLOW_HANDLER in case of overflow when rounding. + */ +#define MPFR_RNDRAW(inexact, dest, srcp, sprec, rnd, sign, OVERFLOW_HANDLER) \ + MPFR_RNDRAW_GEN (inexact, dest, srcp, sprec, rnd, sign, \ + if ((sp[0] & ulp) == 0) \ + { \ + inexact = -sign; \ + goto trunc_doit; \ + } \ + else \ + goto addoneulp; \ + , OVERFLOW_HANDLER) + +/* + * Round mantissa (srcp, sprec) to mpfr_t dest using rounding mode rnd + * assuming dest's sign is sign. + * Execute OVERFLOW_HANDLER in case of overflow when rounding. + * Set inexact to +/- MPFR_EVEN_INEX in case of even rounding. + */ +#define MPFR_RNDRAW_EVEN(inexact, dest, srcp, sprec, rnd, sign, \ + OVERFLOW_HANDLER) \ + MPFR_RNDRAW_GEN (inexact, dest, srcp, sprec, rnd, sign, \ + if ((sp[0] & ulp) == 0) \ + { \ + inexact = -MPFR_EVEN_INEX * sign; \ + goto trunc_doit; \ + } \ + else \ + { \ + inexact = MPFR_EVEN_INEX * sign; \ + goto addoneulp_doit; \ + } \ + , OVERFLOW_HANDLER) + /* Return TRUE if b is non singular and we can round it to precision 'prec' with rounding mode 'rnd', and with error at most 'error' */ #define MPFR_CAN_ROUND(b,err,prec,rnd) \ @@ -1498,6 +1443,13 @@ __MPFR_DECLSPEC int mpfr_const_log2_internal _MPFR_PROTO((mpfr_ptr,mp_rnd_t)); __MPFR_DECLSPEC int mpfr_const_euler_internal _MPFR_PROTO((mpfr_ptr, mp_rnd_t)); __MPFR_DECLSPEC int mpfr_const_catalan_internal _MPFR_PROTO((mpfr_ptr, mp_rnd_t)); + +__MPFR_DECLSPEC void mpfr_init_cache _MPFR_PROTO ((mpfr_cache_t, + int(*)(mpfr_ptr,mpfr_rnd_t))); +__MPFR_DECLSPEC void mpfr_clear_cache _MPFR_PROTO ((mpfr_cache_t)); +__MPFR_DECLSPEC int mpfr_cache _MPFR_PROTO ((mpfr_ptr, mpfr_cache_t, + mpfr_rnd_t)); + __MPFR_DECLSPEC void mpfr_mulhigh_n _MPFR_PROTO ((mp_ptr, mp_srcptr, mp_srcptr, mp_size_t)); diff -Naurd mpfr-2.2.0-p4/mpfr.h mpfr-2.2.0-p5/mpfr.h --- mpfr-2.2.0-p4/mpfr.h 2005-09-29 11:36:36.000000000 +0000 +++ mpfr-2.2.0-p5/mpfr.h 2005-11-23 09:04:29.000000000 +0000 @@ -118,14 +118,6 @@ /* For those who needs a direct access and fast access to the sign field */ #define MPFR_SIGN(x) ((x)->_mpfr_sign) -/* Cache struct */ -struct __gmpfr_cache_s { - mpfr_t x; - int inexact; - int (*func)(mpfr_ptr, mpfr_rnd_t); -}; -typedef struct __gmpfr_cache_s mpfr_cache_t[1]; - /* Stack interface */ typedef enum { MPFR_NAN_KIND = 0, @@ -544,11 +536,6 @@ __MPFR_DECLSPEC int mpfr_sum _MPFR_PROTO ((mpfr_ptr, mpfr_ptr *__gmp_const, unsigned long, mpfr_rnd_t)); -__MPFR_DECLSPEC void mpfr_init_cache _MPFR_PROTO ((mpfr_cache_t, - int(*)(mpfr_ptr,mpfr_rnd_t))); -__MPFR_DECLSPEC void mpfr_clear_cache _MPFR_PROTO ((mpfr_cache_t)); -__MPFR_DECLSPEC int mpfr_cache _MPFR_PROTO ((mpfr_ptr, mpfr_cache_t, - mpfr_rnd_t)); __MPFR_DECLSPEC void mpfr_free_cache _MPFR_PROTO ((void)); __MPFR_DECLSPEC int mpfr_subnormalize _MPFR_PROTO ((mpfr_ptr, int, diff -Naurd mpfr-2.2.0-p4/round_near_x.c mpfr-2.2.0-p5/round_near_x.c --- mpfr-2.2.0-p4/round_near_x.c 2005-08-18 16:35:12.000000000 +0000 +++ mpfr-2.2.0-p5/round_near_x.c 2005-11-23 09:04:29.000000000 +0000 @@ -21,13 +21,13 @@ #include "mpfr-impl.h" -/* Uses MPFR_FAST_COMPUTE_IF_SMALL_INPUT instead (a simple wrapper) */ +/* Use MPFR_FAST_COMPUTE_IF_SMALL_INPUT instead (a simple wrapper) */ /* int mpfr_round_near_x (mpfr_ptr y, mpfr_srcptr x, mp_exp_t err, int dir, mp_rnd_t rnd) Assuming y = o(f(x)) = o(x + g(x)) with |g(x)| < 2^(EXP(x)-error) - If x is small enought, y ~= x. This function checks and does this. + If x is small enough, y ~= x. This function checks and does this. It assumes that f(x) is not representable exactly as a FP number. x must not be a singular value (NAN, INF or ZERO). @@ -42,7 +42,7 @@ Otherwise it returns the ternary flag (It can't return an exact value). */ -/* What "small enought" means? +/* What "small enough" means? We work with the positive values. Assuming err > Prec (y)+1 @@ -50,7 +50,7 @@ i = [ y = o(x)] // i = inexact flag If i == 0 Setting x in y is exact. We have: - y = [XXXXXXXXX[...]]0[...] + error where [..] are optionnal zeros + y = [XXXXXXXXX[...]]0[...] + error where [..] are optional zeros if dirError = ToInf, x < f(x) < x + 2^(EXP(x)-err) since x=y, and ulp (y)/2 > 2^(EXP(x)-err), we have: @@ -172,9 +172,17 @@ sign = MPFR_SIGN (x); MPFR_SET_EXP (y, MPFR_GET_EXP (x)); MPFR_SET_SIGN (y, sign); - MPFR_RNDRAW_EVEN (inexact, y, MPFR_MANT (x), MPFR_PREC (x), rnd, sign, - if (MPFR_UNLIKELY ( ++MPFR_EXP (y) > __gmpfr_emax)) - mpfr_overflow (y, rnd, sign) ); + MPFR_RNDRAW_GEN (inexact, y, MPFR_MANT (x), MPFR_PREC (x), rnd, sign, + if (dir == 0) + { + inexact = -sign; + goto trunc_doit; + } + else + goto addoneulp; + , if (MPFR_UNLIKELY (++MPFR_EXP (y) > __gmpfr_emax)) + mpfr_overflow (y, rnd, sign) + ); /* Fix it in some cases */ MPFR_ASSERTD (!MPFR_IS_NAN (y) && !MPFR_IS_ZERO (y)); @@ -212,15 +220,6 @@ } } } - /* The even rule has been used. But due to error term, we should never - use this rule. That's why we have to fix some wrong rounding */ - else if (inexact == MPFR_EVEN_INEX || inexact == -MPFR_EVEN_INEX) - { - if (inexact*sign > 0 && dir == 0) - goto nexttozero; - else if (inexact*sign < 0 && dir == 1) - goto nexttoinf; - } MPFR_RET (inexact); } diff -Naurd mpfr-2.2.0-p4/tests/thypot.c mpfr-2.2.0-p5/tests/thypot.c --- mpfr-2.2.0-p4/tests/thypot.c 2005-06-02 16:12:04.000000000 +0000 +++ mpfr-2.2.0-p5/tests/thypot.c 2005-11-23 09:04:29.000000000 +0000 @@ -131,7 +131,26 @@ inexact = mpfr_hypot (z, x, y, GMP_RNDN); if (inexact >= 0 || mpfr_cmp (x, z)) { - printf ("Error in test_large_small\n"); + 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) instead of\n "); + mpfr_out_str (stdout, 2, 2, x, GMP_RNDU); + printf ("\n"); exit (1); }