Mercurial > dropbear

comparison libtommath/bn_mp_toom_mul.c @ 293:9d110777f345 contrib-blacklist

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
propagate from branch 'au.asn.ucc.matt.dropbear' (head 7ad1775ed65e75dbece27fe6b65bf1a234db386a) to branch 'au.asn.ucc.matt.dropbear.contrib.blacklist' (head 1d86a4f0a401cc68c2670d821a2f6366c37af143)
author Matt Johnston <matt@ucc.asn.au>
date Fri, 10 Mar 2006 06:31:29 +0000
parents eed26cff980b
children 5ff8218bcee9
comparison
equal deleted inserted replaced
247:c07de41b53d7 293:9d110777f345
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