![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to sub.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpf|
Return to sub.c CVS log | Up to [local] / OpenXM_contrib / gmp / mpf |
1.1 ! maekawa 1: /* mpf_sub -- Subtract two floats.
! 2:
! 3: Copyright (C) 1993, 1994, 1995, 1996 Free Software Foundation, Inc.
! 4:
! 5: This file is part of the GNU MP Library.
! 6:
! 7: The GNU MP Library is free software; you can redistribute it and/or modify
! 8: it under the terms of the GNU Library General Public License as published by
! 9: the Free Software Foundation; either version 2 of the License, or (at your
! 10: option) any later version.
! 11:
! 12: The GNU MP Library is distributed in the hope that it will be useful, but
! 13: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
! 14: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
! 15: License for more details.
! 16:
! 17: You should have received a copy of the GNU Library General Public License
! 18: along with the GNU MP 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:
! 25: void
! 26: #if __STDC__
! 27: mpf_sub (mpf_ptr r, mpf_srcptr u, mpf_srcptr v)
! 28: #else
! 29: mpf_sub (r, u, v)
! 30: mpf_ptr r;
! 31: mpf_srcptr u;
! 32: mpf_srcptr v;
! 33: #endif
! 34: {
! 35: mp_srcptr up, vp;
! 36: mp_ptr rp, tp;
! 37: mp_size_t usize, vsize, rsize;
! 38: mp_size_t prec;
! 39: mp_exp_t exp;
! 40: mp_size_t ediff;
! 41: int negate;
! 42: TMP_DECL (marker);
! 43:
! 44: usize = u->_mp_size;
! 45: vsize = v->_mp_size;
! 46:
! 47: /* Handle special cases that don't work in generic code below. */
! 48: if (usize == 0)
! 49: {
! 50: mpf_neg (r, v);
! 51: return;
! 52: }
! 53: if (vsize == 0)
! 54: {
! 55: mpf_set (r, u);
! 56: return;
! 57: }
! 58:
! 59: /* If signs of U and V are different, perform addition. */
! 60: if ((usize ^ vsize) < 0)
! 61: {
! 62: __mpf_struct v_negated;
! 63: v_negated._mp_size = -vsize;
! 64: v_negated._mp_exp = v->_mp_exp;
! 65: v_negated._mp_d = v->_mp_d;
! 66: mpf_add (r, u, &v_negated);
! 67: return;
! 68: }
! 69:
! 70: TMP_MARK (marker);
! 71:
! 72: /* Signs are now known to be the same. */
! 73: negate = usize < 0;
! 74:
! 75: /* Make U be the operand with the largest exponent. */
! 76: if (u->_mp_exp < v->_mp_exp)
! 77: {
! 78: mpf_srcptr t;
! 79: t = u; u = v; v = t;
! 80: negate ^= 1;
! 81: usize = u->_mp_size;
! 82: vsize = v->_mp_size;
! 83: }
! 84:
! 85: usize = ABS (usize);
! 86: vsize = ABS (vsize);
! 87: up = u->_mp_d;
! 88: vp = v->_mp_d;
! 89: rp = r->_mp_d;
! 90: prec = r->_mp_prec + 1;
! 91: exp = u->_mp_exp;
! 92: ediff = u->_mp_exp - v->_mp_exp;
! 93:
! 94: /* If ediff is 0 or 1, we might have a situation where the operands are
! 95: extremely close. We need to scan the operands from the most significant
! 96: end ignore the initial parts that are equal. */
! 97: if (ediff <= 1)
! 98: {
! 99: if (ediff == 0)
! 100: {
! 101: /* Skip leading limbs in U and V that are equal. */
! 102: if (up[usize - 1] == vp[vsize - 1])
! 103: {
! 104: /* This loop normally exits immediately. Optimize for that. */
! 105: do
! 106: {
! 107: usize--;
! 108: vsize--;
! 109: exp--;
! 110:
! 111: if (usize == 0)
! 112: {
! 113: rsize = vsize;
! 114: tp = (mp_ptr) vp;
! 115: negate ^= 1;
! 116: goto normalize;
! 117: }
! 118: if (vsize == 0)
! 119: {
! 120: rsize = usize;
! 121: tp = (mp_ptr) up;
! 122: goto normalize;
! 123: }
! 124: }
! 125: while (up[usize - 1] == vp[vsize - 1]);
! 126: }
! 127:
! 128: if (up[usize - 1] < vp[vsize - 1])
! 129: {
! 130: /* For simplicity, swap U and V. Note that since the loop above
! 131: wouldn't have exited unless up[usize - 1] and vp[vsize - 1]
! 132: were non-equal, this if-statement catches all cases where U
! 133: is smaller than V. */
! 134: { mp_srcptr tp = up; up = vp; vp = tp; }
! 135: { mp_size_t tsize = usize; usize = vsize; vsize = tsize; }
! 136: negate ^= 1;
! 137: /* negating ediff not necessary since it is 0. */
! 138: }
! 139:
! 140: /* Check for
! 141: x+1 00000000 ...
! 142: x ffffffff ... */
! 143: if (up[usize - 1] != vp[vsize - 1] + 1)
! 144: goto general_case;
! 145: usize--;
! 146: vsize--;
! 147: exp--;
! 148: }
! 149: else /* ediff == 1 */
! 150: {
! 151: /* Check for
! 152: 1 00000000 ...
! 153: 0 ffffffff ... */
! 154:
! 155: if (up[usize - 1] != 1 || vp[vsize - 1] != ~(mp_limb_t) 0
! 156: || (usize >= 2 && up[usize - 2] != 0))
! 157: goto general_case;
! 158:
! 159: usize--;
! 160: exp--;
! 161: }
! 162:
! 163: /* Skip sequences of 00000000/ffffffff */
! 164: while (vsize != 0 && usize != 0 && up[usize - 1] == 0
! 165: && vp[vsize - 1] == ~(mp_limb_t) 0)
! 166: {
! 167: usize--;
! 168: vsize--;
! 169: exp--;
! 170: }
! 171:
! 172: if (usize == 0)
! 173: {
! 174: while (vsize != 0 && vp[vsize - 1] == ~(mp_limb_t) 0)
! 175: {
! 176: vsize--;
! 177: exp--;
! 178: }
! 179: }
! 180:
! 181: if (usize > prec - 1)
! 182: {
! 183: up += usize - (prec - 1);
! 184: usize = prec - 1;
! 185: }
! 186: if (vsize > prec - 1)
! 187: {
! 188: vp += vsize - (prec - 1);
! 189: vsize = prec - 1;
! 190: }
! 191:
! 192: tp = (mp_ptr) TMP_ALLOC (prec * BYTES_PER_MP_LIMB);
! 193: {
! 194: mp_limb_t cy_limb;
! 195: if (vsize == 0)
! 196: {
! 197: mp_size_t size, i;
! 198: size = usize;
! 199: for (i = 0; i < size; i++)
! 200: tp[i] = up[i];
! 201: tp[size] = 1;
! 202: rsize = size + 1;
! 203: exp++;
! 204: goto normalize;
! 205: }
! 206: if (usize == 0)
! 207: {
! 208: mp_size_t size, i;
! 209: size = vsize;
! 210: for (i = 0; i < size; i++)
! 211: tp[i] = ~vp[i];
! 212: cy_limb = 1 - mpn_add_1 (tp, tp, vsize, (mp_limb_t) 1);
! 213: rsize = vsize;
! 214: if (cy_limb == 0)
! 215: {
! 216: tp[rsize] = 1;
! 217: rsize++;
! 218: exp++;
! 219: }
! 220: goto normalize;
! 221: }
! 222: if (usize >= vsize)
! 223: {
! 224: /* uuuu */
! 225: /* vv */
! 226: mp_size_t size;
! 227: size = usize - vsize;
! 228: MPN_COPY (tp, up, size);
! 229: cy_limb = mpn_sub_n (tp + size, up + size, vp, vsize);
! 230: rsize = usize;
! 231: }
! 232: else /* (usize < vsize) */
! 233: {
! 234: /* uuuu */
! 235: /* vvvvvvv */
! 236: mp_size_t size, i;
! 237: size = vsize - usize;
! 238: for (i = 0; i < size; i++)
! 239: tp[i] = ~vp[i];
! 240: cy_limb = mpn_sub_n (tp + size, up, vp + size, usize);
! 241: cy_limb+= mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1);
! 242: cy_limb-= mpn_add_1 (tp, tp, vsize, (mp_limb_t) 1);
! 243: rsize = vsize;
! 244: }
! 245: if (cy_limb == 0)
! 246: {
! 247: tp[rsize] = 1;
! 248: rsize++;
! 249: exp++;
! 250: }
! 251: goto normalize;
! 252: }
! 253: }
! 254:
! 255: general_case:
! 256: /* If U extends beyond PREC, ignore the part that does. */
! 257: if (usize > prec)
! 258: {
! 259: up += usize - prec;
! 260: usize = prec;
! 261: }
! 262:
! 263: /* If V extends beyond PREC, ignore the part that does.
! 264: Note that this may make vsize negative. */
! 265: if (vsize + ediff > prec)
! 266: {
! 267: vp += vsize + ediff - prec;
! 268: vsize = prec - ediff;
! 269: }
! 270:
! 271: /* Allocate temp space for the result. Allocate
! 272: just vsize + ediff later??? */
! 273: tp = (mp_ptr) TMP_ALLOC (prec * BYTES_PER_MP_LIMB);
! 274:
! 275: if (ediff >= prec)
! 276: {
! 277: /* V completely cancelled. */
! 278: if (tp != up)
! 279: MPN_COPY (rp, up, usize);
! 280: rsize = usize;
! 281: }
! 282: else
! 283: {
! 284: /* Locate the least significant non-zero limb in (the needed
! 285: parts of) U and V, to simplify the code below. */
! 286: for (;;)
! 287: {
! 288: if (vsize == 0)
! 289: {
! 290: MPN_COPY (rp, up, usize);
! 291: rsize = usize;
! 292: goto done;
! 293: }
! 294: if (vp[0] != 0)
! 295: break;
! 296: vp++, vsize--;
! 297: }
! 298: for (;;)
! 299: {
! 300: if (usize == 0)
! 301: {
! 302: MPN_COPY (rp, vp, vsize);
! 303: rsize = vsize;
! 304: negate ^= 1;
! 305: goto done;
! 306: }
! 307: if (up[0] != 0)
! 308: break;
! 309: up++, usize--;
! 310: }
! 311:
! 312: /* uuuu | uuuu | uuuu | uuuu | uuuu */
! 313: /* vvvvvvv | vv | vvvvv | v | vv */
! 314:
! 315: if (usize > ediff)
! 316: {
! 317: /* U and V partially overlaps. */
! 318: if (ediff == 0)
! 319: {
! 320: /* Have to compare the leading limbs of u and v
! 321: to determine whether to compute u - v or v - u. */
! 322: if (usize >= vsize)
! 323: {
! 324: /* uuuu */
! 325: /* vv */
! 326: mp_size_t size;
! 327: size = usize - vsize;
! 328: MPN_COPY (tp, up, size);
! 329: mpn_sub_n (tp + size, up + size, vp, vsize);
! 330: rsize = usize;
! 331: }
! 332: else /* (usize < vsize) */
! 333: {
! 334: /* uuuu */
! 335: /* vvvvvvv */
! 336: mp_size_t size, i;
! 337: size = vsize - usize;
! 338: tp[0] = -vp[0];
! 339: for (i = 1; i < size; i++)
! 340: tp[i] = ~vp[i];
! 341: mpn_sub_n (tp + size, up, vp + size, usize);
! 342: mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1);
! 343: rsize = vsize;
! 344: }
! 345: }
! 346: else
! 347: {
! 348: if (vsize + ediff <= usize)
! 349: {
! 350: /* uuuu */
! 351: /* v */
! 352: mp_size_t size;
! 353: size = usize - ediff - vsize;
! 354: MPN_COPY (tp, up, size);
! 355: mpn_sub (tp + size, up + size, usize - size, vp, vsize);
! 356: rsize = usize;
! 357: }
! 358: else
! 359: {
! 360: /* uuuu */
! 361: /* vvvvv */
! 362: mp_size_t size, i;
! 363: size = vsize + ediff - usize;
! 364: tp[0] = -vp[0];
! 365: for (i = 1; i < size; i++)
! 366: tp[i] = ~vp[i];
! 367: mpn_sub (tp + size, up, usize, vp + size, usize - ediff);
! 368: mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1);
! 369: rsize = vsize + ediff;
! 370: }
! 371: }
! 372: }
! 373: else
! 374: {
! 375: /* uuuu */
! 376: /* vv */
! 377: mp_size_t size, i;
! 378: size = vsize + ediff - usize;
! 379: tp[0] = -vp[0];
! 380: for (i = 1; i < vsize; i++)
! 381: tp[i] = ~vp[i];
! 382: for (i = vsize; i < size; i++)
! 383: tp[i] = ~(mp_limb_t) 0;
! 384: mpn_sub_1 (tp + size, up, usize, (mp_limb_t) 1);
! 385: rsize = size + usize;
! 386: }
! 387:
! 388: normalize:
! 389: /* Full normalize. Optimize later. */
! 390: while (rsize != 0 && tp[rsize - 1] == 0)
! 391: {
! 392: rsize--;
! 393: exp--;
! 394: }
! 395: MPN_COPY (rp, tp, rsize);
! 396: }
! 397:
! 398: done:
! 399: r->_mp_size = negate ? -rsize : rsize;
! 400: r->_mp_exp = exp;
! 401: TMP_FREE (marker);
! 402: }