Annotation of OpenXM_contrib/gmp/mpfr/sub1.c, Revision 1.1
1.1 ! ohara 1: /* mpfr_sub1 -- internal function to perform a "real" subtraction
! 2:
! 3: Copyright 2001 Free Software Foundation.
! 4: Contributed by the Spaces project, INRIA Lorraine.
! 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: /* compute sign(b) * (|b| - |c|), with |b| > |c|, diff_exp = EXP(b) - EXP(c)
! 29: Returns 0 iff result is exact,
! 30: a negative value when the result is less than the exact value,
! 31: a positive value otherwise.
! 32: */
! 33:
! 34: int
! 35: mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode,
! 36: int sub)
! 37: {
! 38: int sign;
! 39: mp_exp_unsigned_t diff_exp;
! 40: mp_prec_t cancel, cancel1;
! 41: mp_size_t cancel2, an, bn, cn, cn0;
! 42: mp_limb_t *ap, *bp, *cp;
! 43: mp_limb_t carry, bb, cc, borrow = 0;
! 44: int inexact = 0, shift_b, shift_c, is_exact = 1, down = 0, add_exp = 0;
! 45: int sh, k;
! 46: TMP_DECL(marker);
! 47:
! 48: TMP_MARK(marker);
! 49: ap = MPFR_MANT(a);
! 50: an = 1 + (MPFR_PREC(a) - 1) / BITS_PER_MP_LIMB;
! 51:
! 52: sign = mpfr_cmp2 (b, c, &cancel);
! 53: if (sign == 0)
! 54: {
! 55: if (rnd_mode == GMP_RNDD)
! 56: MPFR_SET_NEG(a);
! 57: else
! 58: MPFR_SET_POS(a);
! 59: MPFR_SET_ZERO(a);
! 60: MPFR_RET(0);
! 61: }
! 62:
! 63: /* If subtraction: sign(a) = sign * sign(b) */
! 64: if (sub && MPFR_SIGN(a) != sign * MPFR_SIGN(b))
! 65: MPFR_CHANGE_SIGN(a);
! 66:
! 67: if (sign < 0) /* swap b and c so that |b| > |c| */
! 68: {
! 69: mpfr_srcptr t;
! 70: t = b; b = c; c = t;
! 71: }
! 72:
! 73: /* If addition: sign(a) = sign of the larger argument in absolute value */
! 74: if (!sub)
! 75: MPFR_SET_SAME_SIGN(a, b);
! 76:
! 77: diff_exp = (mp_exp_unsigned_t) MPFR_EXP(b) - MPFR_EXP(c);
! 78:
! 79: /* reserve a space to store b aligned with the result, i.e. shifted by
! 80: (-cancel) % BITS_PER_MP_LIMB to the right */
! 81: bn = 1 + (MPFR_PREC(b) - 1) / BITS_PER_MP_LIMB;
! 82: shift_b = cancel % BITS_PER_MP_LIMB;
! 83: if (shift_b)
! 84: shift_b = BITS_PER_MP_LIMB - shift_b;
! 85: cancel1 = (cancel + shift_b) / BITS_PER_MP_LIMB;
! 86: /* the high cancel1 limbs from b should not be taken into account */
! 87: if (shift_b == 0)
! 88: bp = MPFR_MANT(b); /* no need of an extra space */
! 89: else
! 90: {
! 91: bp = TMP_ALLOC ((bn + 1) * BYTES_PER_MP_LIMB);
! 92: bp[0] = mpn_rshift (bp + 1, MPFR_MANT(b), bn++, shift_b);
! 93: }
! 94:
! 95: /* reserve a space to store c aligned with the result, i.e. shifted by
! 96: (diff_exp-cancel) % BITS_PER_MP_LIMB to the right */
! 97: cn = 1 + (MPFR_PREC(c) - 1) / BITS_PER_MP_LIMB;
! 98: shift_c = diff_exp - (cancel % BITS_PER_MP_LIMB);
! 99: shift_c = (shift_c + BITS_PER_MP_LIMB) % BITS_PER_MP_LIMB;
! 100: if (shift_c == 0)
! 101: cp = MPFR_MANT(c);
! 102: else
! 103: {
! 104: cp = TMP_ALLOC ((cn + 1) * BYTES_PER_MP_LIMB);
! 105: cp[0] = mpn_rshift (cp + 1, MPFR_MANT(c), cn++, shift_c);
! 106: }
! 107:
! 108: #ifdef DEBUG
! 109: printf("shift_b=%u shift_c=%u\n", shift_b, shift_c);
! 110: #endif
! 111:
! 112: /* ensure ap != bp and ap != cp */
! 113: if (ap == bp)
! 114: {
! 115: bp = (mp_ptr) TMP_ALLOC(bn * BYTES_PER_MP_LIMB);
! 116: MPN_COPY (bp, ap, bn);
! 117: /* ap == cp cannot occur since we would have b=c, which is detected
! 118: in mpfr_add or mpfr_sub */
! 119: }
! 120: else if (ap == cp)
! 121: {
! 122: cp = (mp_ptr) TMP_ALLOC (cn * BYTES_PER_MP_LIMB);
! 123: MPN_COPY(cp, ap, cn);
! 124: }
! 125:
! 126: /* here we have shift_c = (diff_exp - cancel) % BITS_PER_MP_LIMB,
! 127: thus we want cancel2 = ceil((cancel - diff_exp) / BITS_PER_MP_LIMB) */
! 128:
! 129: cancel2 = (long int) (cancel - (diff_exp - shift_c)) / BITS_PER_MP_LIMB;
! 130: /* the high cancel2 limbs from b should not be taken into account */
! 131: #ifdef DEBUG
! 132: printf("cancel=%u cancel1=%u cancel2=%d\n", cancel, cancel1, cancel2);
! 133: #endif
! 134:
! 135: /* ap[an-1] ap[0]
! 136: <----------------+-----------|---->
! 137: <----------PREC(a)----------><-sh->
! 138: cancel1
! 139: limbs bp[bn-cancel1-1]
! 140: <--...-----><----------------+-----------+----------->
! 141: cancel2
! 142: limbs cp[cn-cancel2-1] cancel2 >= 0
! 143: <--...--><----------------+----------------+---------------->
! 144: (-cancel2) cancel2 < 0
! 145: limbs <----------------+---------------->
! 146: */
! 147:
! 148: /* first part: put in ap[0..an-1] the value of high(b) - high(c),
! 149: where high(b) consists of the high an+cancel1 limbs of b,
! 150: and high(c) consists of the high an+cancel2 limbs of c.
! 151: */
! 152:
! 153: /* copy high(b) into a */
! 154: if (an + cancel1 <= bn) /* a: <----------------+-----------|---->
! 155: b: <-----------------------------------------> */
! 156: MPN_COPY (ap, bp + bn - (an + cancel1), an);
! 157: else /* a: <----------------+-----------|---->
! 158: b: <-------------------------> */
! 159: if (cancel1 < bn) /* otherwise b does not overlap with a */
! 160: {
! 161: MPN_ZERO (ap, an + cancel1 - bn);
! 162: MPN_COPY (ap + an + cancel1 - bn, bp, bn - cancel1);
! 163: }
! 164: else
! 165: MPN_ZERO (ap, an);
! 166:
! 167: #ifdef DEBUG
! 168: printf("after copying high(b), a="); mpfr_print_binary(a); putchar('\n');
! 169: #endif
! 170:
! 171: /* subtract high(c) */
! 172: if (an + cancel2 > 0) /* otherwise c does not overlap with a */
! 173: {
! 174: mp_limb_t *ap2;
! 175:
! 176: if (cancel2 >= 0)
! 177: {
! 178: if (an + cancel2 <= cn) /* a: <----------------------------->
! 179: c: <-----------------------------------------> */
! 180: mpn_sub_n (ap, ap, cp + cn - (an + cancel2), an);
! 181: else /* a: <---------------------------->
! 182: c: <-------------------------> */
! 183: {
! 184: ap2 = ap + an + cancel2 - cn;
! 185: if (cn > cancel2)
! 186: mpn_sub_n (ap2, ap2, cp, cn - cancel2);
! 187: }
! 188: }
! 189: else /* cancel2 < 0 */
! 190: {
! 191: if (an + cancel2 <= cn) /* a: <----------------------------->
! 192: c: <-----------------------------> */
! 193: borrow = mpn_sub_n (ap, ap, cp + cn - (an + cancel2), an + cancel2);
! 194: else /* a: <---------------------------->
! 195: c: <----------------> */
! 196: {
! 197: ap2 = ap + an + cancel2 - cn;
! 198: borrow = mpn_sub_n (ap2, ap2, cp, cn);
! 199: }
! 200: ap2 = ap + an + cancel2;
! 201: mpn_sub_1 (ap2, ap2, -cancel2, borrow);
! 202: }
! 203: }
! 204:
! 205: #ifdef DEBUG
! 206: printf("after subtracting high(c), a="); mpfr_print_binary(a); putchar('\n');
! 207: #endif
! 208:
! 209: /* now perform rounding */
! 210: sh = an * BITS_PER_MP_LIMB - MPFR_PREC(a); /* last unused bits from a */
! 211: carry = ap[0] & ((MP_LIMB_T_ONE << sh) - MP_LIMB_T_ONE);
! 212: ap[0] -= carry;
! 213:
! 214: if (rnd_mode == GMP_RNDN)
! 215: {
! 216: if (sh)
! 217: {
! 218: is_exact = (carry == 0);
! 219: /* can decide except when carry = 2^(sh-1) [middle]
! 220: or carry = 0 [truncate, but cannot decide inexact flag] */
! 221: down = (carry < (MP_LIMB_T_ONE << (sh - 1)));
! 222: if (carry > (MP_LIMB_T_ONE << (sh - 1)))
! 223: goto add_one_ulp;
! 224: else if ((0 < carry) && down)
! 225: {
! 226: inexact = -1; /* result if smaller than exact value */
! 227: goto truncate;
! 228: }
! 229: }
! 230: }
! 231: else /* directed rounding: set rnd_mode to RNDZ iff towards zero */
! 232: {
! 233: if (((rnd_mode == GMP_RNDD) && (MPFR_SIGN(a) > 0)) ||
! 234: ((rnd_mode == GMP_RNDU) && (MPFR_SIGN(a) < 0)))
! 235: rnd_mode = GMP_RNDZ;
! 236:
! 237: if (carry)
! 238: {
! 239: if (rnd_mode == GMP_RNDZ)
! 240: {
! 241: inexact = -1;
! 242: goto truncate;
! 243: }
! 244: else /* round away */
! 245: goto add_one_ulp;
! 246: }
! 247: }
! 248:
! 249: /* we have to consider the low (bn - (an+cancel1)) limbs from b,
! 250: and the (cn - (an+cancel2)) limbs from c. */
! 251: bn -= an + cancel1;
! 252: cn0 = cn;
! 253: cn -= (long int) an + cancel2;
! 254: #ifdef DEBUG
! 255: printf("last %d bits from a are %lu, bn=%ld, cn=%ld\n", sh, carry, bn, cn);
! 256: #endif
! 257:
! 258: for (k = 0; (bn > 0) || (cn > 0); k = 1)
! 259: {
! 260: bb = (bn > 0) ? bp[--bn] : 0;
! 261: if ((cn > 0) && (cn-- <= cn0))
! 262: cc = cp[cn];
! 263: else
! 264: cc = 0;
! 265:
! 266: if (down == 0)
! 267: down = (bb < cc);
! 268:
! 269: if ((rnd_mode == GMP_RNDN) && !k && sh == 0)
! 270: {
! 271: mp_limb_t half = GMP_LIMB_HIGHBIT;
! 272:
! 273: is_exact = (bb == cc);
! 274:
! 275: /* add one ulp if bb > cc + half
! 276: truncate if cc - half < bb < cc + half
! 277: sub one ulp if bb < cc - half
! 278: */
! 279:
! 280: if (down)
! 281: {
! 282: if (cc >= half)
! 283: cc -= half;
! 284: else
! 285: bb += half;
! 286: }
! 287: else /* bb >= cc */
! 288: {
! 289: if (cc < half)
! 290: cc += half;
! 291: else
! 292: bb -= half;
! 293: }
! 294: }
! 295:
! 296: #ifdef DEBUG
! 297: printf(" bb=%lu cc=%lu down=%d is_exact=%d\n", bb, cc, down, is_exact);
! 298: #endif
! 299: if (bb < cc)
! 300: {
! 301: if (rnd_mode == GMP_RNDZ)
! 302: goto sub_one_ulp;
! 303: else if (rnd_mode != GMP_RNDN) /* round away */
! 304: {
! 305: inexact = 1;
! 306: goto truncate;
! 307: }
! 308: else /* round to nearest */
! 309: {
! 310: if (is_exact && sh == 0)
! 311: {
! 312: inexact = 0;
! 313: goto truncate;
! 314: }
! 315: else if (down && sh == 0)
! 316: goto sub_one_ulp;
! 317: else
! 318: {
! 319: inexact = (is_exact) ? 1 : -1;
! 320: goto truncate;
! 321: }
! 322: }
! 323: }
! 324: else if (bb > cc)
! 325: {
! 326: if (rnd_mode == GMP_RNDZ)
! 327: {
! 328: inexact = -1;
! 329: goto truncate;
! 330: }
! 331: else if (rnd_mode != GMP_RNDN) /* round away */
! 332: goto add_one_ulp;
! 333: else /* round to nearest */
! 334: {
! 335: if (is_exact)
! 336: {
! 337: inexact = -1;
! 338: goto truncate;
! 339: }
! 340: else if (down)
! 341: {
! 342: inexact = 1;
! 343: goto truncate;
! 344: }
! 345: else
! 346: goto add_one_ulp;
! 347: }
! 348: }
! 349: }
! 350:
! 351: if ((rnd_mode == GMP_RNDN) && !is_exact)
! 352: {
! 353: /* even rounding rule */
! 354: if ((ap[0] >> sh) & 1)
! 355: {
! 356: if (down) goto sub_one_ulp;
! 357: else goto add_one_ulp;
! 358: }
! 359: else
! 360: inexact = (down) ? 1 : -1;
! 361: }
! 362: else
! 363: inexact = 0;
! 364: goto truncate;
! 365:
! 366: sub_one_ulp: /* add one unit in last place to a */
! 367: mpn_sub_1 (ap, ap, an, MP_LIMB_T_ONE << sh);
! 368: inexact = -1;
! 369: goto end_of_sub;
! 370:
! 371: add_one_ulp: /* add one unit in last place to a */
! 372: if (mpn_add_1 (ap, ap, an, MP_LIMB_T_ONE << sh)) /* result is a power of 2 */
! 373: {
! 374: ap[an-1] = GMP_LIMB_HIGHBIT;
! 375: add_exp = 1;
! 376: }
! 377: inexact = 1; /* result larger than exact value */
! 378:
! 379: truncate:
! 380: if ((ap[an-1] >> (BITS_PER_MP_LIMB - 1)) == 0) /* case 1 - epsilon */
! 381: {
! 382: ap[an-1] = GMP_LIMB_HIGHBIT;
! 383: add_exp = 1;
! 384: }
! 385:
! 386: end_of_sub:
! 387: /* we have to set MPFR_EXP(a) to MPFR_EXP(b) - cancel + add_exp, taking
! 388: care of underflows/overflows in that computation, and of the allowed
! 389: exponent range */
! 390: if (cancel)
! 391: {
! 392: mp_exp_t exp_b;
! 393:
! 394: cancel -= add_exp; /* still valid as unsigned long */
! 395: exp_b = MPFR_EXP(b); /* save it in case a equals b */
! 396: MPFR_EXP(a) = MPFR_EXP(b) - cancel;
! 397: if ((MPFR_EXP(a) > exp_b) /* underflow in type mp_exp_t */
! 398: || (MPFR_EXP(a) < __mpfr_emin))
! 399: {
! 400: TMP_FREE(marker);
! 401: return mpfr_set_underflow (a, rnd_mode, MPFR_SIGN(a));
! 402: }
! 403: }
! 404: else /* cancel = 0: MPFR_EXP(a) <- MPFR_EXP(b) + add_exp */
! 405: {
! 406: /* in case cancel = 0, add_exp can still be 1, in case b is just
! 407: below a power of two, c is very small, prec(a) < prec(b),
! 408: and rnd=away or nearest */
! 409: if (add_exp && MPFR_EXP(b) == __mpfr_emax)
! 410: {
! 411: TMP_FREE(marker);
! 412: return mpfr_set_overflow (a, rnd_mode, MPFR_SIGN(a));
! 413: }
! 414: MPFR_EXP(a) = MPFR_EXP(b) + add_exp;
! 415: }
! 416: TMP_FREE(marker);
! 417: #ifdef DEBUG
! 418: printf ("result is a="); mpfr_print_binary(a); putchar('\n');
! 419: #endif
! 420: /* check that result is msb-normalized */
! 421: MPFR_ASSERTN(ap[an-1] > ~ap[an-1]);
! 422: return inexact * MPFR_SIGN(a);
! 423: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to sub1.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr