diff options
| -rw-r--r-- | div.c | 10 | ||||
| -rw-r--r-- | exceptions.c | 10 | ||||
| -rw-r--r-- | get_d.c | 4 | ||||
| -rw-r--r-- | get_str.c | 10 | ||||
| -rw-r--r-- | mpfr-impl.h | 53 | ||||
| -rw-r--r-- | mpfr.h | 40 | ||||
| -rw-r--r-- | mul.c | 11 | ||||
| -rw-r--r-- | round_prec.c | 154 | ||||
| -rw-r--r-- | round_raw_generic.c | 168 | ||||
| -rw-r--r-- | set_str.c | 4 | ||||
| -rw-r--r-- | sub1.c | 24 | ||||
| -rw-r--r-- | tests/Makefile.am | 3 |
12 files changed, 293 insertions, 198 deletions
@@ -354,9 +354,9 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mp_rnd_t rnd_mode) } } - cc = mpfr_round_raw_generic(qp, qp, err, - (MPFR_IS_NEG_SIGN(sign_quotient) ? 1 : 0), - MPFR_PREC(q), rnd_mode, &inex, 1); + cc = mpfr_round_raw_3(qp, qp, err, + (MPFR_IS_NEG_SIGN(sign_quotient) ? 1 : 0), + MPFR_PREC(q), rnd_mode, &inex); qp += qsize - MPFR_LIMB_SIZE(q); /* 0 or 1 */ qsize = MPFR_LIMB_SIZE(q); @@ -387,8 +387,8 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mp_rnd_t rnd_mode) if (k >= 0) /* Remainder is nonzero. */ { - if ((rnd_mode == GMP_RNDD && MPFR_IS_NEG_SIGN(sign_quotient)) - || (rnd_mode == GMP_RNDU && MPFR_IS_POS_SIGN(sign_quotient))) + if (MPFR_IS_RNDUTEST_OR_RNDDNOTTEST(rnd_mode, + MPFR_IS_POS_SIGN(sign_quotient))) /* Rounding to infinity. */ { inex = MPFR_FROM_SIGN_TO_INT( sign_quotient ); diff --git a/exceptions.c b/exceptions.c index fe6c90fcb..d35e6c10d 100644 --- a/exceptions.c +++ b/exceptions.c @@ -137,7 +137,7 @@ mpfr_check_range (mpfr_ptr x, int t, mp_rnd_t rnd_mode) if (rnd_mode == GMP_RNDN && (exp + 1 < __gmpfr_emin || (mpfr_powerof2_raw(x) && - (MPFR_SIGN(x) < 0 ? t <= 0 : t >= 0)))) + (MPFR_IS_NEG(x) ? t <= 0 : t >= 0)))) rnd_mode = GMP_RNDZ; return mpfr_set_underflow(x, rnd_mode, MPFR_SIGN(x)); } @@ -191,10 +191,10 @@ mpfr_set_underflow (mpfr_ptr x, mp_rnd_t rnd_mode, int sign) { int inex; + MPFR_ASSERT_SIGN(sign); MPFR_CLEAR_FLAGS(x); if (rnd_mode == GMP_RNDN - || (rnd_mode == GMP_RNDU && sign > 0) - || (rnd_mode == GMP_RNDD && sign < 0)) + || MPFR_IS_RNDUTEST_OR_RNDDNOTTEST(rnd_mode, sign > 0)) { mpfr_setmin (x, __gmpfr_emin); inex = 1; @@ -216,10 +216,10 @@ mpfr_set_overflow (mpfr_ptr x, mp_rnd_t rnd_mode, int sign) { int inex; + MPFR_ASSERT_SIGN(sign); MPFR_CLEAR_FLAGS(x); if (rnd_mode == GMP_RNDN - || (rnd_mode == GMP_RNDU && sign > 0) - || (rnd_mode == GMP_RNDD && sign < 0)) + || MPFR_IS_RNDUTEST_OR_RNDDNOTTEST(rnd_mode, sign > 0)) { MPFR_SET_INF(x); inex = 1; @@ -202,8 +202,8 @@ mpfr_get_d3 (mpfr_srcptr src, mp_exp_t e, mp_rnd_t rnd_mode) } np = (nbits - 1) / BITS_PER_MP_LIMB; tp = (mp_ptr) (*__gmp_allocate_func)((np+1) * BYTES_PER_MP_LIMB); - carry = mpfr_round_raw (tp, MPFR_MANT(src), MPFR_PREC(src), negative, - nbits, rnd_mode, (int *) 0); + carry = mpfr_round_raw_4 (tp, MPFR_MANT(src), MPFR_PREC(src), negative, + nbits, rnd_mode); if (carry) d = 1.0; else @@ -377,8 +377,8 @@ mpfr_get_str_aux (char *str, mp_exp_t *exp, mp_limb_t *r, mp_size_t n, i0 = (-f) / BITS_PER_MP_LIMB; j0 = (-f) % BITS_PER_MP_LIMB; - ret = mpfr_round_raw_generic (r + i0, r, n * BITS_PER_MP_LIMB, 0, - n * BITS_PER_MP_LIMB + f, rnd, &dir, 0); + ret = mpfr_round_raw (r + i0, r, n * BITS_PER_MP_LIMB, 0, + n * BITS_PER_MP_LIMB + f, rnd, &dir); /* warning: mpfr_round_raw_generic returns 2 or -2 in case of even rounding */ @@ -641,9 +641,9 @@ mpfr_get_str (char *s, mp_exp_t *e, int b, size_t m, mpfr_srcptr x, mp_rnd_t rnd nb = n * BITS_PER_MP_LIMB - prec; /* round xp to the precision prec, and put it into x1 put the carry into x1[n] */ - if ((x1[n] = mpfr_round_raw_generic (x1, xp, MPFR_PREC(x), - MPFR_IS_STRICTNEG(x), - prec, rnd, &inexp, 0))) + if ((x1[n] = mpfr_round_raw (x1, xp, MPFR_PREC(x), + MPFR_IS_STRICTNEG(x), + prec, rnd, &inexp))) { /* overflow when rounding x: x1 = 2^prec */ if (r == pow2) /* prec = m * pow2, diff --git a/mpfr-impl.h b/mpfr-impl.h index f70e8bbd2..602024a23 100644 --- a/mpfr-impl.h +++ b/mpfr-impl.h @@ -50,9 +50,13 @@ MA 02111-1307, USA. */ /* Definition of MPFR_LIMB_HIGHBIT */ #ifdef GMP_LIMB_HIGHBIT -#define MPFR_LIMB_HIGHBIT GMP_LIMB_HIGHBIT +# define MPFR_LIMB_HIGHBIT GMP_LIMB_HIGHBIT #else -#define MPFR_LIMB_HIGHBIT MP_LIMB_T_HIGHBIT +# ifdef MP_LIMB_T_HIGHBIT +# define MPFR_LIMB_HIGHBIT MP_LIMB_T_HIGHBIT +# else +# error "Neither GMP_LIMB_HIGHBIT nor MP_LIMB_T_HIGHBIT defined in GMP" +# endif #endif #if GMP_NAIL_BITS != 0 @@ -251,6 +255,29 @@ long double __gmpfr_longdouble_volatile _MPFR_PROTO ((long double)) ATTRIBUTE_CO # endif #endif +/* We want to test this : + * (rnd == GMP_RNDU && test) || (rnd == RNDD && !test) + * This macro does this test faster*/ +#define MPFR_IS_RNDUTEST_OR_RNDDNOTTEST(rnd, test) \ + (((rnd) + (test)) == GMP_RNDD) + +/* Calcul s = (-a) % BITS_PER_MP_LIMB + * a is unsigned! Check if it works, + * otherwise tries another way to calcul it */ +#define MPFR_UNSIGNED_MINUS_MODULO(s, a) \ + do { \ + if ((UINT_MAX % BITS_PER_MP_LIMB) == (BITS_PER_MP_LIMB-1) \ + && ((-(unsigned) 1)%BITS_PER_MP_LIMB > 0)) \ + (s) = (-(a)) % BITS_PER_MP_LIMB; \ + else \ + {(s) = (a) % BITS_PER_MP_LIMB; \ + if (s) \ + (s) = BITS_PER_MP_LIMB - (s); \ + } \ + MPFR_ASSERTD( (s) >= 0 && (s) < BITS_PER_MP_LIMB); \ + } while (0) + + /* Definition of the special values of the exponent */ /* * Clear flags macros are still defined and should be still used @@ -311,11 +338,11 @@ long double __gmpfr_longdouble_volatile _MPFR_PROTO ((long double)) ATTRIBUTE_CO #define MPFR_CHANGE_SIGN(x) (MPFR_SIGN(x) = -MPFR_SIGN(x)) #define MPFR_SET_SAME_SIGN(x, y) (MPFR_SIGN(x) = MPFR_SIGN(y)) -#define MPFR_SET_OPPOSITE_SIGN(x, y) (MPFR_SIGN(x) = - MPFR_SIGN(y)) -#define MPFR_CHECK_SIGN(s) \ +#define MPFR_SET_OPPOSITE_SIGN(x, y) (MPFR_SIGN(x) = -MPFR_SIGN(y)) +#define MPFR_ASSERT_SIGN(s) \ (MPFR_ASSERTD((s) == MPFR_SIGN_POS || (s) == MPFR_SIGN_NEG)) #define MPFR_SET_SIGN(x, s) \ - (MPFR_CHECK_SIGN(s), MPFR_SIGN(x) = s) + (MPFR_ASSERT_SIGN(s), MPFR_SIGN(x) = s) #define MPFR_IS_POS_SIGN(s1) (s1 > 0) #define MPFR_IS_NEG_SIGN(s1) (s1 < 0) #define MPFR_MULT_SIGN(s1, s2) ((s1) * (s2)) @@ -391,9 +418,6 @@ void mpfr_save_emin_emax _MPFR_PROTO ((void)); void mpfr_restore_emin_emax _MPFR_PROTO ((void)); int mpfr_add1 _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t)); int mpfr_sub1 _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t)); -int mpfr_round_raw_generic _MPFR_PROTO ((mp_limb_t *, mp_limb_t *, - mp_prec_t, int, - mp_prec_t, mp_rnd_t, int *, int)); int mpfr_can_round_raw _MPFR_PROTO ((mp_limb_t *, mp_size_t, int, mp_exp_t, mp_rnd_t, mp_rnd_t, mp_prec_t)); double mpfr_get_d3 _MPFR_PROTO ((mpfr_srcptr, mp_exp_t, mp_rnd_t)); @@ -413,12 +437,17 @@ long mpn_exp _MPFR_PROTO ((mp_limb_t *, mp_exp_t *, int, void mpfr_print_binary _MPFR_PROTO ((mpfr_srcptr)); void mpfr_set_str_binary _MPFR_PROTO ((mpfr_ptr, __gmp_const char *)); -#define mpfr_round_raw(yp, xp, xprec, neg, yprec, r, inexp) \ - mpfr_round_raw_generic((yp), (xp), (xprec), (neg), (yprec), (r), (inexp), 0) +int mpfr_round_raw _MPFR_PROTO ((mp_limb_t *, mp_limb_t *, + mp_prec_t, int, mp_prec_t, mp_rnd_t, int *)); +int mpfr_round_raw_2 _MPFR_PROTO ((mp_limb_t *, mp_limb_t *, + mp_prec_t, int, mp_prec_t, mp_rnd_t)); +int mpfr_round_raw_3 _MPFR_PROTO ((mp_limb_t *, mp_limb_t *, + mp_prec_t, int, mp_prec_t, mp_rnd_t, int *)); +int mpfr_round_raw_4 _MPFR_PROTO ((mp_limb_t *, mp_limb_t *, + mp_prec_t, int, mp_prec_t, mp_rnd_t)); #define mpfr_round_raw2(xp, xn, neg, r, prec) \ - mpfr_round_raw_generic(0, (xp), (xn) * BITS_PER_MP_LIMB, (neg), \ - (prec), (r), 0, 1); + mpfr_round_raw_2(0, (xp), (xn) * BITS_PER_MP_LIMB, (neg), (prec), (r) ) #if defined (__cplusplus) } @@ -24,12 +24,12 @@ MA 02111-1307, USA. */ /* Check if GMP is included, and try to include it (Works with local GMP) */ #ifndef __GMP_H__ -#include <gmp.h> +# include <gmp.h> #endif /* Check if stdio.h is included */ #if defined (EOF) -#define _MPFR_H_HAVE_FILE 1 +# define _MPFR_H_HAVE_FILE 1 #endif /* Definition of rounding modes */ @@ -79,9 +79,14 @@ typedef struct { mp_limb_t *_mpfr_d; } __mpfr_struct; -/* Compatibility with previous versions of MPFR */ +/* Compatibility with previous types of MPFR */ #define mp_rnd_t mpfr_rnd_t #define mp_prec_t mpfr_prec_t +#define __gmp_default_rounding_mode __gmpfr_default_rounding_mode +#define __mpfr_emin __gmpfr_emin +#define __mpfr_emax __gmpfr_emax +#define __mpfr_flags __gmpfr_flags +#define __mpfr_default_fp_bit_precision __gmpfr_default_fp_bit_precision /* The represented number is @@ -103,13 +108,13 @@ typedef __gmp_const __mpfr_struct *mpfr_srcptr; #define MPFR_SIGN(x) (((x)->_mpfr_sign)) +/* size_t is defined by GMP */ + /* Prototypes */ -#ifndef _MPFR_PROTO -# if defined (__STDC__) || defined (__cplusplus) -# define _MPFR_PROTO(x) x -# else -# define _MPFR_PROTO(x) () -# endif +#if defined (__STDC__) || defined (__cplusplus) +# define _MPFR_PROTO(x) x +#else +# define _MPFR_PROTO(x) () #endif #if defined (__cplusplus) @@ -119,6 +124,7 @@ extern "C" { extern unsigned int __gmpfr_flags; extern mp_exp_t __gmpfr_emin; extern mp_exp_t __gmpfr_emax; + mp_exp_t mpfr_get_emin _MPFR_PROTO ((void)); int mpfr_set_emin _MPFR_PROTO ((mp_exp_t)); mp_exp_t mpfr_get_emax _MPFR_PROTO ((void)); @@ -137,7 +143,6 @@ int mpfr_inexflag_p _MPFR_PROTO ((void)); void mpfr_init2 _MPFR_PROTO ((mpfr_ptr, mp_prec_t)); void mpfr_init _MPFR_PROTO ((mpfr_ptr)); int mpfr_prec_round _MPFR_PROTO ((mpfr_ptr, mp_prec_t, mpfr_rnd_t)); -#define mpfr_round_prec(x,r,p) mpfr_prec_round(x,p,r) /* compatibility 2.0.1 */ int mpfr_can_round _MPFR_PROTO ((mpfr_ptr, mp_exp_t, mpfr_rnd_t, mpfr_rnd_t, mp_prec_t)); mp_exp_t mpfr_get_exp _MPFR_PROTO ((mpfr_srcptr)); @@ -166,8 +171,10 @@ void mpfr_nextabove _MPFR_PROTO ((mpfr_ptr)); void mpfr_nextbelow _MPFR_PROTO ((mpfr_ptr)); void mpfr_nexttoward _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr)); int mpfr_set_str _MPFR_PROTO ((mpfr_ptr, __gmp_const char *, int, mpfr_rnd_t)); -int mpfr_init_set_str _MPFR_PROTO ((mpfr_ptr, __gmp_const char *, int, mpfr_rnd_t)); -char* mpfr_get_str _MPFR_PROTO ((char *, mp_exp_t *, int, size_t, mpfr_srcptr, mpfr_rnd_t)); +int mpfr_init_set_str _MPFR_PROTO ((mpfr_ptr, __gmp_const char *, int, + mpfr_rnd_t)); +char* mpfr_get_str _MPFR_PROTO ((char *, mp_exp_t *, int, size_t, mpfr_srcptr, + mpfr_rnd_t)); #ifdef _MPFR_H_HAVE_FILE #define mpfr_inp_str mpfr_inp_str_internal #define mpfr_out_str mpfr_out_str_internal @@ -176,8 +183,8 @@ size_t mpfr_out_str _MPFR_PROTO ((FILE *, int, size_t, mpfr_srcptr, mpfr_rnd_t)) #endif int mpfr_mul _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t)); int mpfr_pow_ui _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr, unsigned long int, mpfr_rnd_t)); - int mpfr_ui_pow_ui _MPFR_PROTO ((mpfr_ptr, unsigned long int, unsigned long int, - mpfr_rnd_t)); +int mpfr_ui_pow_ui _MPFR_PROTO ((mpfr_ptr, unsigned long int, + unsigned long int, mpfr_rnd_t)); int mpfr_div _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t)); int mpfr_agm _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t)); int mpfr_sqrt _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_rnd_t)); @@ -241,7 +248,6 @@ int mpfr_set4 _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_rnd_t, int)); int mpfr_cmp3 _MPFR_PROTO ((mpfr_srcptr, mpfr_srcptr, int)); int mpfr_cmp_d _MPFR_PROTO ((mpfr_srcptr, double)); int mpfr_cmpabs _MPFR_PROTO ((mpfr_srcptr, mpfr_srcptr)); -#define mpfr_cmp_abs mpfr_cmpabs /* keep for compatibility with mpfr-2.0.1 */ int mpfr_nan_p _MPFR_PROTO((mpfr_srcptr)); int mpfr_inf_p _MPFR_PROTO((mpfr_srcptr)); int mpfr_number_p _MPFR_PROTO((mpfr_srcptr)); @@ -307,6 +313,10 @@ int mpfr_sgn _MPFR_PROTO ((mpfr_srcptr)); } #endif +/* Compatibility with 2.0.1 */ +#define mpfr_cmp_abs mpfr_cmpabs +#define mpfr_round_prec(x,r,p) mpfr_prec_round(x,p,r) + /* prevent from using mpfr_get_e{min,max} as lvalues */ #define mpfr_get_emin() (__gmpfr_emin + 0) #define mpfr_get_emax() (__gmpfr_emax + 0) @@ -96,12 +96,13 @@ mpfr_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) * It is usefull iff the precision is big, there is an overflow * and we are doing further mults... Probable ? */ /* - if (bx + cx > __gmpfr_emax + 1) + if (MPFR_UNLIKELY(bx + cx > __gmpfr_emax + 1)) return mpfr_set_overflow (a, rnd_mode, sign_product); - if (bx + cx < __gmpfr_emin - 2) + if (MPFR_UNLIKELY(bx + cx < __gmpfr_emin - 2)) return mpfr_set_underflow (a, rnd_mode == GMP_RNDN ? GMP_RNDZ : rnd_mode, - sign_product); + sign_product); */ + ap = MPFR_MANT(a); bp = MPFR_MANT(b); cp = MPFR_MANT(c); @@ -109,8 +110,8 @@ mpfr_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) aq = MPFR_PREC(a); bq = MPFR_PREC(b); cq = MPFR_PREC(c); - if (MPFR_UNLIKELY(bq+cq < bq)) - MPFR_ASSERTN(0); /* no integer overflow */ + + MPFR_ASSERTD(bq+cq > bq); /* PREC_MAX is /2 so no integer overflow */ an = (aq-1)/BITS_PER_MP_LIMB + 1; /* number of significant limbs of a */ bn = (bq-1)/BITS_PER_MP_LIMB + 1; /* number of significant limbs of b */ diff --git a/round_prec.c b/round_prec.c index d938ffb84..514629143 100644 --- a/round_prec.c +++ b/round_prec.c @@ -25,129 +25,25 @@ MA 02111-1307, USA. */ #include "mpfr.h" #include "mpfr-impl.h" -/* - * If flag = 0, puts in y the value of xp (with precision xprec and - * sign 1 if negative=0, -1 otherwise) rounded to precision yprec and - * direction rnd_mode. Supposes x is not zero nor NaN nor +/- Infinity - * (i.e. *xp != 0). - * - * If inexp != NULL, computes the inexact flag of the rounding. - * (In case of even rounding when rnd = GMP_RNDN, puts 2 or -2 in *inexp.) - * - * If flag = 1, just returns whether one should add 1 or not for rounding. - * - * Note: yprec may be < MPFR_PREC_MIN; in particular, it may be equal - * to 1. In this case, the even rounding is done away from 0, which is - * a natural generalization. Indeed, a number with 1-bit precision can - * be seen as a denormalized number with more precision. - */ - -int -mpfr_round_raw_generic(mp_limb_t *yp, mp_limb_t *xp, mp_prec_t xprec, - int neg, mp_prec_t yprec, mp_rnd_t rnd_mode, - int *inexp, int flag) -{ - mp_size_t xsize, nw; - mp_limb_t himask, lomask; - int rw, carry = 0; - - xsize = (xprec-1)/BITS_PER_MP_LIMB + 1; - - nw = yprec / BITS_PER_MP_LIMB; - rw = yprec & (BITS_PER_MP_LIMB - 1); - - if (flag && !inexp && (rnd_mode==GMP_RNDZ || xprec <= yprec - || (rnd_mode==GMP_RNDU && neg) - || (rnd_mode==GMP_RNDD && neg==0))) - return 0; - - if (rw) - { - nw++; - lomask = ((MP_LIMB_T_ONE << (BITS_PER_MP_LIMB - rw)) - MP_LIMB_T_ONE); - himask = ~lomask; - } - else - { - lomask = -1; - himask = -1; - } - MPFR_ASSERTD(nw >= 1); - - if (xprec <= yprec) - { /* No rounding is necessary. */ - /* if yp=xp, maybe an overlap: MPN_COPY_DECR is ok when src <= dst */ - MPFR_ASSERTD(nw >= xsize); - if (inexp) - *inexp = 0; - if (flag) - return 0; - MPN_COPY_DECR(yp + (nw - xsize), xp, xsize); - MPN_ZERO(yp, nw - xsize); - } - else - { - mp_limb_t sb; - - if ((rnd_mode == GMP_RNDU && neg) || - (rnd_mode == GMP_RNDD && !neg)) - rnd_mode = GMP_RNDZ; - - if (inexp || rnd_mode != GMP_RNDZ) - { - mp_size_t k; - - k = xsize - nw; - if (!rw) - k--; - MPFR_ASSERTD(k >= 0); - sb = xp[k] & lomask; /* First non-significant bits */ - if (rnd_mode == GMP_RNDN) - { - mp_limb_t rbmask = MP_LIMB_T_ONE << (BITS_PER_MP_LIMB - rw - 1); - if (sb & rbmask) /* rounding bit */ - sb &= ~rbmask; /* it is 1, clear it */ - else - rnd_mode = GMP_RNDZ; /* it is 0, behave like rounding to 0 */ - } - while (sb == 0 && k > 0) - sb = xp[--k]; - if (rnd_mode == GMP_RNDN) - { /* rounding to nearest, with rounding bit = 1 */ - if (sb == 0) /* Even rounding. */ - { - sb = xp[xsize - nw] & (himask ^ (himask << 1)); - if (inexp) - *inexp = ((neg != 0) ^ (sb != 0)) - ? MPFR_EVEN_INEX : -MPFR_EVEN_INEX; - } - else /* sb != 0 */ - { - if (inexp) - *inexp = (neg == 0) ? 1 : -1; - } - } - else if (inexp) - *inexp = sb == 0 ? 0 - : (((neg != 0) ^ (rnd_mode != GMP_RNDZ)) ? 1 : -1); - } - else - sb = 0; - - if (flag) - return sb != 0 && rnd_mode != GMP_RNDZ; - - if (sb != 0 && rnd_mode != GMP_RNDZ) - carry = mpn_add_1(yp, xp + xsize - nw, nw, - rw ? MP_LIMB_T_ONE << (BITS_PER_MP_LIMB - rw) : 1); - else - MPN_COPY_INCR(yp, xp + xsize - nw, nw); - - yp[0] &= himask; - } - - return carry; -} +#define mpfr_round_raw_generic mpfr_round_raw +#define flag 0 +#define use_inexp 1 +#include "round_raw_generic.c" + +#define mpfr_round_raw_generic mpfr_round_raw_2 +#define flag 1 +#define use_inexp 0 +#include "round_raw_generic.c" + +#define mpfr_round_raw_generic mpfr_round_raw_3 +#define flag 1 +#define use_inexp 1 +#include "round_raw_generic.c" + +#define mpfr_round_raw_generic mpfr_round_raw_4 +#define flag 0 +#define use_inexp 0 +#include "round_raw_generic.c" int mpfr_prec_round (mpfr_ptr x, mp_prec_t prec, mp_rnd_t rnd_mode) @@ -168,7 +64,7 @@ mpfr_prec_round (mpfr_ptr x, mp_prec_t prec, mp_rnd_t rnd_mode) /* Realloc mantissa */ mp_ptr tmp = (mp_ptr) (*__gmp_reallocate_func) (MPFR_GET_REAL_PTR(x), MPFR_ALLOC_SIZE(ow), MPFR_ALLOC_SIZE(nw)); - MPFR_SET_MANT_PTR(x, tmp); /* mant pre must be set before alloc size */ + MPFR_SET_MANT_PTR(x, tmp); /* mant ptr must be set before alloc size */ MPFR_SET_ALLOC_SIZE(x, nw); /* new number of allocated limbs */ } @@ -252,6 +148,7 @@ mpfr_can_round_raw (mp_limb_t *bp, mp_size_t bn, int neg, mp_exp_t err0, /* can round only in rounding to the nearest and err0 >= prec + 2 */ } + MPFR_ASSERT_SIGN(neg); neg = MPFR_IS_NEG_SIGN(neg); /* if the error is smaller than ulp(b), then anyway it will propagate @@ -261,15 +158,12 @@ mpfr_can_round_raw (mp_limb_t *bp, mp_size_t bn, int neg, mp_exp_t err0, /* warning: if k = m*BITS_PER_MP_LIMB, consider limb m-1 and not m */ k = (err - 1) / BITS_PER_MP_LIMB; - s = err % BITS_PER_MP_LIMB; - if (s) - s = BITS_PER_MP_LIMB - s; + MPFR_UNSIGNED_MINUS_MODULO(s, err); + /* the error corresponds to bit s in limb k, the most significant limb being limb 0 */ k1 = (prec - 1) / BITS_PER_MP_LIMB; - s1 = prec % BITS_PER_MP_LIMB; - if (s1) - s1 = BITS_PER_MP_LIMB - s1; + MPFR_UNSIGNED_MINUS_MODULO(s1, prec); /* the last significant bit is bit s1 in limb k1 */ diff --git a/round_raw_generic.c b/round_raw_generic.c new file mode 100644 index 000000000..a8b4cc758 --- /dev/null +++ b/round_raw_generic.c @@ -0,0 +1,168 @@ +/* mpfr_round_raw_generic -- Generic roundin function + +Copyright 1999, 2000, 2001, 2002, 2003 Free Software Foundation, Inc. + +This file is part of the MPFR Library. + +The 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 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 MPFR Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#ifndef flag +# error "ERROR: flag must be defined (0 / 1)" +#endif +#ifndef use_inexp +# error "ERROR: use_enexp must be defined (0 / 1)" +#endif +#ifndef mpfr_round_raw_generic +# error "ERROR: mpfr_round_raw_generic must be defined" +#endif + +/* + * If flag = 0, puts in y the value of xp (with precision xprec and + * sign 1 if negative=0, -1 otherwise) rounded to precision yprec and + * direction rnd_mode. Supposes x is not zero nor NaN nor +/- Infinity + * (i.e. *xp != 0). + * + * If inexp != NULL, computes the inexact flag of the rounding. + * (In case of even rounding when rnd = GMP_RNDN, puts 2 or -2 in *inexp.) + * + * If flag = 1, just returns whether one should add 1 or not for rounding. + * + * Note: yprec may be < MPFR_PREC_MIN; in particular, it may be equal + * to 1. In this case, the even rounding is done away from 0, which is + * a natural generalization. Indeed, a number with 1-bit precision can + * be seen as a denormalized number with more precision. + */ + +int +mpfr_round_raw_generic(mp_limb_t *yp, mp_limb_t *xp, mp_prec_t xprec, + int neg, mp_prec_t yprec, mp_rnd_t rnd_mode +#if use_inexp != 0 + , int *inexp +#endif + ) +{ + mp_size_t xsize, nw; + mp_limb_t himask, lomask; + mp_limb_t sb; + int rw, carry; +#if use_inexp == 0 + int *inexp; +#endif + + if (use_inexp) + MPFR_ASSERTD(inexp != ((int*) 0)); + MPFR_ASSERTD(neg == 0 || neg == 1); + + if (flag && !use_inexp && + (rnd_mode==GMP_RNDZ || xprec <= yprec + ||MPFR_IS_RNDUTEST_OR_RNDDNOTTEST(rnd_mode, neg))) + return 0; + + xsize = (xprec-1)/BITS_PER_MP_LIMB + 1; + nw = yprec / BITS_PER_MP_LIMB; + rw = yprec & (BITS_PER_MP_LIMB - 1); + + if (MPFR_UNLIKELY(xprec <= yprec)) + { /* No rounding is necessary. */ + /* if yp=xp, maybe an overlap: MPN_COPY_DECR is ok when src <= dst */ + if (MPFR_LIKELY(rw)) + nw++; + MPFR_ASSERTD(nw >= 1); + MPFR_ASSERTD(nw >= xsize); + if (use_inexp) + *inexp = 0; + if (!flag) + { + MPN_COPY_DECR(yp + (nw - xsize), xp, xsize); + MPN_ZERO(yp, nw - xsize); + } + return 0; + } + + /* Rounding is necessary */ + if (MPFR_LIKELY(rw)) + { + nw++; + lomask = ((MP_LIMB_T_ONE << (BITS_PER_MP_LIMB - rw)) - MP_LIMB_T_ONE); + himask = ~lomask; + } + else + { + lomask = -1; + himask = -1; + } + + if (MPFR_IS_RNDUTEST_OR_RNDDNOTTEST(rnd_mode, neg)) + rnd_mode = GMP_RNDZ; + + if (use_inexp || rnd_mode != GMP_RNDZ) + { + mp_size_t k; + + k = xsize - nw; + if (MPFR_UNLIKELY(!rw)) + k--; + MPFR_ASSERTD(k >= 0); + sb = xp[k] & lomask; /* First non-significant bits */ + if (rnd_mode == GMP_RNDN) + { + mp_limb_t rbmask = MP_LIMB_T_ONE << (BITS_PER_MP_LIMB - rw - 1); + if (sb & rbmask) /* rounding bit */ + sb &= ~rbmask; /* it is 1, clear it */ + else + rnd_mode = GMP_RNDZ; /* it is 0, behave like rounding to 0 */ + } + while (sb == 0 && k > 0) + sb = xp[--k]; + if (rnd_mode == GMP_RNDN) + { /* rounding to nearest, with rounding bit = 1 */ + if (sb == 0) /* Even rounding. */ + { + sb = xp[xsize - nw] & (himask ^ (himask << 1)); + if (use_inexp) + *inexp = ((neg != 0) ^ (sb != 0)) + ? MPFR_EVEN_INEX : -MPFR_EVEN_INEX; + } + else /* sb != 0 */ + if (use_inexp) + *inexp = (neg == 0) ? 1 : -1; + } + else if (use_inexp) + *inexp = sb == 0 ? 0 + : (((neg != 0) ^ (rnd_mode != GMP_RNDZ)) ? 1 : -1); + } + else + sb = 0; + + if (flag) + return sb != 0 && rnd_mode != GMP_RNDZ; + + if (sb != 0 && rnd_mode != GMP_RNDZ) + carry = mpn_add_1(yp, xp + xsize - nw, nw, + rw ? MP_LIMB_T_ONE << (BITS_PER_MP_LIMB - rw) : 1); + else + { + MPN_COPY_INCR(yp, xp + xsize - nw, nw); + carry = 0; + } + yp[0] &= himask; + + return carry; +} + +#undef flag +#undef use_inexp +#undef mpfr_round_raw_generic @@ -351,8 +351,8 @@ mpfr_set_str (mpfr_t x, const char *str, int base, mp_rnd_t rnd) /* round y */ - if (mpfr_round_raw_generic (MPFR_MANT(x), result + n, n * BITS_PER_MP_LIMB, - negative, MPFR_PREC(x), rnd, NULL, (int) 0)) + if (mpfr_round_raw_4 (MPFR_MANT(x), result + n, n * BITS_PER_MP_LIMB, + negative, MPFR_PREC(x), rnd )) { /* overflaw when rounding y */ MPFR_MANT(x)[MPFR_LIMB_SIZE(x) - 1] @@ -74,29 +74,21 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) * So sign(a) = - sign(b) */ - MPFR_SET_SAME_SIGN(a,b); if (sign < 0) /* swap b and c so that |b| > |c| */ { mpfr_srcptr t; - MPFR_CHANGE_SIGN(a); + MPFR_SET_OPPOSITE_SIGN(a,b); t = b; b = c; c = t; } + else + MPFR_SET_SAME_SIGN(a,b); diff_exp = (mp_exp_unsigned_t) MPFR_GET_EXP (b) - MPFR_GET_EXP (c); /* reserve a space to store b aligned with the result, i.e. shifted by (-cancel) % BITS_PER_MP_LIMB to the right */ bn = MPFR_LIMB_SIZE(b); - if ((UINT_MAX % BITS_PER_MP_LIMB) == (BITS_PER_MP_LIMB-1)) - shift_b = (-cancel) % BITS_PER_MP_LIMB; - else - { - shift_b = cancel % BITS_PER_MP_LIMB; - if (shift_b) - shift_b = BITS_PER_MP_LIMB - shift_b; - } - MPFR_ASSERTD( shift_b >= 0 && shift_b < BITS_PER_MP_LIMB); - + MPFR_UNSIGNED_MINUS_MODULO(shift_b, cancel); cancel1 = (cancel + shift_b) / BITS_PER_MP_LIMB; /* the high cancel1 limbs from b should not be taken into account */ @@ -117,9 +109,10 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) } /* reserve a space to store c aligned with the result, i.e. shifted by - (diff_exp-cancel) % BITS_PER_MP_LIMB to the right */ + (diff_exp-cancel) % BITS_PER_MP_LIMB to the right */ cn = MPFR_LIMB_SIZE(c); - if ((UINT_MAX % BITS_PER_MP_LIMB) == (BITS_PER_MP_LIMB-1)) + if ((UINT_MAX % BITS_PER_MP_LIMB) == (BITS_PER_MP_LIMB-1) + && ((-(unsigned) 1)%BITS_PER_MP_LIMB > 0)) shift_c = (diff_exp - cancel) % BITS_PER_MP_LIMB; else { @@ -265,8 +258,7 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) } else /* directed rounding: set rnd_mode to RNDZ iff towards zero */ { - if (((rnd_mode == GMP_RNDD) && (MPFR_IS_POS(a))) || - ((rnd_mode == GMP_RNDU) && (MPFR_IS_NEG(a)))) + if (MPFR_IS_RNDUTEST_OR_RNDDNOTTEST(rnd_mode, MPFR_IS_NEG(a))) rnd_mode = GMP_RNDZ; if (carry) diff --git a/tests/Makefile.am b/tests/Makefile.am index 9edf17db1..83645de60 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -1,6 +1,7 @@ AUTOMAKE_OPTIONS = gnu -check_PROGRAMS = tinits reuse tabs tadd tagm tcan_round tcmp tcmp2 tcmp_ui tdiv tdiv_ui teq texp tfrac tget_str tlog tconst_log2 tmul tmul_2exp tmul_ui tout_str tconst_pi tpow trandom tround_prec tset_d tset_f tset_q tset_si tset_str tset_z tsqrt tsqrt_ui tui_div tui_sub tswap ttrunc trint tisnan tget_d tget_d_2exp tatan tcosh tsinh ttanh tacosh tasinh tatanh thyperbolic texp2 tfactorial tsub tasin tconst_euler tcos tsin ttan tsub_ui tset tlog1p texpm1 tlog2 tlog10 tui_pow tpow3 tadd_ui texceptions tfma thypot tacos tgamma tset_ld tcbrt tsin_cos mpf_compat mpfr_compat tzeta tcmp_d terf tnext tcomparisons tsgn +check_PROGRAMS = tinits tsgn tisnan texceptions tset tabs tset_d tset_f tset_q tset_si tset_str tset_z tset_ld tadd tsub tmul tdiv tcmp tcmp2 tcan_round tround_prec tswap reuse tsqrt tnext tcomparisons teq tadd_ui tsub_ui tcmp_ui tdiv_ui tmul_ui tsqrt_ui tui_div tui_sub tcmp_d tmul_2exp tfma tfrac tget_str tout_str tget_d tget_d_2exp tconst_log2 tconst_pi tconst_euler trandom ttrunc trint texp texp2 texpm1 tlog tlog2 tlog10 tlog1p tpow tui_pow tpow3 tcosh tsinh ttanh tacosh tasinh tatanh thyperbolic tasin tacos tcos tatan tsin ttan tsin_cos tagm thypot tfactorial tgamma terf tcbrt mpf_compat mpfr_compat tzeta + EXTRA_DIST = tgeneric.c mpf_compat.h LDADD = libfrtests.a $(MPFR_LIBM) $(top_builddir)/libmpfr.a @LDADD@ |