|
282
|
1 #include <tommath.h> |
|
|
2 #ifdef BN_MP_TOOM_MUL_C |
|
|
3 /* LibTomMath, multiple-precision integer library -- Tom St Denis |
|
|
4 * |
|
|
5 * LibTomMath is a library that provides multiple-precision |
|
|
6 * integer arithmetic as well as number theoretic functionality. |
|
|
7 * |
|
|
8 * The library was designed directly after the MPI library by |
|
|
9 * Michael Fromberger but has been written from scratch with |
|
|
10 * additional optimizations in place. |
|
|
11 * |
|
|
12 * The library is free for all purposes without any express |
|
|
13 * guarantee it works. |
|
|
14 * |
|
386
|
15 * Tom St Denis, [email protected], http://math.libtomcrypt.com |
|
282
|
16 */ |
|
|
17 |
|
|
18 /* multiplication using the Toom-Cook 3-way algorithm |
|
|
19 * |
|
|
20 * Much more complicated than Karatsuba but has a lower |
|
|
21 * asymptotic running time of O(N**1.464). This algorithm is |
|
|
22 * only particularly useful on VERY large inputs |
|
|
23 * (we're talking 1000s of digits here...). |
|
|
24 */ |
|
|
25 int mp_toom_mul(mp_int *a, mp_int *b, mp_int *c) |
|
|
26 { |
|
|
27 mp_int w0, w1, w2, w3, w4, tmp1, tmp2, a0, a1, a2, b0, b1, b2; |
|
|
28 int res, B; |
|
|
29 |
|
|
30 /* init temps */ |
|
|
31 if ((res = mp_init_multi(&w0, &w1, &w2, &w3, &w4, |
|
|
32 &a0, &a1, &a2, &b0, &b1, |
|
|
33 &b2, &tmp1, &tmp2, NULL)) != MP_OKAY) { |
|
|
34 return res; |
|
|
35 } |
|
|
36 |
|
|
37 /* B */ |
|
|
38 B = MIN(a->used, b->used) / 3; |
|
|
39 |
|
|
40 /* a = a2 * B**2 + a1 * B + a0 */ |
|
|
41 if ((res = mp_mod_2d(a, DIGIT_BIT * B, &a0)) != MP_OKAY) { |
|
|
42 goto ERR; |
|
|
43 } |
|
|
44 |
|
|
45 if ((res = mp_copy(a, &a1)) != MP_OKAY) { |
|
|
46 goto ERR; |
|
|
47 } |
|
|
48 mp_rshd(&a1, B); |
|
|
49 mp_mod_2d(&a1, DIGIT_BIT * B, &a1); |
|
|
50 |
|
|
51 if ((res = mp_copy(a, &a2)) != MP_OKAY) { |
|
|
52 goto ERR; |
|
|
53 } |
|
|
54 mp_rshd(&a2, B*2); |
|
|
55 |
|
|
56 /* b = b2 * B**2 + b1 * B + b0 */ |
|
|
57 if ((res = mp_mod_2d(b, DIGIT_BIT * B, &b0)) != MP_OKAY) { |
|
|
58 goto ERR; |
|
|
59 } |
|
|
60 |
|
|
61 if ((res = mp_copy(b, &b1)) != MP_OKAY) { |
|
|
62 goto ERR; |
|
|
63 } |
|
|
64 mp_rshd(&b1, B); |
|
|
65 mp_mod_2d(&b1, DIGIT_BIT * B, &b1); |
|
|
66 |
|
|
67 if ((res = mp_copy(b, &b2)) != MP_OKAY) { |
|
|
68 goto ERR; |
|
|
69 } |
|
|
70 mp_rshd(&b2, B*2); |
|
|
71 |
|
|
72 /* w0 = a0*b0 */ |
|
|
73 if ((res = mp_mul(&a0, &b0, &w0)) != MP_OKAY) { |
|
|
74 goto ERR; |
|
|
75 } |
|
|
76 |
|
|
77 /* w4 = a2 * b2 */ |
|
|
78 if ((res = mp_mul(&a2, &b2, &w4)) != MP_OKAY) { |
|
|
79 goto ERR; |
|
|
80 } |
|
|
81 |
|
|
82 /* w1 = (a2 + 2(a1 + 2a0))(b2 + 2(b1 + 2b0)) */ |
|
|
83 if ((res = mp_mul_2(&a0, &tmp1)) != MP_OKAY) { |
|
|
84 goto ERR; |
|
|
85 } |
|
|
86 if ((res = mp_add(&tmp1, &a1, &tmp1)) != MP_OKAY) { |
|
|
87 goto ERR; |
|
|
88 } |
|
|
89 if ((res = mp_mul_2(&tmp1, &tmp1)) != MP_OKAY) { |
|
|
90 goto ERR; |
|
|
91 } |
|
|
92 if ((res = mp_add(&tmp1, &a2, &tmp1)) != MP_OKAY) { |
|
|
93 goto ERR; |
|
|
94 } |
|
|
95 |
|
|
96 if ((res = mp_mul_2(&b0, &tmp2)) != MP_OKAY) { |
|
|
97 goto ERR; |
|
|
98 } |
|
|
99 if ((res = mp_add(&tmp2, &b1, &tmp2)) != MP_OKAY) { |
|
|
100 goto ERR; |
|
|
101 } |
|
|
102 if ((res = mp_mul_2(&tmp2, &tmp2)) != MP_OKAY) { |
|
|
103 goto ERR; |
|
|
104 } |
|
|
105 if ((res = mp_add(&tmp2, &b2, &tmp2)) != MP_OKAY) { |
|
|
106 goto ERR; |
|
|
107 } |
|
|
108 |
|
|
109 if ((res = mp_mul(&tmp1, &tmp2, &w1)) != MP_OKAY) { |
|
|
110 goto ERR; |
|
|
111 } |
|
|
112 |
|
|
113 /* w3 = (a0 + 2(a1 + 2a2))(b0 + 2(b1 + 2b2)) */ |
|
|
114 if ((res = mp_mul_2(&a2, &tmp1)) != MP_OKAY) { |
|
|
115 goto ERR; |
|
|
116 } |
|
|
117 if ((res = mp_add(&tmp1, &a1, &tmp1)) != MP_OKAY) { |
|
|
118 goto ERR; |
|
|
119 } |
|
|
120 if ((res = mp_mul_2(&tmp1, &tmp1)) != MP_OKAY) { |
|
|
121 goto ERR; |
|
|
122 } |
|
|
123 if ((res = mp_add(&tmp1, &a0, &tmp1)) != MP_OKAY) { |
|
|
124 goto ERR; |
|
|
125 } |
|
|
126 |
|
|
127 if ((res = mp_mul_2(&b2, &tmp2)) != MP_OKAY) { |
|
|
128 goto ERR; |
|
|
129 } |
|
|
130 if ((res = mp_add(&tmp2, &b1, &tmp2)) != MP_OKAY) { |
|
|
131 goto ERR; |
|
|
132 } |
|
|
133 if ((res = mp_mul_2(&tmp2, &tmp2)) != MP_OKAY) { |
|
|
134 goto ERR; |
|
|
135 } |
|
|
136 if ((res = mp_add(&tmp2, &b0, &tmp2)) != MP_OKAY) { |
|
|
137 goto ERR; |
|
|
138 } |
|
|
139 |
|
|
140 if ((res = mp_mul(&tmp1, &tmp2, &w3)) != MP_OKAY) { |
|
|
141 goto ERR; |
|
|
142 } |
|
|
143 |
|
|
144 |
|
|
145 /* w2 = (a2 + a1 + a0)(b2 + b1 + b0) */ |
|
|
146 if ((res = mp_add(&a2, &a1, &tmp1)) != MP_OKAY) { |
|
|
147 goto ERR; |
|
|
148 } |
|
|
149 if ((res = mp_add(&tmp1, &a0, &tmp1)) != MP_OKAY) { |
|
|
150 goto ERR; |
|
|
151 } |
|
|
152 if ((res = mp_add(&b2, &b1, &tmp2)) != MP_OKAY) { |
|
|
153 goto ERR; |
|
|
154 } |
|
|
155 if ((res = mp_add(&tmp2, &b0, &tmp2)) != MP_OKAY) { |
|
|
156 goto ERR; |
|
|
157 } |
|
|
158 if ((res = mp_mul(&tmp1, &tmp2, &w2)) != MP_OKAY) { |
|
|
159 goto ERR; |
|
|
160 } |
|
|
161 |
|
|
162 /* now solve the matrix |
|
|
163 |
|
|
164 0 0 0 0 1 |
|
|
165 1 2 4 8 16 |
|
|
166 1 1 1 1 1 |
|
|
167 16 8 4 2 1 |
|
|
168 1 0 0 0 0 |
|
|
169 |
|
|
170 using 12 subtractions, 4 shifts, |
|
|
171 2 small divisions and 1 small multiplication |
|
|
172 */ |
|
|
173 |
|
|
174 /* r1 - r4 */ |
|
|
175 if ((res = mp_sub(&w1, &w4, &w1)) != MP_OKAY) { |
|
|
176 goto ERR; |
|
|
177 } |
|
|
178 /* r3 - r0 */ |
|
|
179 if ((res = mp_sub(&w3, &w0, &w3)) != MP_OKAY) { |
|
|
180 goto ERR; |
|
|
181 } |
|
|
182 /* r1/2 */ |
|
|
183 if ((res = mp_div_2(&w1, &w1)) != MP_OKAY) { |
|
|
184 goto ERR; |
|
|
185 } |
|
|
186 /* r3/2 */ |
|
|
187 if ((res = mp_div_2(&w3, &w3)) != MP_OKAY) { |
|
|
188 goto ERR; |
|
|
189 } |
|
|
190 /* r2 - r0 - r4 */ |
|
|
191 if ((res = mp_sub(&w2, &w0, &w2)) != MP_OKAY) { |
|
|
192 goto ERR; |
|
|
193 } |
|
|
194 if ((res = mp_sub(&w2, &w4, &w2)) != MP_OKAY) { |
|
|
195 goto ERR; |
|
|
196 } |
|
|
197 /* r1 - r2 */ |
|
|
198 if ((res = mp_sub(&w1, &w2, &w1)) != MP_OKAY) { |
|
|
199 goto ERR; |
|
|
200 } |
|
|
201 /* r3 - r2 */ |
|
|
202 if ((res = mp_sub(&w3, &w2, &w3)) != MP_OKAY) { |
|
|
203 goto ERR; |
|
|
204 } |
|
|
205 /* r1 - 8r0 */ |
|
|
206 if ((res = mp_mul_2d(&w0, 3, &tmp1)) != MP_OKAY) { |
|
|
207 goto ERR; |
|
|
208 } |
|
|
209 if ((res = mp_sub(&w1, &tmp1, &w1)) != MP_OKAY) { |
|
|
210 goto ERR; |
|
|
211 } |
|
|
212 /* r3 - 8r4 */ |
|
|
213 if ((res = mp_mul_2d(&w4, 3, &tmp1)) != MP_OKAY) { |
|
|
214 goto ERR; |
|
|
215 } |
|
|
216 if ((res = mp_sub(&w3, &tmp1, &w3)) != MP_OKAY) { |
|
|
217 goto ERR; |
|
|
218 } |
|
|
219 /* 3r2 - r1 - r3 */ |
|
|
220 if ((res = mp_mul_d(&w2, 3, &w2)) != MP_OKAY) { |
|
|
221 goto ERR; |
|
|
222 } |
|
|
223 if ((res = mp_sub(&w2, &w1, &w2)) != MP_OKAY) { |
|
|
224 goto ERR; |
|
|
225 } |
|
|
226 if ((res = mp_sub(&w2, &w3, &w2)) != MP_OKAY) { |
|
|
227 goto ERR; |
|
|
228 } |
|
|
229 /* r1 - r2 */ |
|
|
230 if ((res = mp_sub(&w1, &w2, &w1)) != MP_OKAY) { |
|
|
231 goto ERR; |
|
|
232 } |
|
|
233 /* r3 - r2 */ |
|
|
234 if ((res = mp_sub(&w3, &w2, &w3)) != MP_OKAY) { |
|
|
235 goto ERR; |
|
|
236 } |
|
|
237 /* r1/3 */ |
|
|
238 if ((res = mp_div_3(&w1, &w1, NULL)) != MP_OKAY) { |
|
|
239 goto ERR; |
|
|
240 } |
|
|
241 /* r3/3 */ |
|
|
242 if ((res = mp_div_3(&w3, &w3, NULL)) != MP_OKAY) { |
|
|
243 goto ERR; |
|
|
244 } |
|
|
245 |
|
|
246 /* at this point shift W[n] by B*n */ |
|
|
247 if ((res = mp_lshd(&w1, 1*B)) != MP_OKAY) { |
|
|
248 goto ERR; |
|
|
249 } |
|
|
250 if ((res = mp_lshd(&w2, 2*B)) != MP_OKAY) { |
|
|
251 goto ERR; |
|
|
252 } |
|
|
253 if ((res = mp_lshd(&w3, 3*B)) != MP_OKAY) { |
|
|
254 goto ERR; |
|
|
255 } |
|
|
256 if ((res = mp_lshd(&w4, 4*B)) != MP_OKAY) { |
|
|
257 goto ERR; |
|
|
258 } |
|
|
259 |
|
|
260 if ((res = mp_add(&w0, &w1, c)) != MP_OKAY) { |
|
|
261 goto ERR; |
|
|
262 } |
|
|
263 if ((res = mp_add(&w2, &w3, &tmp1)) != MP_OKAY) { |
|
|
264 goto ERR; |
|
|
265 } |
|
|
266 if ((res = mp_add(&w4, &tmp1, &tmp1)) != MP_OKAY) { |
|
|
267 goto ERR; |
|
|
268 } |
|
|
269 if ((res = mp_add(&tmp1, c, c)) != MP_OKAY) { |
|
|
270 goto ERR; |
|
|
271 } |
|
|
272 |
|
|
273 ERR: |
|
|
274 mp_clear_multi(&w0, &w1, &w2, &w3, &w4, |
|
|
275 &a0, &a1, &a2, &b0, &b1, |
|
|
276 &b2, &tmp1, &tmp2, NULL); |
|
|
277 return res; |
|
|
278 } |
|
|
279 |
|
|
280 #endif |
|
386
|
281 |
|
|
282 /* $Source: /cvs/libtom/libtommath/bn_mp_toom_mul.c,v $ */ |
|
|
283 /* $Revision: 1.3 $ */ |
|
|
284 /* $Date: 2006/03/31 14:18:44 $ */ |
