1 /* mpz_addmul, mpz_submul -- add or subtract multiple. 2 3 Copyright 2001, 2004, 2005 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 Lesser General Public License as published by 9 the Free Software Foundation; either version 3 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 Lesser General Public 15 License for more details. 16 17 You should have received a copy of the GNU Lesser General Public License 18 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ 19 20 #include "gmp.h" 21 #include "gmp-impl.h" 22 23 24 /* expecting x and y both with non-zero high limbs */ 25 #define mpn_cmp_twosizes_lt(xp,xsize, yp,ysize) \ 26 ((xsize) < (ysize) \ 27 || ((xsize) == (ysize) && mpn_cmp (xp, yp, xsize) < 0)) 28 29 30 /* sub>=0 means an addmul w += x*y, sub<0 means a submul w -= x*y. 31 32 The signs of w, x and y are fully accounted for by each flipping "sub". 33 34 The sign of w is retained for the result, unless the absolute value 35 submul underflows, in which case it flips. */ 36 37 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); 38 #define mpz_aorsmul(w,x,y,sub) __gmpz_aorsmul (REGPARM_3_1 (w, x, y, sub)) 39 40 REGPARM_ATTR (1) static void 41 mpz_aorsmul (mpz_ptr w, mpz_srcptr x, mpz_srcptr y, mp_size_t sub) 42 { 43 mp_size_t xsize, ysize, tsize, wsize, wsize_signed; 44 mp_ptr wp, tp; 45 mp_limb_t c, high; 46 TMP_DECL; 47 48 /* w unaffected if x==0 or y==0 */ 49 xsize = SIZ(x); 50 ysize = SIZ(y); 51 if (xsize == 0 || ysize == 0) 52 return; 53 54 /* make x the bigger of the two */ 55 if (ABS(ysize) > ABS(xsize)) 56 { 57 MPZ_SRCPTR_SWAP (x, y); 58 MP_SIZE_T_SWAP (xsize, ysize); 59 } 60 61 sub ^= ysize; 62 ysize = ABS(ysize); 63 64 /* use mpn_addmul_1/mpn_submul_1 if possible */ 65 if (ysize == 1) 66 { 67 mpz_aorsmul_1 (w, x, PTR(y)[0], sub); 68 return; 69 } 70 71 sub ^= xsize; 72 xsize = ABS(xsize); 73 74 wsize_signed = SIZ(w); 75 sub ^= wsize_signed; 76 wsize = ABS(wsize_signed); 77 78 tsize = xsize + ysize; 79 MPZ_REALLOC (w, MAX (wsize, tsize) + 1); 80 wp = PTR(w); 81 82 if (wsize_signed == 0) 83 { 84 /* Nothing to add to, just set w=x*y. No w==x or w==y overlap here, 85 since we know x,y!=0 but w==0. */ 86 high = mpn_mul (wp, PTR(x),xsize, PTR(y),ysize); 87 tsize -= (high == 0); 88 SIZ(w) = (sub >= 0 ? tsize : -tsize); 89 return; 90 } 91 92 TMP_MARK; 93 tp = TMP_ALLOC_LIMBS (tsize); 94 95 high = mpn_mul (tp, PTR(x),xsize, PTR(y),ysize); 96 tsize -= (high == 0); 97 ASSERT (tp[tsize-1] != 0); 98 if (sub >= 0) 99 { 100 mp_srcptr up = wp; 101 mp_size_t usize = wsize; 102 103 if (usize < tsize) 104 { 105 up = tp; 106 usize = tsize; 107 tp = wp; 108 tsize = wsize; 109 110 wsize = usize; 111 } 112 113 c = mpn_add (wp, up,usize, tp,tsize); 114 wp[wsize] = c; 115 wsize += (c != 0); 116 } 117 else 118 { 119 mp_srcptr up = wp; 120 mp_size_t usize = wsize; 121 122 if (mpn_cmp_twosizes_lt (up,usize, tp,tsize)) 123 { 124 up = tp; 125 usize = tsize; 126 tp = wp; 127 tsize = wsize; 128 129 wsize = usize; 130 wsize_signed = -wsize_signed; 131 } 132 133 ASSERT_NOCARRY (mpn_sub (wp, up,usize, tp,tsize)); 134 wsize = usize; 135 MPN_NORMALIZE (wp, wsize); 136 } 137 138 SIZ(w) = (wsize_signed >= 0 ? wsize : -wsize); 139 140 TMP_FREE; 141 } 142 143 144 void 145 mpz_addmul (mpz_ptr w, mpz_srcptr u, mpz_srcptr v) 146 { 147 mpz_aorsmul (w, u, v, (mp_size_t) 0); 148 } 149 150 void 151 mpz_submul (mpz_ptr w, mpz_srcptr u, mpz_srcptr v) 152 { 153 mpz_aorsmul (w, u, v, (mp_size_t) -1); 154 } 155