1*86d7f5d3SJohn Marino /* mpn_toom33_mul -- Multiply {ap,an} and {p,bn} where an and bn are close in
2*86d7f5d3SJohn Marino size. Or more accurately, bn <= an < (3/2)bn.
3*86d7f5d3SJohn Marino
4*86d7f5d3SJohn Marino Contributed to the GNU project by Torbjorn Granlund.
5*86d7f5d3SJohn Marino Additional improvements by Marco Bodrato.
6*86d7f5d3SJohn Marino
7*86d7f5d3SJohn Marino THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE. IT IS ONLY
8*86d7f5d3SJohn Marino SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST
9*86d7f5d3SJohn Marino GUARANTEED THAT IT WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
10*86d7f5d3SJohn Marino
11*86d7f5d3SJohn Marino Copyright 2006, 2007, 2008, 2010 Free Software Foundation, Inc.
12*86d7f5d3SJohn Marino
13*86d7f5d3SJohn Marino This file is part of the GNU MP Library.
14*86d7f5d3SJohn Marino
15*86d7f5d3SJohn Marino The GNU MP Library is free software; you can redistribute it and/or modify
16*86d7f5d3SJohn Marino it under the terms of the GNU Lesser General Public License as published by
17*86d7f5d3SJohn Marino the Free Software Foundation; either version 3 of the License, or (at your
18*86d7f5d3SJohn Marino option) any later version.
19*86d7f5d3SJohn Marino
20*86d7f5d3SJohn Marino The GNU MP Library is distributed in the hope that it will be useful, but
21*86d7f5d3SJohn Marino WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
22*86d7f5d3SJohn Marino or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
23*86d7f5d3SJohn Marino License for more details.
24*86d7f5d3SJohn Marino
25*86d7f5d3SJohn Marino You should have received a copy of the GNU Lesser General Public License
26*86d7f5d3SJohn Marino along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
27*86d7f5d3SJohn Marino
28*86d7f5d3SJohn Marino
29*86d7f5d3SJohn Marino #include "gmp.h"
30*86d7f5d3SJohn Marino #include "gmp-impl.h"
31*86d7f5d3SJohn Marino
32*86d7f5d3SJohn Marino /* Evaluate in: -1, 0, +1, +2, +inf
33*86d7f5d3SJohn Marino
34*86d7f5d3SJohn Marino <-s--><--n--><--n--><--n-->
35*86d7f5d3SJohn Marino ____ ______ ______ ______
36*86d7f5d3SJohn Marino |_a3_|___a2_|___a1_|___a0_|
37*86d7f5d3SJohn Marino |b3_|___b2_|___b1_|___b0_|
38*86d7f5d3SJohn Marino <-t-><--n--><--n--><--n-->
39*86d7f5d3SJohn Marino
40*86d7f5d3SJohn Marino v0 = a0 * b0 # A(0)*B(0)
41*86d7f5d3SJohn Marino v1 = (a0+ a1+ a2)*(b0+ b1+ b2) # A(1)*B(1) ah <= 2 bh <= 2
42*86d7f5d3SJohn Marino vm1 = (a0- a1+ a2)*(b0- b1+ b2) # A(-1)*B(-1) |ah| <= 1 bh <= 1
43*86d7f5d3SJohn Marino v2 = (a0+2a1+4a2)*(b0+2b1+4b2) # A(2)*B(2) ah <= 6 bh <= 6
44*86d7f5d3SJohn Marino vinf= a2 * b2 # A(inf)*B(inf)
45*86d7f5d3SJohn Marino */
46*86d7f5d3SJohn Marino
47*86d7f5d3SJohn Marino #if TUNE_PROGRAM_BUILD
48*86d7f5d3SJohn Marino #define MAYBE_mul_basecase 1
49*86d7f5d3SJohn Marino #define MAYBE_mul_toom33 1
50*86d7f5d3SJohn Marino #else
51*86d7f5d3SJohn Marino #define MAYBE_mul_basecase \
52*86d7f5d3SJohn Marino (MUL_TOOM33_THRESHOLD < 3 * MUL_TOOM22_THRESHOLD)
53*86d7f5d3SJohn Marino #define MAYBE_mul_toom33 \
54*86d7f5d3SJohn Marino (MUL_TOOM44_THRESHOLD >= 3 * MUL_TOOM33_THRESHOLD)
55*86d7f5d3SJohn Marino #endif
56*86d7f5d3SJohn Marino
57*86d7f5d3SJohn Marino /* FIXME: TOOM33_MUL_N_REC is not quite right for a balanced
58*86d7f5d3SJohn Marino multiplication at the infinity point. We may have
59*86d7f5d3SJohn Marino MAYBE_mul_basecase == 0, and still get s just below
60*86d7f5d3SJohn Marino MUL_TOOM22_THRESHOLD. If MUL_TOOM33_THRESHOLD == 7, we can even get
61*86d7f5d3SJohn Marino s == 1 and mpn_toom22_mul will crash.
62*86d7f5d3SJohn Marino */
63*86d7f5d3SJohn Marino
64*86d7f5d3SJohn Marino #define TOOM33_MUL_N_REC(p, a, b, n, ws) \
65*86d7f5d3SJohn Marino do { \
66*86d7f5d3SJohn Marino if (MAYBE_mul_basecase \
67*86d7f5d3SJohn Marino && BELOW_THRESHOLD (n, MUL_TOOM22_THRESHOLD)) \
68*86d7f5d3SJohn Marino mpn_mul_basecase (p, a, n, b, n); \
69*86d7f5d3SJohn Marino else if (! MAYBE_mul_toom33 \
70*86d7f5d3SJohn Marino || BELOW_THRESHOLD (n, MUL_TOOM33_THRESHOLD)) \
71*86d7f5d3SJohn Marino mpn_toom22_mul (p, a, n, b, n, ws); \
72*86d7f5d3SJohn Marino else \
73*86d7f5d3SJohn Marino mpn_toom33_mul (p, a, n, b, n, ws); \
74*86d7f5d3SJohn Marino } while (0)
75*86d7f5d3SJohn Marino
76*86d7f5d3SJohn Marino void
mpn_toom33_mul(mp_ptr pp,mp_srcptr ap,mp_size_t an,mp_srcptr bp,mp_size_t bn,mp_ptr scratch)77*86d7f5d3SJohn Marino mpn_toom33_mul (mp_ptr pp,
78*86d7f5d3SJohn Marino mp_srcptr ap, mp_size_t an,
79*86d7f5d3SJohn Marino mp_srcptr bp, mp_size_t bn,
80*86d7f5d3SJohn Marino mp_ptr scratch)
81*86d7f5d3SJohn Marino {
82*86d7f5d3SJohn Marino mp_size_t n, s, t;
83*86d7f5d3SJohn Marino int vm1_neg;
84*86d7f5d3SJohn Marino mp_limb_t cy, vinf0;
85*86d7f5d3SJohn Marino mp_ptr gp;
86*86d7f5d3SJohn Marino mp_ptr as1, asm1, as2;
87*86d7f5d3SJohn Marino mp_ptr bs1, bsm1, bs2;
88*86d7f5d3SJohn Marino
89*86d7f5d3SJohn Marino #define a0 ap
90*86d7f5d3SJohn Marino #define a1 (ap + n)
91*86d7f5d3SJohn Marino #define a2 (ap + 2*n)
92*86d7f5d3SJohn Marino #define b0 bp
93*86d7f5d3SJohn Marino #define b1 (bp + n)
94*86d7f5d3SJohn Marino #define b2 (bp + 2*n)
95*86d7f5d3SJohn Marino
96*86d7f5d3SJohn Marino n = (an + 2) / (size_t) 3;
97*86d7f5d3SJohn Marino
98*86d7f5d3SJohn Marino s = an - 2 * n;
99*86d7f5d3SJohn Marino t = bn - 2 * n;
100*86d7f5d3SJohn Marino
101*86d7f5d3SJohn Marino ASSERT (an >= bn);
102*86d7f5d3SJohn Marino
103*86d7f5d3SJohn Marino ASSERT (0 < s && s <= n);
104*86d7f5d3SJohn Marino ASSERT (0 < t && t <= n);
105*86d7f5d3SJohn Marino
106*86d7f5d3SJohn Marino as1 = scratch + 4 * n + 4;
107*86d7f5d3SJohn Marino asm1 = scratch + 2 * n + 2;
108*86d7f5d3SJohn Marino as2 = pp + n + 1;
109*86d7f5d3SJohn Marino
110*86d7f5d3SJohn Marino bs1 = pp;
111*86d7f5d3SJohn Marino bsm1 = scratch + 3 * n + 3; /* we need 4n+4 <= 4n+s+t */
112*86d7f5d3SJohn Marino bs2 = pp + 2 * n + 2;
113*86d7f5d3SJohn Marino
114*86d7f5d3SJohn Marino gp = scratch;
115*86d7f5d3SJohn Marino
116*86d7f5d3SJohn Marino vm1_neg = 0;
117*86d7f5d3SJohn Marino
118*86d7f5d3SJohn Marino /* Compute as1 and asm1. */
119*86d7f5d3SJohn Marino cy = mpn_add (gp, a0, n, a2, s);
120*86d7f5d3SJohn Marino #if HAVE_NATIVE_mpn_add_n_sub_n
121*86d7f5d3SJohn Marino if (cy == 0 && mpn_cmp (gp, a1, n) < 0)
122*86d7f5d3SJohn Marino {
123*86d7f5d3SJohn Marino cy = mpn_add_n_sub_n (as1, asm1, a1, gp, n);
124*86d7f5d3SJohn Marino as1[n] = cy >> 1;
125*86d7f5d3SJohn Marino asm1[n] = 0;
126*86d7f5d3SJohn Marino vm1_neg = 1;
127*86d7f5d3SJohn Marino }
128*86d7f5d3SJohn Marino else
129*86d7f5d3SJohn Marino {
130*86d7f5d3SJohn Marino mp_limb_t cy2;
131*86d7f5d3SJohn Marino cy2 = mpn_add_n_sub_n (as1, asm1, gp, a1, n);
132*86d7f5d3SJohn Marino as1[n] = cy + (cy2 >> 1);
133*86d7f5d3SJohn Marino asm1[n] = cy - (cy2 & 1);
134*86d7f5d3SJohn Marino }
135*86d7f5d3SJohn Marino #else
136*86d7f5d3SJohn Marino as1[n] = cy + mpn_add_n (as1, gp, a1, n);
137*86d7f5d3SJohn Marino if (cy == 0 && mpn_cmp (gp, a1, n) < 0)
138*86d7f5d3SJohn Marino {
139*86d7f5d3SJohn Marino mpn_sub_n (asm1, a1, gp, n);
140*86d7f5d3SJohn Marino asm1[n] = 0;
141*86d7f5d3SJohn Marino vm1_neg = 1;
142*86d7f5d3SJohn Marino }
143*86d7f5d3SJohn Marino else
144*86d7f5d3SJohn Marino {
145*86d7f5d3SJohn Marino cy -= mpn_sub_n (asm1, gp, a1, n);
146*86d7f5d3SJohn Marino asm1[n] = cy;
147*86d7f5d3SJohn Marino }
148*86d7f5d3SJohn Marino #endif
149*86d7f5d3SJohn Marino
150*86d7f5d3SJohn Marino /* Compute as2. */
151*86d7f5d3SJohn Marino #if HAVE_NATIVE_mpn_rsblsh1_n
152*86d7f5d3SJohn Marino cy = mpn_add_n (as2, a2, as1, s);
153*86d7f5d3SJohn Marino if (s != n)
154*86d7f5d3SJohn Marino cy = mpn_add_1 (as2 + s, as1 + s, n - s, cy);
155*86d7f5d3SJohn Marino cy += as1[n];
156*86d7f5d3SJohn Marino cy = 2 * cy + mpn_rsblsh1_n (as2, a0, as2, n);
157*86d7f5d3SJohn Marino #else
158*86d7f5d3SJohn Marino #if HAVE_NATIVE_mpn_addlsh1_n
159*86d7f5d3SJohn Marino cy = mpn_addlsh1_n (as2, a1, a2, s);
160*86d7f5d3SJohn Marino if (s != n)
161*86d7f5d3SJohn Marino cy = mpn_add_1 (as2 + s, a1 + s, n - s, cy);
162*86d7f5d3SJohn Marino cy = 2 * cy + mpn_addlsh1_n (as2, a0, as2, n);
163*86d7f5d3SJohn Marino #else
164*86d7f5d3SJohn Marino cy = mpn_add_n (as2, a2, as1, s);
165*86d7f5d3SJohn Marino if (s != n)
166*86d7f5d3SJohn Marino cy = mpn_add_1 (as2 + s, as1 + s, n - s, cy);
167*86d7f5d3SJohn Marino cy += as1[n];
168*86d7f5d3SJohn Marino cy = 2 * cy + mpn_lshift (as2, as2, n, 1);
169*86d7f5d3SJohn Marino cy -= mpn_sub_n (as2, as2, a0, n);
170*86d7f5d3SJohn Marino #endif
171*86d7f5d3SJohn Marino #endif
172*86d7f5d3SJohn Marino as2[n] = cy;
173*86d7f5d3SJohn Marino
174*86d7f5d3SJohn Marino /* Compute bs1 and bsm1. */
175*86d7f5d3SJohn Marino cy = mpn_add (gp, b0, n, b2, t);
176*86d7f5d3SJohn Marino #if HAVE_NATIVE_mpn_add_n_sub_n
177*86d7f5d3SJohn Marino if (cy == 0 && mpn_cmp (gp, b1, n) < 0)
178*86d7f5d3SJohn Marino {
179*86d7f5d3SJohn Marino cy = mpn_add_n_sub_n (bs1, bsm1, b1, gp, n);
180*86d7f5d3SJohn Marino bs1[n] = cy >> 1;
181*86d7f5d3SJohn Marino bsm1[n] = 0;
182*86d7f5d3SJohn Marino vm1_neg ^= 1;
183*86d7f5d3SJohn Marino }
184*86d7f5d3SJohn Marino else
185*86d7f5d3SJohn Marino {
186*86d7f5d3SJohn Marino mp_limb_t cy2;
187*86d7f5d3SJohn Marino cy2 = mpn_add_n_sub_n (bs1, bsm1, gp, b1, n);
188*86d7f5d3SJohn Marino bs1[n] = cy + (cy2 >> 1);
189*86d7f5d3SJohn Marino bsm1[n] = cy - (cy2 & 1);
190*86d7f5d3SJohn Marino }
191*86d7f5d3SJohn Marino #else
192*86d7f5d3SJohn Marino bs1[n] = cy + mpn_add_n (bs1, gp, b1, n);
193*86d7f5d3SJohn Marino if (cy == 0 && mpn_cmp (gp, b1, n) < 0)
194*86d7f5d3SJohn Marino {
195*86d7f5d3SJohn Marino mpn_sub_n (bsm1, b1, gp, n);
196*86d7f5d3SJohn Marino bsm1[n] = 0;
197*86d7f5d3SJohn Marino vm1_neg ^= 1;
198*86d7f5d3SJohn Marino }
199*86d7f5d3SJohn Marino else
200*86d7f5d3SJohn Marino {
201*86d7f5d3SJohn Marino cy -= mpn_sub_n (bsm1, gp, b1, n);
202*86d7f5d3SJohn Marino bsm1[n] = cy;
203*86d7f5d3SJohn Marino }
204*86d7f5d3SJohn Marino #endif
205*86d7f5d3SJohn Marino
206*86d7f5d3SJohn Marino /* Compute bs2. */
207*86d7f5d3SJohn Marino #if HAVE_NATIVE_mpn_rsblsh1_n
208*86d7f5d3SJohn Marino cy = mpn_add_n (bs2, b2, bs1, t);
209*86d7f5d3SJohn Marino if (t != n)
210*86d7f5d3SJohn Marino cy = mpn_add_1 (bs2 + t, bs1 + t, n - t, cy);
211*86d7f5d3SJohn Marino cy += bs1[n];
212*86d7f5d3SJohn Marino cy = 2 * cy + mpn_rsblsh1_n (bs2, b0, bs2, n);
213*86d7f5d3SJohn Marino #else
214*86d7f5d3SJohn Marino #if HAVE_NATIVE_mpn_addlsh1_n
215*86d7f5d3SJohn Marino cy = mpn_addlsh1_n (bs2, b1, b2, t);
216*86d7f5d3SJohn Marino if (t != n)
217*86d7f5d3SJohn Marino cy = mpn_add_1 (bs2 + t, b1 + t, n - t, cy);
218*86d7f5d3SJohn Marino cy = 2 * cy + mpn_addlsh1_n (bs2, b0, bs2, n);
219*86d7f5d3SJohn Marino #else
220*86d7f5d3SJohn Marino cy = mpn_add_n (bs2, bs1, b2, t);
221*86d7f5d3SJohn Marino if (t != n)
222*86d7f5d3SJohn Marino cy = mpn_add_1 (bs2 + t, bs1 + t, n - t, cy);
223*86d7f5d3SJohn Marino cy += bs1[n];
224*86d7f5d3SJohn Marino cy = 2 * cy + mpn_lshift (bs2, bs2, n, 1);
225*86d7f5d3SJohn Marino cy -= mpn_sub_n (bs2, bs2, b0, n);
226*86d7f5d3SJohn Marino #endif
227*86d7f5d3SJohn Marino #endif
228*86d7f5d3SJohn Marino bs2[n] = cy;
229*86d7f5d3SJohn Marino
230*86d7f5d3SJohn Marino ASSERT (as1[n] <= 2);
231*86d7f5d3SJohn Marino ASSERT (bs1[n] <= 2);
232*86d7f5d3SJohn Marino ASSERT (asm1[n] <= 1);
233*86d7f5d3SJohn Marino ASSERT (bsm1[n] <= 1);
234*86d7f5d3SJohn Marino ASSERT (as2[n] <= 6);
235*86d7f5d3SJohn Marino ASSERT (bs2[n] <= 6);
236*86d7f5d3SJohn Marino
237*86d7f5d3SJohn Marino #define v0 pp /* 2n */
238*86d7f5d3SJohn Marino #define v1 (pp + 2 * n) /* 2n+1 */
239*86d7f5d3SJohn Marino #define vinf (pp + 4 * n) /* s+t */
240*86d7f5d3SJohn Marino #define vm1 scratch /* 2n+1 */
241*86d7f5d3SJohn Marino #define v2 (scratch + 2 * n + 1) /* 2n+2 */
242*86d7f5d3SJohn Marino #define scratch_out (scratch + 5 * n + 5)
243*86d7f5d3SJohn Marino
244*86d7f5d3SJohn Marino /* vm1, 2n+1 limbs */
245*86d7f5d3SJohn Marino #ifdef SMALLER_RECURSION
246*86d7f5d3SJohn Marino TOOM33_MUL_N_REC (vm1, asm1, bsm1, n, scratch_out);
247*86d7f5d3SJohn Marino cy = 0;
248*86d7f5d3SJohn Marino if (asm1[n] != 0)
249*86d7f5d3SJohn Marino cy = bsm1[n] + mpn_add_n (vm1 + n, vm1 + n, bsm1, n);
250*86d7f5d3SJohn Marino if (bsm1[n] != 0)
251*86d7f5d3SJohn Marino cy += mpn_add_n (vm1 + n, vm1 + n, asm1, n);
252*86d7f5d3SJohn Marino vm1[2 * n] = cy;
253*86d7f5d3SJohn Marino #else
254*86d7f5d3SJohn Marino TOOM33_MUL_N_REC (vm1, asm1, bsm1, n + 1, scratch_out);
255*86d7f5d3SJohn Marino #endif
256*86d7f5d3SJohn Marino
257*86d7f5d3SJohn Marino TOOM33_MUL_N_REC (v2, as2, bs2, n + 1, scratch_out); /* v2, 2n+1 limbs */
258*86d7f5d3SJohn Marino
259*86d7f5d3SJohn Marino /* vinf, s+t limbs */
260*86d7f5d3SJohn Marino if (s > t) mpn_mul (vinf, a2, s, b2, t);
261*86d7f5d3SJohn Marino else TOOM33_MUL_N_REC (vinf, a2, b2, s, scratch_out);
262*86d7f5d3SJohn Marino
263*86d7f5d3SJohn Marino vinf0 = vinf[0]; /* v1 overlaps with this */
264*86d7f5d3SJohn Marino
265*86d7f5d3SJohn Marino #ifdef SMALLER_RECURSION
266*86d7f5d3SJohn Marino /* v1, 2n+1 limbs */
267*86d7f5d3SJohn Marino TOOM33_MUL_N_REC (v1, as1, bs1, n, scratch_out);
268*86d7f5d3SJohn Marino if (as1[n] == 1)
269*86d7f5d3SJohn Marino {
270*86d7f5d3SJohn Marino cy = bs1[n] + mpn_add_n (v1 + n, v1 + n, bs1, n);
271*86d7f5d3SJohn Marino }
272*86d7f5d3SJohn Marino else if (as1[n] != 0)
273*86d7f5d3SJohn Marino {
274*86d7f5d3SJohn Marino #if HAVE_NATIVE_mpn_addlsh1_n
275*86d7f5d3SJohn Marino cy = 2 * bs1[n] + mpn_addlsh1_n (v1 + n, v1 + n, bs1, n);
276*86d7f5d3SJohn Marino #else
277*86d7f5d3SJohn Marino cy = 2 * bs1[n] + mpn_addmul_1 (v1 + n, bs1, n, CNST_LIMB(2));
278*86d7f5d3SJohn Marino #endif
279*86d7f5d3SJohn Marino }
280*86d7f5d3SJohn Marino else
281*86d7f5d3SJohn Marino cy = 0;
282*86d7f5d3SJohn Marino if (bs1[n] == 1)
283*86d7f5d3SJohn Marino {
284*86d7f5d3SJohn Marino cy += mpn_add_n (v1 + n, v1 + n, as1, n);
285*86d7f5d3SJohn Marino }
286*86d7f5d3SJohn Marino else if (bs1[n] != 0)
287*86d7f5d3SJohn Marino {
288*86d7f5d3SJohn Marino #if HAVE_NATIVE_mpn_addlsh1_n
289*86d7f5d3SJohn Marino cy += mpn_addlsh1_n (v1 + n, v1 + n, as1, n);
290*86d7f5d3SJohn Marino #else
291*86d7f5d3SJohn Marino cy += mpn_addmul_1 (v1 + n, as1, n, CNST_LIMB(2));
292*86d7f5d3SJohn Marino #endif
293*86d7f5d3SJohn Marino }
294*86d7f5d3SJohn Marino v1[2 * n] = cy;
295*86d7f5d3SJohn Marino #else
296*86d7f5d3SJohn Marino cy = vinf[1];
297*86d7f5d3SJohn Marino TOOM33_MUL_N_REC (v1, as1, bs1, n + 1, scratch_out);
298*86d7f5d3SJohn Marino vinf[1] = cy;
299*86d7f5d3SJohn Marino #endif
300*86d7f5d3SJohn Marino
301*86d7f5d3SJohn Marino TOOM33_MUL_N_REC (v0, ap, bp, n, scratch_out); /* v0, 2n limbs */
302*86d7f5d3SJohn Marino
303*86d7f5d3SJohn Marino mpn_toom_interpolate_5pts (pp, v2, vm1, n, s + t, vm1_neg, vinf0);
304*86d7f5d3SJohn Marino }
305