Coverage Report

Created: 2020-06-30 13:58

/src/botan/src/lib/math/mp/mp_karat.cpp
Line
Count
Source (jump to first uncovered line)
1
/*
2
* Multiplication and Squaring
3
* (C) 1999-2010,2018 Jack Lloyd
4
*     2016 Matthias Gierlings
5
*
6
* Botan is released under the Simplified BSD License (see license.txt)
7
*/
8
9
#include <botan/internal/mp_core.h>
10
#include <botan/internal/mp_asmi.h>
11
#include <botan/internal/ct_utils.h>
12
#include <botan/mem_ops.h>
13
#include <botan/exceptn.h>
14
15
namespace Botan {
16
17
namespace {
18
19
const size_t KARATSUBA_MULTIPLY_THRESHOLD = 32;
20
const size_t KARATSUBA_SQUARE_THRESHOLD = 32;
21
22
/*
23
* Simple O(N^2) Multiplication
24
*/
25
void basecase_mul(word z[], size_t z_size,
26
                  const word x[], size_t x_size,
27
                  const word y[], size_t y_size)
28
7.08M
   {
29
7.08M
   if(z_size < x_size + y_size)
30
0
      throw Invalid_Argument("basecase_mul z_size too small");
31
7.08M
32
7.08M
   const size_t x_size_8 = x_size - (x_size % 8);
33
7.08M
34
7.08M
   clear_mem(z, z_size);
35
7.08M
36
92.8M
   for(size_t i = 0; i != y_size; ++i)
37
85.7M
      {
38
85.7M
      const word y_i = y[i];
39
85.7M
40
85.7M
      word carry = 0;
41
85.7M
42
280M
      for(size_t j = 0; j != x_size_8; j += 8)
43
194M
         carry = word8_madd3(z + i + j, x + j, y_i, carry);
44
85.7M
45
428M
      for(size_t j = x_size_8; j != x_size; ++j)
46
342M
         z[i+j] = word_madd3(x[j], y_i, z[i+j], &carry);
47
85.7M
48
85.7M
      z[x_size+i] = carry;
49
85.7M
      }
50
7.08M
   }
51
52
void basecase_sqr(word z[], size_t z_size,
53
                  const word x[], size_t x_size)
54
1.52M
   {
55
1.52M
   if(z_size < 2*x_size)
56
0
      throw Invalid_Argument("basecase_sqr z_size too small");
57
1.52M
58
1.52M
   const size_t x_size_8 = x_size - (x_size % 8);
59
1.52M
60
1.52M
   clear_mem(z, z_size);
61
1.52M
62
29.0M
   for(size_t i = 0; i != x_size; ++i)
63
27.5M
      {
64
27.5M
      const word x_i = x[i];
65
27.5M
66
27.5M
      word carry = 0;
67
27.5M
68
99.1M
      for(size_t j = 0; j != x_size_8; j += 8)
69
71.6M
         carry = word8_madd3(z + i + j, x + j, x_i, carry);
70
27.5M
71
163M
      for(size_t j = x_size_8; j != x_size; ++j)
72
136M
         z[i+j] = word_madd3(x[j], x_i, z[i+j], &carry);
73
27.5M
74
27.5M
      z[x_size+i] = carry;
75
27.5M
      }
76
1.52M
   }
77
78
/*
79
* Karatsuba Multiplication Operation
80
*/
81
void karatsuba_mul(word z[], const word x[], const word y[], size_t N,
82
                   word workspace[])
83
7.71M
   {
84
7.71M
   if(N < KARATSUBA_MULTIPLY_THRESHOLD || N % 2)
85
5.61M
      {
86
5.61M
      switch(N)
87
5.61M
         {
88
0
         case 6:
89
0
            return bigint_comba_mul6(z, x, y);
90
0
         case 8:
91
0
            return bigint_comba_mul8(z, x, y);
92
0
         case 9:
93
0
            return bigint_comba_mul9(z, x, y);
94
5.27M
         case 16:
95
5.27M
            return bigint_comba_mul16(z, x, y);
96
321k
         case 24:
97
321k
            return bigint_comba_mul24(z, x, y);
98
11.1k
         default:
99
11.1k
            return basecase_mul(z, 2*N, x, N, y, N);
100
2.10M
         }
101
2.10M
      }
102
2.10M
103
2.10M
   const size_t N2 = N / 2;
104
2.10M
105
2.10M
   const word* x0 = x;
106
2.10M
   const word* x1 = x + N2;
107
2.10M
   const word* y0 = y;
108
2.10M
   const word* y1 = y + N2;
109
2.10M
   word* z0 = z;
110
2.10M
   word* z1 = z + N;
111
2.10M
112
2.10M
   word* ws0 = workspace;
113
2.10M
   word* ws1 = workspace + N;
114
2.10M
115
2.10M
   clear_mem(workspace, 2*N);
116
2.10M
117
2.10M
   /*
118
2.10M
   * If either of cmp0 or cmp1 is zero then z0 or z1 resp is zero here,
119
2.10M
   * resulting in a no-op - z0*z1 will be equal to zero so we don't need to do
120
2.10M
   * anything, clear_mem above already set the correct result.
121
2.10M
   *
122
2.10M
   * However we ignore the result of the comparisons and always perform the
123
2.10M
   * subtractions and recursively multiply to avoid the timing channel.
124
2.10M
   */
125
2.10M
126
2.10M
   // First compute (X_lo - X_hi)*(Y_hi - Y_lo)
127
2.10M
   const auto cmp0 = bigint_sub_abs(z0, x0, x1, N2, workspace);
128
2.10M
   const auto cmp1 = bigint_sub_abs(z1, y1, y0, N2, workspace);
129
2.10M
   const auto neg_mask = ~(cmp0 ^ cmp1);
130
2.10M
131
2.10M
   karatsuba_mul(ws0, z0, z1, N2, ws1);
132
2.10M
133
2.10M
   // Compute X_lo * Y_lo
134
2.10M
   karatsuba_mul(z0, x0, y0, N2, ws1);
135
2.10M
136
2.10M
   // Compute X_hi * Y_hi
137
2.10M
   karatsuba_mul(z1, x1, y1, N2, ws1);
138
2.10M
139
2.10M
   const word ws_carry = bigint_add3_nc(ws1, z0, N, z1, N);
140
2.10M
   word z_carry = bigint_add2_nc(z + N2, N, ws1, N);
141
2.10M
142
2.10M
   z_carry += bigint_add2_nc(z + N + N2, N2, &ws_carry, 1);
143
2.10M
   bigint_add2_nc(z + N + N2, N2, &z_carry, 1);
144
2.10M
145
2.10M
   clear_mem(workspace + N, N2);
146
2.10M
147
2.10M
   bigint_cnd_add_or_sub(neg_mask, z + N2, workspace, 2*N-N2);
148
2.10M
   }
149
150
/*
151
* Karatsuba Squaring Operation
152
*/
153
void karatsuba_sqr(word z[], const word x[], size_t N, word workspace[])
154
21.6M
   {
155
21.6M
   if(N < KARATSUBA_SQUARE_THRESHOLD || N % 2)
156
15.8M
      {
157
15.8M
      switch(N)
158
15.8M
         {
159
0
         case 6:
160
0
            return bigint_comba_sqr6(z, x);
161
0
         case 8:
162
0
            return bigint_comba_sqr8(z, x);
163
0
         case 9:
164
0
            return bigint_comba_sqr9(z, x);
165
14.1M
         case 16:
166
14.1M
            return bigint_comba_sqr16(z, x);
167
1.06M
         case 24:
168
1.06M
            return bigint_comba_sqr24(z, x);
169
647k
         default:
170
647k
            return basecase_sqr(z, 2*N, x, N);
171
5.84M
         }
172
5.84M
      }
173
5.84M
174
5.84M
   const size_t N2 = N / 2;
175
5.84M
176
5.84M
   const word* x0 = x;
177
5.84M
   const word* x1 = x + N2;
178
5.84M
   word* z0 = z;
179
5.84M
   word* z1 = z + N;
180
5.84M
181
5.84M
   word* ws0 = workspace;
182
5.84M
   word* ws1 = workspace + N;
183
5.84M
184
5.84M
   clear_mem(workspace, 2*N);
185
5.84M
186
5.84M
   // See comment in karatsuba_mul
187
5.84M
   bigint_sub_abs(z0, x0, x1, N2, workspace);
188
5.84M
   karatsuba_sqr(ws0, z0, N2, ws1);
189
5.84M
190
5.84M
   karatsuba_sqr(z0, x0, N2, ws1);
191
5.84M
   karatsuba_sqr(z1, x1, N2, ws1);
192
5.84M
193
5.84M
   const word ws_carry = bigint_add3_nc(ws1, z0, N, z1, N);
194
5.84M
   word z_carry = bigint_add2_nc(z + N2, N, ws1, N);
195
5.84M
196
5.84M
   z_carry += bigint_add2_nc(z + N + N2, N2, &ws_carry, 1);
197
5.84M
   bigint_add2_nc(z + N + N2, N2, &z_carry, 1);
198
5.84M
199
5.84M
   /*
200
5.84M
   * This is only actually required if cmp (result of bigint_sub_abs) is != 0,
201
5.84M
   * however if cmp==0 then ws0[0:N] == 0 and avoiding the jump hides a
202
5.84M
   * timing channel.
203
5.84M
   */
204
5.84M
   bigint_sub2(z + N2, 2*N-N2, ws0, N);
205
5.84M
   }
206
207
/*
208
* Pick a good size for the Karatsuba multiply
209
*/
210
size_t karatsuba_size(size_t z_size,
211
                      size_t x_size, size_t x_sw,
212
                      size_t y_size, size_t y_sw)
213
1.68M
   {
214
1.68M
   if(x_sw > x_size || x_sw > y_size || y_sw > x_size || y_sw > y_size)
215
418
      return 0;
216
1.68M
217
1.68M
   if(((x_size == x_sw) && (x_size % 2)) ||
218
1.68M
      ((y_size == y_sw) && (y_size % 2)))
219
0
      return 0;
220
1.68M
221
1.68M
   const size_t start = (x_sw > y_sw) ? x_sw : y_sw;
222
1.68M
   const size_t end = (x_size < y_size) ? x_size : y_size;
223
1.68M
224
1.68M
   if(start == end)
225
990k
      {
226
990k
      if(start % 2)
227
0
         return 0;
228
990k
      return start;
229
990k
      }
230
697k
231
975k
   for(size_t j = start; j <= end; ++j)
232
975k
      {
233
975k
      if(j % 2)
234
278k
         continue;
235
697k
236
697k
      if(2*j > z_size)
237
276k
         return 0;
238
421k
239
421k
      if(x_sw <= j && j <= x_size && y_sw <= j && j <= y_size)
240
421k
         {
241
421k
         if(j % 4 == 2 &&
242
421k
            (j+2) <= x_size && (j+2) <= y_size && 2*(j+2) <= z_size)
243
380
            return j+2;
244
420k
         return j;
245
420k
         }
246
421k
      }
247
697k
248
697k
   return 0;
249
697k
   }
250
251
/*
252
* Pick a good size for the Karatsuba squaring
253
*/
254
size_t karatsuba_size(size_t z_size, size_t x_size, size_t x_sw)
255
4.43M
   {
256
4.43M
   if(x_sw == x_size)
257
12.9k
      {
258
12.9k
      if(x_sw % 2)
259
0
         return 0;
260
12.9k
      return x_sw;
261
12.9k
      }
262
4.41M
263
5.23M
   for(size_t j = x_sw; j <= x_size; ++j)
264
5.23M
      {
265
5.23M
      if(j % 2)
266
817k
         continue;
267
4.41M
268
4.41M
      if(2*j > z_size)
269
293k
         return 0;
270
4.12M
271
4.12M
      if(j % 4 == 2 && (j+2) <= x_size && 2*(j+2) <= z_size)
272
80
         return j+2;
273
4.12M
      return j;
274
4.12M
      }
275
4.41M
276
4.41M
   return 0;
277
4.41M
   }
278
279
template<size_t SZ>
280
inline bool sized_for_comba_mul(size_t x_sw, size_t x_size,
281
                                size_t y_sw, size_t y_size,
282
                                size_t z_size)
283
500M
   {
284
500M
   return (x_sw <= SZ && x_size >= SZ &&
285
500M
           y_sw <= SZ && y_size >= SZ &&
286
500M
           z_size >= 2*SZ);
287
500M
   }
mp_karat.cpp:bool Botan::(anonymous namespace)::sized_for_comba_mul<4ul>(unsigned long, unsigned long, unsigned long, unsigned long, unsigned long)
Line
Count
Source
283
173M
   {
284
173M
   return (x_sw <= SZ && x_size >= SZ &&
285
173M
           y_sw <= SZ && y_size >= SZ &&
286
173M
           z_size >= 2*SZ);
287
173M
   }
mp_karat.cpp:bool Botan::(anonymous namespace)::sized_for_comba_mul<6ul>(unsigned long, unsigned long, unsigned long, unsigned long, unsigned long)
Line
Count
Source
283
128M
   {
284
128M
   return (x_sw <= SZ && x_size >= SZ &&
285
128M
           y_sw <= SZ && y_size >= SZ &&
286
128M
           z_size >= 2*SZ);
287
128M
   }
mp_karat.cpp:bool Botan::(anonymous namespace)::sized_for_comba_mul<8ul>(unsigned long, unsigned long, unsigned long, unsigned long, unsigned long)
Line
Count
Source
283
95.1M
   {
284
95.1M
   return (x_sw <= SZ && x_size >= SZ &&
285
95.1M
           y_sw <= SZ && y_size >= SZ &&
286
95.1M
           z_size >= 2*SZ);
287
95.1M
   }
mp_karat.cpp:bool Botan::(anonymous namespace)::sized_for_comba_mul<9ul>(unsigned long, unsigned long, unsigned long, unsigned long, unsigned long)
Line
Count
Source
283
85.9M
   {
284
85.9M
   return (x_sw <= SZ && x_size >= SZ &&
285
85.9M
           y_sw <= SZ && y_size >= SZ &&
286
85.9M
           z_size >= 2*SZ);
287
85.9M
   }
mp_karat.cpp:bool Botan::(anonymous namespace)::sized_for_comba_mul<16ul>(unsigned long, unsigned long, unsigned long, unsigned long, unsigned long)
Line
Count
Source
283
8.91M
   {
284
8.91M
   return (x_sw <= SZ && x_size >= SZ &&
285
8.91M
           y_sw <= SZ && y_size >= SZ &&
286
8.91M
           z_size >= 2*SZ);
287
8.91M
   }
mp_karat.cpp:bool Botan::(anonymous namespace)::sized_for_comba_mul<24ul>(unsigned long, unsigned long, unsigned long, unsigned long, unsigned long)
Line
Count
Source
283
8.67M
   {
284
8.67M
   return (x_sw <= SZ && x_size >= SZ &&
285
8.67M
           y_sw <= SZ && y_size >= SZ &&
286
8.67M
           z_size >= 2*SZ);
287
8.67M
   }
288
289
template<size_t SZ>
290
inline bool sized_for_comba_sqr(size_t x_sw, size_t x_size,
291
                                size_t z_size)
292
454M
   {
293
454M
   return (x_sw <= SZ && x_size >= SZ && z_size >= 2*SZ);
294
454M
   }
mp_karat.cpp:bool Botan::(anonymous namespace)::sized_for_comba_sqr<4ul>(unsigned long, unsigned long, unsigned long)
Line
Count
Source
292
157M
   {
293
157M
   return (x_sw <= SZ && x_size >= SZ && z_size >= 2*SZ);
294
157M
   }
mp_karat.cpp:bool Botan::(anonymous namespace)::sized_for_comba_sqr<6ul>(unsigned long, unsigned long, unsigned long)
Line
Count
Source
292
120M
   {
293
120M
   return (x_sw <= SZ && x_size >= SZ && z_size >= 2*SZ);
294
120M
   }
mp_karat.cpp:bool Botan::(anonymous namespace)::sized_for_comba_sqr<8ul>(unsigned long, unsigned long, unsigned long)
Line
Count
Source
292
86.8M
   {
293
86.8M
   return (x_sw <= SZ && x_size >= SZ && z_size >= 2*SZ);
294
86.8M
   }
mp_karat.cpp:bool Botan::(anonymous namespace)::sized_for_comba_sqr<9ul>(unsigned long, unsigned long, unsigned long)
Line
Count
Source
292
79.2M
   {
293
79.2M
   return (x_sw <= SZ && x_size >= SZ && z_size >= 2*SZ);
294
79.2M
   }
mp_karat.cpp:bool Botan::(anonymous namespace)::sized_for_comba_sqr<16ul>(unsigned long, unsigned long, unsigned long)
Line
Count
Source
292
5.39M
   {
293
5.39M
   return (x_sw <= SZ && x_size >= SZ && z_size >= 2*SZ);
294
5.39M
   }
mp_karat.cpp:bool Botan::(anonymous namespace)::sized_for_comba_sqr<24ul>(unsigned long, unsigned long, unsigned long)
Line
Count
Source
292
5.17M
   {
293
5.17M
   return (x_sw <= SZ && x_size >= SZ && z_size >= 2*SZ);
294
5.17M
   }
295
296
}
297
298
void bigint_mul(word z[], size_t z_size,
299
                const word x[], size_t x_size, size_t x_sw,
300
                const word y[], size_t y_size, size_t y_sw,
301
                word workspace[], size_t ws_size)
302
173M
   {
303
173M
   clear_mem(z, z_size);
304
173M
305
173M
   if(x_sw == 1)
306
186k
      {
307
186k
      bigint_linmul3(z, y, y_sw, x[0]);
308
186k
      }
309
173M
   else if(y_sw == 1)
310
1.62k
      {
311
1.62k
      bigint_linmul3(z, x, x_sw, y[0]);
312
1.62k
      }
313
173M
   else if(sized_for_comba_mul<4>(x_sw, x_size, y_sw, y_size, z_size))
314
45.4M
      {
315
45.4M
      bigint_comba_mul4(z, x, y);
316
45.4M
      }
317
128M
   else if(sized_for_comba_mul<6>(x_sw, x_size, y_sw, y_size, z_size))
318
33.1M
      {
319
33.1M
      bigint_comba_mul6(z, x, y);
320
33.1M
      }
321
95.1M
   else if(sized_for_comba_mul<8>(x_sw, x_size, y_sw, y_size, z_size))
322
9.14M
      {
323
9.14M
      bigint_comba_mul8(z, x, y);
324
9.14M
      }
325
85.9M
   else if(sized_for_comba_mul<9>(x_sw, x_size, y_sw, y_size, z_size))
326
77.0M
      {
327
77.0M
      bigint_comba_mul9(z, x, y);
328
77.0M
      }
329
8.91M
   else if(sized_for_comba_mul<16>(x_sw, x_size, y_sw, y_size, z_size))
330
245k
      {
331
245k
      bigint_comba_mul16(z, x, y);
332
245k
      }
333
8.67M
   else if(sized_for_comba_mul<24>(x_sw, x_size, y_sw, y_size, z_size))
334
191k
      {
335
191k
      bigint_comba_mul24(z, x, y);
336
191k
      }
337
8.47M
   else if(x_sw < KARATSUBA_MULTIPLY_THRESHOLD ||
338
8.47M
           y_sw < KARATSUBA_MULTIPLY_THRESHOLD ||
339
8.47M
           !workspace)
340
6.79M
      {
341
6.79M
      basecase_mul(z, z_size, x, x_sw, y, y_sw);
342
6.79M
      }
343
1.68M
   else
344
1.68M
      {
345
1.68M
      const size_t N = karatsuba_size(z_size, x_size, x_sw, y_size, y_sw);
346
1.68M
347
1.68M
      if(N && z_size >= 2*N && ws_size >= 2*N)
348
1.40M
         karatsuba_mul(z, x, y, N, workspace);
349
278k
      else
350
278k
         basecase_mul(z, z_size, x, x_sw, y, y_sw);
351
1.68M
      }
352
173M
   }
353
354
/*
355
* Squaring Algorithm Dispatcher
356
*/
357
void bigint_sqr(word z[], size_t z_size,
358
                const word x[], size_t x_size, size_t x_sw,
359
                word workspace[], size_t ws_size)
360
157M
   {
361
157M
   clear_mem(z, z_size);
362
157M
363
157M
   BOTAN_ASSERT(z_size/2 >= x_sw, "Output size is sufficient");
364
157M
365
157M
   if(x_sw == 1)
366
192k
      {
367
192k
      bigint_linmul3(z, x, x_sw, x[0]);
368
192k
      }
369
157M
   else if(sized_for_comba_sqr<4>(x_sw, x_size, z_size))
370
37.3M
      {
371
37.3M
      bigint_comba_sqr4(z, x);
372
37.3M
      }
373
120M
   else if(sized_for_comba_sqr<6>(x_sw, x_size, z_size))
374
33.5M
      {
375
33.5M
      bigint_comba_sqr6(z, x);
376
33.5M
      }
377
86.8M
   else if(sized_for_comba_sqr<8>(x_sw, x_size, z_size))
378
7.59M
      {
379
7.59M
      bigint_comba_sqr8(z, x);
380
7.59M
      }
381
79.2M
   else if(sized_for_comba_sqr<9>(x_sw, x_size, z_size))
382
73.8M
      {
383
73.8M
      bigint_comba_sqr9(z, x);
384
73.8M
      }
385
5.39M
   else if(sized_for_comba_sqr<16>(x_sw, x_size, z_size))
386
226k
      {
387
226k
      bigint_comba_sqr16(z, x);
388
226k
      }
389
5.17M
   else if(sized_for_comba_sqr<24>(x_sw, x_size, z_size))
390
158k
      {
391
158k
      bigint_comba_sqr24(z, x);
392
158k
      }
393
5.01M
   else if(x_size < KARATSUBA_SQUARE_THRESHOLD || !workspace)
394
582k
      {
395
582k
      basecase_sqr(z, z_size, x, x_sw);
396
582k
      }
397
4.43M
   else
398
4.43M
      {
399
4.43M
      const size_t N = karatsuba_size(z_size, x_size, x_sw);
400
4.43M
401
4.43M
      if(N && z_size >= 2*N && ws_size >= 2*N)
402
4.13M
         karatsuba_sqr(z, x, N, workspace);
403
293k
      else
404
293k
         basecase_sqr(z, z_size, x, x_sw);
405
4.43M
      }
406
157M
   }
407
408
}