/* mpz_addmul, mpz_submul -- add or subtract multiple. Copyright 2001, 2004, 2005 Free Software Foundation, Inc. This file is part of the GNU MP Library. The GNU MP Library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version. The GNU MP Library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ #include "gmp.h" #include "gmp-impl.h" /* expecting x and y both with non-zero high limbs */ #define mpn_cmp_twosizes_lt(xp,xsize, yp,ysize) \ ((xsize) < (ysize) \ || ((xsize) == (ysize) && mpn_cmp (xp, yp, xsize) < 0)) /* sub>=0 means an addmul w += x*y, sub<0 means a submul w -= x*y. The signs of w, x and y are fully accounted for by each flipping "sub". The sign of w is retained for the result, unless the absolute value submul underflows, in which case it flips. */ static void __gmpz_aorsmul __GMP_PROTO ((REGPARM_3_1 (mpz_ptr w, mpz_srcptr x, mpz_srcptr y, mp_size_t sub))) REGPARM_ATTR (1); #define mpz_aorsmul(w,x,y,sub) __gmpz_aorsmul (REGPARM_3_1 (w, x, y, sub)) REGPARM_ATTR (1) static void mpz_aorsmul (mpz_ptr w, mpz_srcptr x, mpz_srcptr y, mp_size_t sub) { mp_size_t xsize, ysize, tsize, wsize, wsize_signed; mp_ptr wp, tp; mp_limb_t c, high; TMP_DECL; /* w unaffected if x==0 or y==0 */ xsize = SIZ(x); ysize = SIZ(y); if (xsize == 0 || ysize == 0) return; /* make x the bigger of the two */ if (ABS(ysize) > ABS(xsize)) { MPZ_SRCPTR_SWAP (x, y); MP_SIZE_T_SWAP (xsize, ysize); } sub ^= ysize; ysize = ABS(ysize); /* use mpn_addmul_1/mpn_submul_1 if possible */ if (ysize == 1) { mpz_aorsmul_1 (w, x, PTR(y)[0], sub); return; } sub ^= xsize; xsize = ABS(xsize); wsize_signed = SIZ(w); sub ^= wsize_signed; wsize = ABS(wsize_signed); tsize = xsize + ysize; MPZ_REALLOC (w, MAX (wsize, tsize) + 1); wp = PTR(w); if (wsize_signed == 0) { /* Nothing to add to, just set w=x*y. No w==x or w==y overlap here, since we know x,y!=0 but w==0. */ high = mpn_mul (wp, PTR(x),xsize, PTR(y),ysize); tsize -= (high == 0); SIZ(w) = (sub >= 0 ? tsize : -tsize); return; } TMP_MARK; tp = TMP_ALLOC_LIMBS (tsize); high = mpn_mul (tp, PTR(x),xsize, PTR(y),ysize); tsize -= (high == 0); ASSERT (tp[tsize-1] != 0); if (sub >= 0) { mp_srcptr up = wp; mp_size_t usize = wsize; if (usize < tsize) { up = tp; usize = tsize; tp = wp; tsize = wsize; wsize = usize; } c = mpn_add (wp, up,usize, tp,tsize); wp[wsize] = c; wsize += (c != 0); } else { mp_srcptr up = wp; mp_size_t usize = wsize; if (mpn_cmp_twosizes_lt (up,usize, tp,tsize)) { up = tp; usize = tsize; tp = wp; tsize = wsize; wsize = usize; wsize_signed = -wsize_signed; } ASSERT_NOCARRY (mpn_sub (wp, up,usize, tp,tsize)); wsize = usize; MPN_NORMALIZE (wp, wsize); } SIZ(w) = (wsize_signed >= 0 ? wsize : -wsize); TMP_FREE; } void mpz_addmul (mpz_ptr w, mpz_srcptr u, mpz_srcptr v) { mpz_aorsmul (w, u, v, (mp_size_t) 0); } void mpz_submul (mpz_ptr w, mpz_srcptr u, mpz_srcptr v) { mpz_aorsmul (w, u, v, (mp_size_t) -1); }