1 /* mpn_toom44_mul -- Multiply {ap,an} and {bp,bn} where an and bn are close in
2 size. Or more accurately, bn <= an < (4/3)bn.
3
4 Contributed to the GNU project by Torbjorn Granlund and Marco Bodrato.
5
6 THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE. IT IS ONLY
7 SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST
8 GUARANTEED THAT IT WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
9
10 Copyright 2006-2008, 2013 Free Software Foundation, Inc.
11
12 This file is part of the GNU MP Library.
13
14 The GNU MP Library is free software; you can redistribute it and/or modify
15 it under the terms of either:
16
17 * the GNU Lesser General Public License as published by the Free
18 Software Foundation; either version 3 of the License, or (at your
19 option) any later version.
20
21 or
22
23 * the GNU General Public License as published by the Free Software
24 Foundation; either version 2 of the License, or (at your option) any
25 later version.
26
27 or both in parallel, as here.
28
29 The GNU MP Library is distributed in the hope that it will be useful, but
30 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
31 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
32 for more details.
33
34 You should have received copies of the GNU General Public License and the
35 GNU Lesser General Public License along with the GNU MP Library. If not,
36 see https://www.gnu.org/licenses/. */
37
38
39 #include "gmp-impl.h"
40
41 /* Evaluate in: 0, +1, -1, +2, -2, 1/2, +inf
42
43 <-s--><--n--><--n--><--n-->
44 ____ ______ ______ ______
45 |_a3_|___a2_|___a1_|___a0_|
46 |b3_|___b2_|___b1_|___b0_|
47 <-t-><--n--><--n--><--n-->
48
49 v0 = a0 * b0 # A(0)*B(0)
50 v1 = ( a0+ a1+ a2+ a3)*( b0+ b1+ b2+ b3) # A(1)*B(1) ah <= 3 bh <= 3
51 vm1 = ( a0- a1+ a2- a3)*( b0- b1+ b2- b3) # A(-1)*B(-1) |ah| <= 1 |bh| <= 1
52 v2 = ( a0+2a1+4a2+8a3)*( b0+2b1+4b2+8b3) # A(2)*B(2) ah <= 14 bh <= 14
53 vm2 = ( a0-2a1+4a2-8a3)*( b0-2b1+4b2-8b3) # A(-2)*B(-2) |ah| <= 9 |bh| <= 9
54 vh = (8a0+4a1+2a2+ a3)*(8b0+4b1+2b2+ b3) # A(1/2)*B(1/2) ah <= 14 bh <= 14
55 vinf= a3 * b2 # A(inf)*B(inf)
56 */
57
58 #if TUNE_PROGRAM_BUILD
59 #define MAYBE_mul_basecase 1
60 #define MAYBE_mul_toom22 1
61 #define MAYBE_mul_toom44 1
62 #else
63 #define MAYBE_mul_basecase \
64 (MUL_TOOM44_THRESHOLD < 4 * MUL_TOOM22_THRESHOLD)
65 #define MAYBE_mul_toom22 \
66 (MUL_TOOM44_THRESHOLD < 4 * MUL_TOOM33_THRESHOLD)
67 #define MAYBE_mul_toom44 \
68 (MUL_TOOM6H_THRESHOLD >= 4 * MUL_TOOM44_THRESHOLD)
69 #endif
70
71 #define TOOM44_MUL_N_REC(p, a, b, n, ws) \
72 do { \
73 if (MAYBE_mul_basecase \
74 && BELOW_THRESHOLD (n, MUL_TOOM22_THRESHOLD)) \
75 mpn_mul_basecase (p, a, n, b, n); \
76 else if (MAYBE_mul_toom22 \
77 && BELOW_THRESHOLD (n, MUL_TOOM33_THRESHOLD)) \
78 mpn_toom22_mul (p, a, n, b, n, ws); \
79 else if (! MAYBE_mul_toom44 \
80 || BELOW_THRESHOLD (n, MUL_TOOM44_THRESHOLD)) \
81 mpn_toom33_mul (p, a, n, b, n, ws); \
82 else \
83 mpn_toom44_mul (p, a, n, b, n, ws); \
84 } while (0)
85
86 /* Use of scratch space. In the product area, we store
87
88 ___________________
89 |vinf|____|_v1_|_v0_|
90 s+t 2n-1 2n+1 2n
91
92 The other recursive products, vm1, v2, vm2, vh are stored in the
93 scratch area. When computing them, we use the product area for
94 intermediate values.
95
96 Next, we compute v1. We can store the intermediate factors at v0
97 and at vh + 2n + 2.
98
99 Finally, for v0 and vinf, factors are parts of the input operands,
100 and we need scratch space only for the recursive multiplication.
101
102 In all, if S(an) is the scratch need, the needed space is bounded by
103
104 S(an) <= 4 (2*ceil(an/4) + 1) + 1 + S(ceil(an/4) + 1)
105
106 which should give S(n) = 8 n/3 + c log(n) for some constant c.
107 */
108
109 void
110 mpn_toom44_mul (mp_ptr pp,
111 mp_srcptr ap, mp_size_t an,
112 mp_srcptr bp, mp_size_t bn,
113 mp_ptr scratch)
114 {
115 mp_size_t n, s, t;
116 mp_limb_t cy;
117 enum toom7_flags flags;
118
119 #define a0 ap
120 #define a1 (ap + n)
121 #define a2 (ap + 2*n)
122 #define a3 (ap + 3*n)
123 #define b0 bp
124 #define b1 (bp + n)
125 #define b2 (bp + 2*n)
126 #define b3 (bp + 3*n)
127
128 ASSERT (an >= bn);
129
130 n = (an + 3) >> 2;
131
132 s = an - 3 * n;
133 t = bn - 3 * n;
134
135 ASSERT (0 < s && s <= n);
136 ASSERT (0 < t && t <= n);
137 ASSERT (s >= t);
138
139 /* NOTE: The multiplications to v2, vm2, vh and vm1 overwrites the
140 * following limb, so these must be computed in order, and we need a
141 * one limb gap to tp. */
142 #define v0 pp /* 2n */
143 #define v1 (pp + 2 * n) /* 2n+1 */
144 #define vinf (pp + 6 * n) /* s+t */
145 #define v2 scratch /* 2n+1 */
146 #define vm2 (scratch + 2 * n + 1) /* 2n+1 */
147 #define vh (scratch + 4 * n + 2) /* 2n+1 */
148 #define vm1 (scratch + 6 * n + 3) /* 2n+1 */
149 #define tp (scratch + 8*n + 5)
150
151 /* apx and bpx must not overlap with v1 */
152 #define apx pp /* n+1 */
153 #define amx (pp + n + 1) /* n+1 */
154 #define bmx (pp + 2*n + 2) /* n+1 */
155 #define bpx (pp + 4*n + 2) /* n+1 */
156
157 /* Total scratch need: 8*n + 5 + scratch for recursive calls. This
158 gives roughly 32 n/3 + log term. */
159
160 /* Compute apx = a0 + 2 a1 + 4 a2 + 8 a3 and amx = a0 - 2 a1 + 4 a2 - 8 a3. */
161 flags = (enum toom7_flags) (toom7_w1_neg & mpn_toom_eval_dgr3_pm2 (apx, amx, ap, n, s, tp));
162
163 /* Compute bpx = b0 + 2 b1 + 4 b2 + 8 b3 and bmx = b0 - 2 b1 + 4 b2 - 8 b3. */
164 flags = (enum toom7_flags) (flags ^ (toom7_w1_neg & mpn_toom_eval_dgr3_pm2 (bpx, bmx, bp, n, t, tp)));
165
166 TOOM44_MUL_N_REC (v2, apx, bpx, n + 1, tp); /* v2, 2n+1 limbs */
167 TOOM44_MUL_N_REC (vm2, amx, bmx, n + 1, tp); /* vm2, 2n+1 limbs */
168
169 /* Compute apx = 8 a0 + 4 a1 + 2 a2 + a3 = (((2*a0 + a1) * 2 + a2) * 2 + a3 */
170 #if HAVE_NATIVE_mpn_addlsh1_n
171 cy = mpn_addlsh1_n (apx, a1, a0, n);
172 cy = 2*cy + mpn_addlsh1_n (apx, a2, apx, n);
173 if (s < n)
174 {
175 mp_limb_t cy2;
176 cy2 = mpn_addlsh1_n (apx, a3, apx, s);
177 apx[n] = 2*cy + mpn_lshift (apx + s, apx + s, n - s, 1);
178 MPN_INCR_U (apx + s, n+1-s, cy2);
179 }
180 else
181 apx[n] = 2*cy + mpn_addlsh1_n (apx, a3, apx, n);
182 #else
183 cy = mpn_lshift (apx, a0, n, 1);
184 cy += mpn_add_n (apx, apx, a1, n);
185 cy = 2*cy + mpn_lshift (apx, apx, n, 1);
186 cy += mpn_add_n (apx, apx, a2, n);
187 cy = 2*cy + mpn_lshift (apx, apx, n, 1);
188 apx[n] = cy + mpn_add (apx, apx, n, a3, s);
189 #endif
190
191 /* Compute bpx = 8 b0 + 4 b1 + 2 b2 + b3 = (((2*b0 + b1) * 2 + b2) * 2 + b3 */
192 #if HAVE_NATIVE_mpn_addlsh1_n
193 cy = mpn_addlsh1_n (bpx, b1, b0, n);
194 cy = 2*cy + mpn_addlsh1_n (bpx, b2, bpx, n);
195 if (t < n)
196 {
197 mp_limb_t cy2;
198 cy2 = mpn_addlsh1_n (bpx, b3, bpx, t);
199 bpx[n] = 2*cy + mpn_lshift (bpx + t, bpx + t, n - t, 1);
200 MPN_INCR_U (bpx + t, n+1-t, cy2);
201 }
202 else
203 bpx[n] = 2*cy + mpn_addlsh1_n (bpx, b3, bpx, n);
204 #else
205 cy = mpn_lshift (bpx, b0, n, 1);
206 cy += mpn_add_n (bpx, bpx, b1, n);
207 cy = 2*cy + mpn_lshift (bpx, bpx, n, 1);
208 cy += mpn_add_n (bpx, bpx, b2, n);
209 cy = 2*cy + mpn_lshift (bpx, bpx, n, 1);
210 bpx[n] = cy + mpn_add (bpx, bpx, n, b3, t);
211 #endif
212
213 ASSERT (apx[n] < 15);
214 ASSERT (bpx[n] < 15);
215
216 TOOM44_MUL_N_REC (vh, apx, bpx, n + 1, tp); /* vh, 2n+1 limbs */
217
218 /* Compute apx = a0 + a1 + a2 + a3 and amx = a0 - a1 + a2 - a3. */
219 flags = (enum toom7_flags) (flags | (toom7_w3_neg & mpn_toom_eval_dgr3_pm1 (apx, amx, ap, n, s, tp)));
220
221 /* Compute bpx = b0 + b1 + b2 + b3 and bmx = b0 - b1 + b2 - b3. */
222 flags = (enum toom7_flags) (flags ^ (toom7_w3_neg & mpn_toom_eval_dgr3_pm1 (bpx, bmx, bp, n, t, tp)));
223
224 ASSERT (amx[n] <= 1);
225 ASSERT (bmx[n] <= 1);
226
227 vm1 [2 * n] = 0;
228 TOOM44_MUL_N_REC (vm1, amx, bmx, n + (bmx[n] | amx[n]), tp); /* vm1, 2n+1 limbs */
229 /* Clobbers amx, bmx. */
230 TOOM44_MUL_N_REC (v1, apx, bpx, n + 1, tp); /* v1, 2n+1 limbs */
231
232 TOOM44_MUL_N_REC (v0, a0, b0, n, tp);
233 if (s > t)
234 mpn_mul (vinf, a3, s, b3, t);
235 else
236 TOOM44_MUL_N_REC (vinf, a3, b3, s, tp); /* vinf, s+t limbs */
237
238 mpn_toom_interpolate_7pts (pp, n, flags, vm2, vm1, v2, vh, s + t, tp);
239 }