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