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