Annotation of OpenXM_contrib/gmp/mpfr/div.c, Revision 1.1.1.2
1.1 maekawa 1: /* mpfr_div -- divide two floating-point numbers
2:
1.1.1.2 ! ohara 3: Copyright 1999, 2001 Free Software Foundation.
1.1 maekawa 4:
5: This file is part of the MPFR Library.
6:
7: The MPFR Library is free software; you can redistribute it and/or modify
1.1.1.2 ! ohara 8: it under the terms of the GNU Lesser General Public License as published by
! 9: the Free Software Foundation; either version 2.1 of the License, or (at your
1.1 maekawa 10: option) any later version.
11:
12: The MPFR Library is distributed in the hope that it will be useful, but
13: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.1.1.2 ! ohara 14: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
1.1 maekawa 15: License for more details.
16:
1.1.1.2 ! ohara 17: You should have received a copy of the GNU Lesser General Public License
1.1 maekawa 18: along with the MPFR Library; see the file COPYING.LIB. If not, write to
19: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
20: MA 02111-1307, USA. */
21:
22: #include "gmp.h"
23: #include "gmp-impl.h"
24: #include "longlong.h"
1.1.1.2 ! ohara 25: #include "mpfr.h"
! 26: #include "mpfr-impl.h"
1.1 maekawa 27:
1.1.1.2 ! ohara 28: int
! 29: mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mp_rnd_t rnd_mode)
1.1 maekawa 30: {
1.1.1.2 ! ohara 31: mp_srcptr up, vp, bp;
! 32: mp_size_t usize, vsize;
! 33:
! 34: mp_ptr ap, qp, rp;
! 35: mp_size_t asize, bsize, qsize, rsize;
! 36: mp_exp_t qexp;
! 37:
! 38: mp_rnd_t rnd_mode1, rnd_mode2;
! 39:
! 40: mp_size_t err, k;
! 41: mp_limb_t near;
! 42: int inex, sh, can_round, can_round2, sign_quotient;
! 43: unsigned int cc = 0, rw;
! 44:
1.1 maekawa 45: TMP_DECL (marker);
46:
47:
1.1.1.2 ! ohara 48: /**************************************************************************
! 49: * *
! 50: * This part of the code deals with special cases *
! 51: * *
! 52: **************************************************************************/
1.1 maekawa 53:
1.1.1.2 ! ohara 54: if (MPFR_IS_NAN(u) || MPFR_IS_NAN(v))
1.1 maekawa 55: {
1.1.1.2 ! ohara 56: MPFR_SET_NAN(q);
! 57: MPFR_RET_NAN;
1.1 maekawa 58: }
59:
1.1.1.2 ! ohara 60: MPFR_CLEAR_NAN(q);
1.1 maekawa 61:
1.1.1.2 ! ohara 62: sign_quotient = MPFR_SIGN(u) * MPFR_SIGN(v);
! 63: if (MPFR_SIGN(q) != sign_quotient)
! 64: MPFR_CHANGE_SIGN(q);
1.1 maekawa 65:
1.1.1.2 ! ohara 66: if (MPFR_IS_INF(u))
1.1 maekawa 67: {
1.1.1.2 ! ohara 68: if (MPFR_IS_INF(v))
! 69: {
! 70: MPFR_SET_NAN(q);
! 71: MPFR_RET_NAN;
! 72: }
! 73: else
! 74: {
! 75: MPFR_SET_INF(q);
! 76: MPFR_RET(0);
! 77: }
1.1 maekawa 78: }
1.1.1.2 ! ohara 79: else
! 80: if (MPFR_IS_INF(v))
! 81: {
! 82: MPFR_CLEAR_INF(q);
! 83: MPFR_SET_ZERO(q);
! 84: MPFR_RET(0);
! 85: }
1.1 maekawa 86:
1.1.1.2 ! ohara 87: MPFR_CLEAR_INF(q); /* clear Inf flag */
1.1 maekawa 88:
1.1.1.2 ! ohara 89: if (MPFR_IS_ZERO(v))
! 90: {
! 91: if (MPFR_IS_ZERO(u))
! 92: {
! 93: MPFR_SET_NAN(q);
! 94: MPFR_RET_NAN;
! 95: }
! 96: else
! 97: {
! 98: MPFR_SET_INF(q);
! 99: MPFR_RET(0);
! 100: }
! 101: }
1.1 maekawa 102:
1.1.1.2 ! ohara 103: if (MPFR_IS_ZERO(u))
1.1 maekawa 104: {
1.1.1.2 ! ohara 105: MPFR_SET_ZERO(q);
! 106: MPFR_RET(0);
! 107: }
1.1 maekawa 108:
1.1.1.2 ! ohara 109: /**************************************************************************
! 110: * *
! 111: * End of the part concerning special values. *
! 112: * *
! 113: **************************************************************************/
! 114:
! 115:
! 116: up = MPFR_MANT(u);
! 117: vp = MPFR_MANT(v);
! 118:
! 119: TMP_MARK (marker);
! 120: usize = MPFR_ESIZE(u);
! 121: vsize = MPFR_ESIZE(v);
! 122:
! 123: /**************************************************************************
! 124: * *
! 125: * First try to use only part of u, v. If this is not sufficient, *
! 126: * use the full u and v, to avoid long computations eg. in the case *
! 127: * u = v. *
! 128: * *
! 129: **************************************************************************/
1.1 maekawa 130:
1.1.1.2 ! ohara 131: /* The dividend is a, length asize. The divisor is b, length bsize. */
1.1 maekawa 132:
1.1.1.2 ! ohara 133: qsize = (MPFR_PREC(q) + 3)/BITS_PER_MP_LIMB + 1;
! 134: if (vsize < qsize)
! 135: {
! 136: bsize = vsize;
! 137: bp = vp;
! 138: }
! 139: else
! 140: {
! 141: bsize = qsize;
! 142: bp = (mp_srcptr)vp + vsize - qsize;
! 143: }
! 144:
! 145: asize = bsize + qsize;
! 146: ap = (mp_ptr) TMP_ALLOC(asize * BYTES_PER_MP_LIMB);
! 147: if (asize > usize)
! 148: {
! 149: MPN_COPY(ap + asize - usize, up, usize);
! 150: MPN_ZERO(ap, asize - usize);
! 151: }
! 152: else
! 153: MPN_COPY(ap, up + usize - asize, asize);
! 154:
! 155: /* Allocate limbs for quotient and remainder. */
! 156: qp = (mp_ptr) TMP_ALLOC ((qsize + 1) * BYTES_PER_MP_LIMB);
! 157: rp = (mp_ptr) TMP_ALLOC (bsize * BYTES_PER_MP_LIMB);
! 158: rsize = bsize;
1.1 maekawa 159:
1.1.1.2 ! ohara 160: mpn_tdiv_qr(qp, rp, 0, ap, asize, bp, bsize);
! 161:
! 162: /* Estimate number of correct bits. */
! 163:
! 164: err = qsize * BITS_PER_MP_LIMB;
! 165: if (bsize < vsize) err -= 2; else if (asize < usize) err --;
! 166:
! 167: /* We want to check if rounding is possible, but without normalizing
! 168: because we might have to divide again if rounding is impossible, or
! 169: if the result might be exact. We have however to mimic normalization */
! 170:
! 171: if (qp[qsize] != 0) { sh = -1; }
! 172: else { count_leading_zeros(sh, qp[qsize - 1]); }
! 173:
! 174: /*
! 175: To detect asap if the result is inexact, so as to avoid doing the
! 176: division completely, we perform the following check :
! 177:
! 178: - if rnd_mode == GMP_RNDN, and the result is exact, we are unable
! 179: to round simultaneously to zero and to infinity ;
! 180:
! 181: - if rnd_mode == GMP_RNDN, and if we can round to zero with one extra
! 182: bit of precision, we can decide rounding. Hence in that case, check
! 183: as in the case of GMP_RNDN, with one extra bit. Note that in the case
! 184: of close to even rounding we shall do the division completely, but
! 185: this is necessary anyway : we need to know whether this is really
! 186: even rounding or not.
! 187: */
! 188:
! 189: if (rnd_mode == GMP_RNDN)
! 190: {
! 191: rnd_mode1 = GMP_RNDZ;
! 192: near = 1;
! 193: }
! 194: else
! 195: {
! 196: rnd_mode1 = rnd_mode;
! 197: near = 0;
! 198: }
! 199:
! 200: sh += near;
! 201: can_round = mpfr_can_round_raw(qp, qsize + 1, sign_quotient, err + sh +
! 202: BITS_PER_MP_LIMB, GMP_RNDN, rnd_mode1,
! 203: MPFR_PREC(q) + sh + BITS_PER_MP_LIMB);
! 204:
! 205: switch (rnd_mode1)
! 206: {
! 207: case GMP_RNDU : rnd_mode2 = GMP_RNDD; break;
! 208: case GMP_RNDD : rnd_mode2 = GMP_RNDU; break;
! 209: case GMP_RNDZ : rnd_mode2 = sign_quotient == 1 ? GMP_RNDU : GMP_RNDD;
! 210: break;
! 211: default : rnd_mode2 = GMP_RNDZ;
! 212: }
! 213:
! 214: can_round2 = mpfr_can_round_raw(qp, qsize + 1, sign_quotient, err + sh +
! 215: BITS_PER_MP_LIMB, GMP_RNDN, rnd_mode2,
! 216: MPFR_PREC(q) + sh + BITS_PER_MP_LIMB);
! 217:
! 218: sh -= near;
! 219:
! 220: /* If either can_round or can_round2 is 0, either we cannot round or
! 221: the result might be exact. If asize >= usize and bsize >= vsize, we
! 222: can just check this by looking at the remainder. Otherwise, we
! 223: have to correct our first approximation. */
! 224:
! 225: if ((!can_round || !can_round2) && (asize < usize || bsize < vsize))
! 226: {
! 227: int b = 0;
! 228: mp_ptr rem, rem2;
! 229:
! 230: /**************************************************************************
! 231: * *
! 232: * The attempt to use only part of u and v failed. We first compute a *
! 233: * correcting term, then perform the full division. *
! 234: * Put u = uhi + ulo, v = vhi + vlo. We have uhi = vhi * qp + rp, *
! 235: * thus u - qp * v = rp + ulo - qp * vlo, that we shall divide by v. *
! 236: * *
! 237: **************************************************************************/
! 238:
! 239: rsize = qsize + 1 +
! 240: (usize - asize > vsize - bsize
! 241: ? usize - asize
! 242: : vsize - bsize);
! 243:
! 244: /*
! 245: TODO : One operand is probably enough, but then we have to
! 246: perform one further comparison (compute first vlo * q,
! 247: try to substract r, try to substract ulo. Which is best ?
! 248: NB : ulo and r do not overlap. Draw advantage of this
! 249: [eg. HI(vlo*q) = r => compare LO(vlo*q) with b.]
! 250: */
! 251:
! 252: rem = TMP_ALLOC(rsize * BYTES_PER_MP_LIMB);
! 253: rem2 = TMP_ALLOC(rsize * BYTES_PER_MP_LIMB);
! 254:
! 255: rem[rsize - 1] = rem2 [rsize - 1] = 0;
! 256:
! 257: if (bsize < vsize)
! 258: {
! 259: /* Compute vlo * q */
! 260: if (qsize + 1 > vsize - bsize)
! 261: mpn_mul(rem + rsize - vsize - qsize - 1 + bsize,
! 262: qp, qsize + 1, vp, vsize - bsize);
! 263: else
! 264: mpn_mul(rem + rsize - vsize - qsize - 1 + bsize,
! 265: vp, vsize - bsize, qp, qsize + 1);
! 266:
! 267: MPN_ZERO(rem, rsize - vsize - qsize - 1 + bsize);
! 268: }
! 269: else MPN_ZERO(rem, rsize);
! 270:
! 271: /* Compute ulo + r. The two of them do not overlap. */
! 272: MPN_COPY(rem2 + rsize - 1 - qsize, rp, bsize);
1.1 maekawa 273:
1.1.1.2 ! ohara 274: if (qsize + 1 > bsize)
! 275: MPN_ZERO(rem2 + rsize - 1 - qsize + bsize, qsize + 1 - bsize);
! 276:
! 277: if (asize < usize)
1.1 maekawa 278: {
1.1.1.2 ! ohara 279: MPN_COPY(rem2 + rsize - 1 - qsize - usize + asize,
! 280: up, usize - asize);
! 281: MPN_ZERO(rem2, rsize - 1 - qsize - usize + asize);
1.1 maekawa 282: }
1.1.1.2 ! ohara 283: else
! 284: MPN_ZERO(rem2, rsize - 1 - qsize);
! 285:
! 286: b = 0;
! 287: if (mpn_cmp(rem2, rem, rsize) >= 0)
! 288: {
! 289: /* Positive correction is at most 1. */
! 290:
! 291: mpn_sub_n(rem, rem2, rem, rsize);
! 292: if (rem[rsize - 1] != 0 ||
! 293: mpn_cmp(rem + rsize - vsize - 1, vp, vsize) >= 0)
! 294: {
! 295: rem[rsize - 1] -=
! 296: mpn_sub_n(rem + rsize - vsize - 1,
! 297: rem + rsize - vsize - 1,
! 298: vp, vsize);
! 299: qp[qsize] -= mpn_add_1(qp, qp, qsize, 1);
! 300: }
! 301: }
! 302: else
1.1 maekawa 303: {
1.1.1.2 ! ohara 304: /* Negative correction is at most 3 */
! 305: do
! 306: {
! 307: b++;
! 308: rem2[rsize - 1] +=
! 309: mpn_add_n(rem2 + rsize - vsize - 1,
! 310: rem2 + rsize - vsize - 1, vp, vsize);
! 311: }
! 312: while (mpn_cmp(rem2, rem, rsize) < 0);
! 313:
! 314: qp[qsize] -= mpn_sub_1(qp, qp, qsize, b);
! 315: mpn_sub_n(rem, rem2, rem, rsize);
1.1 maekawa 316: }
1.1.1.2 ! ohara 317:
! 318: if (qp[qsize] != 0)
! 319: sh = -1;
! 320: else
! 321: count_leading_zeros(sh, qp[qsize - 1]);
1.1 maekawa 322:
1.1.1.2 ! ohara 323: err = BITS_PER_MP_LIMB * qsize;
! 324: rp = rem;
! 325: }
! 326:
! 327: /**************************************************************************
! 328: * *
! 329: * Final stuff (rounding and so.) *
! 330: * From now on : qp is the quotient [size qsize], rp the remainder *
! 331: * [size rsize]. *
! 332: **************************************************************************/
1.1 maekawa 333:
1.1.1.2 ! ohara 334: qexp = MPFR_EXP(u) - MPFR_EXP(v);
1.1 maekawa 335:
1.1.1.2 ! ohara 336: if (qp[qsize] != 0)
! 337: /* Hack : qp[qsize] is 0, 1 or 2, hence if not 0, = 2^(qp[qsize] - 1). */
! 338: {
! 339: near = mpn_rshift(qp, qp, qsize, qp[qsize]);
! 340: qp[qsize - 1] |= GMP_LIMB_HIGHBIT; qexp += qp[qsize];
! 341: }
! 342: else
! 343: {
! 344: near = 0;
! 345: if (sh != 0)
1.1 maekawa 346: {
1.1.1.2 ! ohara 347: mpn_lshift(qp, qp, qsize, sh);
! 348: qexp -= sh;
1.1 maekawa 349: }
350: }
1.1.1.2 ! ohara 351:
! 352: cc = mpfr_round_raw_generic(qp, qp, err, (sign_quotient == -1 ? 1 : 0),
! 353: MPFR_PREC(q), rnd_mode, &inex, 1);
1.1 maekawa 354:
1.1.1.2 ! ohara 355: qp += qsize - MPFR_ESIZE(q); /* 0 or 1 */
! 356: qsize = MPFR_ESIZE(q);
! 357:
! 358: /*
! 359: At that point, either we were able to round from the beginning,
! 360: and know thus that the result is inexact.
! 361:
! 362: Or we have performed a full division. In that case, we might still
! 363: be wrong if both
! 364: - the remainder is nonzero ;
! 365: - we are rounding to infinity or to nearest (the nasty case of even
! 366: rounding).
! 367: - inex = 0, meaning that the non-significant bits of the quotients are 0,
! 368: except when rounding to nearest (the nasty case of even rounding again).
! 369: */
! 370:
! 371: if (!can_round || !can_round2) /* Lazy case. */
1.1 maekawa 372: {
1.1.1.2 ! ohara 373: if (inex == 0)
! 374: {
! 375: k = rsize - 1;
! 376:
! 377: /* If a bit has been shifted out during normalization, hence
! 378: the remainder is nonzero. */
! 379: if (near == 0)
! 380: while (k >= 0) { if (rp[k]) break; k--; }
! 381:
! 382: if (k >= 0) /* Remainder is nonzero. */
1.1 maekawa 383: {
1.1.1.2 ! ohara 384: if ((rnd_mode == GMP_RNDD && sign_quotient == -1)
! 385: || (rnd_mode == GMP_RNDU && sign_quotient == 1))
! 386: /* Rounding to infinity. */
1.1 maekawa 387: {
1.1.1.2 ! ohara 388: inex = sign_quotient;
! 389: cc = 1;
1.1 maekawa 390: }
1.1.1.2 ! ohara 391: /* rounding to zero. */
! 392: else inex = -sign_quotient;
1.1 maekawa 393: }
1.1.1.2 ! ohara 394: }
! 395: else /* We might have to correct an even rounding if remainder
! 396: is nonzero and if even rounding was towards 0. */
! 397: if (rnd_mode == GMP_RNDN && (inex == MPFR_EVEN_INEX
! 398: || inex == -MPFR_EVEN_INEX))
! 399: {
! 400: k = rsize - 1;
! 401:
! 402: /* If a bit has been shifted out during normalization, hence
! 403: the remainder is nonzero. */
! 404: if (near == 0)
! 405: while (k >= 0)
! 406: {
! 407: if (rp[k])
! 408: break;
! 409: k--;
! 410: }
! 411:
! 412: if (k >= 0) /* In fact the quotient is larger than expected */
! 413: {
! 414: inex = sign_quotient; /* To infinity, finally. */
! 415: cc = 1;
! 416: }
! 417: }
! 418: }
! 419:
! 420: /* Final modification due to rounding */
! 421: if (cc)
! 422: {
! 423: sh = MPFR_PREC(q) & (BITS_PER_MP_LIMB - 1);
! 424: if (sh)
! 425: cc = mpn_add_1 (qp, qp, qsize,
! 426: MP_LIMB_T_ONE << (BITS_PER_MP_LIMB - sh));
! 427: else
! 428: cc = mpn_add_1 (qp, qp, qsize, MP_LIMB_T_ONE);
1.1 maekawa 429:
1.1.1.2 ! ohara 430: if (cc)
! 431: {
! 432: mpn_rshift (qp, qp, qsize, 1);
! 433: qp[qsize-1] |= GMP_LIMB_HIGHBIT;
! 434: qexp++;
! 435: }
! 436: }
1.1 maekawa 437:
1.1.1.2 ! ohara 438: rw = qsize * BITS_PER_MP_LIMB - MPFR_PREC(q);
! 439: MPN_COPY(MPFR_MANT(q), qp, qsize);
1.1 maekawa 440: TMP_FREE (marker);
1.1.1.2 ! ohara 441:
! 442: MPFR_MANT(q)[0] &= ~((MP_LIMB_T_ONE << rw) - MP_LIMB_T_ONE);
! 443: MPFR_EXP(q) = qexp;
! 444:
! 445: MPFR_RET(inex);
1.1 maekawa 446: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to div.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr