Annotation of OpenXM_contrib/gmp/mpz/aorsmul_i.c, Revision 1.1
1.1 ! ohara 1: /* mpz_addmul_ui, mpz_submul_ui - add or subtract small multiple.
! 2:
! 3: THE mpz_aorsmul_1 FUNCTION IN THIS FILE IS FOR INTERNAL USE ONLY AND IS
! 4: ALMOST CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR
! 5: COMPLETELY IN FUTURE GNU MP RELEASES.
! 6:
! 7: Copyright 2001, 2002 Free Software Foundation, Inc.
! 8:
! 9: This file is part of the GNU MP Library.
! 10:
! 11: The GNU MP Library is free software; you can redistribute it and/or modify
! 12: it under the terms of the GNU Lesser General Public License as published by
! 13: the Free Software Foundation; either version 2.1 of the License, or (at your
! 14: option) any later version.
! 15:
! 16: The GNU MP Library is distributed in the hope that it will be useful, but
! 17: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
! 18: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
! 19: License for more details.
! 20:
! 21: You should have received a copy of the GNU Lesser General Public License
! 22: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
! 23: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
! 24: MA 02111-1307, USA. */
! 25:
! 26: #include "gmp.h"
! 27: #include "gmp-impl.h"
! 28:
! 29:
! 30: #if HAVE_NATIVE_mpn_mul_1c
! 31: #define MPN_MUL_1C(cout, dst, src, size, n, cin) \
! 32: do { \
! 33: (cout) = mpn_mul_1c (dst, src, size, n, cin); \
! 34: } while (0)
! 35: #else
! 36: #define MPN_MUL_1C(cout, dst, src, size, n, cin) \
! 37: do { \
! 38: mp_limb_t __cy; \
! 39: __cy = mpn_mul_1 (dst, src, size, n); \
! 40: (cout) = __cy + mpn_add_1 (dst, dst, size, cin); \
! 41: } while (0)
! 42: #endif
! 43:
! 44:
! 45: /* sub>=0 means an addmul w += x*y, sub<0 means a submul w -= x*y.
! 46:
! 47: All that's needed to account for negative w or x is to flip "sub".
! 48:
! 49: The final w will retain its sign, unless an underflow occurs in a submul
! 50: of absolute values, in which case it's flipped.
! 51:
! 52: If x has more limbs than w, then mpn_submul_1 followed by mpn_com_n is
! 53: used. The alternative would be mpn_mul_1 into temporary space followed
! 54: by mpn_sub_n. Avoiding temporary space seem good, and submul+com stands
! 55: a chance of being faster since it involves only one set of carry
! 56: propagations, not two. Note that doing an addmul_1 with a
! 57: twos-complement negative y doesn't work, because it effectively adds an
! 58: extra x * 2^BITS_PER_MP_LIMB. */
! 59:
! 60: void
! 61: mpz_aorsmul_1 (mpz_ptr w, mpz_srcptr x, mp_limb_t y, mp_size_t sub)
! 62: {
! 63: mp_size_t xsize, wsize, wsize_signed, new_wsize, min_size, dsize;
! 64: mp_srcptr xp;
! 65: mp_ptr wp;
! 66: mp_limb_t cy;
! 67:
! 68: /* w unaffected if x==0 or y==0 */
! 69: xsize = SIZ (x);
! 70: if (xsize == 0 || y == 0)
! 71: return;
! 72:
! 73: sub ^= xsize;
! 74: xsize = ABS (xsize);
! 75:
! 76: wsize_signed = SIZ (w);
! 77: if (wsize_signed == 0)
! 78: {
! 79: /* nothing to add to, just set x*y, "sub" gives the sign */
! 80: MPZ_REALLOC (w, xsize+1);
! 81: wp = PTR (w);
! 82: cy = mpn_mul_1 (wp, PTR(x), xsize, y);
! 83: wp[xsize] = cy;
! 84: xsize += (cy != 0);
! 85: SIZ (w) = (sub >= 0 ? xsize : -xsize);
! 86: return;
! 87: }
! 88:
! 89: sub ^= wsize_signed;
! 90: wsize = ABS (wsize_signed);
! 91:
! 92: new_wsize = MAX (wsize, xsize);
! 93: MPZ_REALLOC (w, new_wsize+1);
! 94: wp = PTR (w);
! 95: xp = PTR (x);
! 96: min_size = MIN (wsize, xsize);
! 97:
! 98: if (sub >= 0)
! 99: {
! 100: /* addmul of absolute values */
! 101:
! 102: cy = mpn_addmul_1 (wp, xp, min_size, y);
! 103: wp += min_size;
! 104: xp += min_size;
! 105:
! 106: dsize = xsize - wsize;
! 107: #if HAVE_NATIVE_mpn_mul_1c
! 108: if (dsize > 0)
! 109: cy = mpn_mul_1c (wp, xp, dsize, y, cy);
! 110: else if (dsize < 0)
! 111: {
! 112: dsize = -dsize;
! 113: cy = mpn_add_1 (wp, wp, dsize, cy);
! 114: }
! 115: #else
! 116: if (dsize != 0)
! 117: {
! 118: mp_limb_t cy2;
! 119: if (dsize > 0)
! 120: cy2 = mpn_mul_1 (wp, xp, dsize, y);
! 121: else
! 122: {
! 123: dsize = -dsize;
! 124: cy2 = 0;
! 125: }
! 126: cy = cy2 + mpn_add_1 (wp, wp, dsize, cy);
! 127: }
! 128: #endif
! 129:
! 130: wp[dsize] = cy;
! 131: new_wsize += (cy != 0);
! 132: }
! 133: else
! 134: {
! 135: /* submul of absolute values */
! 136:
! 137: cy = mpn_submul_1 (wp, xp, min_size, y);
! 138: if (wsize >= xsize)
! 139: {
! 140: /* if w bigger than x, then propagate borrow through it */
! 141: if (wsize != xsize)
! 142: cy = mpn_sub_1 (wp+xsize, wp+xsize, wsize-xsize, cy);
! 143:
! 144: if (cy != 0)
! 145: {
! 146: /* Borrow out of w, take twos complement negative to get
! 147: absolute value, flip sign of w. */
! 148: wp[new_wsize] = ~-cy; /* extra limb is 0-cy */
! 149: mpn_com_n (wp, wp, new_wsize);
! 150: new_wsize++;
! 151: MPN_INCR_U (wp, new_wsize, CNST_LIMB(1));
! 152: wsize_signed = -wsize_signed;
! 153: }
! 154: }
! 155: else /* wsize < xsize */
! 156: {
! 157: /* x bigger than w, so want x*y-w. Submul has given w-x*y, so
! 158: take twos complement and use an mpn_mul_1 for the rest. */
! 159:
! 160: mp_limb_t cy2;
! 161:
! 162: /* -(-cy*b^n + w-x*y) = (cy-1)*b^n + ~(w-x*y) + 1 */
! 163: mpn_com_n (wp, wp, wsize);
! 164: cy += mpn_add_1 (wp, wp, wsize, CNST_LIMB(1));
! 165: cy -= 1;
! 166:
! 167: /* If cy-1 == -1 then hold that -1 for latter. mpn_submul_1 never
! 168: returns cy==MP_LIMB_T_MAX so that value always indicates a -1. */
! 169: cy2 = (cy == MP_LIMB_T_MAX);
! 170: cy += cy2;
! 171: MPN_MUL_1C (cy, wp+wsize, xp+wsize, xsize-wsize, y, cy);
! 172: wp[new_wsize] = cy;
! 173: new_wsize += (cy != 0);
! 174:
! 175: /* Apply any -1 from above. The value at wp+wsize is non-zero
! 176: because y!=0 and the high limb of x will be non-zero. */
! 177: if (cy2)
! 178: MPN_DECR_U (wp+wsize, new_wsize-wsize, CNST_LIMB(1));
! 179:
! 180: wsize_signed = -wsize_signed;
! 181: }
! 182:
! 183: /* submul can produce high zero limbs due to cancellation, both when w
! 184: has more limbs or x has more */
! 185: MPN_NORMALIZE (wp, new_wsize);
! 186: }
! 187:
! 188: SIZ (w) = (wsize_signed >= 0 ? new_wsize : -new_wsize);
! 189:
! 190: ASSERT (new_wsize == 0 || PTR(w)[new_wsize-1] != 0);
! 191: }
! 192:
! 193:
! 194: void
! 195: mpz_addmul_ui (mpz_ptr w, mpz_srcptr x, unsigned long y)
! 196: {
! 197: mpz_aorsmul_1 (w, x, (mp_limb_t) y & GMP_NUMB_MASK, (mp_size_t) 0);
! 198: #if GMP_NAIL_BITS != 0
! 199: if (y > GMP_NUMB_MAX && SIZ(x) != 0)
! 200: {
! 201: mpz_t t;
! 202: mp_ptr tp;
! 203: mp_size_t xn;
! 204: TMP_DECL (mark);
! 205: TMP_MARK (mark);
! 206: xn = SIZ (x);
! 207: MPZ_TMP_INIT (t, ABS (xn) + 1);
! 208: tp = PTR (t);
! 209: tp[0] = 0;
! 210: MPN_COPY (tp + 1, PTR(x), ABS (xn));
! 211: SIZ(t) = xn >= 0 ? xn + 1 : xn - 1;
! 212: mpz_aorsmul_1 (w, t, (mp_limb_t) y >> GMP_NUMB_BITS, (mp_size_t) 0);
! 213: TMP_FREE (mark);
! 214: }
! 215: #endif
! 216: }
! 217:
! 218: void
! 219: mpz_submul_ui (mpz_ptr w, mpz_srcptr x, unsigned long y)
! 220: {
! 221: mpz_aorsmul_1 (w, x, (mp_limb_t) y & GMP_NUMB_MASK, (mp_size_t) -1);
! 222: #if GMP_NAIL_BITS != 0
! 223: if (y > GMP_NUMB_MAX && SIZ(x) != 0)
! 224: {
! 225: mpz_t t;
! 226: mp_ptr tp;
! 227: mp_size_t xn;
! 228: TMP_DECL (mark);
! 229: TMP_MARK (mark);
! 230: xn = SIZ (x);
! 231: MPZ_TMP_INIT (t, ABS (xn) + 1);
! 232: tp = PTR (t);
! 233: tp[0] = 0;
! 234: MPN_COPY (tp + 1, PTR(x), ABS (xn));
! 235: SIZ(t) = xn >= 0 ? xn + 1 : xn - 1;
! 236: mpz_aorsmul_1 (w, t, (mp_limb_t) y >> GMP_NUMB_BITS, (mp_size_t) -1);
! 237: TMP_FREE (mark);
! 238: }
! 239: #endif
! 240: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to aorsmul_i.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpz