diff options
| -rw-r--r-- | Makefile.am | 2 | ||||
| -rw-r--r-- | add.c | 13 | ||||
| -rw-r--r-- | add1sp.c | 336 | ||||
| -rw-r--r-- | copysign.c | 6 | ||||
| -rw-r--r-- | div_2ui.c | 2 | ||||
| -rw-r--r-- | mpfr-impl.h | 1 | ||||
| -rw-r--r-- | sqrt.c | 3 | ||||
| -rw-r--r-- | sub.c | 12 | ||||
| -rw-r--r-- | tests/Makefile.am | 2 | ||||
| -rw-r--r-- | tests/tadd1sp.c | 135 | ||||
| -rw-r--r-- | tests/tcopysign.c | 71 | ||||
| -rw-r--r-- | tests/tdim.c | 2 | ||||
| -rw-r--r-- | tests/texp.c | 26 | ||||
| -rw-r--r-- | tests/tminmax.c | 170 | ||||
| -rw-r--r-- | tests/tmul_2exp.c | 33 |
15 files changed, 780 insertions, 34 deletions
diff --git a/Makefile.am b/Makefile.am index 5eb18b750..c6807c517 100644 --- a/Makefile.am +++ b/Makefile.am @@ -7,7 +7,7 @@ include_HEADERS = mpfr.h mpf2mpfr.h lib_LIBRARIES = libmpfr.a -libmpfr_a_SOURCES = mpfr.h mpf2mpfr.h mpfr-impl.h mpfr-test.h log_b2.h exceptions.c save_expo.c extract.c uceil_exp2.c uceil_log2.c ufloor_log2.c add.c add1.c add_one_ulp.c add_ui.c agm.c clear.c cmp.c cmp_abs.c cmp_si.c cmp_ui.c comparisons.c div_2exp.c div_2si.c div_2ui.c div.c div_ui.c dump.c eq.c exp10.c exp2.c exp3.c exp.c frac.c get_d.c get_exp.c get_str.c init.c inp_str.c isinteger.c isinf.c isnan.c isnum.c const_log2.c log.c mul_2exp.c mul_2si.c mul_2ui.c mul.c mul_ui.c neg.c next.c out_str.c const_pi.c pow.c pow_si.c pow_ui.c print_raw.c print_rnd_mode.c random2.c random.c reldiff.c round_prec.c set.c setmax.c setmin.c set_d.c set_dfl_prec.c set_exp.c set_rnd.c set_f.c set_prc_raw.c set_prec.c set_q.c set_si.c set_str.c set_str_raw.c set_ui.c set_z.c sqrt.c sqrt_ui.c sub.c sub1.c sub_one_ulp.c sub_ui.c rint.c ui_div.c ui_sub.c urandomb.c get_z_exp.c swap.c factorial.c cosh.c sinh.c tanh.c acosh.c asinh.c atanh.c atan.c cmp2.c exp_2.c asin.c const_euler.c cos.c sin.c tan.c fma.c hypot.c log1p.c expm1.c log2.c log10.c ui_pow.c ui_pow_ui.c minmax.c dim.c copysign.c gmp_op.c init2.c acos.c sin_cos.c set_nan.c set_inf.c powerof2.c gamma.c set_ld.c get_ld.c cbrt.c volatile.c fits_s.h fits_sshort.c fits_sint.c fits_slong.c fits_u.h fits_ushort.c fits_uint.c fits_ulong.c get_si.c get_ui.c zeta.c cmp_d.c erf.c inits.c inits2.c clears.c sgn.c check.c sub1sp.c version.c mpn_exp.c mpfr-gmp.c mp_clz_tab.c sum.c +libmpfr_a_SOURCES = mpfr.h mpf2mpfr.h mpfr-impl.h mpfr-test.h log_b2.h exceptions.c save_expo.c extract.c uceil_exp2.c uceil_log2.c ufloor_log2.c add.c add1.c add_one_ulp.c add_ui.c agm.c clear.c cmp.c cmp_abs.c cmp_si.c cmp_ui.c comparisons.c div_2exp.c div_2si.c div_2ui.c div.c div_ui.c dump.c eq.c exp10.c exp2.c exp3.c exp.c frac.c get_d.c get_exp.c get_str.c init.c inp_str.c isinteger.c isinf.c isnan.c isnum.c const_log2.c log.c mul_2exp.c mul_2si.c mul_2ui.c mul.c mul_ui.c neg.c next.c out_str.c const_pi.c pow.c pow_si.c pow_ui.c print_raw.c print_rnd_mode.c random2.c random.c reldiff.c round_prec.c set.c setmax.c setmin.c set_d.c set_dfl_prec.c set_exp.c set_rnd.c set_f.c set_prc_raw.c set_prec.c set_q.c set_si.c set_str.c set_str_raw.c set_ui.c set_z.c sqrt.c sqrt_ui.c sub.c sub1.c sub_one_ulp.c sub_ui.c rint.c ui_div.c ui_sub.c urandomb.c get_z_exp.c swap.c factorial.c cosh.c sinh.c tanh.c acosh.c asinh.c atanh.c atan.c cmp2.c exp_2.c asin.c const_euler.c cos.c sin.c tan.c fma.c hypot.c log1p.c expm1.c log2.c log10.c ui_pow.c ui_pow_ui.c minmax.c dim.c copysign.c gmp_op.c init2.c acos.c sin_cos.c set_nan.c set_inf.c powerof2.c gamma.c set_ld.c get_ld.c cbrt.c volatile.c fits_s.h fits_sshort.c fits_sint.c fits_slong.c fits_u.h fits_ushort.c fits_uint.c fits_ulong.c get_si.c get_ui.c zeta.c cmp_d.c erf.c inits.c inits2.c clears.c sgn.c check.c sub1sp.c version.c mpn_exp.c mpfr-gmp.c mp_clz_tab.c sum.c add1sp.c libmpfr_a_LIBADD = @LIBOBJS@ @@ -87,9 +87,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_GET_EXP(b) < MPFR_GET_EXP(c)) - return mpfr_add1(a, c, b, rnd_mode); + if (MPFR_LIKELY(MPFR_PREC(a) == MPFR_PREC(b) + && MPFR_PREC(b) == MPFR_PREC(c))) + if (MPFR_GET_EXP(b) < MPFR_GET_EXP(c)) + return mpfr_add1sp(a, c, b, rnd_mode); + else + return mpfr_add1sp(a, b, c, rnd_mode); else - return mpfr_add1(a, b, c, rnd_mode); + if (MPFR_GET_EXP(b) < MPFR_GET_EXP(c)) + return mpfr_add1(a, c, b, rnd_mode); + else + return mpfr_add1(a, b, c, rnd_mode); } } diff --git a/add1sp.c b/add1sp.c new file mode 100644 index 000000000..da58592e8 --- /dev/null +++ b/add1sp.c @@ -0,0 +1,336 @@ +/* mpfr_add1sp -- internal function to perform a "real" addition + All the op must have the same precision + +Copyright 2004 Free Software Foundation. +Contributed by the Spaces project, INRIA Lorraine. + +This file is part of the MPFR Library. + +The MPFR Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 of the License, or (at your +option) any later version. + +The MPFR Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the MPFR Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include <stdio.h> + +#define MPFR_NEED_LONGLONG_H +#include "mpfr-impl.h" + +/* Define some internal MACROS */ + +/*#define DEBUG*/ +/*#define CHECK_AGAINST_SUB1*/ + +#ifdef DEBUG +# undef DEBUG +# define DEBUG(x) (x) +#else +# define DEBUG(x) /**/ +#endif + +/* compute sign(b) * (|b| + |c|) + Returns 0 iff result is exact, + a negative value when the result is less than the exact value, + a positive value otherwise. +*/ +int +mpfr_add1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) +{ + mp_exp_unsigned_t d; + mp_prec_t p, sh; + mp_size_t n; + mp_limb_t *ap, *cp; + mp_exp_t bx; + mp_limb_t limb; + int inexact; + TMP_DECL(marker); + + TMP_MARK(marker); + + MPFR_ASSERTD(MPFR_PREC(a) == MPFR_PREC(b) && MPFR_PREC(b) == MPFR_PREC(c)); + MPFR_ASSERTD(MPFR_IS_PURE_FP(b) && MPFR_IS_PURE_FP(c)); + MPFR_ASSERTD(MPFR_GET_EXP(b) >= MPFR_GET_EXP(c)); + + /* Read prec and num of limbs */ + p = MPFR_PREC(b); + n = (p+BITS_PER_MP_LIMB-1)/BITS_PER_MP_LIMB; + MPFR_UNSIGNED_MINUS_MODULO(sh, p); + bx = MPFR_GET_EXP(b); + d = (mpfr_uexp_t) (bx - MPFR_GET_EXP(c)); + + DEBUG( printf("New add1sp with diff=%lu\n", d) ); + + if (d == 0) + { + /* d==0 */ + DEBUG( mpfr_print_mant_binary("C= ", MPFR_MANT(c), p) ); + DEBUG( mpfr_print_mant_binary("B= ", MPFR_MANT(b), p) ); + bx++; /* exp + 1 */ + ap = MPFR_MANT(a); + limb = mpn_add_n(ap, MPFR_MANT(b), MPFR_MANT(c), n); + DEBUG( mpfr_print_mant_binary("A= ", ap, p) ); + MPFR_ASSERTD(limb != 0); /* There must be a carry */ + limb = ap[0]; /* Get LSB (In fact, LSW) */ + mpn_rshift(ap, ap, n, 1); /* Shift mantissa A */ + ap[n-1] |= MPFR_LIMB_HIGHBIT; /* Set MSB */ + ap[0] &= ~MPFR_LIMB_MASK(sh); /* Clear LSB bit */ + if ((limb&(MPFR_LIMB_ONE<<sh)) == 0) /* Check exact case */ + { inexact = 0; goto set_exponent; } + /* Zero: Truncate + Nearest: Even Rule => truncate or add 1 + Away: Add 1 */ + if (MPFR_LIKELY(rnd_mode==GMP_RNDN)) + { + if ((ap[0]&(MPFR_LIMB_ONE<<sh))==0) + { inexact = -1; goto set_exponent; } + else + goto add_one_ulp; + } + MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(b)); + if (rnd_mode==GMP_RNDZ) + { inexact = -1; goto set_exponent; } + else + goto add_one_ulp; + } + else if (d >= p) + { + if (d > p) + { + /* d > p : Copy B in A */ + ap = MPFR_MANT(a); + MPN_COPY (ap, MPFR_MANT(b), n); + /* Away: Add 1 + Nearest: Trunc + Zero: Trunc */ + if (MPFR_LIKELY(rnd_mode==GMP_RNDN)) + { inexact = -1; goto set_exponent; } + MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(b)); + if (rnd_mode==GMP_RNDZ) + { inexact = -1; goto set_exponent; } + else + goto add_one_ulp; + } + else + { + /* d==p : Copy B in A */ + ap = MPFR_MANT(a); + MPN_COPY (ap, MPFR_MANT(b), n); + /* Away: Add 1 + Nearest: Even Rule if C is a power of 2, else Add 1 + Zero: Trunc */ + if (MPFR_LIKELY(rnd_mode==GMP_RNDN)) + { + /* Check if C was a power of 2 */ + cp = MPFR_MANT(c); + if (MPFR_UNLIKELY(cp[n-1] == MPFR_LIMB_HIGHBIT)) + { + mp_size_t k = n-1; + do { + k--; + } while (k>=0 && cp[k]==0); + if (MPFR_UNLIKELY(k<0)) + /* Power of 2: Even rule */ + if ((ap[0]&(MPFR_LIMB_ONE<<sh))==0) + { inexact = -1; goto set_exponent; } + } + /* Not a Power of 2 */ + goto add_one_ulp; + } + MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(b)); + if (rnd_mode==GMP_RNDZ) + { inexact = -1; goto set_exponent; } + else + goto add_one_ulp; + } + } + else + { + mp_limb_t mask; + mp_limb_t bcp, bcp1; /* Cp and C'p+1 */ + + /* General case: 1 <= d < p */ + cp = (mp_limb_t*) TMP_ALLOC(n * BYTES_PER_MP_LIMB); + + /* Shift c in temporary allocated place */ + { + mp_exp_unsigned_t dm; + mp_size_t m; + + dm = d % BITS_PER_MP_LIMB; + m = d / BITS_PER_MP_LIMB; + if (MPFR_UNLIKELY(dm == 0)) + { + /* dm = 0 and m > 0: Just copy */ + MPFR_ASSERTD(m!=0); + /* Must use INCR since src & dest may overlap and dest <= src */ + MPN_COPY_INCR(cp, MPFR_MANT(c)+m, n-m); + MPN_ZERO(cp+n-m, m); + } + else if (MPFR_LIKELY(m == 0)) + { + /* dm >=1 and m == 0: just shift */ + MPFR_ASSERTD(dm >= 1); + mpn_rshift(cp, MPFR_MANT(c), n, dm); + } + else + { + /* dm > 0 and m > 0: shift and zero */ + mpn_rshift(cp, MPFR_MANT(c)+m, n-m, dm); + MPN_ZERO(cp+n-m, m); + } + } + + DEBUG( mpfr_print_mant_binary("Before", MPFR_MANT(c), p) ); + DEBUG( mpfr_print_mant_binary("B= ", MPFR_MANT(b), p) ); + DEBUG( mpfr_print_mant_binary("After ", cp, p) ); + + /* Compute bcp=Cp and bcp1=C'p+1 */ + if (MPFR_LIKELY(sh)) + { + /* Try to compute them from C' rather than C */ + bcp = (cp[0] & (MPFR_LIMB_ONE<<(sh-1))) ; + if (MPFR_LIKELY(cp[0]&MPFR_LIMB_MASK(sh-1))) + bcp1 = 1; + else + { + /* We can't compute C'p+1 from C'. Compute it from C */ + mp_limb_t *tp = MPFR_MANT(c); + /* Start from bit x=p-d+sh in mantissa C + (+sh since we have already looked sh bits in C'!) */ + mpfr_prec_t x = p-d+sh-1; + if (MPFR_UNLIKELY(x>p)) + /* We are already looked at all the bits of c, so C'p+1 = 0*/ + bcp1 = 0; + else + { + mp_size_t kx = n-1 - (x / BITS_PER_MP_LIMB); + mpfr_prec_t sx = BITS_PER_MP_LIMB-1-(x%BITS_PER_MP_LIMB); + DEBUG( printf("(First) x=%lu Kx=%ld Sx=%lu\n", x, kx, sx)); + /* Looks at the last bits of limb kx (if sx=0 does nothing)*/ + if (tp[kx] & MPFR_LIMB_MASK(sx)) + bcp1 = -1; + else + { + /*kx += (sx==0);*/ + /*If sx==0, tp[kx] hasn't been checked*/ + do { + kx--; + } while (kx>=0 && tp[kx]==0); + bcp1 = (kx >= 0); + } + } + } + } + else /* sh == 0 */ + { + /* Compute Cp and C'p+1 from C with sh=0 */ + mp_limb_t *tp = MPFR_MANT(c); + /* Start from bit x=p-d in mantissa C */ + mpfr_prec_t x = p-d; + mp_size_t kx = n-1 - (x / BITS_PER_MP_LIMB); + mpfr_prec_t sx = BITS_PER_MP_LIMB-1-(x%BITS_PER_MP_LIMB); + MPFR_ASSERTD(p >= d); + bcp = tp[kx] & (MPFR_LIMB_ONE<<sx); + /* Looks at the last bits of limb kx (If sx=0, does nothing)*/ + if (tp[kx]&MPFR_LIMB_MASK(sx)) + bcp1 = 1; + else + { + do { + kx--; + } while (kx>=0 && tp[kx]==0); + bcp1 = (kx>=0); + } + } + DEBUG( printf("sh=%lu Cp=%lu C'p+1=%lu\n", sh, bcp, bcp1) ); + + /* Clean shifted C' */ + mask = ~MPFR_LIMB_MASK(sh); + cp[0] &= mask; + + /* Add the mantissa c from b in a */ + ap = MPFR_MANT(a); + limb = mpn_add_n (ap, MPFR_MANT(b), cp, n); + DEBUG( mpfr_print_mant_binary("Add= ", ap, p) ); + + /* Check for overflow */ + if (limb) + { + limb = ap[0] & (MPFR_LIMB_ONE<<sh); /* Get LSB */ + mpn_rshift(ap, ap, n, 1); /* Shift mantissa*/ + bx++; /* Fix exponent */ + ap[n-1] |= MPFR_LIMB_HIGHBIT; /* Set MSB */ + ap[0] &= mask; /* Clear LSB bit */ + bcp1 |= bcp; /* Recompute C'p+1 */ + bcp = limb; /* Recompute Cp */ + DEBUG( printf("(Overflow) Cp=%lu C'p+1=%lu\n", bcp, bcp1) ); + DEBUG( mpfr_print_mant_binary("Add= ", ap, p) ); + } + + /* Round: + Zero: Truncate but could be exact. + Away: Add 1 if Cp or C'p+1 !=0 + Nearest: Truncate but could be exact if Cp==0 + Add 1 if C'p+1 !=0, + Even rule else */ + if (MPFR_LIKELY(rnd_mode == GMP_RNDN)) + { + if (bcp == 0) + { inexact = MPFR_LIKELY(bcp1) ? -1 : 0; goto set_exponent; } + else if (MPFR_UNLIKELY(bcp1==0) && (ap[0]&(MPFR_LIMB_ONE<<sh))==0) + { inexact = -1; goto set_exponent; } + else + goto add_one_ulp; + } + MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(b)); + if (rnd_mode == GMP_RNDZ) + { + inexact = MPFR_LIKELY(bcp || bcp1) ? -1 : 0; + goto set_exponent; + } + else + { + if (MPFR_UNLIKELY(bcp==0 && bcp1==0)) + { inexact = 0; goto set_exponent; } + else + goto add_one_ulp; + } + } + MPFR_ASSERTN(0); + + add_one_ulp: + /* add one unit in last place to a */ + DEBUG( printf("AddOneUlp\n") ); + if (MPFR_UNLIKELY( mpn_add_1(ap, ap, n, MPFR_LIMB_ONE<<sh) )) + { + /* Case 100000x0 = 0x1111x1 + 1*/ + DEBUG( printf("Pow of 2\n") ); + bx++; + ap[n-1] = MPFR_LIMB_HIGHBIT; + } + inexact = 1; + + set_exponent: + if (MPFR_UNLIKELY(bx >= __gmpfr_emax)) /* Check for overflow */ + { + DEBUG( printf("Overflow\n") ); + TMP_FREE(marker); + MPFR_SET_SAME_SIGN(a,b); + return mpfr_set_overflow(a, rnd_mode, MPFR_SIGN(a)); + } + MPFR_SET_EXP (a, bx); + MPFR_SET_SAME_SIGN(a,b); + + TMP_FREE(marker); + return inexact*MPFR_INT_SIGN(a); +} diff --git a/copysign.c b/copysign.c index ec96fac7c..a763057f0 100644 --- a/copysign.c +++ b/copysign.c @@ -21,9 +21,9 @@ MA 02111-1307, USA. */ #include "mpfr-impl.h" - /* The computation of z with magnitude of x and sign of y - - z = sign(y) * abs(x) + /* + The computation of z with magnitude of x and sign of y + z = sign(y) * abs(x) */ int @@ -22,7 +22,7 @@ MA 02111-1307, USA. */ #include "mpfr-impl.h" int -mpfr_div_2ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int n, mp_rnd_t rnd_mode) +mpfr_div_2ui (mpfr_ptr y, mpfr_srcptr x, unsigned long n, mp_rnd_t rnd_mode) { int inexact; diff --git a/mpfr-impl.h b/mpfr-impl.h index fc0f9760d..a4f7513d4 100644 --- a/mpfr-impl.h +++ b/mpfr-impl.h @@ -474,6 +474,7 @@ void mpfr_restore_emin_emax _MPFR_PROTO ((void)); int mpfr_add1 _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t)); int mpfr_sub1 _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t)); +int mpfr_add1sp _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t)); int mpfr_sub1sp _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t)); int mpfr_can_round_raw _MPFR_PROTO ((mp_limb_t *, mp_size_t, int, mp_exp_t, mp_rnd_t, mp_rnd_t, mp_prec_t)); @@ -53,8 +53,9 @@ mpfr_sqrt (mpfr_ptr r, mpfr_srcptr u, mp_rnd_t rnd_mode) MPFR_SET_ZERO(r); MPFR_RET(0); /* zero is exact */ } - else if (MPFR_IS_INF(u)) + else { + MPFR_ASSERTD(MPFR_IS_INF(u)); /* sqrt(-Inf) = NAN */ if (MPFR_IS_NEG(u)) { @@ -95,13 +95,21 @@ mpfr_sub (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) rnd_mode = GMP_RNDD; else if (rnd_mode == GMP_RNDD) rnd_mode = GMP_RNDU; - inexact = mpfr_add1(a, c, b, rnd_mode); + if (MPFR_LIKELY(MPFR_PREC(a) == MPFR_PREC(b) + && MPFR_PREC(b) == MPFR_PREC(c))) + inexact = mpfr_add1sp(a, c, b, rnd_mode); + else + inexact = mpfr_add1(a, c, b, rnd_mode); MPFR_CHANGE_SIGN(a); return -inexact; } else { - return mpfr_add1(a, b, c, rnd_mode); + if (MPFR_LIKELY(MPFR_PREC(a) == MPFR_PREC(b) + && MPFR_PREC(b) == MPFR_PREC(c))) + return mpfr_add1sp(a, b, c, rnd_mode); + else + return mpfr_add1(a, b, c, rnd_mode); } } } diff --git a/tests/Makefile.am b/tests/Makefile.am index bb5a722d0..38eb12033 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -1,6 +1,6 @@ AUTOMAKE_OPTIONS = gnu -check_PROGRAMS = tversion tinits tsgn tcheck tisnan texceptions tset tabs tset_d tset_f tset_q tset_si tset_str tset_z tset_ld tadd tsub tmul tdiv tcmp tcmp2 tcan_round tround_prec tswap reuse tsqrt tnext tcomparisons teq tadd_ui tsub_ui tcmp_ui tdiv_ui tmul_ui tsqrt_ui tui_div tui_sub tcmp_d tmul_2exp tfma tfrac tget_str tout_str tget_d tget_d_2exp tconst_log2 tconst_pi tconst_euler trandom ttrunc trint texp texp2 texpm1 tlog tlog2 tlog10 tlog1p tpow tui_pow tpow3 tcosh tsinh ttanh tacosh tasinh tatanh thyperbolic tasin tacos tcos tatan tsin ttan tsin_cos tagm thypot tfactorial tgamma terf tcbrt mpf_compat mpfr_compat tzeta tsub1sp tgmpop tsum tdim +check_PROGRAMS = tversion tinits tsgn tcheck tisnan texceptions tset tabs tset_d tset_f tset_q tset_si tset_str tset_z tset_ld tadd tsub tmul tdiv tcmp tcmp2 tcan_round tround_prec tswap reuse tsqrt tnext tcomparisons teq tadd_ui tsub_ui tcmp_ui tdiv_ui tmul_ui tsqrt_ui tui_div tui_sub tcmp_d tmul_2exp tfma tfrac tget_str tout_str tget_d tget_d_2exp tconst_log2 tconst_pi tconst_euler trandom ttrunc trint texp texp2 texpm1 tlog tlog2 tlog10 tlog1p tpow tui_pow tpow3 tcosh tsinh ttanh tacosh tasinh tatanh thyperbolic tasin tacos tcos tatan tsin ttan tsin_cos tagm thypot tfactorial tgamma terf tcbrt mpf_compat mpfr_compat tzeta tsub1sp tadd1sp tgmpop tsum tdim tminmax tcopysign EXTRA_DIST = tgeneric.c mpf_compat.h diff --git a/tests/tadd1sp.c b/tests/tadd1sp.c new file mode 100644 index 000000000..0782ff9f0 --- /dev/null +++ b/tests/tadd1sp.c @@ -0,0 +1,135 @@ +/* Test file for mpfr_add1sp. + +Copyright 2004 Free Software Foundation. + +This file is part of the MPFR Library. + +The MPFR Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 of the License, or (at your +option) any later version. + +The MPFR Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the MPFR Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include <stdio.h> +#include <stdlib.h> + +#include "mpfr-test.h" + +void check_special(void); +void check_random(mpfr_prec_t p); + +int main(void) +{ + mpfr_prec_t p; + + tests_start_mpfr (); + + check_special (); + for(p = 2 ; p < 200 ; p++) + check_random (p); + + tests_end_mpfr (); + return 0; +} + +#define STD_ERROR \ + {\ + printf("ERROR: for %s and p=%lu and i=%d:\nB=",\ + mpfr_print_rnd_mode(r), p, i);\ + mpfr_print_binary(b);\ + printf("\nC="); mpfr_print_binary(c);\ + printf("\nadd1 : "); mpfr_print_binary(a1);\ + printf("\nadd1sp: "); mpfr_print_binary(a2);\ + putchar('\n');\ + abort();\ + } + +#define STD_ERROR2 \ + {\ + printf("ERROR: Wrong inexact flag for %s and p=%lu and i=%d:\nB=",\ + mpfr_print_rnd_mode(r), p, i);\ + mpfr_print_binary(b);\ + printf("\nC="); mpfr_print_binary(c);\ + printf("\nA="); mpfr_print_binary(a1);\ + printf("\nAdd1: %d. Add1sp: %d\n", \ + inexact1, inexact2); \ + abort();\ + } + +#define SET_PREC(_p) \ + { p = _p; \ + mpfr_set_prec(a1, _p); mpfr_set_prec(a2, _p); \ + mpfr_set_prec(b, _p); mpfr_set_prec(c, _p); \ + } + + +void check_random(mp_prec_t p) +{ + mpfr_t a1,b,c,a2; + mp_rnd_t r; + int i, inexact1, inexact2; + + mpfr_inits2(p, a1,b,c,a2, NULL); + + for(i = 0 ; i < 500 ; i++) + { + mpfr_random(b); + mpfr_random(c); + if (MPFR_GET_EXP(b) < MPFR_GET_EXP(c)) + mpfr_swap(b, c); + if (MPFR_IS_PURE_FP(b) && MPFR_IS_PURE_FP(c)) + for(r = 0 ; r < GMP_RND_MAX ; r++) + { + inexact1 = mpfr_add1(a1, b, c, r); + inexact2 = mpfr_add1sp(a2, b, c, r); + if (mpfr_cmp(a1, a2)) + STD_ERROR; + if (inexact1 != inexact2) + STD_ERROR2; + } + } + + mpfr_clears(a1,a2,b,c,NULL); +} + +void check_special(void) +{ + mpfr_t a1,a2,b,c; + mp_rnd_t r; + mpfr_prec_t p; + int i = -1, inexact1, inexact2; + + mpfr_inits(a1,a2,b,c,NULL); + + for(r = 0 ; r < GMP_RND_MAX ; r++) + { + SET_PREC(53); + mpfr_set_str1 (b, "1@100"); + mpfr_set_str1 (c, "1@1"); + inexact1 = mpfr_add1(a1, b, c, r); + inexact2 = mpfr_add1sp(a2, b, c, r); + if (mpfr_cmp(a1, a2)) + STD_ERROR; + if (inexact1 != inexact2) + STD_ERROR2; + mpfr_set_str_binary (b, "1E53"); + mpfr_set_str_binary (c, "1E0"); + inexact1 = mpfr_add1(a1, b, c, r); + inexact2 = mpfr_add1sp(a2, b, c, r); + if (mpfr_cmp(a1, a2)) + STD_ERROR; + if (inexact1 != inexact2) + STD_ERROR2; + } + + mpfr_clears(a1,a2,b,c,NULL); +} diff --git a/tests/tcopysign.c b/tests/tcopysign.c new file mode 100644 index 000000000..7c7799d25 --- /dev/null +++ b/tests/tcopysign.c @@ -0,0 +1,71 @@ +/* Test file for mpfr_copysign. + +Copyright 2004 Free Software Foundation, Inc. + +This file is part of the MPFR Library. + +The MPFR Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 of the License, or (at your +option) any later version. + +The MPFR Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the MPFR Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include <stdio.h> +#include <stdlib.h> + +#include "mpfr-test.h" + +int +main (void) +{ + mpfr_t x, y, z; + + tests_start_mpfr (); + + mpfr_init (x); + mpfr_init (y); + mpfr_init (z); + + /* case y=NaN */ + mpfr_set_nan (y); + mpfr_set_ui (x, 1250, GMP_RNDN); + mpfr_copysign (z, x, y, GMP_RNDN); + if (!mpfr_nan_p (z)) + { + printf ("Error in mpfr_copysign (NaN)\n"); + exit (1); + } + /* case y!=NaN */ + mpfr_set_ui (y, 123, GMP_RNDN); + mpfr_set_ui (x, 1250, GMP_RNDN); + mpfr_copysign (z, x, y, GMP_RNDN); + if (mpfr_cmp_ui (z, 1250)) + { + printf ("Error in mpfr_copysign (1250)\n"); + exit (1); + } + mpfr_set_si (y, -17, GMP_RNDN); + mpfr_set_ui (x, 42, GMP_RNDN); + mpfr_copysign (z, x, y, GMP_RNDN); + if (mpfr_cmp_si (z, -42)) + { + printf ("Error in mpfr_copysign (-42)\n"); + exit (1); + } + + mpfr_clear (x); + mpfr_clear (y); + mpfr_clear (z); + + tests_end_mpfr (); + return 0; +} diff --git a/tests/tdim.c b/tests/tdim.c index cccca8aeb..8749d0c7a 100644 --- a/tests/tdim.c +++ b/tests/tdim.c @@ -71,7 +71,7 @@ main (void) mpfr_dim (z, x, y, GMP_RNDN); if (mpfr_cmp_ui (z, 0) || mpfr_sgn (z) < 0) { - printf ("Error in mpfr_dim (+Inf, 0)\n"); + printf ("Error in mpfr_dim (+Inf, +Inf)\n"); exit (1); } diff --git a/tests/texp.c b/tests/texp.c index ba2f9c193..d574a8c55 100644 --- a/tests/texp.c +++ b/tests/texp.c @@ -247,6 +247,30 @@ check_inexact (void) mpfr_clear (y); } +static void +check_exp10(void) +{ + mpfr_t x; + int inexact; + + mpfr_init2 (x, 200); + mpfr_set_ui(x, 4, GMP_RNDN); + + inexact = mpfr_exp10 (x, x, GMP_RNDN); + if (mpfr_cmp_ui(x, 10*10*10*10)) + { + printf ("exp10: Wrong returned value\n"); + exit (1); + } + if (inexact != 0) + { + printf ("exp10: Wrong inexact flag\n"); + exit (1); + } + + mpfr_clear (x); +} + int main (int argc, char *argv[]) { @@ -302,7 +326,7 @@ main (int argc, char *argv[]) check3("5.30015757134837031117e+02", GMP_RNDD, "1.5237672861171573939e230"); check3("5.16239362447650933063e+02", GMP_RNDZ, "1.5845518406744492105e224"); check3("6.00812634798592370977e-01", GMP_RNDN, "1.823600119339019443"); - + check_exp10 (); done: tests_end_mpfr (); return 0; diff --git a/tests/tminmax.c b/tests/tminmax.c new file mode 100644 index 000000000..2592c420e --- /dev/null +++ b/tests/tminmax.c @@ -0,0 +1,170 @@ +/* Test file for mpfr_min & mpfr_max. + +Copyright 2004 Free Software Foundation, Inc. + +This file is part of the MPFR Library. + +The MPFR Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 of the License, or (at your +option) any later version. + +The MPFR Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the MPFR Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include <stdio.h> +#include <stdlib.h> + +#include "mpfr-test.h" + +int +main (void) +{ + mpfr_t x, y, z; + + tests_start_mpfr (); + + mpfr_init (x); + mpfr_init (y); + mpfr_init (z); + + /* case x=NaN && y=NAN */ + mpfr_set_nan (x); + mpfr_set_nan (y); + mpfr_min (z, x, y, GMP_RNDN); + if (!mpfr_nan_p (z)) + { + printf ("Error in mpfr_min (NaN, NaN)\n"); + exit (1); + } + mpfr_max (z, x, y, GMP_RNDN); + if (!mpfr_nan_p (z)) + { + printf ("Error in mpfr_max (NaN, NaN)\n"); + exit (1); + } + /* case x=NaN */ + mpfr_set_nan (x); + mpfr_set_ui (y, 0, GMP_RNDN); + mpfr_min (z, x, y, GMP_RNDN); + if (mpfr_cmp_ui (z, 0)) + { + printf ("Error in mpfr_min (NaN, 0)\n"); + exit (1); + } + mpfr_min (z, y, x, GMP_RNDN); + if (mpfr_cmp_ui (z, 0)) + { + printf ("Error in mpfr_min (0, NaN)\n"); + exit (1); + } + mpfr_max (z, x, y, GMP_RNDN); + if (mpfr_cmp_ui (z, 0)) + { + printf ("Error in mpfr_max (NaN, 0)\n"); + exit (1); + } + mpfr_max (z, y, x, GMP_RNDN); + if (mpfr_cmp_ui (z, 0)) + { + printf ("Error in mpfr_max (0, NaN)\n"); + exit (1); + } + /* Case x=0+ and x=0- */ + mpfr_set_ui (x, 0, GMP_RNDN); + mpfr_set_ui (y, 0, GMP_RNDN); MPFR_SET_NEG(y); + mpfr_max (z, x, y, GMP_RNDN); + if (!MPFR_IS_ZERO(z) || MPFR_IS_NEG(z)) + { + printf ("Error in mpfr_max (0+, 0-)\n"); + exit (1); + } + mpfr_min (z, x, y, GMP_RNDN); + if (!MPFR_IS_ZERO(z) || MPFR_IS_POS(z)) + { + printf ("Error in mpfr_min (0+, 0-)\n"); + exit (1); + } + /* Case x=0- and y=0+ */ + mpfr_set_ui (x, 0, GMP_RNDN); MPFR_SET_NEG(x); + mpfr_set_ui (y, 0, GMP_RNDN); + mpfr_max (z, x, y, GMP_RNDN); + if (!MPFR_IS_ZERO(z) || MPFR_IS_NEG(z)) + { + printf ("Error in mpfr_max (0+, 0-)\n"); + exit (1); + } + mpfr_min (z, x, y, GMP_RNDN); + if (!MPFR_IS_ZERO(z) || MPFR_IS_POS(z)) + { + printf ("Error in mpfr_min (0+, 0-)\n"); + exit (1); + } + + /* case x=+Inf */ + mpfr_set_inf (x, 1); + mpfr_set_si (y, -12, GMP_RNDN); + mpfr_min (z, x, y, GMP_RNDN); + if ( mpfr_cmp_si (z, -12) ) + { + printf ("Error in mpfr_min (+Inf, -12)\n"); + exit (1); + } + mpfr_max (z, x, y, GMP_RNDN); + if ( !MPFR_IS_INF(z) || MPFR_IS_NEG(z) ) + { + printf ("Error in mpfr_max (+Inf, 12)\n"); + exit (1); + } + /* case x=-Inf */ + mpfr_set_inf (x, -1); + mpfr_set_ui (y, 12, GMP_RNDN); + mpfr_max (z, x, y, GMP_RNDN); + if ( mpfr_cmp_ui (z, 12) ) + { + printf ("Error in mpfr_max (-Inf, 12)\n"); + exit (1); + } + mpfr_min (z, x, y, GMP_RNDN); + if ( !MPFR_IS_INF(z) || MPFR_IS_POS(z) ) + { + printf ("Error in mpfr_min (-Inf, 12)\n"); + exit (1); + } + + /* case x=17 and y=42 */ + mpfr_set_ui (x, 17, GMP_RNDN); + mpfr_set_ui (y, 42, GMP_RNDN); + mpfr_max (z, x, y, GMP_RNDN); + if ( mpfr_cmp_ui (z, 42) ) + { + printf ("Error in mpfr_max (17, 42)\n"); + exit (1); + } + mpfr_max (z, y, x, GMP_RNDN); + if ( mpfr_cmp_ui (z, 42) ) + { + printf ("Error in mpfr_max (42, 17)\n"); + exit (1); + } + mpfr_min (z, y, x, GMP_RNDN); + if ( mpfr_cmp_ui (z, 17) ) + { + printf ("Error in mpfr_min (42, 17)\n"); + exit (1); + } + + mpfr_clear (x); + mpfr_clear (y); + mpfr_clear (z); + + tests_end_mpfr (); + return 0; +} diff --git a/tests/tmul_2exp.c b/tests/tmul_2exp.c index f3229a8b4..ebe6d49f5 100644 --- a/tests/tmul_2exp.c +++ b/tests/tmul_2exp.c @@ -83,27 +83,20 @@ main (int argc, char *argv[]) exit(-1); } } - - /* for (k = 0; k < 100000; k++) + mpfr_set_ui(w, 1, GMP_RNDN); + mpfr_mul_2ui(w, w, (unsigned long) MPFR_EXP_MAX+12, GMP_RNDN); + if (!mpfr_inf_p(w)) { - x = DBL_RAND (); - mpfr_set_d (w, x, 0); - mpfr_mul_2exp (w, w, 10, GMP_RNDZ); - if (x != (z = mpfr_get_d1 (w)/1024)) - { - printf ("%f != %f\n", x, z); - exit (1); - } - - mpfr_set_d(w, x, 0); - mpfr_div_2exp(w, w, 10, GMP_RNDZ); - if (x != (z = mpfr_get_d1 (w)*1024)) - { - printf ("%f != %f\n", x, z); - mpfr_clear (w); - exit (1); - } - }*/ + printf("Overflow error!\n"); + exit(1); + } + mpfr_set_ui(w, 1, GMP_RNDN); + mpfr_div_2ui(w, w, (unsigned long) MPFR_EXP_MAX+12, GMP_RNDN); + if (mpfr_cmp_ui(w, 0)) + { + printf("Underflow error!\n"); + exit(1); + } mpfr_clears(w,z,NULL); |