Mercurial > dropbear
comparison libtommath/bn_mp_toom_mul.c @ 284:eed26cff980b
propagate from branch 'au.asn.ucc.matt.ltm.dropbear' (head 6c790cad5a7fa866ad062cb3a0c279f7ba788583)
to branch 'au.asn.ucc.matt.dropbear' (head fff0894a0399405a9410ea1c6d118f342cf2aa64)
author | Matt Johnston <matt@ucc.asn.au> |
---|---|
date | Wed, 08 Mar 2006 13:23:49 +0000 |
parents | |
children | 5ff8218bcee9 |
comparison
equal
deleted
inserted
replaced
283:bd240aa12ba7 | 284:eed26cff980b |
---|---|
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 * | |
15 * Tom St Denis, [email protected], http://math.libtomcrypt.org | |
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 |