Botan 1.10.17
mp_karat.cpp
Go to the documentation of this file.
1/*
2* Karatsuba Multiplication/Squaring
3* (C) 1999-2010 Jack Lloyd
4*
5* Distributed under the terms of the Botan license
6*/
7
8#include <botan/internal/mp_core.h>
9#include <botan/internal/mp_asmi.h>
10#include <botan/internal/assert.h>
11#include <botan/mem_ops.h>
12
13namespace Botan {
14
15namespace {
16
17/*
18* Karatsuba Multiplication Operation
19*/
20void karatsuba_mul(word z[], const word x[], const word y[], size_t N,
21 word workspace[])
22 {
23 if(N < BOTAN_KARAT_MUL_THRESHOLD || N % 2)
24 {
25 if(N == 6)
26 return bigint_comba_mul6(z, x, y);
27 else if(N == 8)
28 return bigint_comba_mul8(z, x, y);
29 else if(N == 16)
30 return bigint_comba_mul16(z, x, y);
31 else
32 return bigint_simple_mul(z, x, N, y, N);
33 }
34
35 const size_t N2 = N / 2;
36
37 const word* x0 = x;
38 const word* x1 = x + N2;
39 const word* y0 = y;
40 const word* y1 = y + N2;
41 word* z0 = z;
42 word* z1 = z + N;
43
44 const s32bit cmp0 = bigint_cmp(x0, N2, x1, N2);
45 const s32bit cmp1 = bigint_cmp(y1, N2, y0, N2);
46
47 clear_mem(workspace, 2*N);
48
49 //if(cmp0 && cmp1)
50 {
51 if(cmp0 > 0)
52 bigint_sub3(z0, x0, N2, x1, N2);
53 else
54 bigint_sub3(z0, x1, N2, x0, N2);
55
56 if(cmp1 > 0)
57 bigint_sub3(z1, y1, N2, y0, N2);
58 else
59 bigint_sub3(z1, y0, N2, y1, N2);
60
61 karatsuba_mul(workspace, z0, z1, N2, workspace+N);
62 }
63
64 karatsuba_mul(z0, x0, y0, N2, workspace+N);
65 karatsuba_mul(z1, x1, y1, N2, workspace+N);
66
67 const size_t blocks_of_8 = N - (N % 8);
68
69 word ws_carry = 0;
70
71 for(size_t j = 0; j != blocks_of_8; j += 8)
72 ws_carry = word8_add3(workspace + N + j, z0 + j, z1 + j, ws_carry);
73
74 for(size_t j = blocks_of_8; j != N; ++j)
75 workspace[N + j] = word_add(z0[j], z1[j], &ws_carry);
76
77 word z_carry = 0;
78
79 for(size_t j = 0; j != blocks_of_8; j += 8)
80 z_carry = word8_add2(z + N2 + j, workspace + N + j, z_carry);
81
82 for(size_t j = blocks_of_8; j != N; ++j)
83 z[N2 + j] = word_add(z[N2 + j], workspace[N + j], &z_carry);
84
85 z[N + N2] = word_add(z[N + N2], ws_carry, &z_carry);
86
87 if(z_carry)
88 for(size_t j = 1; j != N2; ++j)
89 if(++z[N + N2 + j])
90 break;
91
92 if((cmp0 == cmp1) || (cmp0 == 0) || (cmp1 == 0))
93 bigint_add2(z + N2, 2*N-N2, workspace, N);
94 else
95 bigint_sub2(z + N2, 2*N-N2, workspace, N);
96 }
97
98/*
99* Karatsuba Squaring Operation
100*/
101void karatsuba_sqr(word z[], const word x[], size_t N, word workspace[])
102 {
103 if(N < BOTAN_KARAT_SQR_THRESHOLD || N % 2)
104 {
105 if(N == 6)
106 return bigint_comba_sqr6(z, x);
107 else if(N == 8)
108 return bigint_comba_sqr8(z, x);
109 else if(N == 16)
110 return bigint_comba_sqr16(z, x);
111 else
112 return bigint_simple_sqr(z, x, N);
113 }
114
115 const size_t N2 = N / 2;
116
117 const word* x0 = x;
118 const word* x1 = x + N2;
119 word* z0 = z;
120 word* z1 = z + N;
121
122 const s32bit cmp = bigint_cmp(x0, N2, x1, N2);
123
124 clear_mem(workspace, 2*N);
125
126 //if(cmp)
127 {
128 if(cmp > 0)
129 bigint_sub3(z0, x0, N2, x1, N2);
130 else
131 bigint_sub3(z0, x1, N2, x0, N2);
132
133 karatsuba_sqr(workspace, z0, N2, workspace+N);
134 }
135
136 karatsuba_sqr(z0, x0, N2, workspace+N);
137 karatsuba_sqr(z1, x1, N2, workspace+N);
138
139 const size_t blocks_of_8 = N - (N % 8);
140
141 word ws_carry = 0;
142
143 for(size_t j = 0; j != blocks_of_8; j += 8)
144 ws_carry = word8_add3(workspace + N + j, z0 + j, z1 + j, ws_carry);
145
146 for(size_t j = blocks_of_8; j != N; ++j)
147 workspace[N + j] = word_add(z0[j], z1[j], &ws_carry);
148
149 word z_carry = 0;
150
151 for(size_t j = 0; j != blocks_of_8; j += 8)
152 z_carry = word8_add2(z + N2 + j, workspace + N + j, z_carry);
153
154 for(size_t j = blocks_of_8; j != N; ++j)
155 z[N2 + j] = word_add(z[N2 + j], workspace[N + j], &z_carry);
156
157 z[N + N2] = word_add(z[N + N2], ws_carry, &z_carry);
158
159 if(z_carry)
160 for(size_t j = 1; j != N2; ++j)
161 if(++z[N + N2 + j])
162 break;
163
164 /*
165 * This is only actually required if cmp is != 0, however
166 * if cmp==0 then workspace[0:N] == 0 and avoiding the jump
167 * hides a timing channel.
168 */
169 bigint_sub2(z + N2, 2*N-N2, workspace, N);
170 }
171
172/*
173* Pick a good size for the Karatsuba multiply
174*/
175size_t karatsuba_size(size_t z_size,
176 size_t x_size, size_t x_sw,
177 size_t y_size, size_t y_sw)
178 {
179 if(x_sw > x_size || x_sw > y_size || y_sw > x_size || y_sw > y_size)
180 return 0;
181
182 if(((x_size == x_sw) && (x_size % 2)) ||
183 ((y_size == y_sw) && (y_size % 2)))
184 return 0;
185
186 const size_t start = (x_sw > y_sw) ? x_sw : y_sw;
187 const size_t end = (x_size < y_size) ? x_size : y_size;
188
189 if(start == end)
190 {
191 if(start % 2)
192 return 0;
193 return start;
194 }
195
196 for(size_t j = start; j <= end; ++j)
197 {
198 if(j % 2)
199 continue;
200
201 if(2*j > z_size)
202 return 0;
203
204 if(x_sw <= j && j <= x_size && y_sw <= j && j <= y_size)
205 {
206 if(j % 4 == 2 &&
207 (j+2) <= x_size && (j+2) <= y_size && 2*(j+2) <= z_size)
208 return j+2;
209 return j;
210 }
211 }
212
213 return 0;
214 }
215
216/*
217* Pick a good size for the Karatsuba squaring
218*/
219size_t karatsuba_size(size_t z_size, size_t x_size, size_t x_sw)
220 {
221 if(x_sw == x_size)
222 {
223 if(x_sw % 2)
224 return 0;
225 return x_sw;
226 }
227
228 for(size_t j = x_sw; j <= x_size; ++j)
229 {
230 if(j % 2)
231 continue;
232
233 if(2*j > z_size)
234 return 0;
235
236 if(j % 4 == 2 && (j+2) <= x_size && 2*(j+2) <= z_size)
237 return j+2;
238 return j;
239 }
240
241 return 0;
242 }
243
244}
245
246/*
247* Multiplication Algorithm Dispatcher
248*/
249void bigint_mul(word z[], size_t z_size, word workspace[],
250 const word x[], size_t x_size, size_t x_sw,
251 const word y[], size_t y_size, size_t y_sw)
252 {
253 BOTAN_ASSERT(z_size > x_sw && z_size > y_sw && z_size - x_sw >= y_sw, "Sufficient output size");
254
255 if(x_sw == 1)
256 {
257 bigint_linmul3(z, y, y_sw, x[0]);
258 }
259 else if(y_sw == 1)
260 {
261 bigint_linmul3(z, x, x_sw, y[0]);
262 }
263 else if(x_sw <= 4 && x_size >= 4 &&
264 y_sw <= 4 && y_size >= 4 && z_size >= 8)
265 {
266 bigint_comba_mul4(z, x, y);
267 }
268 else if(x_sw <= 6 && x_size >= 6 &&
269 y_sw <= 6 && y_size >= 6 && z_size >= 12)
270 {
271 bigint_comba_mul6(z, x, y);
272 }
273 else if(x_sw <= 8 && x_size >= 8 &&
274 y_sw <= 8 && y_size >= 8 && z_size >= 16)
275 {
276 bigint_comba_mul8(z, x, y);
277 }
278 else if(x_sw <= 16 && x_size >= 16 &&
279 y_sw <= 16 && y_size >= 16 && z_size >= 32)
280 {
281 bigint_comba_mul16(z, x, y);
282 }
283 else if(x_sw < BOTAN_KARAT_MUL_THRESHOLD ||
284 y_sw < BOTAN_KARAT_MUL_THRESHOLD ||
285 !workspace)
286 {
287 bigint_simple_mul(z, x, x_sw, y, y_sw);
288 }
289 else
290 {
291 const size_t N = karatsuba_size(z_size, x_size, x_sw, y_size, y_sw);
292
293 if(N)
294 {
295 clear_mem(workspace, 2*N);
296 karatsuba_mul(z, x, y, N, workspace);
297 }
298 else
299 bigint_simple_mul(z, x, x_sw, y, y_sw);
300 }
301 }
302
303/*
304* Squaring Algorithm Dispatcher
305*/
306void bigint_sqr(word z[], size_t z_size, word workspace[],
307 const word x[], size_t x_size, size_t x_sw)
308 {
309 BOTAN_ASSERT(z_size/2 >= x_sw, "Sufficient output size");
310
311 if(x_sw == 1)
312 {
313 bigint_linmul3(z, x, x_sw, x[0]);
314 }
315 else if(x_sw <= 4 && x_size >= 4 && z_size >= 8)
316 {
317 bigint_comba_sqr4(z, x);
318 }
319 else if(x_sw <= 6 && x_size >= 6 && z_size >= 12)
320 {
321 bigint_comba_sqr6(z, x);
322 }
323 else if(x_sw <= 8 && x_size >= 8 && z_size >= 16)
324 {
325 bigint_comba_sqr8(z, x);
326 }
327 else if(x_sw <= 16 && x_size >= 16 && z_size >= 32)
328 {
329 bigint_comba_sqr16(z, x);
330 }
331 else if(x_size < BOTAN_KARAT_SQR_THRESHOLD || !workspace)
332 {
333 bigint_simple_sqr(z, x, x_sw);
334 }
335 else
336 {
337 const size_t N = karatsuba_size(z_size, x_size, x_sw);
338
339 if(N)
340 {
341 clear_mem(workspace, 2*N);
342 karatsuba_sqr(z, x, N, workspace);
343 }
344 else
345 bigint_simple_sqr(z, x, x_sw);
346 }
347 }
348
349}
#define BOTAN_ASSERT(expr, msg)
Definition assert.h:19
void bigint_linmul3(word z[], const word x[], size_t x_size, word y)
Definition mp_asm.cpp:238
word bigint_sub2(word x[], size_t x_size, const word y[], size_t y_size)
Definition mp_asm.cpp:158
void bigint_comba_mul4(word z[8], const word x[4], const word y[4])
Definition mp_comba.cpp:51
void bigint_sqr(word z[], size_t z_size, word workspace[], const word x[], size_t x_size, size_t x_sw)
Definition mp_karat.cpp:306
signed int s32bit
Definition types.h:37
void bigint_simple_mul(word z[], const word x[], size_t x_size, const word y[], size_t y_size)
Definition mp_mulop.cpp:20
void bigint_comba_sqr6(word z[12], const word x[6])
Definition mp_comba.cpp:90
void bigint_mul(word z[], size_t z_size, word workspace[], const word x[], size_t x_size, size_t x_sw, const word y[], size_t y_size, size_t y_sw)
Definition mp_karat.cpp:249
void bigint_comba_sqr16(word z[32], const word x[16])
Definition mp_comba.cpp:387
s32bit bigint_cmp(const word x[], size_t x_size, const word y[], size_t y_size)
Definition mp_misc.cpp:41
void bigint_simple_sqr(word z[], const word x[], size_t x_size)
Definition mp_mulop.cpp:54
word word8_add2(word x[8], const word y[8], word carry)
Definition mp_asmi.h:33
void bigint_comba_mul8(word z[16], const word x[8], const word y[8])
Definition mp_comba.cpp:284
void bigint_comba_sqr4(word z[8], const word x[4])
Definition mp_comba.cpp:18
word word8_add3(word z[8], const word x[8], const word y[8], word carry)
Definition mp_asmi.h:49
void bigint_add2(word x[], size_t x_size, const word y[], size_t y_size)
Definition mp_asm.cpp:139
word word_add(word x, word y, word *carry)
Definition mp_asmi.h:21
void bigint_comba_mul16(word z[32], const word x[16], const word y[16])
Definition mp_comba.cpp:594
void bigint_comba_mul6(word z[12], const word x[6], const word y[6])
Definition mp_comba.cpp:142
void clear_mem(T *ptr, size_t n)
Definition mem_ops.h:32
word bigint_sub3(word z[], const word x[], size_t x_size, const word y[], size_t y_size)
Definition mp_asm.cpp:198
void bigint_comba_sqr8(word z[16], const word x[8])
Definition mp_comba.cpp:209