Annotation of OpenXM_contrib/gmp/mpfr/round.c, Revision 1.1
1.1 ! maekawa 1: /* mpfr_round_raw2, mpfr_round_raw, mpfr_round, mpfr_can_round,
! 2: mpfr_can_round_raw -- various rounding functions
! 3:
! 4: Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
! 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 Library General Public License as published by
! 10: the Free Software Foundation; either version 2 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 Library General Public
! 16: License for more details.
! 17:
! 18: You should have received a copy of the GNU Library 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 <stdio.h>
! 24: #include "gmp.h"
! 25: #include "gmp-impl.h"
! 26: #include "mpfr.h"
! 27:
! 28: #ifdef Exp
! 29: #include "longlong.h"
! 30: #endif
! 31:
! 32: /* returns 0 if round(sign*xp[0..xn-1], prec, rnd) =
! 33: round(sign*xp[0..xn-1], prec, GMP_RNDZ), 1 otherwise,
! 34: where sign=1 if neg=0, sign=-1 otherwise.
! 35:
! 36: Does *not* modify anything.
! 37: */
! 38: int
! 39: #if __STDC__
! 40: mpfr_round_raw2(mp_limb_t *xp, unsigned long xn,
! 41: char neg, char rnd, unsigned long prec)
! 42: #else
! 43: mpfr_round_raw2(xp, xn, neg, rnd, prec)
! 44: mp_limb_t *xp;
! 45: unsigned long xn;
! 46: char neg;
! 47: char rnd;
! 48: unsigned long prec;
! 49: #endif
! 50: {
! 51: unsigned long nw; long wd; char rw; short l; mp_limb_t mask;
! 52:
! 53: nw = prec / BITS_PER_MP_LIMB; rw = prec & (BITS_PER_MP_LIMB - 1);
! 54: if (rw) nw++;
! 55: if (rnd==GMP_RNDZ || xn<nw || (rnd==GMP_RNDU && neg)
! 56: || (rnd==GMP_RNDD && neg==0)) return 0;
! 57:
! 58: mask = ~((((mp_limb_t)1)<<(BITS_PER_MP_LIMB - rw)) - 1);
! 59: switch (rnd)
! 60: {
! 61: case GMP_RNDU:
! 62: case GMP_RNDD:
! 63: if (xp[xn - nw] & ~mask) return 1;
! 64: for (l = nw + 1;l <= xn; l++)
! 65: if (xp[xn - l]) break;
! 66: return (l <= xn);
! 67:
! 68: case GMP_RNDN:
! 69: /* First check if we are just halfway between two representable numbers */
! 70: wd = xn - nw;
! 71: if (!rw)
! 72: {
! 73: if (!wd) /* all bits are significative */ return 0;
! 74: wd--;
! 75: if (xp[wd] == ((mp_limb_t)1 << (BITS_PER_MP_LIMB - 1)))
! 76: {
! 77: do wd--; while (wd > 0 && !xp[wd]);
! 78: if (!wd) { return 1; } else return xp[xn - nw] & 1;
! 79: }
! 80:
! 81: return xp[wd]>>(BITS_PER_MP_LIMB - 1);
! 82: }
! 83: else
! 84: if (rw + 1 < BITS_PER_MP_LIMB)
! 85: {
! 86: if ((xp[wd] & (~mask)) == (((mp_limb_t)1) << (BITS_PER_MP_LIMB - rw - 1)))
! 87: do { wd--; } while (wd >= 0 && !xp[wd]);
! 88: else return ((xp[wd]>>(BITS_PER_MP_LIMB - rw - 1)) & 1);
! 89:
! 90: /* first limb was in the middle, and others down to wd+1 were 0 */
! 91: if (wd>=0) return 1;
! 92: else
! 93: return ((xp[xn - nw] & mask) >> (BITS_PER_MP_LIMB - rw)) & 1;
! 94: }
! 95: else
! 96: /* Modified PZ, 27 May 1999:
! 97: rw, i.e. the number of bits stored in xp[xn-nw], is
! 98: BITS_PER_MP_LIMB-1, i.e. there is exactly one non significant bit.
! 99: We are just halfway iff xp[wd] has its low significant bit
! 100: set and all limbs xp[0]...xp[wd-1] are zero */
! 101: {
! 102: if (xp[wd] & 1)
! 103: do wd--; while (wd >= 0 && !xp[wd]);
! 104: return ((wd<0) ? xp[xn-nw]>>1 : xp[xn-nw]) & 1;
! 105: }
! 106: default: return 0;
! 107: }
! 108: }
! 109:
! 110: /* puts in y the value of xp (with precision xprec and sign 1 if negative=0,
! 111: -1 otherwise) rounded to precision yprec and direction RND_MODE
! 112: Supposes x is not zero nor NaN nor +/- Infinity (i.e. *xp != 0).
! 113: */
! 114: int
! 115: #if __STDC__
! 116: mpfr_round_raw(mp_limb_t *y, mp_limb_t *xp, unsigned long xprec, char negative,
! 117: unsigned long yprec, char RND_MODE)
! 118: #else
! 119: mpfr_round_raw(y, xp, xprec, negative, yprec, RND_MODE)
! 120: mp_limb_t *y;
! 121: mp_limb_t *xp;
! 122: unsigned long xprec;
! 123: char negative;
! 124: unsigned long yprec;
! 125: char RND_MODE;
! 126: #endif
! 127: {
! 128: unsigned long nw, xsize; mp_limb_t mask;
! 129: char rw, xrw, carry = 0;
! 130:
! 131: xsize = (xprec-1)/BITS_PER_MP_LIMB + 1;
! 132: xrw = xprec % BITS_PER_MP_LIMB; if (xrw == 0) { xrw = BITS_PER_MP_LIMB; }
! 133:
! 134: #ifdef Exp
! 135: count_leading_zeros(flag, xp[xsize-1]);
! 136: yprec += flag;
! 137: #endif
! 138:
! 139: nw = yprec / BITS_PER_MP_LIMB; rw = yprec & (BITS_PER_MP_LIMB - 1);
! 140: if (rw) nw++;
! 141: /* number of words needed to represent x */
! 142:
! 143: mask = ~((((mp_limb_t)1)<<(BITS_PER_MP_LIMB - rw)) - (mp_limb_t)1);
! 144:
! 145: /* precision is larger than the size of x. No rounding is necessary.
! 146: Just add zeroes at the end */
! 147: if (xsize < nw) {
! 148: MPN_COPY(y + nw - xsize, xp, xsize);
! 149: MPN_ZERO(y, nw - xsize); /* PZ 27 May 99 */
! 150: y[0] &= mask;
! 151: return 0;
! 152: }
! 153:
! 154: if (mpfr_round_raw2(xp, xsize, negative, RND_MODE, yprec))
! 155: carry = mpn_add_1(y, xp + xsize - nw, nw,
! 156: ((mp_limb_t)1) << (BITS_PER_MP_LIMB - rw));
! 157: else MPN_COPY(y, xp + xsize - nw, nw);
! 158:
! 159: y[0] &= mask;
! 160: return carry;
! 161: }
! 162:
! 163: void
! 164: #if __STDC__
! 165: mpfr_round(mpfr_t x, char RND_MODE, unsigned long prec)
! 166: #else
! 167: mpfr_round(x, RND_MODE, prec)
! 168: mpfr_t x;
! 169: char RND_MODE;
! 170: unsigned long prec;
! 171: #endif
! 172: {
! 173: mp_limb_t *tmp; int carry; unsigned long nw;
! 174: TMP_DECL(marker);
! 175:
! 176: nw = prec / BITS_PER_MP_LIMB;
! 177: if (prec & (BITS_PER_MP_LIMB - 1)) nw++;
! 178: TMP_MARK(marker);
! 179: tmp = TMP_ALLOC (nw * BYTES_PER_MP_LIMB);
! 180: carry = mpfr_round_raw(tmp, MANT(x), PREC(x), (SIGN(x)<0), prec, RND_MODE);
! 181:
! 182: if (carry)
! 183: {
! 184: mpn_rshift(tmp, tmp, nw, 1);
! 185: tmp [nw-1] |= (((mp_limb_t)1) << (BITS_PER_MP_LIMB - 1));
! 186: EXP(x)++;
! 187: }
! 188:
! 189: if (SIGN(x)<0) { SIZE(x)=nw; CHANGE_SIGN(x); } else SIZE(x)=nw;
! 190: PREC(x) = prec;
! 191: MPN_COPY(MANT(x), tmp, nw);
! 192: TMP_FREE(marker);
! 193: }
! 194:
! 195: /* hypotheses : BITS_PER_MP_LIMB est une puissance de 2
! 196: dans un test, la premiere partie du && est evaluee d'abord */
! 197:
! 198:
! 199: /* assuming b is an approximation of x in direction rnd1
! 200: with error at most 2^(EXP(b)-err), returns 1 if one is
! 201: able to round exactly x to precision prec with direction rnd2,
! 202: and 0 otherwise.
! 203:
! 204: Side effects: none.
! 205: */
! 206:
! 207: int
! 208: #if __STDC__
! 209: mpfr_can_round(mpfr_t b, unsigned long err, unsigned char rnd1,
! 210: unsigned char rnd2, unsigned long prec)
! 211: #else
! 212: mpfr_can_round(b, err, rnd1, rnd2, prec)
! 213: mpfr_t b;
! 214: unsigned long err;
! 215: unsigned char rnd1;
! 216: unsigned char rnd2;
! 217: unsigned long prec;
! 218: #endif
! 219: {
! 220: return mpfr_can_round_raw(MANT(b), (PREC(b) - 1)/BITS_PER_MP_LIMB + 1,
! 221: SIGN(b), err, rnd1, rnd2, prec);
! 222: }
! 223:
! 224: int
! 225: #if __STDC__
! 226: mpfr_can_round_raw(mp_limb_t *bp, unsigned long bn, int neg,
! 227: unsigned long err, unsigned char rnd1, unsigned char rnd2,
! 228: unsigned long prec)
! 229: #else
! 230: mpfr_can_round_raw(bp, bn, neg, err, rnd1, rnd2, prec)
! 231: mp_limb_t *bp;
! 232: unsigned long bn;
! 233: int neg;
! 234: unsigned long err;
! 235: unsigned char rnd1;
! 236: unsigned char rnd2;
! 237: unsigned long prec;
! 238: #endif
! 239: {
! 240: int k, k1, l, l1, tn; mp_limb_t cc, cc2, *tmp;
! 241: TMP_DECL(marker);
! 242:
! 243: if (err<=prec) return 0;
! 244: neg = (neg > 0 ? 0 : 1);
! 245:
! 246: /* warning: if k = m*BITS_PER_MP_LIMB, consider limb m-1 and not m */
! 247: k = (err-1)/BITS_PER_MP_LIMB;
! 248: l = err % BITS_PER_MP_LIMB; if (l) l = BITS_PER_MP_LIMB-l;
! 249: /* the error corresponds to bit l in limb k */
! 250: k1 = (prec-1)/BITS_PER_MP_LIMB;
! 251: l1 = prec%BITS_PER_MP_LIMB; if (l1) l1 = BITS_PER_MP_LIMB-l1;
! 252:
! 253: /* the last significant bit is bit l1 in limb k1 */
! 254:
! 255: /* don't need to consider the k1 most significant limbs */
! 256: k -= k1; bn -= k1; prec -= k1*BITS_PER_MP_LIMB; k1=0;
! 257:
! 258: if (rnd1==GMP_RNDU) { if (neg) rnd1=GMP_RNDZ; }
! 259: if (rnd1==GMP_RNDD) { if (neg) rnd1=GMP_RNDU; else rnd1=GMP_RNDZ; }
! 260:
! 261: /* in the sequel, RNDU = towards infinity, RNDZ = towards zero */
! 262:
! 263: TMP_MARK(marker);
! 264: tn = bn;
! 265: k++; /* since we work with k+1 everywhere */
! 266:
! 267: switch (rnd1) {
! 268:
! 269: case GMP_RNDZ: /* b <= x <= b+2^(EXP(b)-err) */
! 270: tmp = TMP_ALLOC(tn*BYTES_PER_MP_LIMB);
! 271: cc = (bp[bn-1]>>l1) & 1;
! 272: cc ^= mpfr_round_raw2(bp, bn, neg, rnd2, prec);
! 273:
! 274: /* now round b+2^(EXP(b)-err) */
! 275: cc2 = mpn_add_1(tmp+bn-k, bp+bn-k, k, (mp_limb_t)1<<l);
! 276: /* if carry, then all bits up to err were to 1, and we can round only
! 277: if cc==0 and mpfr_round_raw2 returns 0 below */
! 278: if (cc2 && cc) { TMP_FREE(marker); return 0; }
! 279: cc2 = (tmp[bn-1]>>l1) & 1; /* gives 0 when carry */
! 280: cc2 ^= mpfr_round_raw2(tmp, bn, neg, rnd2, prec);
! 281:
! 282: TMP_FREE(marker);
! 283: return (cc == cc2);
! 284:
! 285: case GMP_RNDU: /* b-2^(EXP(b)-err) <= x <= b */
! 286: tmp = TMP_ALLOC(tn*BYTES_PER_MP_LIMB);
! 287: /* first round b */
! 288: cc = (bp[bn-1]>>l1) & 1;
! 289: cc ^= mpfr_round_raw2(bp, bn, neg, rnd2, prec);
! 290:
! 291: /* now round b-2^(EXP(b)-err) */
! 292: cc2 = mpn_sub_1(tmp+bn-k, bp+bn-k, k, (mp_limb_t)1<<l);
! 293: /* if borrow, then all bits up to err were to 0, and we can round only
! 294: if cc==0 and mpfr_round_raw2 returns 1 below */
! 295: if (cc2 && cc) { TMP_FREE(marker); return 0; }
! 296: cc2 = (tmp[bn-1]>>l1) & 1; /* gives 1 when carry */
! 297: cc2 ^= mpfr_round_raw2(tmp, bn, neg, rnd2, prec);
! 298:
! 299: TMP_FREE(marker);
! 300: return (cc == cc2);
! 301:
! 302: case GMP_RNDN: /* b-2^(EXP(b)-err-1) <= x <= b+2^(EXP(b)-err-1) */
! 303: if (l==0) tn++;
! 304: tmp = TMP_ALLOC(tn*BYTES_PER_MP_LIMB);
! 305:
! 306: /* this case is the same than GMP_RNDZ, except we first have to
! 307: subtract 2^(EXP(b)-err-1) from b */
! 308:
! 309: if (l) {
! 310: l--; /* tn=bn */
! 311: mpn_sub_1(tmp+tn-k, bp+bn-k, k, (mp_limb_t)1<<l);
! 312: }
! 313: else {
! 314: MPN_COPY(tmp+1, bp, bn); *tmp=0; /* extra limb to add or subtract 1 */
! 315: k++; l=BITS_PER_MP_LIMB-1;
! 316: mpn_sub_1(tmp+tn-k, tmp+tn-k, k, (mp_limb_t)1<<l);
! 317: }
! 318:
! 319: /* round b-2^(EXP(b)-err-1) */
! 320: /* we can disregard borrow, since we start from tmp in 2nd case too */
! 321: cc = (tmp[tn-1]>>l1) & 1;
! 322: cc ^= mpfr_round_raw2(tmp, tn, neg, rnd2, prec);
! 323:
! 324: if (l==BITS_PER_MP_LIMB-1) { l=0; k--; } else l++;
! 325:
! 326: /* round b+2^(EXP(b)-err-1) = b-2^(EXP(b)-err-1) + 2^(EXP(b)-err) */
! 327: cc2 = mpn_add_1(tmp+tn-k, tmp+tn-k, k, (mp_limb_t)1<<l);
! 328: /* if carry, then all bits up to err were to 1, and we can round only
! 329: if cc==0 and mpfr_round_raw2 returns 0 below */
! 330: if (cc2 && cc) { TMP_FREE(marker); return 0; }
! 331: cc2 = (tmp[tn-1]>>l1) & 1; /* gives 0 when carry */
! 332: cc2 ^= mpfr_round_raw2(tmp, tn, neg, rnd2, prec);
! 333:
! 334: TMP_FREE(marker);
! 335: return (cc == cc2);
! 336: }
! 337: return 0;
! 338: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to round.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr