Line data Source code
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, tomstdenis@gmail.com, http://libtom.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 0 : 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 0 : if ((res = mp_init_multi(&w0, &w1, &w2, &w3, &w4,
32 : &a0, &a1, &a2, &b0, &b1,
33 : &b2, &tmp1, &tmp2, NULL)) != MP_OKAY) {
34 0 : return res;
35 : }
36 :
37 : /* B */
38 0 : B = MIN(a->used, b->used) / 3;
39 :
40 : /* a = a2 * B**2 + a1 * B + a0 */
41 0 : if ((res = mp_mod_2d(a, DIGIT_BIT * B, &a0)) != MP_OKAY) {
42 0 : goto ERR;
43 : }
44 :
45 0 : if ((res = mp_copy(a, &a1)) != MP_OKAY) {
46 0 : goto ERR;
47 : }
48 0 : mp_rshd(&a1, B);
49 0 : mp_mod_2d(&a1, DIGIT_BIT * B, &a1);
50 :
51 0 : if ((res = mp_copy(a, &a2)) != MP_OKAY) {
52 0 : goto ERR;
53 : }
54 0 : mp_rshd(&a2, B*2);
55 :
56 : /* b = b2 * B**2 + b1 * B + b0 */
57 0 : if ((res = mp_mod_2d(b, DIGIT_BIT * B, &b0)) != MP_OKAY) {
58 0 : goto ERR;
59 : }
60 :
61 0 : if ((res = mp_copy(b, &b1)) != MP_OKAY) {
62 0 : goto ERR;
63 : }
64 0 : mp_rshd(&b1, B);
65 0 : mp_mod_2d(&b1, DIGIT_BIT * B, &b1);
66 :
67 0 : if ((res = mp_copy(b, &b2)) != MP_OKAY) {
68 0 : goto ERR;
69 : }
70 0 : mp_rshd(&b2, B*2);
71 :
72 : /* w0 = a0*b0 */
73 0 : if ((res = mp_mul(&a0, &b0, &w0)) != MP_OKAY) {
74 0 : goto ERR;
75 : }
76 :
77 : /* w4 = a2 * b2 */
78 0 : if ((res = mp_mul(&a2, &b2, &w4)) != MP_OKAY) {
79 0 : goto ERR;
80 : }
81 :
82 : /* w1 = (a2 + 2(a1 + 2a0))(b2 + 2(b1 + 2b0)) */
83 0 : if ((res = mp_mul_2(&a0, &tmp1)) != MP_OKAY) {
84 0 : goto ERR;
85 : }
86 0 : if ((res = mp_add(&tmp1, &a1, &tmp1)) != MP_OKAY) {
87 0 : goto ERR;
88 : }
89 0 : if ((res = mp_mul_2(&tmp1, &tmp1)) != MP_OKAY) {
90 0 : goto ERR;
91 : }
92 0 : if ((res = mp_add(&tmp1, &a2, &tmp1)) != MP_OKAY) {
93 0 : goto ERR;
94 : }
95 :
96 0 : if ((res = mp_mul_2(&b0, &tmp2)) != MP_OKAY) {
97 0 : goto ERR;
98 : }
99 0 : if ((res = mp_add(&tmp2, &b1, &tmp2)) != MP_OKAY) {
100 0 : goto ERR;
101 : }
102 0 : if ((res = mp_mul_2(&tmp2, &tmp2)) != MP_OKAY) {
103 0 : goto ERR;
104 : }
105 0 : if ((res = mp_add(&tmp2, &b2, &tmp2)) != MP_OKAY) {
106 0 : goto ERR;
107 : }
108 :
109 0 : if ((res = mp_mul(&tmp1, &tmp2, &w1)) != MP_OKAY) {
110 0 : goto ERR;
111 : }
112 :
113 : /* w3 = (a0 + 2(a1 + 2a2))(b0 + 2(b1 + 2b2)) */
114 0 : if ((res = mp_mul_2(&a2, &tmp1)) != MP_OKAY) {
115 0 : goto ERR;
116 : }
117 0 : if ((res = mp_add(&tmp1, &a1, &tmp1)) != MP_OKAY) {
118 0 : goto ERR;
119 : }
120 0 : if ((res = mp_mul_2(&tmp1, &tmp1)) != MP_OKAY) {
121 0 : goto ERR;
122 : }
123 0 : if ((res = mp_add(&tmp1, &a0, &tmp1)) != MP_OKAY) {
124 0 : goto ERR;
125 : }
126 :
127 0 : if ((res = mp_mul_2(&b2, &tmp2)) != MP_OKAY) {
128 0 : goto ERR;
129 : }
130 0 : if ((res = mp_add(&tmp2, &b1, &tmp2)) != MP_OKAY) {
131 0 : goto ERR;
132 : }
133 0 : if ((res = mp_mul_2(&tmp2, &tmp2)) != MP_OKAY) {
134 0 : goto ERR;
135 : }
136 0 : if ((res = mp_add(&tmp2, &b0, &tmp2)) != MP_OKAY) {
137 0 : goto ERR;
138 : }
139 :
140 0 : if ((res = mp_mul(&tmp1, &tmp2, &w3)) != MP_OKAY) {
141 0 : goto ERR;
142 : }
143 :
144 :
145 : /* w2 = (a2 + a1 + a0)(b2 + b1 + b0) */
146 0 : if ((res = mp_add(&a2, &a1, &tmp1)) != MP_OKAY) {
147 0 : goto ERR;
148 : }
149 0 : if ((res = mp_add(&tmp1, &a0, &tmp1)) != MP_OKAY) {
150 0 : goto ERR;
151 : }
152 0 : if ((res = mp_add(&b2, &b1, &tmp2)) != MP_OKAY) {
153 0 : goto ERR;
154 : }
155 0 : if ((res = mp_add(&tmp2, &b0, &tmp2)) != MP_OKAY) {
156 0 : goto ERR;
157 : }
158 0 : if ((res = mp_mul(&tmp1, &tmp2, &w2)) != MP_OKAY) {
159 0 : 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 0 : if ((res = mp_sub(&w1, &w4, &w1)) != MP_OKAY) {
176 0 : goto ERR;
177 : }
178 : /* r3 - r0 */
179 0 : if ((res = mp_sub(&w3, &w0, &w3)) != MP_OKAY) {
180 0 : goto ERR;
181 : }
182 : /* r1/2 */
183 0 : if ((res = mp_div_2(&w1, &w1)) != MP_OKAY) {
184 0 : goto ERR;
185 : }
186 : /* r3/2 */
187 0 : if ((res = mp_div_2(&w3, &w3)) != MP_OKAY) {
188 0 : goto ERR;
189 : }
190 : /* r2 - r0 - r4 */
191 0 : if ((res = mp_sub(&w2, &w0, &w2)) != MP_OKAY) {
192 0 : goto ERR;
193 : }
194 0 : if ((res = mp_sub(&w2, &w4, &w2)) != MP_OKAY) {
195 0 : goto ERR;
196 : }
197 : /* r1 - r2 */
198 0 : if ((res = mp_sub(&w1, &w2, &w1)) != MP_OKAY) {
199 0 : goto ERR;
200 : }
201 : /* r3 - r2 */
202 0 : if ((res = mp_sub(&w3, &w2, &w3)) != MP_OKAY) {
203 0 : goto ERR;
204 : }
205 : /* r1 - 8r0 */
206 0 : if ((res = mp_mul_2d(&w0, 3, &tmp1)) != MP_OKAY) {
207 0 : goto ERR;
208 : }
209 0 : if ((res = mp_sub(&w1, &tmp1, &w1)) != MP_OKAY) {
210 0 : goto ERR;
211 : }
212 : /* r3 - 8r4 */
213 0 : if ((res = mp_mul_2d(&w4, 3, &tmp1)) != MP_OKAY) {
214 0 : goto ERR;
215 : }
216 0 : if ((res = mp_sub(&w3, &tmp1, &w3)) != MP_OKAY) {
217 0 : goto ERR;
218 : }
219 : /* 3r2 - r1 - r3 */
220 0 : if ((res = mp_mul_d(&w2, 3, &w2)) != MP_OKAY) {
221 0 : goto ERR;
222 : }
223 0 : if ((res = mp_sub(&w2, &w1, &w2)) != MP_OKAY) {
224 0 : goto ERR;
225 : }
226 0 : if ((res = mp_sub(&w2, &w3, &w2)) != MP_OKAY) {
227 0 : goto ERR;
228 : }
229 : /* r1 - r2 */
230 0 : if ((res = mp_sub(&w1, &w2, &w1)) != MP_OKAY) {
231 0 : goto ERR;
232 : }
233 : /* r3 - r2 */
234 0 : if ((res = mp_sub(&w3, &w2, &w3)) != MP_OKAY) {
235 0 : goto ERR;
236 : }
237 : /* r1/3 */
238 0 : if ((res = mp_div_3(&w1, &w1, NULL)) != MP_OKAY) {
239 0 : goto ERR;
240 : }
241 : /* r3/3 */
242 0 : if ((res = mp_div_3(&w3, &w3, NULL)) != MP_OKAY) {
243 0 : goto ERR;
244 : }
245 :
246 : /* at this point shift W[n] by B*n */
247 0 : if ((res = mp_lshd(&w1, 1*B)) != MP_OKAY) {
248 0 : goto ERR;
249 : }
250 0 : if ((res = mp_lshd(&w2, 2*B)) != MP_OKAY) {
251 0 : goto ERR;
252 : }
253 0 : if ((res = mp_lshd(&w3, 3*B)) != MP_OKAY) {
254 0 : goto ERR;
255 : }
256 0 : if ((res = mp_lshd(&w4, 4*B)) != MP_OKAY) {
257 0 : goto ERR;
258 : }
259 :
260 0 : if ((res = mp_add(&w0, &w1, c)) != MP_OKAY) {
261 0 : goto ERR;
262 : }
263 0 : if ((res = mp_add(&w2, &w3, &tmp1)) != MP_OKAY) {
264 0 : goto ERR;
265 : }
266 0 : if ((res = mp_add(&w4, &tmp1, &tmp1)) != MP_OKAY) {
267 0 : goto ERR;
268 : }
269 0 : if ((res = mp_add(&tmp1, c, c)) != MP_OKAY) {
270 0 : goto ERR;
271 : }
272 :
273 0 : ERR:
274 0 : mp_clear_multi(&w0, &w1, &w2, &w3, &w4,
275 : &a0, &a1, &a2, &b0, &b1,
276 : &b2, &tmp1, &tmp2, NULL);
277 0 : return res;
278 : }
279 :
280 : #endif
281 :
282 : /* $Source: /cvs/libtom/libtommath/bn_mp_toom_mul.c,v $ */
283 : /* $Revision: 1.4 $ */
284 : /* $Date: 2006/12/28 01:25:13 $ */
|
