![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to ui_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 ui_sub.c CVS log | Up to [local] / OpenXM_contrib / gmp / mpf |
1.1 ! maekawa 1: /* mpf_ui_sub -- Subtract a float from an unsigned long int.
! 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_ui_sub (mpf_ptr r, unsigned long int u, mpf_srcptr v)
! 28: #else
! 29: mpf_ui_sub (r, u, v)
! 30: mpf_ptr r;
! 31: unsigned long int 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 uexp;
! 40: mp_size_t ediff;
! 41: int negate;
! 42: mp_limb_t ulimb;
! 43: TMP_DECL (marker);
! 44:
! 45: vsize = v->_mp_size;
! 46:
! 47: /* Handle special cases that don't work in generic code below. */
! 48: if (u == 0)
! 49: {
! 50: mpf_neg (r, v);
! 51: return;
! 52: }
! 53: if (vsize == 0)
! 54: {
! 55: mpf_set_ui (r, u);
! 56: return;
! 57: }
! 58:
! 59: /* If signs of U and V are different, perform addition. */
! 60: if (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_ui (r, &v_negated, u);
! 67: return;
! 68: }
! 69:
! 70: TMP_MARK (marker);
! 71:
! 72: /* Signs are now known to be the same. */
! 73:
! 74: ulimb = u;
! 75: /* Make U be the operand with the largest exponent. */
! 76: if (1 < v->_mp_exp)
! 77: {
! 78: negate = 1;
! 79: usize = ABS (vsize);
! 80: vsize = 1;
! 81: up = v->_mp_d;
! 82: vp = &ulimb;
! 83: rp = r->_mp_d;
! 84: prec = r->_mp_prec + 1;
! 85: uexp = v->_mp_exp;
! 86: ediff = uexp - 1;
! 87: }
! 88: else
! 89: {
! 90: negate = 0;
! 91: usize = 1;
! 92: vsize = ABS (vsize);
! 93: up = &ulimb;
! 94: vp = v->_mp_d;
! 95: rp = r->_mp_d;
! 96: prec = r->_mp_prec;
! 97: uexp = 1;
! 98: ediff = 1 - v->_mp_exp;
! 99: }
! 100:
! 101: /* Ignore leading limbs in U and V that are equal. Doing
! 102: this helps increase the precision of the result. */
! 103: if (ediff == 0)
! 104: {
! 105: /* This loop normally exits immediately. Optimize for that. */
! 106: for (;;)
! 107: {
! 108: usize--;
! 109: vsize--;
! 110: if (up[usize] != vp[vsize])
! 111: break;
! 112: uexp--;
! 113: if (usize == 0)
! 114: goto Lu0;
! 115: if (vsize == 0)
! 116: goto Lv0;
! 117: }
! 118: usize++;
! 119: vsize++;
! 120: /* Note that either operand (but not both operands) might now have
! 121: leading zero limbs. It matters only that U is unnormalized if
! 122: vsize is now zero, and vice versa. And it is only in that case
! 123: that we have to adjust uexp. */
! 124: if (vsize == 0)
! 125: Lv0:
! 126: while (usize != 0 && up[usize - 1] == 0)
! 127: usize--, uexp--;
! 128: if (usize == 0)
! 129: Lu0:
! 130: while (vsize != 0 && vp[vsize - 1] == 0)
! 131: vsize--, uexp--;
! 132: }
! 133:
! 134: /* If U extends beyond PREC, ignore the part that does. */
! 135: if (usize > prec)
! 136: {
! 137: up += usize - prec;
! 138: usize = prec;
! 139: }
! 140:
! 141: /* If V extends beyond PREC, ignore the part that does.
! 142: Note that this may make vsize negative. */
! 143: if (vsize + ediff > prec)
! 144: {
! 145: vp += vsize + ediff - prec;
! 146: vsize = prec - ediff;
! 147: }
! 148:
! 149: /* Allocate temp space for the result. Allocate
! 150: just vsize + ediff later??? */
! 151: tp = (mp_ptr) TMP_ALLOC (prec * BYTES_PER_MP_LIMB);
! 152:
! 153: if (ediff >= prec)
! 154: {
! 155: /* V completely cancelled. */
! 156: if (tp != up)
! 157: MPN_COPY (rp, up, usize);
! 158: rsize = usize;
! 159: }
! 160: else
! 161: {
! 162: /* Locate the least significant non-zero limb in (the needed
! 163: parts of) U and V, to simplify the code below. */
! 164: for (;;)
! 165: {
! 166: if (vsize == 0)
! 167: {
! 168: MPN_COPY (rp, up, usize);
! 169: rsize = usize;
! 170: goto done;
! 171: }
! 172: if (vp[0] != 0)
! 173: break;
! 174: vp++, vsize--;
! 175: }
! 176: for (;;)
! 177: {
! 178: if (usize == 0)
! 179: {
! 180: MPN_COPY (rp, vp, vsize);
! 181: rsize = vsize;
! 182: negate ^= 1;
! 183: goto done;
! 184: }
! 185: if (up[0] != 0)
! 186: break;
! 187: up++, usize--;
! 188: }
! 189:
! 190: /* uuuu | uuuu | uuuu | uuuu | uuuu */
! 191: /* vvvvvvv | vv | vvvvv | v | vv */
! 192:
! 193: if (usize > ediff)
! 194: {
! 195: /* U and V partially overlaps. */
! 196: if (ediff == 0)
! 197: {
! 198: /* Have to compare the leading limbs of u and v
! 199: to determine whether to compute u - v or v - u. */
! 200: if (usize > vsize)
! 201: {
! 202: /* uuuu */
! 203: /* vv */
! 204: int cmp;
! 205: cmp = mpn_cmp (up + usize - vsize, vp, vsize);
! 206: if (cmp >= 0)
! 207: {
! 208: mp_size_t size;
! 209: size = usize - vsize;
! 210: MPN_COPY (tp, up, size);
! 211: mpn_sub_n (tp + size, up + size, vp, vsize);
! 212: rsize = usize;
! 213: }
! 214: else
! 215: {
! 216: /* vv */ /* Swap U and V. */
! 217: /* uuuu */
! 218: mp_size_t size, i;
! 219: size = usize - vsize;
! 220: tp[0] = -up[0];
! 221: for (i = 1; i < size; i++)
! 222: tp[i] = ~up[i];
! 223: mpn_sub_n (tp + size, vp, up + size, vsize);
! 224: mpn_sub_1 (tp + size, tp + size, vsize, (mp_limb_t) 1);
! 225: negate ^= 1;
! 226: rsize = usize;
! 227: }
! 228: }
! 229: else if (usize < vsize)
! 230: {
! 231: /* uuuu */
! 232: /* vvvvvvv */
! 233: int cmp;
! 234: cmp = mpn_cmp (up, vp + vsize - usize, usize);
! 235: if (cmp > 0)
! 236: {
! 237: mp_size_t size, i;
! 238: size = vsize - usize;
! 239: tp[0] = -vp[0];
! 240: for (i = 1; i < size; i++)
! 241: tp[i] = ~vp[i];
! 242: mpn_sub_n (tp + size, up, vp + size, usize);
! 243: mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1);
! 244: rsize = vsize;
! 245: }
! 246: else
! 247: {
! 248: /* vvvvvvv */ /* Swap U and V. */
! 249: /* uuuu */
! 250: /* This is the only place we can get 0.0. */
! 251: mp_size_t size;
! 252: size = vsize - usize;
! 253: MPN_COPY (tp, vp, size);
! 254: mpn_sub_n (tp + size, vp + size, up, usize);
! 255: negate ^= 1;
! 256: rsize = vsize;
! 257: }
! 258: }
! 259: else
! 260: {
! 261: /* uuuu */
! 262: /* vvvv */
! 263: int cmp;
! 264: cmp = mpn_cmp (up, vp + vsize - usize, usize);
! 265: if (cmp > 0)
! 266: {
! 267: mpn_sub_n (tp, up, vp, usize);
! 268: rsize = usize;
! 269: }
! 270: else
! 271: {
! 272: mpn_sub_n (tp, vp, up, usize);
! 273: negate ^= 1;
! 274: rsize = usize;
! 275: /* can give zero */
! 276: }
! 277: }
! 278: }
! 279: else
! 280: {
! 281: if (vsize + ediff <= usize)
! 282: {
! 283: /* uuuu */
! 284: /* v */
! 285: mp_size_t size;
! 286: size = usize - ediff - vsize;
! 287: MPN_COPY (tp, up, size);
! 288: mpn_sub (tp + size, up + size, usize - size, vp, vsize);
! 289: rsize = usize;
! 290: }
! 291: else
! 292: {
! 293: /* uuuu */
! 294: /* vvvvv */
! 295: mp_size_t size, i;
! 296: size = vsize + ediff - usize;
! 297: tp[0] = -vp[0];
! 298: for (i = 1; i < size; i++)
! 299: tp[i] = ~vp[i];
! 300: mpn_sub (tp + size, up, usize, vp + size, usize - ediff);
! 301: mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1);
! 302: rsize = vsize + ediff;
! 303: }
! 304: }
! 305: }
! 306: else
! 307: {
! 308: /* uuuu */
! 309: /* vv */
! 310: mp_size_t size, i;
! 311: size = vsize + ediff - usize;
! 312: tp[0] = -vp[0];
! 313: for (i = 1; i < vsize; i++)
! 314: tp[i] = ~vp[i];
! 315: for (i = vsize; i < size; i++)
! 316: tp[i] = ~(mp_limb_t) 0;
! 317: mpn_sub_1 (tp + size, up, usize, (mp_limb_t) 1);
! 318: rsize = size + usize;
! 319: }
! 320:
! 321: /* Full normalize. Optimize later. */
! 322: while (rsize != 0 && tp[rsize - 1] == 0)
! 323: {
! 324: rsize--;
! 325: uexp--;
! 326: }
! 327: MPN_COPY (rp, tp, rsize);
! 328: }
! 329:
! 330: done:
! 331: r->_mp_size = negate ? -rsize : rsize;
! 332: r->_mp_exp = uexp;
! 333: TMP_FREE (marker);
! 334: }