Annotation of OpenXM_contrib/gmp/mpfr/round_prec.c, Revision 1.1
1.1 ! ohara 1: /* mpfr_round_raw_generic, mpfr_round_raw2, mpfr_round_raw, mpfr_round_prec,
! 2: mpfr_can_round, mpfr_can_round_raw -- various rounding functions
! 3:
! 4: Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
! 5:
! 6: This file is part of the MPFR Library.
! 7:
! 8: The MPFR Library is free software; you can redistribute it and/or modify
! 9: it under the terms of the GNU Lesser General Public License as published by
! 10: the Free Software Foundation; either version 2.1 of the License, or (at your
! 11: option) any later version.
! 12:
! 13: The MPFR Library is distributed in the hope that it will be useful, but
! 14: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
! 15: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
! 16: License for more details.
! 17:
! 18: You should have received a copy of the GNU Lesser General Public License
! 19: along with the MPFR Library; see the file COPYING.LIB. If not, write to
! 20: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
! 21: MA 02111-1307, USA. */
! 22:
! 23: #include "gmp.h"
! 24: #include "gmp-impl.h"
! 25: #include "mpfr.h"
! 26: #include "mpfr-impl.h"
! 27:
! 28: #if (BITS_PER_MP_LIMB & (BITS_PER_MP_LIMB - 1))
! 29: #error "BITS_PER_MP_LIMB must be a power of 2"
! 30: #endif
! 31:
! 32: /*
! 33: * If flag = 0, puts in y the value of xp (with precision xprec and
! 34: * sign 1 if negative=0, -1 otherwise) rounded to precision yprec and
! 35: * direction rnd_mode. Supposes x is not zero nor NaN nor +/- Infinity
! 36: * (i.e. *xp != 0). If inexp != NULL, computes the inexact flag of the
! 37: * rounding.
! 38: *
! 39: * In case of even rounding when rnd = GMP_RNDN, returns 2 or -2.
! 40: *
! 41: * If flag = 1, just returns whether one should add 1 or not for rounding.
! 42: *
! 43: * Note: yprec may be < MPFR_PREC_MIN; in particular, it may be equal
! 44: * to 1. In this case, the even rounding is done away from 0, which is
! 45: * a natural generalization. Indeed, a number with 1-bit precision can
! 46: * be seen as a denormalized number with more precision.
! 47: */
! 48:
! 49: int
! 50: mpfr_round_raw_generic(mp_limb_t *yp, mp_limb_t *xp, mp_prec_t xprec,
! 51: int neg, mp_prec_t yprec, mp_rnd_t rnd_mode,
! 52: int *inexp, int flag)
! 53: {
! 54: mp_size_t xsize, nw;
! 55: mp_limb_t himask, lomask;
! 56: int rw, carry = 0;
! 57:
! 58: xsize = (xprec-1)/BITS_PER_MP_LIMB + 1;
! 59:
! 60: nw = yprec / BITS_PER_MP_LIMB;
! 61: rw = yprec & (BITS_PER_MP_LIMB - 1);
! 62:
! 63: if (flag && !inexp && (rnd_mode==GMP_RNDZ || xprec <= yprec
! 64: || (rnd_mode==GMP_RNDU && neg)
! 65: || (rnd_mode==GMP_RNDD && neg==0)))
! 66: return 0;
! 67:
! 68: if (rw)
! 69: {
! 70: nw++;
! 71: lomask = ((MP_LIMB_T_ONE << (BITS_PER_MP_LIMB - rw)) - MP_LIMB_T_ONE);
! 72: himask = ~lomask;
! 73: }
! 74: else
! 75: {
! 76: lomask = -1;
! 77: himask = -1;
! 78: }
! 79: MPFR_ASSERTN(nw >= 1);
! 80:
! 81: if (xprec <= yprec)
! 82: { /* No rounding is necessary. */
! 83: /* if yp=xp, maybe an overlap: MPN_COPY_DECR is ok when src <= dst */
! 84: MPFR_ASSERTN(nw >= xsize);
! 85: if (inexp)
! 86: *inexp = 0;
! 87: if (flag)
! 88: return 0;
! 89: MPN_COPY_DECR(yp + (nw - xsize), xp, xsize);
! 90: MPN_ZERO(yp, nw - xsize);
! 91: }
! 92: else
! 93: {
! 94: mp_limb_t sb;
! 95:
! 96: if ((rnd_mode == GMP_RNDU && neg) ||
! 97: (rnd_mode == GMP_RNDD && !neg))
! 98: rnd_mode = GMP_RNDZ;
! 99:
! 100: if (inexp || rnd_mode != GMP_RNDZ)
! 101: {
! 102: mp_size_t k;
! 103:
! 104: k = xsize - nw;
! 105: if (!rw)
! 106: k--;
! 107: MPFR_ASSERTN(k >= 0);
! 108: sb = xp[k] & lomask; /* First non-significant bits */
! 109: if (rnd_mode == GMP_RNDN)
! 110: {
! 111: mp_limb_t rbmask = MP_LIMB_T_ONE << (BITS_PER_MP_LIMB - rw - 1);
! 112: if (sb & rbmask) /* rounding bit */
! 113: sb &= ~rbmask; /* it is 1, clear it */
! 114: else
! 115: rnd_mode = GMP_RNDZ; /* it is 0, behave like rounding to 0 */
! 116: }
! 117: while (sb == 0 && k > 0)
! 118: sb = xp[--k];
! 119: if (rnd_mode == GMP_RNDN)
! 120: { /* rounding to nearest, with rounding bit = 1 */
! 121: if (sb == 0) /* Even rounding. */
! 122: {
! 123: sb = xp[xsize - nw] & (himask ^ (himask << 1));
! 124: if (inexp)
! 125: *inexp = ((neg != 0) ^ (sb != 0))
! 126: ? MPFR_EVEN_INEX : -MPFR_EVEN_INEX;
! 127: }
! 128: else /* sb != 0 */
! 129: {
! 130: if (inexp)
! 131: *inexp = (neg == 0) ? 1 : -1;
! 132: }
! 133: }
! 134: else if (inexp)
! 135: *inexp = sb == 0 ? 0
! 136: : (((neg != 0) ^ (rnd_mode != GMP_RNDZ)) ? 1 : -1);
! 137: }
! 138: else
! 139: sb = 0;
! 140:
! 141: if (flag)
! 142: return sb != 0 && rnd_mode != GMP_RNDZ;
! 143:
! 144: if (sb != 0 && rnd_mode != GMP_RNDZ)
! 145: carry = mpn_add_1(yp, xp + xsize - nw, nw,
! 146: rw ? MP_LIMB_T_ONE << (BITS_PER_MP_LIMB - rw) : 1);
! 147: else
! 148: MPN_COPY_INCR(yp, xp + xsize - nw, nw);
! 149:
! 150: yp[0] &= himask;
! 151: }
! 152:
! 153: return carry;
! 154: }
! 155:
! 156: int
! 157: mpfr_round_prec (mpfr_ptr x, mp_rnd_t rnd_mode, mp_prec_t prec)
! 158: {
! 159: mp_limb_t *tmp, *xp;
! 160: int carry, inexact;
! 161: mp_prec_t nw;
! 162: TMP_DECL(marker);
! 163:
! 164: MPFR_ASSERTN(prec >= MPFR_PREC_MIN && prec <= MPFR_PREC_MAX);
! 165:
! 166: nw = 1 + (prec - 1) / BITS_PER_MP_LIMB; /* needed allocated limbs */
! 167:
! 168: /* check if x has enough allocated space for the mantissa */
! 169: if (nw > MPFR_ABSSIZE(x))
! 170: {
! 171: MPFR_MANT(x) = (mp_ptr) (*__gmp_reallocate_func)
! 172: (MPFR_MANT(x), (size_t) MPFR_ABSSIZE(x) * BYTES_PER_MP_LIMB,
! 173: (size_t) nw * BYTES_PER_MP_LIMB);
! 174: MPFR_SET_ABSSIZE(x, nw); /* new number of allocated limbs */
! 175: }
! 176:
! 177: if (MPFR_IS_NAN(x))
! 178: MPFR_RET_NAN;
! 179:
! 180: if (MPFR_IS_INF(x))
! 181: return 0; /* infinity is exact */
! 182:
! 183: /* x is a real number */
! 184:
! 185: TMP_MARK(marker);
! 186: tmp = TMP_ALLOC (nw * BYTES_PER_MP_LIMB);
! 187: xp = MPFR_MANT(x);
! 188: carry = mpfr_round_raw (tmp, xp, MPFR_PREC(x), MPFR_SIGN(x) < 0,
! 189: prec, rnd_mode, &inexact);
! 190: MPFR_PREC(x) = prec;
! 191:
! 192: if (carry)
! 193: {
! 194: mp_exp_t exp = MPFR_EXP(x);
! 195:
! 196: if (exp == __mpfr_emax)
! 197: (void) mpfr_set_overflow(x, rnd_mode, MPFR_SIGN(x));
! 198: else
! 199: {
! 200: MPFR_EXP(x)++;
! 201: xp[nw - 1] = GMP_LIMB_HIGHBIT;
! 202: if (nw - 1 > 0)
! 203: MPN_ZERO(xp, nw - 1);
! 204: }
! 205: }
! 206: else
! 207: MPN_COPY(xp, tmp, nw);
! 208:
! 209: TMP_FREE(marker);
! 210: return inexact;
! 211: }
! 212:
! 213: /* assumption: BITS_PER_MP_LIMB is a power of 2 */
! 214:
! 215: /* assuming b is an approximation of x in direction rnd1
! 216: with error at most 2^(MPFR_EXP(b)-err), returns 1 if one is
! 217: able to round exactly x to precision prec with direction rnd2,
! 218: and 0 otherwise.
! 219:
! 220: Side effects: none.
! 221: */
! 222:
! 223: int
! 224: mpfr_can_round (mpfr_ptr b, mp_exp_t err, mp_rnd_t rnd1,
! 225: mp_rnd_t rnd2, mp_prec_t prec)
! 226: {
! 227: return MPFR_IS_ZERO(b) ? 0 : /* we cannot round if b=0 */
! 228: mpfr_can_round_raw (MPFR_MANT(b),
! 229: (MPFR_PREC(b) - 1)/BITS_PER_MP_LIMB + 1,
! 230: MPFR_SIGN(b), err, rnd1, rnd2, prec);
! 231: }
! 232:
! 233: int
! 234: mpfr_can_round_raw (mp_limb_t *bp, mp_size_t bn, int neg, mp_exp_t err0,
! 235: mp_rnd_t rnd1, mp_rnd_t rnd2, mp_prec_t prec)
! 236: {
! 237: mp_prec_t err;
! 238: mp_size_t k, k1, tn;
! 239: int s, s1;
! 240: mp_limb_t cc, cc2;
! 241: mp_limb_t *tmp;
! 242: TMP_DECL(marker);
! 243:
! 244: if (err0 < 0 || (mp_exp_unsigned_t) err0 <= prec)
! 245: return 0; /* can't round */
! 246:
! 247: neg = neg <= 0;
! 248:
! 249: /* if the error is smaller than ulp(b), then anyway it will propagate
! 250: up to ulp(b) */
! 251: err = ((mp_exp_unsigned_t) err0 > (mp_prec_t) bn * BITS_PER_MP_LIMB) ?
! 252: (mp_prec_t) bn * BITS_PER_MP_LIMB : err0;
! 253:
! 254: /* warning: if k = m*BITS_PER_MP_LIMB, consider limb m-1 and not m */
! 255: k = (err - 1) / BITS_PER_MP_LIMB;
! 256: s = err % BITS_PER_MP_LIMB;
! 257: if (s)
! 258: s = BITS_PER_MP_LIMB - s;
! 259: /* the error corresponds to bit s in limb k, the most significant limb
! 260: being limb 0 */
! 261: k1 = (prec - 1) / BITS_PER_MP_LIMB;
! 262: s1 = prec % BITS_PER_MP_LIMB;
! 263: if (s1)
! 264: s1 = BITS_PER_MP_LIMB - s1;
! 265:
! 266: /* the last significant bit is bit s1 in limb k1 */
! 267:
! 268: /* don't need to consider the k1 most significant limbs */
! 269: k -= k1;
! 270: bn -= k1;
! 271: prec -= (mp_prec_t) k1 * BITS_PER_MP_LIMB;
! 272: /* if when adding or subtracting (1 << s) in bp[bn-1-k], it does not
! 273: change bp[bn-1] >> s1, then we can round */
! 274:
! 275: if (rnd1 == GMP_RNDU)
! 276: if (neg)
! 277: rnd1 = GMP_RNDZ;
! 278:
! 279: if (rnd1 == GMP_RNDD)
! 280: rnd1 = neg ? GMP_RNDU : GMP_RNDZ;
! 281:
! 282: /* in the sequel, RNDU = towards infinity, RNDZ = towards zero */
! 283:
! 284: TMP_MARK(marker);
! 285: tn = bn;
! 286: k++; /* since we work with k+1 everywhere */
! 287: tmp = TMP_ALLOC(tn * BYTES_PER_MP_LIMB);
! 288: if (bn > k)
! 289: MPN_COPY (tmp, bp, bn - k);
! 290:
! 291: if (rnd1 != GMP_RNDN)
! 292: { /* GMP_RNDZ or GMP_RNDU */
! 293: cc = (bp[bn - 1] >> s1) & 1;
! 294: cc ^= mpfr_round_raw2(bp, bn, neg, rnd2, prec);
! 295:
! 296: /* now round b +/- 2^(MPFR_EXP(b)-err) */
! 297: cc2 = rnd1 == GMP_RNDZ ?
! 298: mpn_add_1 (tmp + bn - k, bp + bn - k, k, MP_LIMB_T_ONE << s) :
! 299: mpn_sub_1 (tmp + bn - k, bp + bn - k, k, MP_LIMB_T_ONE << s);
! 300: }
! 301: else
! 302: { /* GMP_RNDN */
! 303: /* first round b+2^(MPFR_EXP(b)-err) */
! 304: cc = mpn_add_1 (tmp + bn - k, bp + bn - k, k, MP_LIMB_T_ONE << s);
! 305: cc = (tmp[bn - 1] >> s1) & 1; /* gives 0 when cc=1 */
! 306: cc ^= mpfr_round_raw2 (tmp, bn, neg, rnd2, prec);
! 307:
! 308: /* now round b-2^(MPFR_EXP(b)-err) */
! 309: cc2 = mpn_sub_1 (tmp + bn - k, bp + bn - k, k, MP_LIMB_T_ONE << s);
! 310: }
! 311:
! 312: if (cc2 && cc)
! 313: {
! 314: TMP_FREE(marker);
! 315: return 0;
! 316: }
! 317:
! 318: cc2 = (tmp[bn - 1] >> s1) & 1;
! 319: cc2 ^= mpfr_round_raw2 (tmp, bn, neg, rnd2, prec);
! 320:
! 321: TMP_FREE(marker);
! 322: return cc == cc2;
! 323: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to round_prec.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr