1 /* $NetBSD: bn_mp_toom_sqr.c,v 1.1.1.2 2014/04/24 12:45:31 pettai Exp $ */
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 */
21 int
mp_toom_sqr(mp_int * a,mp_int * b)22 mp_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
219 ERR:
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