diff options
| -rw-r--r-- | BUGS | 8 | ||||
| -rw-r--r-- | acos.c | 12 | ||||
| -rw-r--r-- | acosh.c | 48 | ||||
| -rw-r--r-- | add.c | 13 | ||||
| -rw-r--r-- | add1.c | 20 | ||||
| -rw-r--r-- | add_one_ulp.c | 6 | ||||
| -rw-r--r-- | add_ui.c | 4 | ||||
| -rw-r--r-- | asin.c | 12 | ||||
| -rw-r--r-- | asinh.c | 45 | ||||
| -rw-r--r-- | atan.c | 14 | ||||
| -rw-r--r-- | atanh.c | 53 | ||||
| -rw-r--r-- | cbrt.c | 4 | ||||
| -rw-r--r-- | cmp.c | 6 | ||||
| -rw-r--r-- | cmp2.c | 8 | ||||
| -rw-r--r-- | cmp_abs.c | 6 | ||||
| -rw-r--r-- | cmp_si.c | 4 | ||||
| -rw-r--r-- | cmp_ui.c | 4 | ||||
| -rw-r--r-- | const_log2.c | 4 | ||||
| -rw-r--r-- | const_pi.c | 230 | ||||
| -rw-r--r-- | cos.c | 15 | ||||
| -rw-r--r-- | div.c | 6 | ||||
| -rw-r--r-- | div_2si.c | 12 | ||||
| -rw-r--r-- | div_2ui.c | 27 | ||||
| -rw-r--r-- | div_ui.c | 8 | ||||
| -rw-r--r-- | eq.c | 4 | ||||
| -rw-r--r-- | exceptions.c | 4 | ||||
| -rw-r--r-- | exp.c | 4 | ||||
| -rw-r--r-- | exp2.c | 2 | ||||
| -rw-r--r-- | exp3.c | 6 | ||||
| -rw-r--r-- | exp_2.c | 125 | ||||
| -rw-r--r-- | expm1.c | 38 | ||||
| -rw-r--r-- | frac.c | 4 | ||||
| -rw-r--r-- | generic.c | 4 | ||||
| -rw-r--r-- | get_d.c | 6 | ||||
| -rw-r--r-- | get_exp.c | 4 | ||||
| -rw-r--r-- | get_ld.c | 6 | ||||
| -rw-r--r-- | get_si.c | 2 | ||||
| -rw-r--r-- | get_str.c | 12 | ||||
| -rw-r--r-- | get_ui.c | 2 | ||||
| -rw-r--r-- | get_z_exp.c | 4 | ||||
| -rw-r--r-- | hypot.c | 6 | ||||
| -rw-r--r-- | isinteger.c | 2 | ||||
| -rw-r--r-- | log.c | 8 | ||||
| -rw-r--r-- | log1p.c | 6 | ||||
| -rw-r--r-- | log2.c | 9 | ||||
| -rw-r--r-- | mpfr-impl.h | 28 | ||||
| -rw-r--r-- | mul.c | 8 | ||||
| -rw-r--r-- | mul_2si.c | 12 | ||||
| -rw-r--r-- | mul_2ui.c | 15 | ||||
| -rw-r--r-- | mul_ui.c | 7 | ||||
| -rw-r--r-- | next.c | 10 | ||||
| -rw-r--r-- | pow.c | 8 | ||||
| -rw-r--r-- | print_raw.c | 4 | ||||
| -rw-r--r-- | random.c | 4 | ||||
| -rw-r--r-- | random2.c | 4 | ||||
| -rw-r--r-- | rint.c | 10 | ||||
| -rw-r--r-- | round_prec.c | 6 | ||||
| -rw-r--r-- | set.c | 8 | ||||
| -rw-r--r-- | set_d.c | 6 | ||||
| -rw-r--r-- | set_exp.c | 4 | ||||
| -rw-r--r-- | set_f.c | 4 | ||||
| -rw-r--r-- | set_si.c | 9 | ||||
| -rw-r--r-- | set_str.c | 2 | ||||
| -rw-r--r-- | set_str_raw.c | 10 | ||||
| -rw-r--r-- | set_ui.c | 9 | ||||
| -rw-r--r-- | set_z.c | 4 | ||||
| -rw-r--r-- | setmax.c | 4 | ||||
| -rw-r--r-- | setmin.c | 4 | ||||
| -rw-r--r-- | sin.c | 17 | ||||
| -rw-r--r-- | sin_cos.c | 8 | ||||
| -rw-r--r-- | sinh.c | 8 | ||||
| -rw-r--r-- | sqrt.c | 12 | ||||
| -rw-r--r-- | sqrt_ui.c | 4 | ||||
| -rw-r--r-- | sub.c | 13 | ||||
| -rw-r--r-- | sub1.c | 15 | ||||
| -rw-r--r-- | sub_one_ulp.c | 6 | ||||
| -rw-r--r-- | sub_ui.c | 4 | ||||
| -rw-r--r-- | swap.c | 11 | ||||
| -rw-r--r-- | tan.c | 5 | ||||
| -rw-r--r-- | tanh.c | 6 | ||||
| -rw-r--r-- | ui_div.c | 5 | ||||
| -rw-r--r-- | ui_pow.c | 39 | ||||
| -rw-r--r-- | ui_sub.c | 4 | ||||
| -rw-r--r-- | urandomb.c | 4 | ||||
| -rw-r--r-- | zeta.c | 6 |
85 files changed, 628 insertions, 546 deletions
@@ -59,3 +59,11 @@ Potential bugs: * mpfr_hypot may fail for x very large, y very small and a very large target precision. + +* Possible bugs related to the exponents. Compile with -DMPFR_EXP_CHECK + and make check to see the potential problems (but in this case, don't + be surprised if nothing is working!). This doesn't mean that the + results may be incorrect though (a further analysis should be done + and the code should be updated to avoid the failed assertions). As a + consequence, MPFR_EXP_CHECK is not currently defined by default even + when WANT_ASSERT is defined. @@ -1,6 +1,6 @@ /* mpfr_acos -- arc-cosinus of a floating-point number -Copyright 2001, 2002 Free Software Foundation. +Copyright 2001, 2002, 2003 Free Software Foundation. This file is part of the MPFR Library, and was contributed by Mathieu Dutour. @@ -32,7 +32,7 @@ mpfr_acos (mpfr_ptr acos, mpfr_srcptr x, mp_rnd_t rnd_mode) mpfr_t xp; mpfr_t arcc; - int signe, suplement; + int signe, supplement; mpfr_t tmp; int Prec; @@ -79,7 +79,7 @@ mpfr_acos (mpfr_ptr acos, mpfr_srcptr x, mp_rnd_t rnd_mode) { mpfr_clear (xp); mpfr_const_pi (acos, rnd_mode); - MPFR_EXP(acos)--; + MPFR_SET_EXP (acos, MPFR_GET_EXP (acos) - 1); return 1; /* inexact */ } @@ -87,15 +87,15 @@ mpfr_acos (mpfr_ptr acos, mpfr_srcptr x, mp_rnd_t rnd_mode) mpfr_ui_sub (xp, 1, xp, GMP_RNDD); if (signe > 0) - suplement = 2 - 2 * MPFR_EXP(xp); + supplement = 2 - 2 * MPFR_GET_EXP (xp); else - suplement = 2 - MPFR_EXP(xp); + supplement = 2 - MPFR_GET_EXP(xp); realprec = prec_acos + 10; while (!good) { - Prec = realprec+suplement; + Prec = realprec + supplement; /* Initialisation */ mpfr_init2 (tmp, Prec); @@ -1,6 +1,6 @@ /* mpfr_acosh -- inverse hyperbolic cosine -Copyright 2001, 2002 Free Software Foundation. +Copyright 2001, 2002, 2003 Free Software Foundation. This file is part of the MPFR Library. @@ -85,28 +85,30 @@ mpfr_acosh (mpfr_ptr y, mpfr_srcptr x , mp_rnd_t rnd_mode) mpfr_save_emin_emax (); /* First computation of acosh */ - do { - - /* reactualisation of the precision */ - mpfr_set_prec (t, Nt); - mpfr_set_prec (te, Nt); - mpfr_set_prec (ti, Nt); - - /* compute acosh */ - mpfr_mul (te, x, x, GMP_RNDD); /* x^2 */ - mpfr_sub_ui (ti, te, 1, GMP_RNDD); /* x^2-1 */ - mpfr_sqrt (t, ti, GMP_RNDN); /* sqrt(x^2-1) */ - mpfr_add (t, t, x, GMP_RNDN); /* sqrt(x^2-1)+x */ - mpfr_log (t, t, GMP_RNDN); /* ln(sqrt(x^2-1)+x)*/ - - /* estimation of the error -- see algorithms.ps */ - err = Nt - (-1 + 2 * MAX(2 + MAX(2 - MPFR_EXP(t), - 1 + MPFR_EXP(te) - MPFR_EXP(ti) - MPFR_EXP(t)), 0)); - - /* actualisation of the precision */ - Nt += 10; - - } while ((err < 0) || !mpfr_can_round (t, err, GMP_RNDN, rnd_mode, Ny)); + do + { + /* reactualisation of the precision */ + mpfr_set_prec (t, Nt); + mpfr_set_prec (te, Nt); + mpfr_set_prec (ti, Nt); + + /* compute acosh */ + mpfr_mul (te, x, x, GMP_RNDD); /* x^2 */ + mpfr_sub_ui (ti, te, 1, GMP_RNDD); /* x^2-1 */ + mpfr_sqrt (t, ti, GMP_RNDN); /* sqrt(x^2-1) */ + mpfr_add (t, t, x, GMP_RNDN); /* sqrt(x^2-1)+x */ + mpfr_log (t, t, GMP_RNDN); /* ln(sqrt(x^2-1)+x)*/ + + /* error estimate -- see algorithms.ps */ + err = Nt - (-1 + 2 * MAX(2 + MAX(2 - MPFR_GET_EXP (t), + 1 + MPFR_GET_EXP (te) + - MPFR_GET_EXP (ti) + - MPFR_GET_EXP (t)), 0)); + + /* actualisation of the precision */ + Nt += 10; + } + while ((err < 0) || !mpfr_can_round (t, err, GMP_RNDN, rnd_mode, Ny)); inexact = mpfr_set (y, t, rnd_mode); @@ -1,6 +1,6 @@ /* mpfr_add -- add two floating-point numbers -Copyright 1999, 2000, 2001, 2002 Free Software Foundation. +Copyright 1999, 2000, 2001, 2002, 2003 Free Software Foundation. Contributed by the Spaces project, INRIA Lorraine. This file is part of the MPFR Library. @@ -91,15 +91,16 @@ mpfr_add (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) } else { /* signs are equal, it's an addition */ - if (MPFR_EXP(b) < MPFR_EXP(c)) + mp_exp_t eb, ec; + eb = MPFR_GET_EXP (b); + ec = MPFR_GET_EXP (c); + if (eb < ec) { - return mpfr_add1(a, c, b, rnd_mode, - (mp_exp_unsigned_t) MPFR_EXP(c) - MPFR_EXP(b)); + return mpfr_add1(a, c, b, rnd_mode, (mp_exp_unsigned_t) ec - eb); } else { - return mpfr_add1(a, b, c, rnd_mode, - (mp_exp_unsigned_t) MPFR_EXP(b) - MPFR_EXP(c)); + return mpfr_add1(a, b, c, rnd_mode, (mp_exp_unsigned_t) eb - ec); } } } @@ -1,6 +1,6 @@ /* mpfr_add1 -- internal function to perform a "real" addition -Copyright 1999, 2000, 2001, 2002 Free Software Foundation. +Copyright 1999, 2000, 2001, 2002, 2003 Free Software Foundation. Contributed by the Spaces project, INRIA Lorraine. This file is part of the MPFR Library. @@ -71,7 +71,7 @@ mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, bn = (bq-1)/BITS_PER_MP_LIMB + 1; /* number of significant limbs of b */ cn = (cq-1)/BITS_PER_MP_LIMB + 1; /* number of significant limbs of c */ - MPFR_EXP(a) = MPFR_EXP(b); + MPFR_SET_EXP (a, MPFR_GET_EXP (b)); MPFR_SET_SAME_SIGN(a, b); /* @@ -159,13 +159,13 @@ mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, if (cc) /* carry */ { - mp_exp_t exp = MPFR_EXP(a); + mp_exp_t exp = MPFR_GET_EXP (a); if (exp == __gmpfr_emax) { inex = mpfr_set_overflow(a, rnd_mode, MPFR_SIGN(a)); goto end_of_add; } - MPFR_EXP(a)++; + MPFR_SET_EXP (a, exp + 1); rb = (ap[0] >> sh) & 1; /* LSB(a) --> rounding bit after the shift */ if (sh) { @@ -301,13 +301,13 @@ mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, && (rb < 0 || (rb ^= 1) == 0) && mpn_add_1(ap, ap, an, MP_LIMB_T_ONE << sh)) { - mp_exp_t exp = MPFR_EXP(a); + mp_exp_t exp = MPFR_GET_EXP (a); if (exp == __gmpfr_emax) { inex = mpfr_set_overflow(a, rnd_mode, MPFR_SIGN(a)); goto end_of_add; } - MPFR_EXP(a)++; + MPFR_SET_EXP (a, exp + 1); ap[an-1] = MPFR_LIMB_HIGHBIT; rb = 0; } @@ -355,13 +355,13 @@ mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, rb ^= 1; if (rb == 0 && mpn_add_1(ap, ap, an, MP_LIMB_T_ONE << sh)) { - mp_exp_t exp = MPFR_EXP(a); + mp_exp_t exp = MPFR_GET_EXP (a); if (exp == __gmpfr_emax) { inex = mpfr_set_overflow(a, rnd_mode, MPFR_SIGN(a)); goto end_of_add; } - MPFR_EXP(a)++; + MPFR_SET_EXP (a, exp + 1); ap[an-1] = MPFR_LIMB_HIGHBIT; } } /* bb < cc */ @@ -532,12 +532,12 @@ mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, add_one_ulp: /* add one unit in last place to a */ if (mpn_add_1(ap, ap, an, MP_LIMB_T_ONE << sh)) { - mp_exp_t exp = MPFR_EXP(a); + mp_exp_t exp = MPFR_GET_EXP (a); if (exp == __gmpfr_emax) inex = mpfr_set_overflow(a, rnd_mode, MPFR_SIGN(a)); else { - MPFR_EXP(a)++; + MPFR_SET_EXP (a, exp + 1); ap[an-1] = MPFR_LIMB_HIGHBIT; } } diff --git a/add_one_ulp.c b/add_one_ulp.c index b831c254f..f1199a535 100644 --- a/add_one_ulp.c +++ b/add_one_ulp.c @@ -1,6 +1,6 @@ /* mpfr_add_one_ulp -- add one unit in last place -Copyright 1999, 2001, 2002 Free Software Foundation, Inc. +Copyright 1999, 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -43,12 +43,12 @@ mpfr_add_one_ulp (mpfr_ptr x, mp_rnd_t rnd_mode) xp = MPFR_MANT(x); if (mpn_add_1 (xp, xp, xn, MP_LIMB_T_ONE << sh)) /* got 1.0000... */ { - mp_exp_t exp = MPFR_EXP(x); + mp_exp_t exp = MPFR_GET_EXP (x); if (exp == __gmpfr_emax) return mpfr_set_overflow(x, rnd_mode, MPFR_SIGN(x)); else { - MPFR_EXP(x)++; + MPFR_SET_EXP (x, exp + 1); xp[xn-1] = MPFR_LIMB_HIGHBIT; } } @@ -1,6 +1,6 @@ /* mpfr_add_ui -- add a floating-point number with a machine integer -Copyright 2000, 2001, 2002 Free Software Foundation. +Copyright 2000, 2001, 2002, 2003 Free Software Foundation. This file is part of the MPFR Library. @@ -39,7 +39,7 @@ mpfr_add_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mp_rnd_t rnd_mode) MPFR_ASSERTN(u == (mp_limb_t) u); count_leading_zeros(cnt, (mp_limb_t) u); *up = (mp_limb_t) u << cnt; - MPFR_EXP(uu) = BITS_PER_MP_LIMB - cnt; + MPFR_SET_EXP (uu, BITS_PER_MP_LIMB - cnt); /* Optimization note: Exponent operations may be removed if mpfr_add works even when uu is out-of-range. */ @@ -1,6 +1,6 @@ /* mpfr_asin -- arc-sinus of a floating-point number -Copyright 2001, 2002 Free Software Foundation. +Copyright 2001, 2002, 2003 Free Software Foundation. This file is part of the MPFR Library, and was contributed by Mathieu Dutour. @@ -32,7 +32,7 @@ mpfr_asin (mpfr_ptr asin, mpfr_srcptr x, mp_rnd_t rnd_mode) mpfr_t xp; mpfr_t arcs; - int signe, suplement; + int signe, supplement; mpfr_t tmp; int Prec; @@ -78,7 +78,7 @@ mpfr_asin (mpfr_ptr asin, mpfr_srcptr x, mp_rnd_t rnd_mode) mpfr_const_pi (asin, rnd_mode); mpfr_neg (asin, asin, rnd_mode); } - MPFR_EXP(asin)--; + MPFR_SET_EXP (asin, MPFR_GET_EXP (asin) - 1); mpfr_clear (xp); return 1; /* inexact */ } @@ -93,15 +93,15 @@ mpfr_asin (mpfr_ptr asin, mpfr_srcptr x, mp_rnd_t rnd_mode) prec_asin = MPFR_PREC(asin); mpfr_ui_sub (xp, 1, xp, GMP_RNDD); - suplement = 2 - MPFR_EXP(xp); + supplement = 2 - MPFR_GET_EXP (xp); #ifdef DEBUG - printf("suplement=%d\n", suplement); + printf("supplement=%d\n", supplement); #endif realprec = prec_asin + 10; while (!good) { - estimated_delta = 1 + suplement; + estimated_delta = 1 + supplement; Prec = realprec+estimated_delta; /* Initialisation */ @@ -1,6 +1,6 @@ /* mpfr_asinh -- inverse hyperbolic sine -Copyright 2001, 2002 Free Software Foundation. +Copyright 2001, 2002, 2003 Free Software Foundation. This file is part of the MPFR Library. @@ -82,27 +82,28 @@ mpfr_asinh (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) mpfr_save_emin_emax (); /* First computation of asinh */ - do { - - /* reactualisation of the precision */ - mpfr_set_prec (t, Nt); - mpfr_set_prec (te, Nt); - mpfr_set_prec (ti, Nt); - - /* compute asinh */ - mpfr_mul (te, x, x, GMP_RNDD); /* x^2 */ - mpfr_add_ui (ti, te, 1, GMP_RNDD); /* x^2+1 */ - mpfr_sqrt (t, ti, GMP_RNDN); /* sqrt(x^2+1) */ - ((neg) ? mpfr_sub : mpfr_add) (t, t, x, GMP_RNDN); /* sqrt(x^2+1)+x */ - mpfr_log (t, t, GMP_RNDN); /* ln(sqrt(x^2+1)+x)*/ - - /* estimation of the error -- see algorithms.ps */ - err = Nt - (MAX(3 - MPFR_EXP(t), 0) + 1); - - /* actualisation of the precision */ - Nt += 10; - - } while ((err < 0) || (!mpfr_can_round (t, err, GMP_RNDN, rnd_mode, Ny) || (MPFR_IS_ZERO(t)))); + do + { + /* reactualisation of the precision */ + mpfr_set_prec (t, Nt); + mpfr_set_prec (te, Nt); + mpfr_set_prec (ti, Nt); + + /* compute asinh */ + mpfr_mul (te, x, x, GMP_RNDD); /* x^2 */ + mpfr_add_ui (ti, te, 1, GMP_RNDD); /* x^2+1 */ + mpfr_sqrt (t, ti, GMP_RNDN); /* sqrt(x^2+1) */ + (neg ? mpfr_sub : mpfr_add) (t, t, x, GMP_RNDN); /* sqrt(x^2+1)+x */ + mpfr_log (t, t, GMP_RNDN); /* ln(sqrt(x^2+1)+x)*/ + + /* error estimate -- see algorithms.ps */ + err = Nt - (MAX(3 - MPFR_GET_EXP (t), 0) + 1); + + /* actualisation of the precision */ + Nt += 10; + } + while ((err < 0) || (!mpfr_can_round (t, err, GMP_RNDN, rnd_mode, Ny) + || MPFR_IS_ZERO(t))); mpfr_restore_emin_emax (); @@ -1,6 +1,6 @@ /* mpfr_atan -- arc-tangent of a floating-point number -Copyright 2001, 2002 Free Software Foundation. +Copyright 2001, 2002, 2003 Free Software Foundation. This file is part of the MPFR Library, and was contributed by Mathieu Dutour. @@ -57,7 +57,7 @@ mpfr_atan (mpfr_ptr arctangent, mpfr_srcptr x, mp_rnd_t rnd_mode) mpfr_t xp; mpfr_t arctgt; - int comparaison, signe, suplement; + int comparaison, signe, supplement; mpfr_t t_arctan; int i; @@ -103,7 +103,7 @@ mpfr_atan (mpfr_ptr arctangent, mpfr_srcptr x, mp_rnd_t rnd_mode) rnd_mode = GMP_RNDU; mpfr_const_pi (arctangent, rnd_mode); } - MPFR_EXP(arctangent)--; + MPFR_SET_EXP (arctangent, MPFR_GET_EXP (arctangent) - 1); return 1; /* inexact */ } @@ -136,9 +136,9 @@ mpfr_atan (mpfr_ptr arctangent, mpfr_srcptr x, mp_rnd_t rnd_mode) return 1; /* inexact result */ } if (comparaison > 0) - suplement = 2; + supplement = 2; else - suplement = 2-MPFR_EXP(xp); + supplement = 2 - MPFR_GET_EXP (xp); prec_x = __gmpfr_ceil_log2 ((double) MPFR_PREC(x) / BITS_PER_MP_LIMB); logn = __gmpfr_ceil_log2 ((double) prec_x); @@ -149,8 +149,8 @@ mpfr_atan (mpfr_ptr arctangent, mpfr_srcptr x, mp_rnd_t rnd_mode) while (!good){ - N0 = __gmpfr_ceil_log2((double) realprec + suplement + CST); - estimated_delta = 1 + suplement + __gmpfr_ceil_log2((double) (3*N0-2)); + N0 = __gmpfr_ceil_log2((double) realprec + supplement + CST); + estimated_delta = 1 + supplement + __gmpfr_ceil_log2((double) (3*N0-2)); Prec = realprec+estimated_delta; /* Initialisation */ @@ -1,6 +1,6 @@ /* mpfr_atanh -- Inverse Hyperbolic Tangente of Unsigned Integer Number -Copyright 2001, 2002 Free Software Foundation. +Copyright 2001, 2002, 2003 Free Software Foundation. This file is part of the MPFR Library. @@ -100,28 +100,30 @@ mpfr_atanh (mpfr_ptr y, mpfr_srcptr xt , mp_rnd_t rnd_mode) mpfr_init(ti); /* First computation of cosh */ - do { - - /* reactualisation of the precision */ - mpfr_set_prec(t,Nt); - mpfr_set_prec(te,Nt); - mpfr_set_prec(ti,Nt); - - /* compute atanh */ - mpfr_ui_sub(te,1,x,GMP_RNDU); /* (1-xt)*/ - mpfr_add_ui(ti,x,1,GMP_RNDD); /* (xt+1)*/ - mpfr_div(te,ti,te,GMP_RNDN); /* (1+xt)/(1-xt)*/ - mpfr_log(te,te,GMP_RNDN); /* ln((1+xt)/(1-xt))*/ - mpfr_div_2ui(t,te,1,GMP_RNDN); /* (1/2)*ln((1+xt)/(1-xt))*/ - - /* estimation of the error see- algorithms.ps*/ - /* err=Nt-__gmpfr_ceil_log2(1+5*pow(2,1-MPFR_EXP(t)));*/ - err=Nt-(MAX(4-MPFR_EXP(t),0)+1); - - /* actualisation of the precision */ - Nt += 10; - - } while ((err < 0) || (!mpfr_can_round(t,err,GMP_RNDN,rnd_mode,Ny) || (MPFR_IS_ZERO(t)))); + do + { + /* reactualisation of the precision */ + mpfr_set_prec(t,Nt); + mpfr_set_prec(te,Nt); + mpfr_set_prec(ti,Nt); + + /* compute atanh */ + mpfr_ui_sub(te,1,x,GMP_RNDU); /* (1-xt)*/ + mpfr_add_ui(ti,x,1,GMP_RNDD); /* (xt+1)*/ + mpfr_div(te,ti,te,GMP_RNDN); /* (1+xt)/(1-xt)*/ + mpfr_log(te,te,GMP_RNDN); /* ln((1+xt)/(1-xt))*/ + mpfr_div_2ui(t,te,1,GMP_RNDN); /* (1/2)*ln((1+xt)/(1-xt))*/ + + /* error estimate see- algorithms.ps*/ + /* err=Nt-__gmpfr_ceil_log2(1+5*pow(2,1-MPFR_EXP(t)));*/ + err = Nt - (MAX (4 - MPFR_GET_EXP (t), 0) + 1); + + /* actualisation of the precision */ + Nt += 10; + + } + while ((err < 0) || + (!mpfr_can_round(t,err,GMP_RNDN,rnd_mode,Ny) || MPFR_IS_ZERO(t))); if(flag_neg) MPFR_CHANGE_SIGN(t); @@ -135,8 +137,3 @@ mpfr_atanh (mpfr_ptr y, mpfr_srcptr xt , mp_rnd_t rnd_mode) mpfr_clear(x); return inexact; } - - - - - @@ -1,6 +1,6 @@ /* mpfr_cbrt -- cube root function. -Copyright 2002 Free Software Foundation. +Copyright 2002, 2003 Free Software Foundation. Contributed by the Spaces project, INRIA Lorraine. This file is part of the MPFR Library. @@ -134,7 +134,7 @@ mpfr_cbrt (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) is not changed; or inexact=0, and inexact is set only when rnd_mode=GMP_RNDN and bit (n+1) from m is 1 */ inexact += mpfr_set_z (y, m, GMP_RNDN); - MPFR_EXP(y) += e / 3; + MPFR_SET_EXP (y, MPFR_GET_EXP (y) + e / 3); if (sign_x < 0) { @@ -1,6 +1,6 @@ /* mpfr_cmp -- compare two floating-point numbers -Copyright 1999, 2001 Free Software Foundation. +Copyright 1999, 2001, 2003 Free Software Foundation. This file is part of the MPFR Library. @@ -59,8 +59,8 @@ mpfr_cmp3 (mpfr_srcptr b, mpfr_srcptr c, int s) /* now signs are equal */ - be = MPFR_EXP(b); - ce = MPFR_EXP(c); + be = MPFR_GET_EXP (b); + ce = MPFR_GET_EXP (c); if (be > ce) return s; if (be < ce) @@ -1,6 +1,6 @@ /* mpfr_cmp2 -- exponent shift when subtracting two 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. @@ -65,10 +65,10 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mp_prec_t *cancel) return 1; } - if (MPFR_EXP(b) >= MPFR_EXP(c)) + if (MPFR_GET_EXP (b) >= MPFR_GET_EXP (c)) { sign = 1; - diff_exp = (mp_exp_unsigned_t) MPFR_EXP(b) - MPFR_EXP(c); + diff_exp = (mp_exp_unsigned_t) MPFR_GET_EXP (b) - MPFR_GET_EXP (c); bp = MPFR_MANT(b); cp = MPFR_MANT(c); @@ -131,7 +131,7 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mp_prec_t *cancel) else /* MPFR_EXP(b) < MPFR_EXP(c) */ { sign = -1; - diff_exp = (mp_exp_unsigned_t) MPFR_EXP(c) - MPFR_EXP(b); + diff_exp = (mp_exp_unsigned_t) MPFR_GET_EXP (c) - MPFR_GET_EXP (b); bp = MPFR_MANT(c); cp = MPFR_MANT(b); @@ -1,6 +1,6 @@ /* mpfr_cmpabs -- compare the absolute values of two FP numbers -Copyright 1999, 2001, 2002 Free Software Foundation, Inc. +Copyright 1999, 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -42,8 +42,8 @@ mpfr_cmpabs (mpfr_srcptr b, mpfr_srcptr c) if (MPFR_IS_ZERO(c)) return 1; - be = MPFR_EXP(b); - ce = MPFR_EXP(c); + be = MPFR_GET_EXP (b); + ce = MPFR_GET_EXP (c); if (be > ce) return 1; if (be < ce) @@ -1,7 +1,7 @@ /* mpfr_cmp_si_2exp -- compare a floating-point number with a signed machine integer multiplied by a power of 2 -Copyright 1999, 2001, 2002 Free Software Foundation, Inc. +Copyright 1999, 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -60,7 +60,7 @@ mpfr_cmp_si_2exp (mpfr_srcptr b, long int i, mp_exp_t f) /* ai must be representable in a mp_limb_t */ MPFR_ASSERTN(ai == (mp_limb_t) ai); - e = MPFR_EXP(b); /* 2^(e-1) <= b < 2^e */ + e = MPFR_GET_EXP (b); /* 2^(e-1) <= b < 2^e */ if (e <= f) return -si; if (f < MPFR_EMAX_MAX - BITS_PER_MP_LIMB && @@ -1,7 +1,7 @@ /* mpfr_cmp_ui_2exp -- compare a floating-point number with an unsigned machine integer multiplied by a power of 2 -Copyright 1999, 2001, 2002 Free Software Foundation, Inc. +Copyright 1999, 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -58,7 +58,7 @@ mpfr_cmp_ui_2exp (mpfr_srcptr b, unsigned long int i, mp_exp_t f) /* i must be representable in a mp_limb_t */ MPFR_ASSERTN(i == (mp_limb_t) i); - e = MPFR_EXP(b); /* 2^(e-1) <= b < 2^e */ + e = MPFR_GET_EXP (b); /* 2^(e-1) <= b < 2^e */ if (e <= f) return -1; if (f < MPFR_EMAX_MAX - BITS_PER_MP_LIMB && diff --git a/const_log2.c b/const_log2.c index df7f29e59..c63ab717a 100644 --- a/const_log2.c +++ b/const_log2.c @@ -1,6 +1,6 @@ /* mpfr_const_log2 -- compute natural logarithm of 2 -Copyright 1999, 2001, 2002 Free Software Foundation, Inc. +Copyright 1999, 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -167,7 +167,7 @@ mpfr_const_log2 (mpfr_ptr x, mp_rnd_t rnd_mode) } mpfr_set_z (x, s, rnd_mode); - MPFR_EXP(x) -= N; + MPFR_SET_EXP (x, MPFR_GET_EXP (x) - N); mpz_clear (s); mpz_clear (t); mpz_clear (u); diff --git a/const_pi.c b/const_pi.c index f50581790..d57c201e6 100644 --- a/const_pi.c +++ b/const_pi.c @@ -1,6 +1,6 @@ /* mpfr_const_pi -- compute Pi -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. @@ -53,67 +53,70 @@ mpfr_pi_machin3 (mpfr_ptr mylog, mp_rnd_t rnd_mode) mpz_t cst; MPFR_CLEAR_FLAGS(mylog); - logn = __gmpfr_ceil_log2 ((double) MPFR_PREC(mylog)); + logn = __gmpfr_ceil_log2 ((double) MPFR_PREC(mylog)); prec_x = prec_i_want + logn + 5; mpz_init(cst); - while (!good){ - prec = __gmpfr_ceil_log2 ((double) prec_x); - - mpfr_init2(tmp1, prec_x); - mpfr_init2(tmp2, prec_x); - mpfr_init2(tmp3, prec_x); - mpfr_init2(tmp4, prec_x); - mpfr_init2(tmp5, prec_x); - mpfr_init2(tmp6, prec_x); - mpfr_init2(result, prec_x); - mpz_set_si(cst, -1); - - mpfr_aux_pi(tmp1, cst, 268*268, prec - 4); - mpfr_div_ui(tmp1, tmp1, 268, GMP_RNDD); - mpfr_mul_ui(tmp1, tmp1, 61, GMP_RNDD); - - mpfr_aux_pi(tmp2, cst, 343*343, prec - 4); - mpfr_div_ui(tmp2, tmp2, 343, GMP_RNDD); - mpfr_mul_ui(tmp2, tmp2, 122, GMP_RNDD); - - mpfr_aux_pi(tmp3, cst, 557*557, prec - 4); - mpfr_div_ui(tmp3, tmp3, 557, GMP_RNDD); - mpfr_mul_ui(tmp3, tmp3, 115, GMP_RNDD); - - mpfr_aux_pi(tmp4, cst, 1068*1068, prec - 4); - mpfr_div_ui(tmp4, tmp4, 1068, GMP_RNDD); - mpfr_mul_ui(tmp4, tmp4, 32, GMP_RNDD); - - mpfr_aux_pi(tmp5, cst, 3458*3458, prec - 4); - mpfr_div_ui(tmp5, tmp5, 3458, GMP_RNDD); - mpfr_mul_ui(tmp5, tmp5, 83, GMP_RNDD); - - mpfr_aux_pi(tmp6, cst, 27493*27493, prec - 4); - mpfr_div_ui(tmp6, tmp6, 27493, GMP_RNDD); - mpfr_mul_ui(tmp6, tmp6, 44, GMP_RNDD); - - mpfr_add(result, tmp1, tmp2, GMP_RNDD); - mpfr_add(result, result, tmp3, GMP_RNDD); - mpfr_sub(result, result, tmp4, GMP_RNDD); - mpfr_add(result, result, tmp5, GMP_RNDD); - mpfr_add(result, result, tmp6, GMP_RNDD); - mpfr_mul_2ui(result, result, 2, GMP_RNDD); - mpfr_clear(tmp1); - mpfr_clear(tmp2); - mpfr_clear(tmp3); - mpfr_clear(tmp4); - mpfr_clear(tmp5); - mpfr_clear(tmp6); - if (mpfr_can_round(result, prec_x - 5, GMP_RNDD, rnd_mode, prec_i_want)){ - mpfr_set(mylog, result, rnd_mode); - mpfr_clear(result); - good = 1; - } else + while (!good) { - mpfr_clear(result); - prec_x += logn; - } - } + prec = __gmpfr_ceil_log2 ((double) prec_x); + + mpfr_init2(tmp1, prec_x); + mpfr_init2(tmp2, prec_x); + mpfr_init2(tmp3, prec_x); + mpfr_init2(tmp4, prec_x); + mpfr_init2(tmp5, prec_x); + mpfr_init2(tmp6, prec_x); + mpfr_init2(result, prec_x); + mpz_set_si(cst, -1); + + mpfr_aux_pi(tmp1, cst, 268*268, prec - 4); + mpfr_div_ui(tmp1, tmp1, 268, GMP_RNDD); + mpfr_mul_ui(tmp1, tmp1, 61, GMP_RNDD); + + mpfr_aux_pi(tmp2, cst, 343*343, prec - 4); + mpfr_div_ui(tmp2, tmp2, 343, GMP_RNDD); + mpfr_mul_ui(tmp2, tmp2, 122, GMP_RNDD); + + mpfr_aux_pi(tmp3, cst, 557*557, prec - 4); + mpfr_div_ui(tmp3, tmp3, 557, GMP_RNDD); + mpfr_mul_ui(tmp3, tmp3, 115, GMP_RNDD); + + mpfr_aux_pi(tmp4, cst, 1068*1068, prec - 4); + mpfr_div_ui(tmp4, tmp4, 1068, GMP_RNDD); + mpfr_mul_ui(tmp4, tmp4, 32, GMP_RNDD); + + mpfr_aux_pi(tmp5, cst, 3458*3458, prec - 4); + mpfr_div_ui(tmp5, tmp5, 3458, GMP_RNDD); + mpfr_mul_ui(tmp5, tmp5, 83, GMP_RNDD); + + mpfr_aux_pi(tmp6, cst, 27493*27493, prec - 4); + mpfr_div_ui(tmp6, tmp6, 27493, GMP_RNDD); + mpfr_mul_ui(tmp6, tmp6, 44, GMP_RNDD); + + mpfr_add(result, tmp1, tmp2, GMP_RNDD); + mpfr_add(result, result, tmp3, GMP_RNDD); + mpfr_sub(result, result, tmp4, GMP_RNDD); + mpfr_add(result, result, tmp5, GMP_RNDD); + mpfr_add(result, result, tmp6, GMP_RNDD); + mpfr_mul_2ui(result, result, 2, GMP_RNDD); + mpfr_clear(tmp1); + mpfr_clear(tmp2); + mpfr_clear(tmp3); + mpfr_clear(tmp4); + mpfr_clear(tmp5); + mpfr_clear(tmp6); + if (mpfr_can_round(result, prec_x - 5, GMP_RNDD, rnd_mode, prec_i_want)) + { + mpfr_set(mylog, result, rnd_mode); + mpfr_clear(result); + good = 1; + } + else + { + mpfr_clear(result); + prec_x += logn; + } + } mpz_clear(cst); return 0; } @@ -168,52 +171,73 @@ mpfr_const_pi (mpfr_ptr x, mp_rnd_t rnd_mode) mpfr_set(x, __mpfr_const_pi, rnd_mode); return; } - if (prec < 20000){ - /* need to recompute */ - N=1; - do { - oldN = N; - N = (prec+3)/4 + __gmpfr_ceil_log2((double) N + 1.0); - } while (N != oldN); - mpz_init(pi); mpz_init(num); mpz_init(den); mpz_init(d3); mpz_init(d2); - mpz_init(tmp); - mpz_set_ui(pi, 0); - mpz_set_ui(num, 16); /* num(-1) */ - mpz_set_ui(den, 21); /* den(-1) */ - mpz_set_si(d3, -2454); - mpz_set_ui(d2, 14736); - /* invariants: num=120*n^2+151*n+47, den=512*n^4+1024*n^3+712*n^2+194*n+15 - d3 = 2048*n^3+400*n-6, d2 = 6144*n^2-6144*n+2448 - */ - for (n=0; n<N; n++) { - /* num(n)-num(n-1) = 240*n+31 */ - mpz_add_ui(num, num, 240*n+31); /* no overflow up to MPFR_PREC=71M */ - /* d2(n) - d2(n-1) = 12288*(n-1) */ - if (n>0) mpz_add_ui(d2, d2, 12288*(n-1)); - else mpz_sub_ui(d2, d2, 12288); - /* d3(n) - d3(n-1) = d2 */ - mpz_add(d3, d3, d2); - /* den(n)-den(n-1) = 2048*n^3 + 400n - 6 = d3 */ - mpz_add(den, den, d3); - mpz_mul_2exp(tmp, num, 4*(N-n)); - mpz_fdiv_q(tmp, tmp, den); - mpz_add(pi, pi, tmp); - } - mpfr_set_z(x, pi, rnd_mode); - mpfr_init2(y, mpfr_get_prec(x)); - mpz_add_ui(pi, pi, N+1); - mpfr_set_z(y, pi, rnd_mode); - if (mpfr_cmp(x, y) != 0) { - fprintf(stderr, "does not converge\n"); exit(1); + if (prec < 20000) + { + /* need to recompute */ + N=1; + do + { + oldN = N; + N = (prec+3)/4 + __gmpfr_ceil_log2((double) N + 1.0); + } + while (N != oldN); + mpz_init(pi); + mpz_init(num); + mpz_init(den); + mpz_init(d3); + mpz_init(d2); + mpz_init(tmp); + mpz_set_ui(pi, 0); + mpz_set_ui(num, 16); /* num(-1) */ + mpz_set_ui(den, 21); /* den(-1) */ + mpz_set_si(d3, -2454); + mpz_set_ui(d2, 14736); + /* invariants: num=120*n^2+151*n+47, + den=512*n^4+1024*n^3+712*n^2+194*n+15 + d3 = 2048*n^3+400*n-6, d2 = 6144*n^2-6144*n+2448 + */ + for (n=0; n<N; n++) + { + /* num(n)-num(n-1) = 240*n+31 */ + mpz_add_ui(num, num, 240*n+31); /* no overflow up to MPFR_PREC=71M */ + /* d2(n) - d2(n-1) = 12288*(n-1) */ + if (n>0) + mpz_add_ui(d2, d2, 12288*(n-1)); + else + mpz_sub_ui(d2, d2, 12288); + /* d3(n) - d3(n-1) = d2 */ + mpz_add(d3, d3, d2); + /* den(n)-den(n-1) = 2048*n^3 + 400n - 6 = d3 */ + mpz_add(den, den, d3); + mpz_mul_2exp(tmp, num, 4*(N-n)); + mpz_fdiv_q(tmp, tmp, den); + mpz_add(pi, pi, tmp); + } + mpfr_set_z(x, pi, rnd_mode); + mpfr_init2(y, mpfr_get_prec(x)); + mpz_add_ui(pi, pi, N+1); + mpfr_set_z(y, pi, rnd_mode); + if (mpfr_cmp(x, y) != 0) + { + fprintf(stderr, "does not converge\n"); + exit(1); + } + MPFR_SET_EXP (x, MPFR_GET_EXP(x) - 4*N); + mpz_clear(pi); + mpz_clear(num); + mpz_clear(den); + mpz_clear(d3); + mpz_clear(d2); + mpz_clear(tmp); + mpfr_clear(y); } - MPFR_EXP(x) -= 4*N; - mpz_clear(pi); mpz_clear(num); mpz_clear(den); mpz_clear(d3); mpz_clear(d2); - mpz_clear(tmp); mpfr_clear(y); - } else - mpfr_pi_machin3(x, rnd_mode); + else + mpfr_pi_machin3(x, rnd_mode); /* store computed value */ - if (__gmpfr_const_pi_prec==0) mpfr_init2(__mpfr_const_pi, prec); - else mpfr_set_prec(__mpfr_const_pi, prec); + if (__gmpfr_const_pi_prec==0) + mpfr_init2(__mpfr_const_pi, prec); + else + mpfr_set_prec(__mpfr_const_pi, prec); mpfr_set(__mpfr_const_pi, x, rnd_mode); __gmpfr_const_pi_prec=prec; __mpfr_const_pi_rnd=rnd_mode; @@ -1,6 +1,6 @@ /* mpfr_cos -- cosine of a floating-point number -Copyright 2001, 2002 Free Software Foundation. +Copyright 2001, 2002, 2003 Free Software Foundation. This file is part of the MPFR Library. @@ -59,7 +59,7 @@ mpfr_cos (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) mpfr_mul (r, x, x, GMP_RNDU); /* err <= 1 ulp */ /* we need that |r| < 1 for mpfr_cos2_aux, i.e. up(x^2)/2^(2K) < 1 */ - K = K0 + MAX(MPFR_EXP(r), 0); + K = K0 + MAX (MPFR_GET_EXP (r), 0); mpfr_div_2ui (r, r, 2 * K, GMP_RNDN); /* r = (x/2^K)^2, err <= 1 ulp */ @@ -77,7 +77,8 @@ mpfr_cos (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) for (k = 2 * K, l = 2 * l + 1; l > 1; k++, l = (l + 1) >> 1); /* now the error is bounded by 2^(k-m) = 2^(EXP(s)-err) */ - l = mpfr_can_round (s, MPFR_EXP(s) + m - k, GMP_RNDN, rnd_mode, precy); + l = mpfr_can_round (s, MPFR_GET_EXP (s) + m - k, + GMP_RNDN, rnd_mode, precy); if (l == 0) { @@ -108,12 +109,12 @@ mpfr_cos2_aux (mpfr_ptr s, mpfr_srcptr r) long int prec_t, m = MPFR_PREC(s); mpfr_t t; - MPFR_ASSERTN (MPFR_EXP(r) <= 0); + MPFR_ASSERTN (MPFR_GET_EXP (r) <= 0); mpfr_init2 (t, m); mpfr_set_ui (t, 1, GMP_RNDN); mpfr_set_ui(s, 1, GMP_RNDN); - for (l = 1; MPFR_EXP(t) + m >= 0; l++) + for (l = 1; MPFR_GET_EXP (t) + m >= 0; l++) { mpfr_mul (t, t, r, GMP_RNDU); /* err <= (3l-1) ulp */ mpfr_div_ui (t, t, (2*l-1)*(2*l), GMP_RNDU); /* err <= 3l ulp */ @@ -121,13 +122,13 @@ mpfr_cos2_aux (mpfr_ptr s, mpfr_srcptr r) mpfr_add (s, s, t, GMP_RNDD); else mpfr_sub (s, s, t, GMP_RNDD); - MPFR_ASSERTN (MPFR_EXP(s) == 0); /* check 1/2 <= s < 1 */ + MPFR_ASSERTN (MPFR_GET_EXP (s) == 0); /* check 1/2 <= s < 1 */ /* err(s) <= l * 2^(-m) */ if (3 * l > (1 << b)) b++; /* now 3l <= 2^b, we want 3l*ulp(t) <= 2^(-m) i.e. b+EXP(t)-PREC(t) <= -m */ - prec_t = m + MPFR_EXP(t) + b; + prec_t = m + MPFR_GET_EXP (t) + b; if (prec_t >= MPFR_PREC_MIN) mpfr_round_prec (t, GMP_RNDN, prec_t); } @@ -1,6 +1,6 @@ /* mpfr_div -- divide two floating-point numbers -Copyright 1999, 2001, 2002 Free Software Foundation. +Copyright 1999, 2001, 2002, 2003 Free Software Foundation. This file is part of the MPFR Library. @@ -331,7 +331,7 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mp_rnd_t rnd_mode) * [size rsize]. * **************************************************************************/ - qexp = MPFR_EXP(u) - MPFR_EXP(v); + qexp = MPFR_GET_EXP (u) - MPFR_GET_EXP (v); if (qp[qsize] != 0) /* Hack : qp[qsize] is 0, 1 or 2, hence if not 0, = 2^(qp[qsize] - 1). */ @@ -440,7 +440,7 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mp_rnd_t rnd_mode) TMP_FREE (marker); MPFR_MANT(q)[0] &= ~((MP_LIMB_T_ONE << rw) - MP_LIMB_T_ONE); - MPFR_EXP(q) = qexp; + MPFR_SET_EXP (q, qexp); MPFR_RET(inex); } @@ -1,6 +1,6 @@ /* mpfr_div_2si -- divide a floating-point number by a power of two -Copyright 1999, 2001, 2002 Free Software Foundation, Inc. +Copyright 1999, 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -33,22 +33,22 @@ mpfr_div_2si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd_mode) if (MPFR_IS_FP(y) && MPFR_NOTZERO(y)) { + mp_exp_t exp = MPFR_GET_EXP (y); if (n > 0 && (__gmpfr_emin > MPFR_EMAX_MAX - n || - MPFR_EXP(y) < __gmpfr_emin + n)) + exp < __gmpfr_emin + n)) { if (rnd_mode == GMP_RNDN && (__gmpfr_emin > MPFR_EMAX_MAX - (n - 1) || - MPFR_EXP(y) < __gmpfr_emin + (n - 1) || - mpfr_powerof2_raw (y))) + exp < __gmpfr_emin + (n - 1) || mpfr_powerof2_raw (y))) rnd_mode = GMP_RNDZ; return mpfr_set_underflow (y, rnd_mode, MPFR_SIGN(y)); } if (n < 0 && (__gmpfr_emax < MPFR_EMIN_MIN - n || - MPFR_EXP(y) > __gmpfr_emax + n)) + exp > __gmpfr_emax + n)) return mpfr_set_overflow (y, rnd_mode, MPFR_SIGN(y)); - MPFR_EXP(y) -= n; + MPFR_SET_EXP (y, exp - n); } return inexact; @@ -1,6 +1,6 @@ /* mpfr_div_2ui -- divide a floating-point number by a power of two -Copyright 1999, 2001, 2002 Free Software Foundation, Inc. +Copyright 1999, 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -48,18 +48,21 @@ mpfr_div_2ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int n, mp_rnd_t rnd_mode) /* MPFR_EMAX_MAX - (long) n is signed and doesn't lead to an integer overflow; the first test useful so that the real test can't lead to an integer overflow. */ - if (__gmpfr_emin > MPFR_EMAX_MAX - (long) n || - MPFR_EXP(y) < __gmpfr_emin + (long) n) - { - if (rnd_mode == GMP_RNDN && - (__gmpfr_emin > MPFR_EMAX_MAX - (long) (n - 1) || - MPFR_EXP(y) < __gmpfr_emin + (long) (n - 1) || - mpfr_powerof2_raw (y))) - rnd_mode = GMP_RNDZ; - return mpfr_set_underflow (y, rnd_mode, MPFR_SIGN(y)); - } + { + mp_exp_t exp = MPFR_GET_EXP (y); + if (__gmpfr_emin > MPFR_EMAX_MAX - (long) n || + exp < __gmpfr_emin + (long) n) + { + if (rnd_mode == GMP_RNDN && + (__gmpfr_emin > MPFR_EMAX_MAX - (long) (n - 1) || + exp < __gmpfr_emin + (long) (n - 1) || + mpfr_powerof2_raw (y))) + rnd_mode = GMP_RNDZ; + return mpfr_set_underflow (y, rnd_mode, MPFR_SIGN(y)); + } - MPFR_EXP(y) -= (long) n; + MPFR_SET_EXP(y, exp - (long) n); + } } return inexact; @@ -1,6 +1,6 @@ /* mpfr_div_ui -- divide a floating-point number by a machine integer -Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc. +Copyright 1999, 2000, 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -80,7 +80,7 @@ mpfr_div_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mp_rnd_t rnd_mode) xp = MPFR_MANT(x); yp = MPFR_MANT(y); - MPFR_EXP(y) = MPFR_EXP(x); + MPFR_SET_EXP (y, MPFR_GET_EXP (x)); dif = yn + 1 - xn; @@ -113,7 +113,7 @@ mpfr_div_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mp_rnd_t rnd_mode) { tmp--; sh = 0; - MPFR_EXP(y) -= BITS_PER_MP_LIMB; + MPFR_SET_EXP (y, MPFR_GET_EXP (y) - BITS_PER_MP_LIMB); } /* now we have yn limbs starting from tmp[1], with tmp[yn]<>0 */ @@ -126,7 +126,7 @@ mpfr_div_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mp_rnd_t rnd_mode) yp[0] += tmp[0] >> (BITS_PER_MP_LIMB - sh); middle = middle || ((tmp[0] << sh) != 0); inexact = inexact || ((tmp[0] << sh) != 0); - MPFR_EXP(y) -= sh; + MPFR_SET_EXP (y, MPFR_GET_EXP (y) - sh); } else MPN_COPY(yp, tmp + 1, yn); @@ -63,8 +63,8 @@ mpfr_eq (mpfr_srcptr u, mpfr_srcptr v, unsigned long int n_bits) return MPFR_IS_ZERO(u) && MPFR_IS_ZERO(v); } - uexp = MPFR_EXP(u); - vexp = MPFR_EXP(v); + uexp = MPFR_GET_EXP (u); + vexp = MPFR_GET_EXP (v); /* 2. Are the exponents different? */ if (uexp != vexp) diff --git a/exceptions.c b/exceptions.c index c2bf923a7..db62a01ee 100644 --- a/exceptions.c +++ b/exceptions.c @@ -1,6 +1,6 @@ /* Exception flags and utilities. -Copyright 2001, 2002 Free Software Foundation. +Copyright 2001, 2002, 2003 Free Software Foundation. This file is part of the MPFR Library. @@ -124,7 +124,7 @@ mpfr_check_range (mpfr_ptr x, int t, mp_rnd_t rnd_mode) { if (MPFR_IS_FP(x) && MPFR_NOTZERO(x)) { /* x is a non-zero FP */ - mp_exp_t exp = MPFR_EXP(x); + mp_exp_t exp = MPFR_EXP (x); /* Do not use MPFR_GET_EXP */ if (exp < __gmpfr_emin) { /* The following test is necessary because in the rounding to the @@ -1,6 +1,6 @@ /* mpfr_exp -- exponential of a floating-point number -Copyright 1999, 2000, 2001, 2002 Free Software Foundation. +Copyright 1999, 2000, 2001, 2002, 2003 Free Software Foundation. Contributed by the Spaces project. This file is part of the MPFR Library. @@ -66,7 +66,7 @@ mpfr_exp (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) if (MPFR_IS_ZERO(x)) return mpfr_set_ui (y, 1, GMP_RNDN); - expx = MPFR_EXP(x); + expx = MPFR_GET_EXP (x); precy = MPFR_PREC(y); /* result is +Inf when exp(x) >= 2^(__gmpfr_emax), i.e. @@ -125,7 +125,7 @@ mpfr_exp2 (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) mpfr_exp (t, te, GMP_RNDN); /* exp(x*ln(2))*/ /* estimate of the error -- see pow function in algorithms.ps*/ - err = Nt - (MPFR_EXP(te) + 2); + err = Nt - (MPFR_GET_EXP (te) + 2); /* actualisation of the precision */ Nt += __gmpfr_isqrt (Nt) + 10; @@ -107,7 +107,7 @@ mpfr_exp_rational (mpfr_ptr y, mpz_srcptr p, int r, int m) mpz_tdiv_q(S[0], S[0], P[0]); mpfr_set_z(y,S[0], GMP_RNDD); - MPFR_EXP(y) += expo; + MPFR_SET_EXP (y, MPFR_GET_EXP (y) + expo); mpfr_div_2ui(y, y, r*(i-1),GMP_RNDN); for (i=0;i<=m;i++) { mpz_clear(P[i]); mpz_clear(S[i]); mpz_clear(ptoj[i]); } @@ -143,7 +143,7 @@ mpfr_exp3 (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) if (prec_x < 0) prec_x = 0; logn = __gmpfr_ceil_log2 ((double) prec_x + MPFR_PREC(y)); if (logn < 2) logn = 2; - ttt = MPFR_EXP(x); + ttt = MPFR_GET_EXP (x); mpfr_init2(x_copy,MPFR_PREC(x)); mpfr_set(x_copy,x,GMP_RNDD); /* we shift to get a number less than 1 */ @@ -151,7 +151,7 @@ mpfr_exp3 (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) { shift_x = ttt; mpfr_div_2ui(x_copy, x, ttt, GMP_RNDN); - ttt = MPFR_EXP(x_copy); + ttt = MPFR_GET_EXP (x_copy); } realprec = MPFR_PREC(y)+logn; mpz_init (uk); @@ -196,7 +196,8 @@ mpfr_exp_2 (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) exps <<= 1; exps += mpz_normalize (ss, ss, q); } - mpfr_set_z (s, ss, GMP_RNDN); MPFR_EXP(s) += exps; + mpfr_set_z (s, ss, GMP_RNDN); + MPFR_SET_EXP(s, MPFR_GET_EXP (s) + exps); TMP_FREE(marker); /* don't need ss anymore */ if (n>0) mpfr_mul_2ui(s, s, n, GMP_RNDU); @@ -301,10 +302,11 @@ mpfr_exp2_aux2 (mpz_t s, mpfr_srcptr r, int q, int *exps) TMP_DECL(marker); /* estimate value of l */ - l = q / (-MPFR_EXP(r)); + l = q / (- MPFR_GET_EXP (r)); m = (int) __gmpfr_isqrt (l); /* we access R[2], thus we need m >= 2 */ - if (m < 2) m = 2; + if (m < 2) + m = 2; TMP_MARK(marker); R = (mpz_t*) TMP_ALLOC((m+1)*sizeof(mpz_t)); /* R[i] stands for r^i */ expR = (int*) TMP_ALLOC((m+1)*sizeof(int)); /* exponent for R[i] */ @@ -313,68 +315,85 @@ mpfr_exp2_aux2 (mpz_t s, mpfr_srcptr r, int q, int *exps) MY_INIT_MPZ(rr, sizer+2); MY_INIT_MPZ(t, 2*sizer); /* double size for products */ mpz_set_ui(s, 0); *exps = 1-q; /* initialize s to zero, 1 ulp = 2^(1-q) */ - for (i=0;i<=m;i++) MY_INIT_MPZ(R[i], sizer+2); + for (i=0;i<=m;i++) + MY_INIT_MPZ(R[i], sizer+2); expR[1] = mpfr_get_z_exp(R[1], r); /* exact operation: no error */ expR[1] = mpz_normalize2(R[1], R[1], expR[1], 1-q); /* error <= 1 ulp */ mpz_mul(t, R[1], R[1]); /* err(t) <= 2 ulps */ mpz_div_2exp(R[2], t, q-1); /* err(R[2]) <= 3 ulps */ expR[2] = 1-q; - for (i=3;i<=m;i++) { - mpz_mul(t, R[i-1], R[1]); /* err(t) <= 2*i-2 */ - mpz_div_2exp(R[i], t, q-1); /* err(R[i]) <= 2*i-1 ulps */ - expR[i] = 1-q; - } - mpz_set_ui(R[0], 1); mpz_mul_2exp(R[0], R[0], q-1); expR[0]=1-q; /* R[0]=1 */ - mpz_set_ui(rr, 1); expr=0; /* rr contains r^l/l! */ + for (i=3;i<=m;i++) + { + mpz_mul(t, R[i-1], R[1]); /* err(t) <= 2*i-2 */ + mpz_div_2exp(R[i], t, q-1); /* err(R[i]) <= 2*i-1 ulps */ + expR[i] = 1-q; + } + mpz_set_ui(R[0], 1); + mpz_mul_2exp(R[0], R[0], q-1); + expR[0]=1-q; /* R[0]=1 */ + mpz_set_ui(rr, 1); + expr=0; /* rr contains r^l/l! */ /* by induction: err(rr) <= 2*l ulps */ l = 0; ql = q; /* precision used for current giant step */ - do { - /* all R[i] must have exponent 1-ql */ - if (l) for (i=0;i<m;i++) { - expR[i] = mpz_normalize2(R[i], R[i], expR[i], 1-ql); - } - /* the absolute error on R[i]*rr is still 2*i-1 ulps */ - expt = mpz_normalize2(t, R[m-1], expR[m-1], 1-ql); - /* err(t) <= 2*m-1 ulps */ - /* computes t = 1 + r/(l+1) + ... + r^(m-1)*l!/(l+m-1)! - using Horner's scheme */ - for (i=m-2;i>=0;i--) { - mpz_div_ui(t, t, l+i+1); /* err(t) += 1 ulp */ - mpz_add(t, t, R[i]); - } - /* now err(t) <= (3m-2) ulps */ - - /* now multiplies t by r^l/l! and adds to s */ - mpz_mul(t, t, rr); expt += expr; - expt = mpz_normalize2(t, t, expt, *exps); - /* err(t) <= (3m-1) + err_rr(l) <= (3m-2) + 2*l */ + do + { + /* all R[i] must have exponent 1-ql */ + if (l) + for (i=0;i<m;i++) + { + expR[i] = mpz_normalize2(R[i], R[i], expR[i], 1-ql); + } + /* the absolute error on R[i]*rr is still 2*i-1 ulps */ + expt = mpz_normalize2(t, R[m-1], expR[m-1], 1-ql); + /* err(t) <= 2*m-1 ulps */ + /* computes t = 1 + r/(l+1) + ... + r^(m-1)*l!/(l+m-1)! + using Horner's scheme */ + for (i=m-2;i>=0;i--) + { + mpz_div_ui(t, t, l+i+1); /* err(t) += 1 ulp */ + mpz_add(t, t, R[i]); + } + /* now err(t) <= (3m-2) ulps */ + + /* now multiplies t by r^l/l! and adds to s */ + mpz_mul(t, t, rr); + expt += expr; + expt = mpz_normalize2(t, t, expt, *exps); + /* err(t) <= (3m-1) + err_rr(l) <= (3m-2) + 2*l */ #ifdef DEBUG - if (expt != *exps) { - fprintf(stderr, "Error: expt != exps %d %d\n", expt, *exps); exit(1); - } + if (expt != *exps) + { + fprintf(stderr, "Error: expt != exps %d %d\n", expt, *exps); + exit(1); + } #endif - mpz_add(s, s, t); /* no error here */ - - /* updates rr, the multiplication of the factors l+i could be done - using binary splitting too, but it is not sure it would save much */ - mpz_mul(t, rr, R[m]); /* err(t) <= err(rr) + 2m-1 */ - expr += expR[m]; - mpz_set_ui (tmp, 1); - for (i=1, c=1; i<=m; i++) { - if (l+i > ~((unsigned long int) 0)/c) { - mpz_mul_ui(tmp, tmp, c); - c = l+i; - } - else c *= (unsigned long int) l+i; + mpz_add(s, s, t); /* no error here */ + + /* updates rr, the multiplication of the factors l+i could be done + using binary splitting too, but it is not sure it would save much */ + mpz_mul(t, rr, R[m]); /* err(t) <= err(rr) + 2m-1 */ + expr += expR[m]; + mpz_set_ui (tmp, 1); + for (i=1, c=1; i<=m; i++) + { + if (l+i > ~((unsigned long int) 0)/c) + { + mpz_mul_ui(tmp, tmp, c); + c = l+i; + } + else + c *= (unsigned long int) l+i; + } + if (c != 1) + mpz_mul_ui (tmp, tmp, c); /* tmp is exact */ + mpz_fdiv_q(t, t, tmp); /* err(t) <= err(rr) + 2m */ + expr += mpz_normalize(rr, t, ql); /* err_rr(l+1) <= err_rr(l) + 2m+1 */ + ql = q - *exps - mpz_sizeinbase(s, 2) + expr + mpz_sizeinbase(rr, 2); + l+=m; } - if (c != 1) mpz_mul_ui (tmp, tmp, c); /* tmp is exact */ - mpz_fdiv_q(t, t, tmp); /* err(t) <= err(rr) + 2m */ - expr += mpz_normalize(rr, t, ql); /* err_rr(l+1) <= err_rr(l) + 2m+1 */ - ql = q - *exps - mpz_sizeinbase(s, 2) + expr + mpz_sizeinbase(rr, 2); - l+=m; - } while (expr+mpz_sizeinbase(rr, 2) > -q); + while (expr+mpz_sizeinbase(rr, 2) > -q); TMP_FREE(marker); mpz_clear(tmp); @@ -1,6 +1,6 @@ /* mpfr_expm1 -- Compute exp(x)-1 -Copyright 2001, 2002 Free Software Foundation. +Copyright 2001, 2002, 2003 Free Software Foundation. This file is part of the MPFR Library. @@ -86,29 +86,31 @@ mpfr_expm1 (mpfr_ptr y, mpfr_srcptr x , mp_rnd_t rnd_mode) mpfr_init(te); /* First computation of cosh */ - do { + do + { - /* reactualisation of the precision */ - mpfr_set_prec(t,Nt); - mpfr_set_prec(te,Nt); - - /* compute expm1 */ - mpfr_exp(te,x,GMP_RNDN); /* exp(x)*/ - mpfr_sub_ui(t,te,1,GMP_RNDN); /* exp(x)-1 */ + /* reactualisation of the precision */ + mpfr_set_prec(t,Nt); + mpfr_set_prec(te,Nt); - /* estimation of the error */ - /*err=Nt-(__gmpfr_ceil_log2(1+pow(2,MPFR_EXP(te)-MPFR_EXP(t))));*/ - err=Nt-(MAX(MPFR_EXP(te)-MPFR_EXP(t),0)+1); + /* compute expm1 */ + mpfr_exp(te,x,GMP_RNDN); /* exp(x)*/ + mpfr_sub_ui(t,te,1,GMP_RNDN); /* exp(x)-1 */ - /* actualisation of the precision */ - Nt += 10; + /* error estimate */ + /*err=Nt-(__gmpfr_ceil_log2(1+pow(2,MPFR_EXP(te)-MPFR_EXP(t))));*/ + err = Nt - (MAX (MPFR_GET_EXP (te) - MPFR_GET_EXP (t), 0) + 1); - } while ((err <0) || !mpfr_can_round(t,err,GMP_RNDN,rnd_mode,Ny)); + /* actualisation of the precision */ + Nt += 10; + } + while ((err <0) || !mpfr_can_round(t,err,GMP_RNDN,rnd_mode,Ny)); - inexact = mpfr_set(y,t,rnd_mode); + inexact = mpfr_set(y,t,rnd_mode); - mpfr_clear(t); - mpfr_clear(te); + mpfr_clear(t); + mpfr_clear(te); } + return inexact; } @@ -50,7 +50,7 @@ mpfr_frac (mpfr_ptr r, mpfr_srcptr u, mp_rnd_t rnd_mode) MPFR_RET(0); /* zero is exact */ } - ue = MPFR_EXP(u); + ue = MPFR_GET_EXP (u); if (ue <= 0) /* |u| < 1 */ return mpfr_set(r, u, rnd_mode); @@ -94,7 +94,7 @@ mpfr_frac (mpfr_ptr r, mpfr_srcptr u, mp_rnd_t rnd_mode) of u into account, because of the mpn_lshift below. */ MPFR_CLEAR_FLAGS(t); MPFR_SET_SAME_SIGN(t, u); - MPFR_EXP(t) = re; + MPFR_SET_EXP (t, re); /* Put the fractional part of u into t */ tn = (MPFR_PREC(t) - 1) / BITS_PER_MP_LIMB; @@ -1,6 +1,6 @@ /* generic file for evaluation of hypergeometric series using binary splitting -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. @@ -178,7 +178,7 @@ GENERIC (mpfr_ptr y, mpz_srcptr p, int r, int m) mpz_tdiv_q(S[0], S[0], P[0]); mpfr_set_z(y, S[0], GMP_RNDD); - MPFR_EXP(y) += expo; + MPFR_SET_EXP (y, MPFR_GET_EXP (y) + expo); #ifdef R_IS_RATIONAL /* exact division */ @@ -186,19 +186,19 @@ double mpfr_get_d (mpfr_srcptr src, mp_rnd_t rnd_mode) { return mpfr_get_d3 (src, MPFR_IS_FP(src) && MPFR_NOTZERO(src) ? - MPFR_EXP(src) : 0, rnd_mode); + MPFR_GET_EXP (src) : 0, rnd_mode); } double mpfr_get_d1 (mpfr_srcptr src) { return mpfr_get_d3 (src, MPFR_IS_FP(src) && MPFR_NOTZERO(src) ? - MPFR_EXP(src) : 0, __gmpfr_default_rounding_mode); + MPFR_GET_EXP (src) : 0, __gmpfr_default_rounding_mode); } double mpfr_get_d_2exp (mp_exp_t *exp, mpfr_srcptr src, mp_rnd_t rnd_mode) { - *exp = MPFR_EXP (src); + *exp = MPFR_GET_EXP (src); return mpfr_get_d3 (src, 0, rnd_mode); } @@ -1,6 +1,6 @@ /* mpfr_get_exp - get the exponent of a floating-point number -Copyright 2002 Free Software Foundation. +Copyright 2002, 2003 Free Software Foundation. This file is part of the MPFR Library. @@ -28,5 +28,5 @@ mp_exp_t mpfr_get_exp (mpfr_srcptr x) { MPFR_ASSERTN(MPFR_IS_FP(x) && MPFR_NOTZERO(x)); - return MPFR_EXP(x); + return MPFR_EXP(x); /* do not use MPFR_GET_EXP of course... */ } @@ -74,16 +74,16 @@ mpfr_get_ld (mpfr_srcptr x, mp_rnd_t rnd_mode) problems) */ mpfr_init2 (y, MPFR_LDBL_MANT_DIG); mpfr_set (y, x, rnd_mode); - e = MPFR_EXP(y); + e = MPFR_GET_EXP (y); if (e > 1023) { sh = e - 1023; - MPFR_EXP(y) = 1023; + MPFR_SET_EXP (y, 1023); } else if (e < -1021) { sh = e + 1021; - MPFR_EXP(y) = -1021; + MPFR_SET_EXP (y, -1021); } else { @@ -47,7 +47,7 @@ mpfr_get_si (mpfr_srcptr f, mp_rnd_t rnd) more limbs */ /* now the result is in the most significant limb of x */ - exp = MPFR_EXP(x); /* since |x| >= 1, exp >= 1 */ + exp = MPFR_GET_EXP (x); /* since |x| >= 1, exp >= 1 */ n = MPFR_ABSSIZE(x); s = MPFR_MANT(x)[n - 1] >> (GMP_NUMB_BITS - exp); @@ -620,8 +620,8 @@ mpfr_get_str (char *s, mp_exp_t *e, int b, size_t m, mpfr_srcptr x, mp_rnd_t rnd pow2 = BITS_PER_MP_LIMB - pow2 - 1; /* base = 2^pow2 */ /* set MPFR_EXP(x) = f*pow2 + r, 1 <= r <= pow2 */ - f = (MPFR_EXP(x) - 1) / pow2; - r = MPFR_EXP(x) - f * pow2; + f = (MPFR_GET_EXP (x) - 1) / pow2; + r = MPFR_GET_EXP (x) - f * pow2; if (r <= 0) { f --; @@ -675,7 +675,7 @@ mpfr_get_str (char *s, mp_exp_t *e, int b, size_t m, mpfr_srcptr x, mp_rnd_t rnd return (s0); } - g = mpfr_get_str_compute_g (b, MPFR_EXP(x) - 1); + g = mpfr_get_str_compute_g (b, MPFR_GET_EXP (x) - 1); exact = 1; prec = (mp_exp_t) _mpfr_ceil ((double) m / log_b2[b-2]) + 1; exp = ((mp_exp_t) m < g) ? g - (mp_exp_t) m : (mp_exp_t) m - g; @@ -705,7 +705,7 @@ mpfr_get_str (char *s, mp_exp_t *e, int b, size_t m, mpfr_srcptr x, mp_rnd_t rnd exact = mpn_scan1 (xp, 0) >= (nx - n) * BITS_PER_MP_LIMB; err = !exact; MPN_COPY2 (a, n, xp, nx); - exp_a = MPFR_EXP(x) - n * BITS_PER_MP_LIMB; + exp_a = MPFR_GET_EXP (x) - n * BITS_PER_MP_LIMB; } else if ((mp_exp_t) m > g) /* we have to multiply x by b^exp */ { @@ -730,7 +730,7 @@ mpfr_get_str (char *s, mp_exp_t *e, int b, size_t m, mpfr_srcptr x, mp_rnd_t rnd /* result = a * x */ result = (mp_limb_t*) TMP_ALLOC ((n + nx1) * sizeof (mp_limb_t)); mpn_mul (result, a, n, x1, nx1); - exp_a += MPFR_EXP (x); + exp_a += MPFR_GET_EXP (x); if (mpn_scan1 (result, 0) < (nx1 * BITS_PER_MP_LIMB)) exact = 0; @@ -763,7 +763,7 @@ mpfr_get_str (char *s, mp_exp_t *e, int b, size_t m, mpfr_srcptr x, mp_rnd_t rnd /* result = x / a */ mpn_tdiv_qr (result, reste, 0, x1, 2 * n, a, n); - exp_a = MPFR_EXP (x) - exp_a - 2 * n * BITS_PER_MP_LIMB; + exp_a = MPFR_GET_EXP (x) - exp_a - 2 * n * BITS_PER_MP_LIMB; /* test if division was exact */ if (exact) @@ -47,7 +47,7 @@ mpfr_get_ui (mpfr_srcptr f, mp_rnd_t rnd) more limbs */ /* now the result is in the most significant limb of x */ - exp = MPFR_EXP(x); /* since |x| >= 1, exp >= 1 */ + exp = MPFR_GET_EXP (x); /* since |x| >= 1, exp >= 1 */ n = MPFR_ABSSIZE(x); s = MPFR_MANT(x)[n - 1] >> (GMP_NUMB_BITS - exp); diff --git a/get_z_exp.c b/get_z_exp.c index 5777b3594..5eaa3c62f 100644 --- a/get_z_exp.c +++ b/get_z_exp.c @@ -67,8 +67,8 @@ mpfr_get_z_exp (mpz_ptr z, mpfr_srcptr f) /* This always fails for very small "f", ie. when MPFR_EXP(f) is equal to or only just above MPFR_EMIN_MIN. - MPFR_ASSERTN((mp_exp_unsigned_t) MPFR_EXP(f) - MPFR_EMIN_MIN + MPFR_ASSERTN((mp_exp_unsigned_t) MPFR_GET_EXP (f) - MPFR_EMIN_MIN >= (mp_exp_unsigned_t) MPFR_PREC(f)); */ - return MPFR_EXP(f) - MPFR_PREC(f); + return MPFR_GET_EXP (f) - MPFR_PREC (f); } @@ -1,6 +1,6 @@ /* mpfr_hypot -- Euclidean distance -Copyright 2001, 2002 Free Software Foundation, Inc. +Copyright 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -78,8 +78,8 @@ mpfr_hypot (mpfr_ptr z, mpfr_srcptr x , mpfr_srcptr y , mp_rnd_t rnd_mode) /* now |x| >= |y| */ - Ex = MPFR_EXP(x); - Ey = MPFR_EXP(y); + Ex = MPFR_GET_EXP (x); + Ey = MPFR_GET_EXP (y); diff_exp = (mp_exp_unsigned_t) Ex - Ey; Nz = MPFR_PREC(z); /* Precision of output variable */ diff --git a/isinteger.c b/isinteger.c index ab2732159..720f4d021 100644 --- a/isinteger.c +++ b/isinteger.c @@ -38,7 +38,7 @@ mpfr_integer_p (mpfr_srcptr x) if (MPFR_IS_ZERO(x)) return 1; - expo = MPFR_EXP(x); + expo = MPFR_GET_EXP (x); if (expo <= 0) return 0; @@ -1,6 +1,6 @@ /* mpfr_log -- natural logarithm of a floating-point number -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. @@ -121,7 +121,7 @@ mpfr_log (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) printf("p=%d\n", p); #endif /* Calculus of m (depends on p) */ - m = (p + 1) / 2 - MPFR_EXP(a) + 1; + m = (p + 1) / 2 - MPFR_GET_EXP (a) + 1; /* All the mpfr_t needed have a precision of p */ TMP_MARK(marker); @@ -145,9 +145,9 @@ mpfr_log (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) mpfr_div (tmp2, cst, tmp1, GMP_RNDN); /* pi/2*AG(1,4/s), err<=5ulps */ mpfr_const_log2 (cst, GMP_RNDN); /* compute log(2), err<=1ulp */ mpfr_mul(tmp1,cst,mm,GMP_RNDN); /* I compute m*log(2), err<=2ulps */ - cancel = MPFR_EXP(tmp2); + cancel = MPFR_GET_EXP (tmp2); mpfr_sub(cst,tmp2,tmp1,GMP_RNDN); /* log(a), err<=7ulps+cancel */ - cancel -= MPFR_EXP(cst); + cancel -= MPFR_GET_EXP (cst); #ifdef DEBUG printf("canceled bits=%d\n", cancel); printf("approx="); mpfr_print_binary(cst); putchar('\n'); @@ -1,6 +1,6 @@ /* mpfr_log1p -- Compute log(1+x) -Copyright 2001, 2002 Free Software Foundation, Inc. +Copyright 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -113,8 +113,8 @@ mpfr_log1p (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) mpfr_log (t, t, GMP_RNDN); /* log(1+x)*/ /* estimation of the error */ - /*err=Nt-(__gmpfr_ceil_log2(1+pow(2,1-MPFR_EXP(t))));*/ - err=Nt-(MAX(1-MPFR_EXP(t),0)+1); + /*err=Nt-(__gmpfr_ceil_log2(1+pow(2,1-MPFR_GET_EXP(t))));*/ + err = Nt - (MAX (1 - MPFR_GET_EXP (t), 0) + 1); /* actualisation of the precision */ Nt += 10; @@ -1,6 +1,6 @@ /* mpfr_log2 -- log base 2 -Copyright 2001, 2002 Free Software Foundation, Inc. +Copyright 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -86,11 +86,10 @@ mpfr_log2 (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) } /* If a is integer, log2(a) is exact*/ - if (mpfr_cmp_ui_2exp(a,1,MPFR_EXP(a)-1) == 0) - return mpfr_set_si(r,MPFR_EXP(a)-1,rnd_mode); + if (mpfr_cmp_ui_2exp (a, 1, MPFR_GET_EXP (a) - 1) == 0) + return mpfr_set_si(r, MPFR_GET_EXP (a) - 1, rnd_mode); - - /* General case */ + /* General case */ { /* Declaration of the intermediary variable */ mpfr_t t, tt; diff --git a/mpfr-impl.h b/mpfr-impl.h index 0da56412f..284b42322 100644 --- a/mpfr-impl.h +++ b/mpfr-impl.h @@ -84,12 +84,17 @@ typedef unsigned long int mp_size_unsigned_t; #define MPFR_ASSERTD(expr) ((void) 0) #endif -#if WANT_ASSERT -#define MPFR_GET_EXP(x) mpfr_get_exp(x) -#define MPFR_SET_EXP(x, exp) MPFR_ASSERTD (!mpfr_set_exp ((x), (exp))) +/* Invalid exponent value (to track bugs...) */ +#define MPFR_EXP_INVALID ((mp_exp_t) 1 << 30) + +#if MPFR_EXP_CHECK +#define MPFR_GET_EXP(x) mpfr_get_exp (x) +#define MPFR_SET_EXP(x, exp) MPFR_ASSERTN (!mpfr_set_exp ((x), (exp))) +#define MPFR_SET_INVALID_EXP(x) ((void) (MPFR_EXP (x) = MPFR_EXP_INVALID)) #else -#define MPFR_GET_EXP(x) MPFR_EXP(x) -#define MPFR_SET_EXP(x, exp) ((void) (MPFR_EXP(x) = (exp))) +#define MPFR_GET_EXP(x) MPFR_EXP (x) +#define MPFR_SET_EXP(x, exp) ((void) (MPFR_EXP (x) = (exp))) +#define MPFR_SET_INVALID_EXP(x) ((void) 0) #endif /* Definition of constants */ @@ -181,11 +186,15 @@ long double __gmpfr_longdouble_volatile __GMP_PROTO ((long double)) ATTRIBUTE_CO #define MPFR_CLEAR_FLAGS(x) \ (((x) -> _mpfr_size &= ~((mp_size_unsigned_t) 3 << 29))) #define MPFR_IS_NAN(x) (((x)->_mpfr_size >> 30) & 1) -#define MPFR_SET_NAN(x) ((x)->_mpfr_size |= ((mp_size_unsigned_t) 1 << 30)) +#define MPFR_SET_NAN(x) \ + (MPFR_SET_INVALID_EXP(x), \ + (x)->_mpfr_size |= ((mp_size_unsigned_t) 1 << 30)) #define MPFR_CLEAR_NAN(x) \ (((x) -> _mpfr_size &= ~((mp_size_unsigned_t) 1 << 30))) #define MPFR_IS_INF(x) (((x)->_mpfr_size >> 29) & 1) -#define MPFR_SET_INF(x) ((x)->_mpfr_size |= ((mp_size_unsigned_t) 1 << 29)) +#define MPFR_SET_INF(x) \ + (MPFR_SET_INVALID_EXP(x), \ + (x)->_mpfr_size |= ((mp_size_unsigned_t) 1 << 29)) #define MPFR_CLEAR_INF(x) ((x)->_mpfr_size &= ~((mp_size_unsigned_t) 1 << 29)) #define MPFR_IS_FP(x) ((((x) -> _mpfr_size >> 29) & 3) == 0) #define MPFR_ABSSIZE(x) \ @@ -209,7 +218,8 @@ long double __gmpfr_longdouble_volatile __GMP_PROTO ((long double)) ATTRIBUTE_CO #define MPFR_IS_ZERO(x) \ (MPFR_MANT(x)[(MPFR_PREC(x)-1)/BITS_PER_MP_LIMB] == (mp_limb_t) 0) #define MPFR_SET_ZERO(x) \ - (MPFR_MANT(x)[(MPFR_PREC(x)-1)/BITS_PER_MP_LIMB] = (mp_limb_t) 0) + (MPFR_SET_INVALID_EXP(x), \ + MPFR_MANT(x)[(MPFR_PREC(x)-1)/BITS_PER_MP_LIMB] = (mp_limb_t) 0) #define MPFR_ESIZE(x) \ ((MPFR_PREC((x)) - 1) / BITS_PER_MP_LIMB + 1) #define MPFR_EVEN_INEX 2 @@ -229,7 +239,7 @@ long double __gmpfr_longdouble_volatile __GMP_PROTO ((long double)) ATTRIBUTE_CO MPFR_PREC(x) = (p), \ MPFR_MANT(x) = (xp), \ MPFR_SIZE(x) = (s), \ - MPFR_SET_EXP((x), 0)) + MPFR_SET_INVALID_EXP(x)) /* same when xp is already allocated */ #define MPFR_INIT1(xp, x, p, s) \ (MPFR_PREC(x) = (p), MPFR_MANT(x) = (xp), MPFR_SIZE(x) = (s)) @@ -1,6 +1,6 @@ /* mpfr_mul -- multiply two floating-point numbers -Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc. +Copyright 1999, 2000, 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -88,8 +88,8 @@ mpfr_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) MPFR_RET(0); /* 0 * 0 is exact */ } - bx = MPFR_EXP(b); - cx = MPFR_EXP(c); + bx = MPFR_GET_EXP (b); + 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)); @@ -155,7 +155,7 @@ mpfr_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) rnd_mode = GMP_RNDZ; return mpfr_set_underflow (a, rnd_mode, sign_product); } - MPFR_EXP(a) = ax; + MPFR_SET_EXP (a, ax); if (MPFR_SIGN(a) != sign_product) MPFR_CHANGE_SIGN(a); @@ -1,6 +1,6 @@ /* mpfr_mul_2si -- multiply a floating-point number by a power of two -Copyright 1999, 2001, 2002 Free Software Foundation, Inc. +Copyright 1999, 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -33,22 +33,22 @@ mpfr_mul_2si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd_mode) if (MPFR_IS_FP(y) && MPFR_NOTZERO(y)) { + mp_exp_t exp = MPFR_GET_EXP (y); if (n > 0 && (__gmpfr_emax < MPFR_EMIN_MIN + n || - MPFR_EXP(y) > __gmpfr_emax - n)) + exp > __gmpfr_emax - n)) return mpfr_set_overflow (y, rnd_mode, MPFR_SIGN(y)); if (n < 0 && (__gmpfr_emin > MPFR_EMAX_MAX + n || - MPFR_EXP(y) < __gmpfr_emin - n)) + exp < __gmpfr_emin - n)) { if (rnd_mode == GMP_RNDN && (__gmpfr_emin > MPFR_EMAX_MAX + (n + 1) || - MPFR_EXP(y) < __gmpfr_emin - (n + 1) || - mpfr_powerof2_raw (y))) + exp < __gmpfr_emin - (n + 1) || mpfr_powerof2_raw (y))) rnd_mode = GMP_RNDZ; return mpfr_set_underflow (y, rnd_mode, MPFR_SIGN(y)); } - MPFR_EXP(y) += n; + MPFR_SET_EXP (y, exp + n); } return inexact; @@ -1,6 +1,6 @@ /* mpfr_mul_2ui -- multiply a floating-point number by a power of two -Copyright 1999, 2001 Free Software Foundation, Inc. +Copyright 1999, 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -48,11 +48,14 @@ mpfr_mul_2ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int n, mp_rnd_t rnd_mode) /* MPFR_EMIN_MIN + (long) n is signed and doesn't lead to an overflow; the first test useful so that the real test can't lead to an overflow. */ - if (__gmpfr_emax < MPFR_EMIN_MIN + (long) n || - MPFR_EXP(y) > __gmpfr_emax - (long) n) - return mpfr_set_overflow (y, rnd_mode, MPFR_SIGN(y)); - - MPFR_EXP(y) += (long) n; + { + mp_exp_t exp = MPFR_GET_EXP (y); + if (__gmpfr_emax < MPFR_EMIN_MIN + (long) n || + exp > __gmpfr_emax - (long) n) + return mpfr_set_overflow (y, rnd_mode, MPFR_SIGN(y)); + + MPFR_SET_EXP (y, exp + (long) n); + } } return inexact; @@ -1,6 +1,6 @@ /* mpfr_mul_ui -- multiply a floating-point number by a machine integer -Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc. +Copyright 1999, 2000, 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -112,10 +112,11 @@ mpfr_mul_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mp_rnd_t rnd_mode) TMP_FREE(marker); - if (__gmpfr_emax < MPFR_EMAX_MIN + cnt || MPFR_EXP(x) > __gmpfr_emax - cnt) + if (__gmpfr_emax < MPFR_EMAX_MIN + cnt || + MPFR_GET_EXP (x) > __gmpfr_emax - cnt) return mpfr_set_overflow(y, rnd_mode, MPFR_SIGN(x)); - MPFR_EXP(y) = MPFR_EXP(x) + cnt; + MPFR_SET_EXP (y, MPFR_GET_EXP (x) + cnt); MPFR_SET_SAME_SIGN(y, x); return inexact; @@ -1,7 +1,7 @@ /* mpfr_nextabove, mpfr_nextbelow, mpfr_nexttoward -- next representable floating-point number -Copyright 1999, 2001, 2002 Free Software Foundation. +Copyright 1999, 2001, 2002, 2003 Free Software Foundation. Contributed by the Spaces project, INRIA Lorraine. This file is part of the MPFR Library. @@ -53,13 +53,13 @@ mpfr_nexttozero (mpfr_ptr x) mpn_sub_1 (xp, xp, xn, MP_LIMB_T_ONE << sh); if (xp[xn-1] >> (BITS_PER_MP_LIMB - 1) == 0) { /* was an exact power of two: not normalized any more */ - mp_exp_t exp = MPFR_EXP(x); + mp_exp_t exp = MPFR_GET_EXP (x); if (exp == __gmpfr_emin) MPFR_SET_ZERO(x); else { mp_size_t i; - MPFR_EXP(x)--; + MPFR_SET_EXP (x, exp - 1); xp[0] = MP_LIMB_T_MAX << sh; for (i = 1; i < xn; i++) xp[i] = MP_LIMB_T_MAX; @@ -87,12 +87,12 @@ mpfr_nexttoinf (mpfr_ptr x) xp = MPFR_MANT(x); if (mpn_add_1 (xp, xp, xn, MP_LIMB_T_ONE << sh)) /* got 1.0000... */ { - mp_exp_t exp = MPFR_EXP(x); + mp_exp_t exp = MPFR_GET_EXP (x); if (exp == __gmpfr_emax) MPFR_SET_INF(x); else { - MPFR_EXP(x)++; + MPFR_SET_EXP (x, exp + 1); xp[xn-1] = MPFR_LIMB_HIGHBIT; } } @@ -41,11 +41,11 @@ mpfr_pow_is_exact (mpfr_srcptr x, mpfr_srcptr y) return 0; if (mpfr_sgn (y) < 0) - return mpfr_cmp_si_2exp (x, MPFR_SIGN(x), MPFR_EXP(x) - 1) == 0; + return mpfr_cmp_si_2exp (x, MPFR_SIGN(x), MPFR_GET_EXP (x) - 1) == 0; /* compute d such that y = c*2^d with c odd integer */ ysize = 1 + (MPFR_PREC(y) - 1) / BITS_PER_MP_LIMB; - d = MPFR_EXP(y) - ysize * BITS_PER_MP_LIMB; + d = MPFR_GET_EXP (y) - ysize * BITS_PER_MP_LIMB; /* since y is not zero, necessarily one of the mantissa limbs is not zero, thus we can simply loop until we find a non zero limb */ yp = MPFR_MANT(y); @@ -98,7 +98,7 @@ is_odd (mpfr_srcptr y) if (MPFR_IS_ZERO(y)) return 0; - expo = MPFR_EXP(y); + expo = MPFR_GET_EXP (y); if (expo <= 0) return 0; /* |y| < 1 and not 0 */ @@ -329,7 +329,7 @@ mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode) /* otherwise MPFR_EXP(te) below doesn't exist */ /* estimate of the error -- see pow function in algorithms.ps */ - err = Nt - (MPFR_EXP(te) + 3); + err = Nt - (MPFR_GET_EXP (te) + 3); /* actualisation of the precision */ Nt += 10; diff --git a/print_raw.c b/print_raw.c index 4622ff764..a04413987 100644 --- a/print_raw.c +++ b/print_raw.c @@ -1,7 +1,7 @@ /* mpfr_print_binary -- print the internal binary representation of a floating-point number -Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc. +Copyright 1999, 2000, 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -34,7 +34,7 @@ mpfr_get_str_raw (char *digit_ptr, mpfr_srcptr x) mp_limb_t *mx, wd, t; long ex, sx, k, l, p; mx = MPFR_MANT(x); - ex = MPFR_EXP(x); + ex = MPFR_GET_EXP (x); p = MPFR_PREC(x); if (MPFR_SIGN(x) < 0) { *digit_ptr = '-'; digit_ptr++; } @@ -1,6 +1,6 @@ /* mpfr_random -- generate a random floating-point number -Copyright 1999, 2001, 2002 Free Software Foundation, Inc. +Copyright 1999, 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -50,7 +50,7 @@ mpfr_random (mpfr_ptr x) count_leading_zeros (cnt, xp[xn - 1]); if (cnt != 0) mpn_lshift (xp, xp, xn, cnt); - MPFR_EXP(x) = -cnt; + MPFR_SET_EXP (x, -cnt); MPFR_SET_POS(x); cnt = (mp_prec_t) xn * BITS_PER_MP_LIMB - prec; @@ -2,7 +2,7 @@ long runs of consecutive ones and zeros in the binary representation. Intended for testing of other MP routines. -Copyright 1999, 2001, 2002 Free Software Foundation, Inc. +Copyright 1999, 2001, 2002, 2003 Free Software Foundation, Inc. (Copied from the GNU MP Library.) This file is part of the MPFR Library. @@ -52,7 +52,7 @@ mpfr_random2 (mpfr_ptr x, mp_size_t size, mp_exp_t exp) mpn_lshift (xp, xp, xn, cnt); mpn_random ((mp_limb_t*) yp, 1); - MPFR_EXP(x) = ABS ((mp_exp_t) yp[0] % (2 * exp + 1)) - exp; + MPFR_SET_EXP (x, ABS ((mp_exp_t) yp[0] % (2 * exp + 1)) - exp); cnt = xn * BITS_PER_MP_LIMB - MPFR_PREC(x); /* cnt is the number of non significant bits in the low limb */ @@ -1,6 +1,6 @@ /* mpfr_rint -- Round to an integer. -Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc. +Copyright 1999, 2000, 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -54,7 +54,7 @@ mpfr_rint (mpfr_ptr r, mpfr_srcptr u, mp_rnd_t rnd_mode) } sign = MPFR_SIGN(u); - exp = MPFR_EXP(u); + exp = MPFR_GET_EXP (u); /* Single out the case where |u| < 1. */ if (exp <= 0) /* 0 < |u| < 1 */ @@ -70,7 +70,7 @@ mpfr_rint (mpfr_ptr r, mpfr_srcptr u, mp_rnd_t rnd_mode) rm = (MPFR_PREC(r) - 1) / BITS_PER_MP_LIMB; rp[rm] = MPFR_LIMB_HIGHBIT; MPN_ZERO(rp, rm); - MPFR_EXP(r) = 1; /* |r| = 1 */ + MPFR_SET_EXP (r, 1); /* |r| = 1 */ MPFR_RET(sign > 0 ? 2 : -2); } else @@ -100,7 +100,7 @@ mpfr_rint (mpfr_ptr r, mpfr_srcptr u, mp_rnd_t rnd_mode) rn = MPFR_ESIZE(r); sh = (mp_prec_t) rn * BITS_PER_MP_LIMB - MPFR_PREC(r); - MPFR_EXP(r) = exp; + MPFR_SET_EXP (r, exp); if ((exp - 1) / BITS_PER_MP_LIMB >= un) { @@ -217,7 +217,7 @@ mpfr_rint (mpfr_ptr r, mpfr_srcptr u, mp_rnd_t rnd_mode) uflags : -uflags; else { - MPFR_EXP(r)++; + MPFR_SET_EXP(r, MPFR_GET_EXP (r) + 1); rp[rn-1] = MPFR_LIMB_HIGHBIT; } } diff --git a/round_prec.c b/round_prec.c index cfd066b7b..86b74e95e 100644 --- a/round_prec.c +++ b/round_prec.c @@ -1,7 +1,7 @@ /* mpfr_round_raw_generic, mpfr_round_raw2, mpfr_round_raw, mpfr_round_prec, mpfr_can_round, mpfr_can_round_raw -- various rounding functions -Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc. +Copyright 1999, 2000, 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -191,13 +191,13 @@ mpfr_round_prec (mpfr_ptr x, mp_rnd_t rnd_mode, mp_prec_t prec) if (carry) { - mp_exp_t exp = MPFR_EXP(x); + mp_exp_t exp = MPFR_GET_EXP (x); if (exp == __gmpfr_emax) (void) mpfr_set_overflow(x, rnd_mode, MPFR_SIGN(x)); else { - MPFR_EXP(x)++; + MPFR_SET_EXP (x, exp + 1); xp[nw - 1] = MPFR_LIMB_HIGHBIT; if (nw - 1 > 0) MPN_ZERO(xp, nw - 1); @@ -1,6 +1,6 @@ /* mpfr_set -- copy of a floating-point number -Copyright 1999, 2001, 2002 Free Software Foundation. +Copyright 1999, 2001, 2002, 2003 Free Software Foundation. This file is part of the MPFR Library. @@ -62,16 +62,16 @@ mpfr_set4 (mpfr_ptr a, mpfr_srcptr b, mp_rnd_t rnd_mode, int signb) carry = mpfr_round_raw(ap, MPFR_MANT(b), MPFR_PREC(b), (signb < 0), aq, rnd_mode, &inex); - MPFR_EXP(a) = MPFR_EXP(b); + MPFR_SET_EXP (a, MPFR_GET_EXP (b)); if (carry) { - mp_exp_t exp = MPFR_EXP(a); + mp_exp_t exp = MPFR_GET_EXP (a); if (exp == __gmpfr_emax) return mpfr_set_overflow(a, rnd_mode, signb); - MPFR_EXP(a)++; + MPFR_SET_EXP(a, exp + 1); ap[(MPFR_PREC(a)-1)/BITS_PER_MP_LIMB] = MPFR_LIMB_HIGHBIT; } } @@ -1,7 +1,7 @@ /* mpfr_set_d -- convert a machine double precision float to a multiple precision floating-point number -Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc. +Copyright 1999, 2000, 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -203,7 +203,7 @@ mpfr_set_d (mpfr_ptr r, double d, mp_rnd_t rnd_mode) signd = (d < 0) ? -1 : 1; d = ABS (d); - MPFR_EXP(tmp) = __mpfr_extract_double (tmpmant, d); + MPFR_SET_EXP(tmp, __mpfr_extract_double (tmpmant, d)); #ifndef NDEBUG /* Failed assertion if the stored value is 0 (e.g., if the exponent range @@ -233,7 +233,7 @@ mpfr_set_d (mpfr_ptr r, double d, mp_rnd_t rnd_mode) if (k) MPN_ZERO (tmpmant, k); - MPFR_EXP(tmp) -= cnt + k * BITS_PER_MP_LIMB; + MPFR_SET_EXP (tmp, MPFR_GET_EXP (tmp) - (cnt + k * BITS_PER_MP_LIMB)); /* tmp is exact since PREC(tmp)=53 */ inexact = mpfr_set4 (r, tmp, rnd_mode, signd); @@ -1,6 +1,6 @@ /* mpfr_set_exp - set the exponent of a floating-point number -Copyright 2002 Free Software Foundation. +Copyright 2002, 2003 Free Software Foundation. This file is part of the MPFR Library. @@ -29,7 +29,7 @@ mpfr_set_exp (mpfr_ptr x, mp_exp_t exponent) { if (exponent >= __gmpfr_emin && exponent <= __gmpfr_emax) { - MPFR_EXP(x) = exponent; + MPFR_EXP(x) = exponent; /* do not use MPFR_SET_EXP of course... */ return 0; } else @@ -1,6 +1,6 @@ /* mpfr_set_f -- set a MPFR number from a GNU MPF number -Copyright 1999, 2000, 2001 Free Software Foundation, Inc. +Copyright 1999, 2000, 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -77,7 +77,7 @@ mpfr_set_f (mpfr_ptr y, mpf_srcptr x, mp_rnd_t rnd_mode) inexact = 0; } - MPFR_EXP(y) = EXP(x) * BITS_PER_MP_LIMB - cnt; + MPFR_SET_EXP(y, EXP(x) * BITS_PER_MP_LIMB - cnt); return inexact; } @@ -1,6 +1,6 @@ /* mpfr_set_si -- set a MPFR number from a machine signed integer -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. @@ -54,7 +54,8 @@ mpfr_set_si (mpfr_ptr x, long i, mp_rnd_t rnd_mode) if ((i < 0) ^ (MPFR_SIGN(x) < 0)) MPFR_CHANGE_SIGN(x); - MPFR_EXP(x) = nbits = BITS_PER_MP_LIMB - cnt; + nbits = BITS_PER_MP_LIMB - cnt; + MPFR_SET_EXP (x, nbits); inex = mpfr_check_range(x, 0, rnd_mode); if (inex) return inex; /* underflow or overflow */ @@ -68,12 +69,12 @@ mpfr_set_si (mpfr_ptr x, long i, mp_rnd_t rnd_mode) rnd_mode, &inex); if (carry) { - mp_exp_t exp = MPFR_EXP(x); + mp_exp_t exp = MPFR_GET_EXP (x); if (exp == __gmpfr_emax) return mpfr_set_overflow(x, rnd_mode, (i < 0 ? -1 : 1)); - MPFR_EXP(x)++; + MPFR_SET_EXP (x, exp + 1); xp[xn] = MPFR_LIMB_HIGHBIT; } } @@ -361,7 +361,7 @@ mpfr_set_str (mpfr_t x, __gmp_const char *str, int base, mp_rnd_t rnd) TMP_FREE(marker); - MPFR_EXP(x) = exp_y + n * BITS_PER_MP_LIMB; + MPFR_SET_EXP (x, exp_y + n * BITS_PER_MP_LIMB); sign_and_flags: MPFR_CLEAR_FLAGS(x); diff --git a/set_str_raw.c b/set_str_raw.c index 9b455920b..b18a390d3 100644 --- a/set_str_raw.c +++ b/set_str_raw.c @@ -125,12 +125,14 @@ mpfr_set_str_raw (mpfr_ptr x, __gmp_const char *str) xp[xsize - k] <<= (BITS_PER_MP_LIMB - j); } - for (; k <= xsize; k++) { xp[xsize - k] = 0; } + for (; k <= xsize; k++) + xp[xsize - k] = 0; - count_leading_zeros(cnt, xp[xsize - 1]); - if (cnt) mpn_lshift(xp, xp, xsize, cnt); + count_leading_zeros(cnt, xp[xsize - 1]); + if (cnt) + mpn_lshift(xp, xp, xsize, cnt); - MPFR_EXP(x) = expn - cnt; + MPFR_SET_EXP (x, expn - cnt); if (MPFR_ISNEG(x) != negative) MPFR_CHANGE_SIGN(x); } @@ -1,6 +1,6 @@ /* mpfr_set_ui -- set a MPFR number from a machine unsigned integer -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. @@ -48,7 +48,8 @@ mpfr_set_ui (mpfr_ptr x, unsigned long i, mp_rnd_t rnd_mode) /* don't forget to put zero in lower limbs */ MPN_ZERO(xp, xn); - MPFR_EXP(x) = nbits = BITS_PER_MP_LIMB - cnt; + nbits = BITS_PER_MP_LIMB - cnt; + MPFR_SET_EXP (x, nbits); inex = mpfr_check_range(x, 0, rnd_mode); if (inex) return inex; /* underflow or overflow */ @@ -61,12 +62,12 @@ mpfr_set_ui (mpfr_ptr x, unsigned long i, mp_rnd_t rnd_mode) rnd_mode, &inex); if (carry) { - mp_exp_t exp = MPFR_EXP(x); + mp_exp_t exp = MPFR_GET_EXP (x); if (exp == __gmpfr_emax) return mpfr_set_overflow(x, rnd_mode, 1); - MPFR_EXP(x)++; + MPFR_SET_EXP (x, exp + 1); xp[xn] = MPFR_LIMB_HIGHBIT; } } @@ -1,6 +1,6 @@ /* mpfr_set_z -- set a floating-point number from a multiple-precision integer -Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc. +Copyright 1999, 2000, 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -178,6 +178,6 @@ mpfr_set_z (mpfr_ptr f, mpz_srcptr z, mp_rnd_t rnd_mode) rnd_mode = GMP_RNDZ; return mpfr_set_underflow(f, rnd_mode, sign_z); } - MPFR_EXP(f) = exp; + MPFR_SET_EXP (f, exp); MPFR_RET(inex); } @@ -1,6 +1,6 @@ /* mpfr_setmax -- maximum representable floating-point number (raw version) -Copyright 2002 Free Software Foundation. +Copyright 2002, 2003 Free Software Foundation. Contributed by the Spaces project, INRIA Lorraine. This file is part of the MPFR Library. @@ -34,7 +34,7 @@ mpfr_setmax (mpfr_ptr x, mp_exp_t e) int sh; mp_limb_t *xp; - MPFR_EXP(x) = e; + MPFR_SET_EXP (x, e); xn = 1 + (MPFR_PREC(x) - 1) / BITS_PER_MP_LIMB; sh = (mp_prec_t) xn * BITS_PER_MP_LIMB - MPFR_PREC(x); xp = MPFR_MANT(x); @@ -1,6 +1,6 @@ /* mpfr_setmin -- minimum representable floating-point number (raw version) -Copyright 2002 Free Software Foundation. +Copyright 2002, 2003 Free Software Foundation. Contributed by the Spaces project, INRIA Lorraine. This file is part of the MPFR Library. @@ -33,7 +33,7 @@ mpfr_setmin (mpfr_ptr x, mp_exp_t e) mp_size_t xn; mp_limb_t *xp; - MPFR_EXP(x) = e; + MPFR_SET_EXP (x, e); xn = (MPFR_PREC(x) - 1) / BITS_PER_MP_LIMB; xp = MPFR_MANT(x); xp[xn] = MPFR_LIMB_HIGHBIT; @@ -1,6 +1,6 @@ /* mpfr_sin -- sine of a floating-point number -Copyright 2001, 2002 Free Software Foundation, Inc. +Copyright 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -51,14 +51,14 @@ mpfr_sin_sign (mpfr_srcptr x) the result is an integer */ mpfr_const_pi (c, GMP_RNDN); /* err <= ulp(c) = 2^(2-m) */ mpfr_div (k, x, c, GMP_RNDN); - MPFR_EXP(k) --; /* x/(2Pi) = 1/2*(x/Pi) */ + MPFR_SET_EXP(k, MPFR_GET_EXP (k) - 1); /* x/(2Pi) = 1/2*(x/Pi) */ mpfr_rint (k, k, GMP_RNDN); if (MPFR_NOTZERO(k)) { - K = MPFR_EXP(k); /* k is an integer, thus K >= 1 */ + K = MPFR_GET_EXP (k); /* k is an integer, thus K >= 1 */ mpfr_mul (k, k, c, GMP_RNDN); /* err <= 2^(K+3-m) */ - MPFR_EXP(k) ++; + MPFR_SET_EXP (k, MPFR_GET_EXP (k) + 1); mpfr_sub (k, x, k, GMP_RNDN); /* err<=2^(4-m)+2^(K+3-m)<=2^(K+4-m) */ y = k; } @@ -74,7 +74,7 @@ mpfr_sin_sign (mpfr_srcptr x) y = k; } } - while (MPFR_IS_ZERO(y) || (MPFR_EXP(y) < K + 5 - (mp_exp_t) m)); + while (MPFR_IS_ZERO (y) || (MPFR_GET_EXP (y) < K + 5 - (mp_exp_t) m)); sign = MPFR_SIGN(y); @@ -105,7 +105,8 @@ mpfr_sin (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) } precy = MPFR_PREC(y); - m = precy + __gmpfr_ceil_log2 ((double) precy) + MAX(0,MPFR_EXP(x)) + 13; + m = precy + __gmpfr_ceil_log2 ((double) precy) + + MAX (0, MPFR_GET_EXP (x)) + 13; sign = mpfr_sin_sign (x); @@ -116,13 +117,13 @@ mpfr_sin (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) mpfr_cos (c, x, GMP_RNDZ); mpfr_mul (c, c, c, GMP_RNDU); mpfr_ui_sub (c, 1, c, GMP_RNDN); - e = 2 + (-MPFR_EXP(c)) / 2; + e = 2 + (- MPFR_GET_EXP (c)) / 2; mpfr_sqrt (c, c, GMP_RNDN); if (sign < 0) mpfr_neg (c, c, GMP_RNDN); /* the absolute error on c is at most 2^(e-m) = 2^(EXP(c)-err) */ - e = MPFR_EXP(c) + m - e; + e = MPFR_GET_EXP (c) + m - e; ok = (e >= 0) && mpfr_can_round (c, e, GMP_RNDN, rnd_mode, precy); if (ok == 0) @@ -1,6 +1,6 @@ /* mpfr_sin_cos -- sine and cosine of a floating-point number -Copyright 2002 Free Software Foundation, Inc. +Copyright 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -48,7 +48,7 @@ mpfr_sin_cos (mpfr_ptr y, mpfr_ptr z, mpfr_srcptr x, mp_rnd_t rnd_mode) } prec = MAX(MPFR_PREC(y), MPFR_PREC(z)); - m = prec + __gmpfr_ceil_log2 ((double) prec) + ABS(MPFR_EXP(x)) + 13; + m = prec + __gmpfr_ceil_log2 ((double) prec) + ABS (MPFR_GET_EXP (x)) + 13; mpfr_init2 (c, m); mpfr_init2 (k, m); @@ -72,13 +72,13 @@ mpfr_sin_cos (mpfr_ptr y, mpfr_ptr z, mpfr_srcptr x, mp_rnd_t rnd_mode) inexact = mpfr_set (z, c, rnd_mode); mpfr_mul (c, c, c, GMP_RNDU); mpfr_ui_sub (c, 1, c, GMP_RNDN); - e = 2 + (-MPFR_EXP(c)) / 2; + e = 2 + (- MPFR_GET_EXP (c)) / 2; mpfr_sqrt (c, c, GMP_RNDN); if (neg) mpfr_neg (c, c, GMP_RNDN); /* the absolute error on c is at most 2^(e-m) = 2^(EXP(c)-err) */ - e = MPFR_EXP(c) + m - e; + e = MPFR_GET_EXP (c) + m - e; ok = (e >= 0) && mpfr_can_round (c, e, GMP_RNDN, rnd_mode, MPFR_PREC(y)); } @@ -1,6 +1,6 @@ /* mpfr_sinh -- hyperbolic sine -Copyright 2001, 2002 Free Software Foundation, Inc. +Copyright 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -116,9 +116,9 @@ mpfr_sinh (mpfr_ptr y, mpfr_srcptr xt, mp_rnd_t rnd_mode) else { /* calculation of the error */ - d = MPFR_EXP(te) - MPFR_EXP(t) + 2; - - /* estimation of the error */ + d = MPFR_GET_EXP (te) - MPFR_GET_EXP (t) + 2; + + /* error estimate */ /* err = Nt-(__gmpfr_ceil_log2(1+pow(2,d)));*/ err = Nt - (MAX(d,0) + 1); } @@ -1,6 +1,6 @@ /* mpfr_sqrt -- square root of a floating-point number -Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc. +Copyright 1999, 2000, 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -92,7 +92,7 @@ mpfr_sqrt (mpfr_ptr r, mpfr_srcptr u, mp_rnd_t rnd_mode) /* Compare the mantissas */ - odd_exp_u = (unsigned int) MPFR_EXP(u) & 1; + odd_exp_u = (unsigned int) MPFR_GET_EXP (u) & 1; MPFR_ASSERTN(MPFR_PREC(r) <= MPFR_INTPREC_MAX - 3); rrsize = (MPFR_PREC(r) + 2 + odd_exp_u) / BITS_PER_MP_LIMB + 1; MPFR_ASSERTN(rrsize <= MP_SIZE_T_MAX/2); @@ -124,9 +124,9 @@ mpfr_sqrt (mpfr_ptr r, mpfr_srcptr u, mp_rnd_t rnd_mode) } } - MPFR_EXP(r) = MPFR_EXP(u) != MPFR_EMAX_MAX - ? (MPFR_EXP(u) + odd_exp_u) / 2 - : (MPFR_EMAX_MAX - 1) / 2 + 1; + MPFR_SET_EXP(r, MPFR_GET_EXP(u) != MPFR_EMAX_MAX + ? (MPFR_GET_EXP(u) + odd_exp_u) / 2 + : (MPFR_EMAX_MAX - 1) / 2 + 1); do { @@ -264,7 +264,7 @@ mpfr_sqrt (mpfr_ptr r, mpfr_srcptr u, mp_rnd_t rnd_mode) /* Is a shift necessary here? Isn't the result 1.0000...? */ mpn_rshift (rp, rp, rrsize, 1); rp[rrsize-1] |= MPFR_LIMB_HIGHBIT; - MPFR_EXP(r)++; + MPFR_SET_EXP (r, MPFR_GET_EXP (r) + 1); } fin: @@ -1,6 +1,6 @@ /* mpfr_sqrt_ui -- square root of a machine integer -Copyright 2000, 2001, 2002 Free Software Foundation, Inc. +Copyright 2000, 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -39,7 +39,7 @@ mpfr_sqrt_ui (mpfr_ptr r, unsigned long u, mp_rnd_t rnd_mode) MPFR_ASSERTN(u == (mp_limb_t) u); count_leading_zeros (cnt, (mp_limb_t) u); *up = (mp_limb_t) u << cnt; - MPFR_EXP(uu) = BITS_PER_MP_LIMB - cnt; + MPFR_SET_EXP (uu, BITS_PER_MP_LIMB - cnt); mpfr_save_emin_emax(); inex = mpfr_sqrt(r, uu, rnd_mode); @@ -1,6 +1,6 @@ /* mpfr_sub -- subtract two floating-point numbers -Copyright 2001, 2002 Free Software Foundation. +Copyright 2001, 2002, 2003 Free Software Foundation. Contributed by the Spaces project, INRIA Lorraine. This file is part of the MPFR Library. @@ -90,22 +90,23 @@ mpfr_sub (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) } else { /* signs differ, it's an addition */ - if (MPFR_EXP(b) < MPFR_EXP(c)) + mp_exp_t eb, ec; + eb = MPFR_GET_EXP (b); + ec = MPFR_GET_EXP (c); + if (eb < ec) { /* exchange rounding modes towards +/- infinity */ int inexact; if (rnd_mode == GMP_RNDU) rnd_mode = GMP_RNDD; else if (rnd_mode == GMP_RNDD) rnd_mode = GMP_RNDU; - inexact = mpfr_add1(a, c, b, rnd_mode, - (mp_exp_unsigned_t) MPFR_EXP(c) - MPFR_EXP(b)); + inexact = mpfr_add1(a, c, b, rnd_mode, (mp_exp_unsigned_t) ec - eb); MPFR_CHANGE_SIGN(a); return -inexact; } else { - return mpfr_add1(a, b, c, rnd_mode, - (mp_exp_unsigned_t) MPFR_EXP(b) - MPFR_EXP(c)); + return mpfr_add1(a, b, c, rnd_mode, (mp_exp_unsigned_t) eb - ec); } } } @@ -1,6 +1,6 @@ /* mpfr_sub1 -- internal function to perform a "real" subtraction -Copyright 2001, 2002 Free Software Foundation. +Copyright 2001, 2002, 2003 Free Software Foundation. Contributed by the Spaces project, INRIA Lorraine. This file is part of the MPFR Library. @@ -74,7 +74,7 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode, if (!sub) MPFR_SET_SAME_SIGN(a, b); - diff_exp = (mp_exp_unsigned_t) MPFR_EXP(b) - MPFR_EXP(c); + 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 */ @@ -393,7 +393,7 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode, mp_exp_t exp_a; cancel -= add_exp; /* still valid as unsigned long */ - exp_a = MPFR_EXP(b) - cancel; + exp_a = MPFR_GET_EXP (b) - cancel; if (exp_a < __gmpfr_emin) { TMP_FREE(marker); @@ -403,19 +403,22 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode, rnd_mode = GMP_RNDZ; return mpfr_set_underflow (a, rnd_mode, MPFR_SIGN(a)); } - MPFR_EXP(a) = exp_a; + MPFR_SET_EXP (a, exp_a); } else /* cancel = 0: MPFR_EXP(a) <- MPFR_EXP(b) + add_exp */ { /* in case cancel = 0, add_exp can still be 1, in case b is just below a power of two, c is very small, prec(a) < prec(b), and rnd=away or nearest */ - if (add_exp && MPFR_EXP(b) == __gmpfr_emax) + mp_exp_t exp_b; + + exp_b = MPFR_GET_EXP (b); + if (add_exp && exp_b == __gmpfr_emax) { TMP_FREE(marker); return mpfr_set_overflow (a, rnd_mode, MPFR_SIGN(a)); } - MPFR_EXP(a) = MPFR_EXP(b) + add_exp; + MPFR_SET_EXP (a, exp_b + add_exp); } TMP_FREE(marker); #ifdef DEBUG diff --git a/sub_one_ulp.c b/sub_one_ulp.c index a14b35867..47b1e4bc8 100644 --- a/sub_one_ulp.c +++ b/sub_one_ulp.c @@ -1,6 +1,6 @@ /* mpfr_sub_one_ulp -- subtract one unit in last place -Copyright 1999, 2001, 2002 Free Software Foundation, Inc. +Copyright 1999, 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -46,7 +46,7 @@ mpfr_sub_one_ulp(mpfr_ptr x, mp_rnd_t rnd_mode) mpn_sub_1 (xp, xp, xn, MP_LIMB_T_ONE << sh); if (xp[xn-1] >> (BITS_PER_MP_LIMB - 1) == 0) { /* was an exact power of two: not normalized any more */ - mp_exp_t exp = MPFR_EXP(x); + mp_exp_t exp = MPFR_GET_EXP (x); /* Note: In case of underflow and rounding to the nearest mode, x won't be changed. Beware of infinite loops! */ if (exp == __gmpfr_emin) @@ -54,7 +54,7 @@ mpfr_sub_one_ulp(mpfr_ptr x, mp_rnd_t rnd_mode) else { mp_size_t i; - MPFR_EXP(x)--; + MPFR_SET_EXP (x, exp - 1); xp[0] = (sh + 1 == BITS_PER_MP_LIMB) ? 0 : MP_LIMB_T_MAX << (sh + 1); for (i = 1; i < xn; i++) xp[i] = MP_LIMB_T_MAX; @@ -1,6 +1,6 @@ /* mpfr_sub_ui -- subtract a floating-point number and a machine integer -Copyright 2000, 2001, 2002 Free Software Foundation, Inc. +Copyright 2000, 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -39,7 +39,7 @@ mpfr_sub_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mp_rnd_t rnd_mode) MPFR_ASSERTN(u == (mp_limb_t) u); count_leading_zeros(cnt, (mp_limb_t) u); *up = (mp_limb_t) u << cnt; - MPFR_EXP(uu) = BITS_PER_MP_LIMB - cnt; + MPFR_SET_EXP (uu, BITS_PER_MP_LIMB - cnt); /* Optimization note: Exponent operations may be removed if mpfr_sub works even when uu is out-of-range. */ @@ -1,6 +1,6 @@ /* mpfr_swap (U, V) -- Swap U and V. -Copyright 2000, 2001 Free Software Foundation, Inc. +Copyright 2000, 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -20,6 +20,7 @@ the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ #include "gmp.h" +#include "gmp-impl.h" #include "mpfr.h" #include "mpfr-impl.h" @@ -41,10 +42,10 @@ mpfr_swap (mpfr_ptr u, mpfr_ptr v) MPFR_SIZE(v) = usize; MPFR_SIZE(u) = vsize; - uexp = MPFR_EXP(u); - vexp = MPFR_EXP(v); - MPFR_EXP(v) = uexp; - MPFR_EXP(u) = vexp; + uexp = MPFR_GET_EXP (u); + vexp = MPFR_GET_EXP (v); + MPFR_SET_EXP (v, uexp); + MPFR_SET_EXP (u, vexp); up = MPFR_MANT(u); vp = MPFR_MANT(v); @@ -1,6 +1,6 @@ /* mpfr_tan -- tangent of a floating-point number -Copyright 2001, 2002 Free Software Foundation, Inc. +Copyright 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -46,7 +46,8 @@ mpfr_tan (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) } precy = MPFR_PREC(y); - m = precy + __gmpfr_ceil_log2 ((double) precy) + ABS(MPFR_EXP(x)) + 13; + m = precy + __gmpfr_ceil_log2 ((double) precy) + + ABS (MPFR_GET_EXP (x)) + 13; mpfr_init2 (s, m); mpfr_init2 (c, m); @@ -1,6 +1,6 @@ /* mpfr_tanh -- hyperbolic tangent -Copyright 2001, 2002 Free Software Foundation. +Copyright 2001, 2002, 2003 Free Software Foundation. This file is part of the MPFR Library. @@ -123,8 +123,8 @@ mpfr_tanh (y, xt, rnd_mode) /* calculation of the error*/ - d = MPFR_EXP(te)-MPFR_EXP(t); - + d = MPFR_GET_EXP (te) - MPFR_GET_EXP (t); + /* estimation of the error */ /*err = Nt-(__gmpfr_ceil_log2(7+pow(2,d+1)));*/ err = Nt-(MAX(d+1,3)+1); @@ -1,6 +1,6 @@ /* mpfr_ui_div -- divide a machine integer by a floating-point number -Copyright 2000, 2001, 2002 Free Software Foundation, Inc. +Copyright 2000, 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -58,8 +58,7 @@ mpfr_ui_div (mpfr_ptr y, unsigned long int u, mpfr_srcptr x, mp_rnd_t rnd_mode) MPFR_ASSERTN(u == (mp_limb_t) u); count_leading_zeros(cnt, (mp_limb_t) u); *up = (mp_limb_t) u << cnt; - MPFR_EXP(uu) = BITS_PER_MP_LIMB-cnt; - + MPFR_SET_EXP (uu, BITS_PER_MP_LIMB - cnt); return mpfr_div (y, uu, x, rnd_mode); } else /* u = 0 */ @@ -1,6 +1,6 @@ /* mpfr_ui_pow -- power of n function n^x -Copyright 2001, 2002 Free Software Foundation, Inc. +Copyright 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -48,7 +48,7 @@ mpfr_ui_pow_is_exact (unsigned long int x, mpfr_srcptr y) /* compute d such that y = c*2^d with c odd integer */ ysize = 1 + (MPFR_PREC(y) - 1) / BITS_PER_MP_LIMB; - d = MPFR_EXP(y) - ysize * BITS_PER_MP_LIMB; + d = MPFR_GET_EXP (y) - ysize * BITS_PER_MP_LIMB; /* since y is not zero, necessarily one of the mantissa limbs is not zero, thus we can simply loop until we find a non zero limb */ yp = MPFR_MANT(y); @@ -128,7 +128,7 @@ mpfr_ui_pow (mpfr_ptr y, unsigned long int n, mpfr_srcptr x, mp_rnd_t rnd_mode) /* General case */ { - /* Declaration of the intermediary variable */ + /* Declaration of the intermediary variable */ mpfr_t t, te, ti; /* Declaration of the size variable */ @@ -148,26 +148,27 @@ mpfr_ui_pow (mpfr_ptr y, unsigned long int n, mpfr_srcptr x, mp_rnd_t rnd_mode) mpfr_init2 (ti, sizeof(unsigned long int) * 8); /* 8 = CHAR_BIT */ mpfr_init2 (te, MPFR_PREC_MIN); - do { - - /* reactualisation of the precision */ - mpfr_set_prec (t, Nt); - mpfr_set_prec (te, Nt); + do + { + /* reactualisation of the precision */ + mpfr_set_prec (t, Nt); + mpfr_set_prec (te, Nt); - /* compute exp(x*ln(n))*/ - mpfr_set_ui (ti, n, GMP_RNDN); /* ti <- n*/ - mpfr_log (t, ti, GMP_RNDU); /* ln(n) */ - mpfr_mul (te, x, t, GMP_RNDU); /* x*ln(n) */ - mpfr_exp (t, te, GMP_RNDN); /* exp(x*ln(n))*/ + /* compute exp(x*ln(n))*/ + mpfr_set_ui (ti, n, GMP_RNDN); /* ti <- n*/ + mpfr_log (t, ti, GMP_RNDU); /* ln(n) */ + mpfr_mul (te, x, t, GMP_RNDU); /* x*ln(n) */ + mpfr_exp (t, te, GMP_RNDN); /* exp(x*ln(n))*/ - /* estimation of the error -- see pow function in algorithms.ps*/ - err = Nt - (MPFR_EXP(te) + 3); + /* error estimate -- see pow function in algorithms.ps */ + err = Nt - (MPFR_GET_EXP (te) + 3); - /* actualisation of the precision */ - Nt += 10; + /* actualisation of the precision */ + Nt += 10; + } + while (inexact && + (err < 0 || !mpfr_can_round (t, err, GMP_RNDN, rnd_mode, Ny))); - } while (inexact && (err < 0 || !mpfr_can_round (t, err, GMP_RNDN, rnd_mode, Ny))); - inexact = mpfr_set (y, t, rnd_mode); mpfr_clear (t); @@ -1,6 +1,6 @@ /* mpfr_ui_sub -- subtract a floating-point number from an integer -Copyright 2000, 2001, 2002 Free Software Foundation, Inc. +Copyright 2000, 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -54,7 +54,7 @@ mpfr_ui_sub (mpfr_ptr y, unsigned long int u, mpfr_srcptr x, mp_rnd_t rnd_mode) MPFR_ASSERTN(u == (mp_limb_t) u); count_leading_zeros (cnt, (mp_limb_t) u); *up = (mp_limb_t) u << cnt; - MPFR_EXP(uu) = BITS_PER_MP_LIMB - cnt; + MPFR_SET_EXP (uu, BITS_PER_MP_LIMB - cnt); return mpfr_sub (y, uu, x, rnd_mode); } else diff --git a/urandomb.c b/urandomb.c index a81ea31ed..2def51763 100644 --- a/urandomb.c +++ b/urandomb.c @@ -3,7 +3,7 @@ using STATE as the random state previously initialized by a call to gmp_randinit(). -Copyright 2000, 2001, 2002 Free Software Foundation, Inc. +Copyright 2000, 2001, 2002, 2003 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -68,6 +68,6 @@ mpfr_urandomb (mpfr_ptr rop, gmp_randstate_t rstate) rp[0] &= ~((MP_LIMB_T_ONE << cnt) - 1); } - MPFR_EXP (rop) = exp; + MPFR_SET_EXP (rop, exp); MPFR_SET_POS (rop); } @@ -189,7 +189,7 @@ mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode) /* Principal loop: we compute, in z_pre, an approximation of Zeta(s), that we send to mpfr_can_round */ mpfr_sub_ui (s1, s, 1, GMP_RNDN); - if (MPFR_EXP(s1) <= -(d-3)/2) + if (MPFR_GET_EXP (s1) <= -(d-3)/2) /* Branch 1: when s-1 is very small, one uses the approximation Zeta(s)=1/(s-1)+gamma, where gamma is Euler's constant */ @@ -356,7 +356,7 @@ mpfr_zeta (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode) precs = mpfr_get_prec (s); /* Precision precs1 needed to represent 1 - s, and s + 2, without any truncation */ - precs1 = precs + 2 + MAX(0, - abs(MPFR_EXP(s))); + precs1 = precs + 2 + MAX (0, - abs (MPFR_GET_EXP (s))); sd = mpfr_get_d (s, GMP_RNDN); /* Precision prec1 is the precision on elementary computations; it ensures a final precision prec1 - add for zeta(s) */ @@ -367,7 +367,7 @@ mpfr_zeta (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode) add = __gmpfr_ceil_log2 (c * c * c * (13.0 + m1)); prec1 = precz + add; /* Note that prec1 is still incremented by 10 at the first entry of loop below */ prec1 = MAX(prec1, precs1); - if (MPFR_SIGN(s) != -1 && MPFR_EXP(s) >= 0) /* Case s >= 1/2 */ + if (MPFR_SIGN (s) != -1 && MPFR_GET_EXP (s) >= 0) /* Case s >= 1/2 */ mpfr_zeta_pos (z, s, rnd_mode); else { |