diff options
author | pelissip <pelissip@280ebfd0-de03-0410-8827-d642c229c3f4> | 2003-11-05 16:18:16 +0000 |
---|---|---|
committer | pelissip <pelissip@280ebfd0-de03-0410-8827-d642c229c3f4> | 2003-11-05 16:18:16 +0000 |
commit | d406c7c3f07cacb170533149e9319c334062792b (patch) | |
tree | da17b87f0205cc16561cf4d73f95ffc893e0385d | |
parent | 3e4ced245d9e603cfc806e64ab56b9b2b7223287 (diff) | |
download | mpfr-d406c7c3f07cacb170533149e9319c334062792b.tar.gz |
Fix some bugs (Use MPFR_ASSERT(1) instead of MPFR_ASSERT(0))
Optimize swap.c and copysign.c.
git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@2536 280ebfd0-de03-0410-8827-d642c229c3f4
-rw-r--r-- | acos.c | 2 | ||||
-rw-r--r-- | acosh.c | 3 | ||||
-rw-r--r-- | add.c | 2 | ||||
-rw-r--r-- | add_one_ulp.c | 7 | ||||
-rw-r--r-- | agm.c | 5 | ||||
-rw-r--r-- | asin.c | 2 | ||||
-rw-r--r-- | asinh.c | 3 | ||||
-rw-r--r-- | atan.c | 3 | ||||
-rw-r--r-- | atanh.c | 2 | ||||
-rw-r--r-- | cbrt.c | 2 | ||||
-rw-r--r-- | cmp.c | 2 | ||||
-rw-r--r-- | cmp_abs.c | 2 | ||||
-rw-r--r-- | cmp_ui.c | 2 | ||||
-rw-r--r-- | copysign.c | 4 | ||||
-rw-r--r-- | cos.c | 2 | ||||
-rw-r--r-- | cosh.c | 3 | ||||
-rw-r--r-- | div.c | 2 | ||||
-rw-r--r-- | div_ui.c | 4 | ||||
-rw-r--r-- | erf.c | 7 | ||||
-rw-r--r-- | exp.c | 2 | ||||
-rw-r--r-- | exp2.c | 2 | ||||
-rw-r--r-- | expm1.c | 2 | ||||
-rw-r--r-- | fma.c | 2 | ||||
-rw-r--r-- | gamma.c | 2 | ||||
-rw-r--r-- | gammaPiAGMformula.c | 18 | ||||
-rw-r--r-- | get_z_exp.c | 6 | ||||
-rw-r--r-- | hypot.c | 2 | ||||
-rw-r--r-- | log.c | 2 | ||||
-rw-r--r-- | log10.c | 2 | ||||
-rw-r--r-- | log1p.c | 2 | ||||
-rw-r--r-- | log2.c | 2 | ||||
-rw-r--r-- | mpfr-impl.h | 70 | ||||
-rw-r--r-- | mpfr.h | 3 | ||||
-rw-r--r-- | mul.c | 6 | ||||
-rw-r--r-- | mul_ui.c | 2 | ||||
-rw-r--r-- | pow.c | 2 | ||||
-rw-r--r-- | pow_si.c | 4 | ||||
-rw-r--r-- | pow_ui.c | 5 | ||||
-rw-r--r-- | rint.c | 2 | ||||
-rw-r--r-- | set.c | 2 | ||||
-rw-r--r-- | sin.c | 2 | ||||
-rw-r--r-- | sin_cos.c | 2 | ||||
-rw-r--r-- | sinh.c | 3 | ||||
-rw-r--r-- | sub.c | 2 | ||||
-rw-r--r-- | sub1.c | 121 | ||||
-rw-r--r-- | swap.c | 30 | ||||
-rw-r--r-- | tan.c | 2 | ||||
-rw-r--r-- | tanh.c | 3 | ||||
-rw-r--r-- | ui_div.c | 4 | ||||
-rw-r--r-- | ui_sub.c | 8 | ||||
-rw-r--r-- | zeta.c | 7 |
51 files changed, 220 insertions, 163 deletions
@@ -55,7 +55,7 @@ mpfr_acos (mpfr_ptr acos, mpfr_srcptr x, mp_rnd_t rnd_mode) MPFR_SET_EXP (acos, MPFR_GET_EXP (acos) - 1); return 1; /* inexact */ } - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } MPFR_CLEAR_FLAGS(x); @@ -51,7 +51,8 @@ mpfr_acosh (mpfr_ptr y, mpfr_srcptr x , mp_rnd_t rnd_mode) MPFR_SET_POS(y); MPFR_RET(0); } - MPFR_ASSERTN(1); + else + MPFR_ASSERTN(0); } comp = mpfr_cmp_ui (x, 1); if (MPFR_UNLIKELY( comp < 0 )) @@ -73,7 +73,7 @@ mpfr_add (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) return mpfr_set(a, b, rnd_mode); } /* Should never reach here */ - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } MPFR_ASSERTD(MPFR_IS_PURE_FP(b) && MPFR_IS_PURE_FP(c)); diff --git a/add_one_ulp.c b/add_one_ulp.c index 2ffe731e9..b56248467 100644 --- a/add_one_ulp.c +++ b/add_one_ulp.c @@ -36,9 +36,10 @@ mpfr_add_one_ulp (mpfr_ptr x, mp_rnd_t rnd_mode) { if (MPFR_IS_NAN(x)) MPFR_RET_NAN; - if (MPFR_IS_INF(x) || MPFR_IS_ZERO(x)) - return 0; - MPFR_ASSERTN(1); + else if (MPFR_IS_INF(x) || MPFR_IS_ZERO(x)) + MPFR_RET(0); + else + MPFR_ASSERTN(0); } xn = 1 + (MPFR_PREC(x) - 1) / BITS_PER_MP_LIMB; @@ -1,6 +1,6 @@ /* mpfr_agm -- arithmetic-geometric mean of two floating-point numbers -Copyright 1999, 2000, 2001, 2002 Free Software Foundation. +Copyright 1999, 2000, 2001, 2002, 2003 Free Software Foundation. This file is part of the MPFR Library. @@ -63,7 +63,8 @@ mpfr_agm (mpfr_ptr r, mpfr_srcptr op2, mpfr_srcptr op1, mp_rnd_t rnd_mode) MPFR_SET_ZERO(r); MPFR_RET(0); /* exact */ } - MPFR_ASSERTN(1); + else + MPFR_ASSERTN(0); } MPFR_CLEAR_FLAGS(r); @@ -54,7 +54,7 @@ mpfr_asin (mpfr_ptr asin, mpfr_srcptr x, mp_rnd_t rnd_mode) mpfr_set_ui (asin, 0, GMP_RNDN); MPFR_RET(0); /* exact result */ } - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } MPFR_CLEAR_FLAGS(asin); @@ -56,7 +56,8 @@ mpfr_asinh (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) MPFR_SET_SAME_SIGN(y, x); MPFR_RET(0); } - MPFR_ASSERTN(1); + else + MPFR_ASSERTN(0); } MPFR_CLEAR_FLAGS(y); @@ -113,7 +113,8 @@ mpfr_atan (mpfr_ptr arctangent, mpfr_srcptr x, mp_rnd_t rnd_mode) mpfr_set_ui (arctangent, 0, GMP_RNDN); return 0; /* exact result */ } - MPFR_ASSERTN(1); + else + MPFR_ASSERTN(0); } MPFR_CLEAR_FLAGS(arctangent); @@ -57,7 +57,7 @@ mpfr_atanh (mpfr_ptr y, mpfr_srcptr xt , mp_rnd_t rnd_mode) MPFR_RET(0); } else - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } /* Useless due to final mpfr_set MPFR_CLEAR_FLAGS(y);*/ @@ -76,7 +76,7 @@ mpfr_cbrt (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) return 0; } else - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } /* Useless due to mpz_init MPFR_CLEAR_FLAGS(y);*/ @@ -55,7 +55,7 @@ mpfr_cmp3 (mpfr_srcptr b, mpfr_srcptr c, int s) return MPFR_SIGN(b); else /* Should never reach here */ - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } /* b and c are real numbers */ if (s != MPFR_SIGN(b)) @@ -45,7 +45,7 @@ mpfr_cmpabs (mpfr_srcptr b, mpfr_srcptr c) else if (MPFR_IS_ZERO (b)) return -1; else - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } be = MPFR_GET_EXP (b); @@ -42,7 +42,7 @@ mpfr_cmp_ui_2exp (mpfr_srcptr b, unsigned long int i, mp_exp_t f) else if (MPFR_IS_ZERO(b)) return i != 0 ? -1 : 0; else - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } if (MPFR_IS_NEG(b)) diff --git a/copysign.c b/copysign.c index 4de5e3314..f38b44f1f 100644 --- a/copysign.c +++ b/copysign.c @@ -1,6 +1,6 @@ /* mpfr_copysign -- Produce a value with the magnitude of x and sign of y -Copyright 2001, 2002 Free Software Foundation. +Copyright 2001, 2002, 2003 Free Software Foundation. This file is part of the MPFR Library. @@ -32,7 +32,7 @@ MA 02111-1307, USA. */ int mpfr_copysign (mpfr_ptr z, mpfr_srcptr x ,mpfr_srcptr y , mp_rnd_t rnd_mode) { - if (MPFR_IS_NAN(y)) + if (MFR_UNLIKELY( MPFR_IS_NAN(y))) { MPFR_SET_NAN(z); MPFR_RET_NAN; @@ -43,7 +43,7 @@ mpfr_cos (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) else if (MPFR_IS_ZERO(x)) return mpfr_set_ui (y, 1, GMP_RNDN); /* Never reach this */ - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } precy = MPFR_PREC(y); @@ -55,7 +55,8 @@ mpfr_cosh (mpfr_ptr y, mpfr_srcptr xt , mp_rnd_t rnd_mode) else if (MPFR_IS_ZERO(xt)) return mpfr_set_ui(y, 1, rnd_mode); /* cosh(0) = 1 */ /* Should never reach this code */ - MPFR_ASSERTN(1); + else + MPFR_ASSERTN(0); } mpfr_init2(x,Nxt); @@ -97,7 +97,7 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mp_rnd_t rnd_mode) MPFR_RET(0); } /* Never reach this !*/ - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } MPFR_CLEAR_FLAGS(q); @@ -63,7 +63,7 @@ mpfr_div_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mp_rnd_t rnd_mode) } } else - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } if (MPFR_UNLIKELY(u == 0)) @@ -188,6 +188,6 @@ mpfr_div_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mp_rnd_t rnd_mode) MPFR_RET(-MPFR_INT_SIGN(x)); } } - MPFR_ASSERTN(1); /* should never go here */ + MPFR_ASSERTN(0); /* should never go here */ return 0; /* To avoid warning*/ } @@ -52,11 +52,12 @@ mpfr_erf (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) MPFR_SET_NAN(y); MPFR_RET_NAN; } - if (MPFR_IS_INF(x)) /* erf(+inf) = +1, erf(-inf) = -1 */ + else if (MPFR_IS_INF(x)) /* erf(+inf) = +1, erf(-inf) = -1 */ return mpfr_set_si (y, MPFR_FROM_SIGN_TO_INT(sign_x), GMP_RNDN); - if (MPFR_IS_ZERO(x)) /* erf(+0) = +0, erf(-0) = -0 */ + else if (MPFR_IS_ZERO(x)) /* erf(+0) = +0, erf(-0) = -0 */ return mpfr_set (y, x, GMP_RNDN); /* should keep the sign of x */ - MPFR_ASSERTN(1); + else + MPFR_ASSERTN(0); } /* now x is neither NaN, Inf nor 0 */ @@ -56,6 +56,8 @@ mpfr_exp (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) } else if (MPFR_IS_ZERO(x)) return mpfr_set_ui (y, 1, GMP_RNDN); + else + MPFR_ASSERTN(0); } MPFR_CLEAR_FLAGS(y); @@ -56,7 +56,7 @@ mpfr_exp2 (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) else if (MPFR_IS_ZERO(x)) return mpfr_set_ui (y, 1, rnd_mode); else - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } /* Useless due to mpfr_set MPFR_CLEAR_FLAGS(y);*/ @@ -60,7 +60,7 @@ mpfr_expm1 (mpfr_ptr y, mpfr_srcptr x , mp_rnd_t rnd_mode) MPFR_RET(0); } else - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } /* Useless due to mpfr_set MPFR_CLEAR_FLAGS(y);*/ @@ -100,7 +100,7 @@ mpfr_fma (mpfr_ptr s, mpfr_srcptr x, mpfr_srcptr y, mpfr_srcptr z, else if (MPFR_IS_ZERO(z)) return mpfr_mul (s, x, y, rnd_mode); else - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } /* Useless since it is done by mpfr_add * MPFR_CLEAR_FLAGS(s); */ @@ -88,7 +88,7 @@ mpfr_gamma (mpfr_ptr gamma, mpfr_srcptr x, mp_rnd_t rnd_mode) return 0; /* exact */ } else - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } MPFR_CLEAR_FLAGS(gamma); diff --git a/gammaPiAGMformula.c b/gammaPiAGMformula.c index c1edad36a..9f64750c5 100644 --- a/gammaPiAGMformula.c +++ b/gammaPiAGMformula.c @@ -37,15 +37,7 @@ int mpfr_gamma _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr, mp_rnd_t)); #define zCST 0.26 /* zCST=1/(2*ln(2*pi)) */ -int -#if __STDC__ -mpfr_gamma (mpfr_ptr gamma, mpfr_srcptr x, mp_rnd_t rnd_mode) -#else -mpfr_gamma (gamma, x, rnd_mode) - mpfr_ptr gamma; - mpfr_srcptr x; - mp_rnd_t rnd_mode; -#endif +int mpfr_gamma (mpfr_ptr gamma, mpfr_srcptr x, mp_rnd_t rnd_mode) { mpfr_t xp; mpfr_t product; @@ -77,19 +69,21 @@ mpfr_gamma (gamma, x, rnd_mode) if (MPFR_IS_NAN(x)) { MPFR_SET_NAN(gamma); + /* FIXME: MPFR_RET_NAN ? */ return 1; } - if (MPFR_ISZERO(x)) + else if (MPFR_ISZERO(x)) { MPFR_SET_INF(gamma); return 1; } - if (MPFR_IS_INF(x)) + else if (MPFR_IS_INF(x)) { MPFR_SET_INF(gamma); return 1; } - MPFR_ASSERTN(1); + else + MPFR_ASSERTN(0); } /* Set x_p=x if x> 1 else set x_p=2-x */ diff --git a/get_z_exp.c b/get_z_exp.c index 763c5a9b0..80fc3df1b 100644 --- a/get_z_exp.c +++ b/get_z_exp.c @@ -51,7 +51,7 @@ mpfr_get_z_exp (mpz_ptr z, mpfr_srcptr f) return __gmpfr_emin; } - fn = 1 + (MPFR_PREC(f) - 1) / BITS_PER_MP_LIMB; + fn = MPFR_LIMB_SIZE(f); /* check whether allocated space for z is enough */ if (ALLOC(z) < fn) @@ -66,9 +66,9 @@ mpfr_get_z_exp (mpz_ptr z, mpfr_srcptr f) SIZ(z) = MPFR_SIGN(f) < 0 ? -fn : fn; /* Test if the result is representable. Later, we could choose - to return MP_EXP_T_MIN if it isn't, or perhaps MP_EXP_T_MAX + to return MPFR_EXP_MIN if it isn't, or perhaps MPFR_EXP_MAX to signal an error. The mantissa would still be meaningful. */ - MPFR_ASSERTN((mp_exp_unsigned_t) MPFR_GET_EXP (f) - MP_EXP_T_MIN + MPFR_ASSERTN((mp_exp_unsigned_t) MPFR_GET_EXP (f) - MPFR_EXP_MIN >= (mp_exp_unsigned_t) MPFR_PREC(f)); return MPFR_GET_EXP (f) - MPFR_PREC (f); @@ -62,7 +62,7 @@ mpfr_hypot (mpfr_ptr z, mpfr_srcptr x , mpfr_srcptr y , mp_rnd_t rnd_mode) else if (MPFR_IS_ZERO(y)) return mpfr_abs (z, x, rnd_mode); else - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } MPFR_CLEAR_FLAGS(z); @@ -83,7 +83,7 @@ mpfr_log (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) MPFR_RET(0); /* log(0) is an exact -infinity */ } else - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } /* If a is negative, the result is NaN */ @@ -67,7 +67,7 @@ mpfr_log10 (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) MPFR_RET(0); /* log10(0) is an exact -infinity */ } else - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } /* If a is negative, the result is NaN */ @@ -63,7 +63,7 @@ mpfr_log1p (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) MPFR_RET(0); } else - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } comp = mpfr_cmp_si(x,-1); @@ -66,7 +66,7 @@ mpfr_log2 (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) MPFR_RET(0); /* log2(0) is an exact -infinity */ } else - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } /* If a is negative, the result is NaN */ diff --git a/mpfr-impl.h b/mpfr-impl.h index 8165e6aa1..2fcc54703 100644 --- a/mpfr-impl.h +++ b/mpfr-impl.h @@ -57,39 +57,39 @@ MA 02111-1307, USA. */ #endif #if GMP_NAIL_BITS != 0 -#error "MPFR doesn't support nonzero values of GMP_NAIL_BITS" +# error "MPFR doesn't support nonzero values of GMP_NAIL_BITS" #endif #if (BITS_PER_MP_LIMB & (BITS_PER_MP_LIMB - 1)) -#error "BITS_PER_MP_LIMB must be a power of 2" +# error "BITS_PER_MP_LIMB must be a power of 2" #endif /* Test if X (positive) is a power of 2 */ #define IS_POW2(X) (((X) & ((X) - 1)) == 0) #define NOT_POW2(X) (((X) & ((X) - 1)) != 0) -/* Set Exponent absolute limits (ie they are invalid exponent values */ -#ifdef MPFR_EXP_FORMAT_INT +/* Defined limits and unsigned types of exponent */ +#ifdef MPFR_EXP_FORMAT_INT + typedef unsigned int mp_exp_unsigned_t; # define MPFR_EXP_MAX (INT_MAX) # define MPFR_EXP_MIN (INT_MIN) #else + typedef unsigned long int mp_exp_unsigned_t; # define MPFR_EXP_MAX (LONG_MAX) # define MPFR_EXP_MIN (LONG_MIN) #endif -/* This unsigned type must correspond to the signed one defined in gmp.h */ -#ifdef MPFR_EXP_FORMAT_INT -typedef unsigned int mp_exp_unsigned_t; -typedef unsigned int mp_size_unsigned_t; +/* Defined limits and unsigned types of size */ +#ifdef MPFR_SIZE_FORMAT_INT + typedef unsigned int mp_size_unsigned_t; +# define MPFR_SIZE_MAX (INT_MAX) #else -typedef unsigned long int mp_exp_unsigned_t; -typedef unsigned long int mp_size_unsigned_t; + typedef unsigned long int mp_size_unsigned_t; +# define MPFR_SIZE_MAX (LONG_MAX) #endif -#define MP_EXP_T_MAX MPFR_EXP_MAX -/* FIXME: Is this really portable? */ -/*#define MP_EXP_T_MAX ((mp_exp_t) ((~ (mp_exp_unsigned_t) 0) >> 1))*/ -#define MP_EXP_T_MIN (-MP_EXP_T_MAX-1) +/*#define MP_EXP_T_MAX ((mp_exp_t) ((~ (mp_exp_unsigned_t) 0) >> 1)) +#define MP_EXP_T_MIN (-MP_EXP_T_MAX-1)*/ #define MP_LIMB_T_ONE ((mp_limb_t) 1) @@ -181,16 +181,16 @@ typedef union ieee_double_extract Ieee_double_extract; /* we only require that LDBL_MANT_DIG is a bound on the mantissa length of the "long double" type */ #ifndef LDBL_MANT_DIG -#define LDBL_MANT_DIG 113 /* works also if long double == quad */ +# define LDBL_MANT_DIG 113 /* works also if long double == quad */ #endif /* Various i386 systems have been seen with incorrect LDBL constants in float.h (notes in set_ld.c), so force the value we know is right for IEEE extended. */ #if HAVE_LDOUBLE_IEEE_EXT_LITTLE -#define MPFR_LDBL_MANT_DIG 64 +# define MPFR_LDBL_MANT_DIG 64 #else -#define MPFR_LDBL_MANT_DIG LDBL_MANT_DIG +# define MPFR_LDBL_MANT_DIG LDBL_MANT_DIG #endif /* LONGDOUBLE_NAN_ACTION executes the code "action" if x is a NaN. */ @@ -200,7 +200,7 @@ typedef union ieee_double_extract Ieee_double_extract; happen only after other comparisons, not sure what's really going on. In any case we can pick apart the bytes to identify a NaN. */ #if HAVE_LDOUBLE_IEEE_QUAD_BIG -#define LONGDOUBLE_NAN_ACTION(x, action) \ +# define LONGDOUBLE_NAN_ACTION(x, action) \ do { \ union { \ long double ld; \ @@ -224,26 +224,26 @@ typedef union ieee_double_extract Ieee_double_extract; "volatile" here stops "cc" on mips64-sgi-irix6.5 from optimizing away x!=x. */ #ifndef LONGDOUBLE_NAN_ACTION -#define LONGDOUBLE_NAN_ACTION(x, action) \ +# define LONGDOUBLE_NAN_ACTION(x, action) \ do { \ volatile long double __x = LONGDOUBLE_VOLATILE (x); \ if ((x) != __x) \ { action; } \ } while (0) -#define WANT_LONGDOUBLE_VOLATILE 1 +# define WANT_LONGDOUBLE_VOLATILE 1 #endif /* If we don't have a proper "volatile" then volatile is #defined to empty, in this case call through an external function to stop the compiler optimizing anything. */ #if WANT_LONGDOUBLE_VOLATILE -#ifdef volatile +# ifdef volatile long double __gmpfr_longdouble_volatile _MPFR_PROTO ((long double)) ATTRIBUTE_CONST; -#define LONGDOUBLE_VOLATILE(x) (__gmpfr_longdouble_volatile (x)) -#define WANT_GMPFR_LONGDOUBLE_VOLATILE 1 -#else -#define LONGDOUBLE_VOLATILE(x) (x) -#endif +# define LONGDOUBLE_VOLATILE(x) (__gmpfr_longdouble_volatile (x)) +# define WANT_GMPFR_LONGDOUBLE_VOLATILE 1 +# else +# define LONGDOUBLE_VOLATILE(x) (x) +# endif #endif /* Definition of the special values of the exponent */ @@ -266,12 +266,17 @@ long double __gmpfr_longdouble_volatile _MPFR_PROTO ((long double)) ATTRIBUTE_CO /* Old ESIZE */ #define MPFR_LIMB_SIZE(x) ((MPFR_PREC((x)) - 1) / BITS_PER_MP_LIMB + 1) -#define MPFR_EXP_ZERO (MPFR_EXP_MIN+0xB) -#define MPFR_EXP_NAN (MPFR_EXP_MIN+0xD) -#define MPFR_EXP_INF (MPFR_EXP_MIN+0xE) +/* Enum special value of exponent. It is an enum for debug reason (gdb) + But ANSI C restricts enumerator values to range of `int'... +typedef enum { + MPFR_EXP_ZERO = MPFR_EXP_MIN+0xB, + MPFR_EXP_NAN = MPFR_EXP_MIN+0xD, + MPFR_EXP_INF = MPFR_EXP_MIN+0xE + } mpfr_special_exp_t; */ -#define MPFR_SIGN_POS (1) -#define MPFR_SIGN_NEG (-1) +#define MPFR_EXP_ZERO (MPFR_EXP_MIN+1) +#define MPFR_EXP_NAN (MPFR_EXP_MIN+2) +#define MPFR_EXP_INF (MPFR_EXP_MIN+3) #define MPFR_CLEAR_FLAGS(x) /*#define MPFR_CLEAR_NAN(x)*/ @@ -294,6 +299,9 @@ long double __gmpfr_longdouble_volatile _MPFR_PROTO ((long double)) ATTRIBUTE_CO /* TODO: Redo all the macros dealing with the signs */ +#define MPFR_SIGN_POS (1) +#define MPFR_SIGN_NEG (-1) + /* Without zeros */ #define MPFR_ISNONNEG(x) (MPFR_NOTZERO((x)) && MPFR_SIGN(x) > 0) #define MPFR_ISNEG(x) (MPFR_NOTZERO((x)) && MPFR_SIGN(x) < 0) @@ -50,6 +50,7 @@ typedef enum { /* Define precision, exponent, sign */ #if __GMP_MP_SIZE_T_INT == 1 #define MPFR_EXP_FORMAT_INT +#define MPFR_SIZE_FORMAT_INT typedef unsigned int mpfr_prec_t; typedef unsigned int mpfr_size_t; typedef int mpfr_exp_t; @@ -173,7 +174,6 @@ 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)); #ifdef _MPFR_H_HAVE_FILE - /* They are only accessible if you include stdio.h first */ #define mpfr_inp_str mpfr_inp_str_internal #define mpfr_out_str mpfr_out_str_internal size_t mpfr_inp_str _MPFR_PROTO ((mpfr_ptr, FILE *, int, mpfr_rnd_t)); @@ -345,7 +345,6 @@ int mpfr_sgn _MPFR_PROTO ((mpfr_srcptr)); #define mpfr_set(a,b,r) mpfr_set4(a,b,r,MPFR_SIGN(b)) #define mpfr_abs(a,b,r) mpfr_set4(a,b,r,1) #define mpfr_cmp(b, c) mpfr_cmp3(b, c, 1) -/*#define mpfr_sgn(x) mpfr_cmp_ui(x,0)*/ #define mpfr_mul_2exp(y,x,n,r) mpfr_mul_2ui((y),(x),(n),(r)) #define mpfr_div_2exp(y,x,n,r) mpfr_div_2ui((y),(x),(n),(r)) @@ -79,7 +79,7 @@ mpfr_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) MPFR_RET(0); /* 0 * 0 is exact */ } /* Should never reach here */ - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } MPFR_CLEAR_FLAGS(a); sign_product = MPFR_MULT_SIGN( MPFR_SIGN(b) , MPFR_SIGN(c) ); @@ -88,8 +88,8 @@ mpfr_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) cx = MPFR_GET_EXP (c); /* Note: the exponent of the exact result will be e = bx + cx + ec with ec in {-1,0,1} and the following assumes that e is representable. */ - MPFR_ASSERTN(MPFR_EMAX_MAX <= (MP_EXP_T_MAX >> 1)); - MPFR_ASSERTN(MPFR_EMIN_MIN >= -(MP_EXP_T_MAX >> 1)); + MPFR_ASSERTN(MPFR_EMAX_MAX <= (MPFR_EXP_MAX >> 1)); + MPFR_ASSERTN(MPFR_EMIN_MIN >= -(MPFR_EXP_MAX >> 1)); if (bx + cx > __gmpfr_emax + 1) return mpfr_set_overflow (a, rnd_mode, sign_product); if (bx + cx < __gmpfr_emin - 2) @@ -62,7 +62,7 @@ mpfr_mul_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mp_rnd_t rnd_mode) MPFR_RET(0); /* zero is exact */ } else - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } else if (MPFR_UNLIKELY(u <= 1)) { @@ -244,7 +244,7 @@ mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode) MPFR_RET(0); } /* Should never reach this code */ - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); return 0; } @@ -38,7 +38,6 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd_mode) return mpfr_pow_ui (y, x, n, rnd_mode); else { - MPFR_CLEAR_FLAGS(y); if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) )) { if (MPFR_IS_NAN(x)) @@ -65,8 +64,9 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd_mode) MPFR_RET(0); } else - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } + MPFR_CLEAR_FLAGS(y); /* detect exact powers: x^(-n) is exact iff x is a power of 2 */ if (mpfr_cmp_si_2exp (x, MPFR_SIGN(x), MPFR_EXP(x) - 1) == 0) @@ -64,7 +64,7 @@ mpfr_pow_ui (mpfr_ptr x, mpfr_srcptr y, unsigned long int n, mp_rnd_t rnd) MPFR_RET(0); } else - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } else if (MPFR_UNLIKELY( n <= 1)) { @@ -75,7 +75,8 @@ mpfr_pow_ui (mpfr_ptr x, mpfr_srcptr y, unsigned long int n, mp_rnd_t rnd) /* y^1 = y */ return mpfr_set(x, y, rnd); } - + /* MPFR_CLEAR_FLAGS useless due to mpfr_set */ + mpfr_save_emin_emax (); mpfr_init (res); @@ -54,7 +54,7 @@ mpfr_rint (mpfr_ptr r, mpfr_srcptr u, mp_rnd_t rnd_mode) MPFR_SET_ZERO(r); MPFR_RET(0); /* zero is exact */ } - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } MPFR_SET_SAME_SIGN(r, u); @@ -50,7 +50,7 @@ mpfr_set4 (mpfr_ptr a, mpfr_srcptr b, mp_rnd_t rnd_mode, int signb) } else /* Should never reach this code */ - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } else if (MPFR_LIKELY(MPFR_PREC(b) == MPFR_PREC(a))) { @@ -104,7 +104,7 @@ mpfr_sin (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) MPFR_SET_SAME_SIGN(y, x); MPFR_RET(0); } - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } precy = MPFR_PREC(y); @@ -47,7 +47,7 @@ mpfr_sin_cos (mpfr_ptr y, mpfr_ptr z, mpfr_srcptr x, mp_rnd_t rnd_mode) MPFR_RET(0); } else - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } /* MPFR_CLEAR_FLAGS is useless since we use mpfr_set to set y and z */ @@ -58,7 +58,8 @@ mpfr_sinh (mpfr_ptr y, mpfr_srcptr xt, mp_rnd_t rnd_mode) MPFR_RET(0); } /* Should never reach this */ - MPFR_ASSERTN(1); + else + MPFR_ASSERTN(0); } mpfr_init2 (x, Nxt); @@ -75,7 +75,7 @@ mpfr_sub (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) return mpfr_set (a, b, rnd_mode); } /* Should never reach here */ - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } MPFR_CLEAR_FLAGS(a); MPFR_ASSERTD(MPFR_IS_PURE_FP(b) && MPFR_IS_PURE_FP(c)); @@ -46,7 +46,7 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) TMP_MARK(marker); ap = MPFR_MANT(a); - an = 1 + (MPFR_PREC(a) - 1) / BITS_PER_MP_LIMB; + an = MPFR_LIMB_SIZE(a); sign = mpfr_cmp2 (b, c, &cancel); if (MPFR_UNLIKELY(sign == 0)) @@ -86,14 +86,26 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) /* reserve a space to store b aligned with the result, i.e. shifted by (-cancel) % BITS_PER_MP_LIMB to the right */ - bn = 1 + (MPFR_PREC(b) - 1) / BITS_PER_MP_LIMB; - shift_b = cancel % BITS_PER_MP_LIMB; - if (shift_b) - shift_b = BITS_PER_MP_LIMB - shift_b; - cancel1 = (cancel + shift_b) / BITS_PER_MP_LIMB; + bn = MPFR_LIMB_SIZE(b); + /* shift_b = cancel % BITS_PER_MP_LIMB; + if (shift_b) + shift_b = BITS_PER_MP_LIMB - shift_b; + I think (-cancel)%BITS_PER_MP_LIMB is ok. I add an assert anyway */ + shift_b = (-cancel) % BITS_PER_MP_LIMB; + cancel1 = (cancel + shift_b) / BITS_PER_MP_LIMB; + MPFR_ASSERTD( shift_b >= 0 && shift_b < BITS_PER_MP_LIMB); + /* the high cancel1 limbs from b should not be taken into account */ - if (shift_b == 0) - bp = MPFR_MANT(b); /* no need of an extra space */ + if (MPFR_UNLIKELY(shift_b == 0)) + { + bp = MPFR_MANT(b); /* no need of an extra space */ + /* Ensure ap != bp */ + if (ap == bp) + { + bp = (mp_ptr) TMP_ALLOC(bn * BYTES_PER_MP_LIMB); + MPN_COPY (bp, ap, bn); + } + } else { bp = TMP_ALLOC ((bn + 1) * BYTES_PER_MP_LIMB); @@ -102,12 +114,23 @@ 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 */ - cn = 1 + (MPFR_PREC(c) - 1) / BITS_PER_MP_LIMB; - shift_c = diff_exp - (cancel % BITS_PER_MP_LIMB); - shift_c = (shift_c + BITS_PER_MP_LIMB) % BITS_PER_MP_LIMB; - if (shift_c == 0) - cp = MPFR_MANT(c); - else + cn = MPFR_LIMB_SIZE(c); + /* shift_c = diff_exp - (cancel % BITS_PER_MP_LIMB); + shift_c = (shift_c + BITS_PER_MP_LIMB) % BITS_PER_MP_LIMB;*/ + shift_c = (diff_exp - cancel) % BITS_PER_MP_LIMB; + MPFR_ASSERTD( shift_c >= 0 && shift_c < BITS_PER_MP_LIMB); + + if (MPFR_UNLIKELY(shift_c == 0)) + { + cp = MPFR_MANT(c); + /* Ensure ap != cp */ + if (ap == cp) + { + cp = (mp_ptr) TMP_ALLOC (cn * BYTES_PER_MP_LIMB); + MPN_COPY(cp, ap, cn); + } + } + else { cp = TMP_ALLOC ((cn + 1) * BYTES_PER_MP_LIMB); cp[0] = mpn_rshift (cp + 1, MPFR_MANT(c), cn++, shift_c); @@ -117,22 +140,8 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) printf("shift_b=%u shift_c=%u\n", shift_b, shift_c); #endif - /* ensure ap != bp and ap != cp */ - if (ap == bp) - { - bp = (mp_ptr) TMP_ALLOC(bn * BYTES_PER_MP_LIMB); - MPN_COPY (bp, ap, bn); - /* ap == cp cannot occur since we would have b = c, and this case - has already been processed (in mpfr_add or mpfr_sub if there is - a special value, else earlier in this function: sign == 0). */ - } - else if (ap == cp) - { - cp = (mp_ptr) TMP_ALLOC (cn * BYTES_PER_MP_LIMB); - MPN_COPY(cp, ap, cn); - } - MPFR_ASSERTN (ap != cp); - MPFR_ASSERTN (bp != cp); + MPFR_ASSERTD (ap != cp); + MPFR_ASSERTD (bp != cp); /* here we have shift_c = (diff_exp - cancel) % BITS_PER_MP_LIMB, thus we want cancel2 = ceil((cancel - diff_exp) / BITS_PER_MP_LIMB) */ @@ -162,10 +171,12 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) */ /* copy high(b) into a */ - if (an + cancel1 <= bn) /* a: <----------------+-----------|----> - b: <-----------------------------------------> */ + if (MPFR_LIKELY(an + cancel1 <= bn)) + /* a: <----------------+-----------|----> + b: <-----------------------------------------> */ MPN_COPY (ap, bp + bn - (an + cancel1), an); - else /* a: <----------------+-----------|----> + else + /* a: <----------------+-----------|----> b: <-------------------------> */ if (cancel1 < bn) /* otherwise b does not overlap with a */ { @@ -186,11 +197,13 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) if (cancel2 >= 0) { - if (an + cancel2 <= cn) /* a: <-----------------------------> - c: <-----------------------------------------> */ + if (an + cancel2 <= cn) + /* a: <-----------------------------> + c: <-----------------------------------------> */ mpn_sub_n (ap, ap, cp + cn - (an + cancel2), an); - else /* a: <----------------------------> - c: <-------------------------> */ + else + /* a: <----------------------------> + c: <-------------------------> */ { ap2 = ap + an + cancel2 - cn; if (cn > cancel2) @@ -199,11 +212,13 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) } else /* cancel2 < 0 */ { - if (an + cancel2 <= cn) /* a: <-----------------------------> - c: <-----------------------------> */ - borrow = mpn_sub_n (ap, ap, cp + cn - (an + cancel2), an + cancel2); - else /* a: <----------------------------> - c: <----------------> */ + if (an + cancel2 <= cn) + /* a: <-----------------------------> + c: <-----------------------------> */ + borrow = mpn_sub_n (ap, ap, cp + cn - (an + cancel2), an + cancel2); + else + /* a: <----------------------------> + c: <----------------> */ { ap2 = ap + an + cancel2 - cn; borrow = mpn_sub_n (ap2, ap2, cp, cn); @@ -225,7 +240,7 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) if (rnd_mode == GMP_RNDN) { - if (sh) + if (MPFR_LIKELY(sh)) { is_exact = (carry == 0); /* can decide except when carry = 2^(sh-1) [middle] @@ -242,8 +257,8 @@ 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_SIGN(a) > 0)) || - ((rnd_mode == GMP_RNDU) && (MPFR_SIGN(a) < 0))) + if (((rnd_mode == GMP_RNDD) && (MPFR_IS_POS(a))) || + ((rnd_mode == GMP_RNDU) && (MPFR_IS_NEG(a)))) rnd_mode = GMP_RNDZ; if (carry) @@ -365,8 +380,10 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) /* even rounding rule */ if ((ap[0] >> sh) & 1) { - if (down) goto sub_one_ulp; - else goto add_one_ulp; + if (down) + goto sub_one_ulp; + else + goto add_one_ulp; } else inexact = (down) ? 1 : -1; @@ -399,13 +416,13 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) /* we have to set MPFR_EXP(a) to MPFR_EXP(b) - cancel + add_exp, taking care of underflows/overflows in that computation, and of the allowed exponent range */ - if (cancel) + if (MPFR_LIKELY(cancel)) { mp_exp_t exp_a; cancel -= add_exp; /* still valid as unsigned long */ exp_a = MPFR_GET_EXP (b) - cancel; - if (exp_a < __gmpfr_emin) + if (MPFR_UNLIKELY(exp_a < __gmpfr_emin)) { TMP_FREE(marker); if (rnd_mode == GMP_RNDN && @@ -424,7 +441,7 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) mp_exp_t exp_b; exp_b = MPFR_GET_EXP (b); - if (add_exp && exp_b == __gmpfr_emax) + if (MPFR_UNLIKELY(add_exp && exp_b == __gmpfr_emax)) { TMP_FREE(marker); return mpfr_set_overflow (a, rnd_mode, MPFR_SIGN(a)); @@ -436,6 +453,6 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) printf ("result is a="); mpfr_print_binary(a); putchar('\n'); #endif /* check that result is msb-normalized */ - MPFR_ASSERTN(ap[an-1] > ~ap[an-1]); - return inexact * MPFR_SIGN(a); + MPFR_ASSERTD(ap[an-1] > ~ap[an-1]); + return inexact * MPFR_INT_SIGN(a); } @@ -22,6 +22,7 @@ MA 02111-1307, USA. */ #include <string.h> #include "gmp.h" #include "mpfr.h" +#include "mpfr-impl.h" /* * We now use memcpy instead of copying the structure field by field. @@ -34,9 +35,34 @@ MA 02111-1307, USA. */ void mpfr_swap (mpfr_ptr u, mpfr_ptr v) { - mpfr_t tmp; + /*mpfr_t tmp; memcpy (tmp, u, sizeof (mpfr_t)); memcpy (u, v, sizeof (mpfr_t)); - memcpy (v, tmp, sizeof (mpfr_t)); + memcpy (v, tmp, sizeof (mpfr_t));*/ + + mpfr_prec_t p1, p2; + mpfr_sign_t s1, s2; + mpfr_exp_t e1, e2; + mp_limb_t *m1, *m2; + + p1 = MPFR_PREC(u); + p2 = MPFR_PREC(v); + MPFR_PREC(v) = p1; + MPFR_PREC(u) = p2; + + s1 = MPFR_SIGN(u); + s2 = MPFR_SIGN(v); + MPFR_SIGN(v) = s1; + MPFR_SIGN(u) = s2; + + e1 = MPFR_EXP(u); + e2 = MPFR_EXP(v); + MPFR_EXP(v) = e1; + MPFR_EXP(u) = e2; + + m1 = MPFR_MANT(u); + m2 = MPFR_MANT(v); + MPFR_MANT(v) = m1; + MPFR_MANT(u) = m2; } @@ -46,7 +46,7 @@ mpfr_tan (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) MPFR_RET(0); } /* Should never reach this point */ - MPFR_ASSERTN(1); + MPFR_ASSERTN(0); } precy = MPFR_PREC(y); @@ -61,7 +61,8 @@ mpfr_tanh (mpfr_ptr y, mpfr_srcptr xt , mp_rnd_t rnd_mode) MPFR_RET(0); } /* Should never reach this point */ - MPFR_ASSERTN(1); + else + MPFR_ASSERTN(0); } mpfr_init2(x,Nxt); @@ -62,8 +62,8 @@ mpfr_ui_div (mpfr_ptr y, unsigned long int u, mpfr_srcptr x, mp_rnd_t rnd_mode) MPFR_RET_NAN; } } - MPFR_ASSERTN(1); - return 0; /* To avoid warning */ + else + MPFR_ASSERTN(0); } else if (u) { @@ -39,19 +39,19 @@ mpfr_ui_sub (mpfr_ptr y, unsigned long int u, mpfr_srcptr x, mp_rnd_t rnd_mode) MPFR_SET_NAN(y); MPFR_RET_NAN; } - if (MPFR_IS_INF(x)) + else if (MPFR_IS_INF(x)) { /* u - Inf = -Inf and u - -Inf = +Inf */ MPFR_SET_INF(y); MPFR_SET_OPPOSITE_SIGN(y,x); MPFR_RET(0); /* +/-infinity is exact */ } - if (MPFR_IS_ZERO(x)) + else if (MPFR_IS_ZERO(x)) /* u - 0 = u */ return mpfr_set_ui(y, u, rnd_mode); /* Should never reach this code */ - MPFR_ASSERTN(1); - return 0; /* To avoid a warning. It isn't reached */ + else + MPFR_ASSERTN(0); } else if (u) { @@ -340,21 +340,22 @@ mpfr_zeta (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode) MPFR_SET_NAN (z); MPFR_RET_NAN; } - if (MPFR_IS_INF(s)) + else if (MPFR_IS_INF(s)) { if (MPFR_SIGN(s) > 0) return mpfr_set_ui (z, 1, GMP_RNDN); /* Zeta(+Inf) = 1 */ MPFR_SET_NAN (z); /* Zeta(-Inf) = NaN */ MPFR_RET_NAN; } - if (MPFR_IS_ZERO(s)) + else if (MPFR_IS_ZERO(s)) { mpfr_set_ui (z, 1, rnd_mode); mpfr_div_2ui (z, z, 1, rnd_mode); MPFR_CHANGE_SIGN(z); MPFR_RET(0); } - MPFR_ASSERTN(1); + else + MPFR_ASSERTN(0); } MPFR_CLEAR_FLAGS(z); |