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