Planet
navi homePPSaboutscreenshotsdownloaddevelopmentforum

source: downloads/tcl8.5.2/libtommath/bn_mp_toom_mul.c @ 35

Last change on this file since 35 was 25, checked in by landauf, 18 years ago

added tcl to libs
File size: 7.0 KB
RevLine 
[25]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://math.libtomcrypt.com
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*/
25int 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     
273ERR:
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
281
282/* $Source: /cvsroot/tcl/libtommath/bn_mp_toom_mul.c,v $ */
283/* $Revision: 1.1.1.4 $ */
284/* $Date: 2006/12/01 00:08:11 $ */
Note: See TracBrowser for help on using the repository browser.