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