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