Coverage Report

Created: 2020-05-23 13:54

/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
6.22M
   {
29
6.22M
   if(z_size < x_size + y_size)
30
0
      throw Invalid_Argument("basecase_mul z_size too small");
31
6.22M
32
6.22M
   const size_t x_size_8 = x_size - (x_size % 8);
33
6.22M
34
6.22M
   clear_mem(z, z_size);
35
6.22M
36
78.4M
   for(size_t i = 0; i != y_size; ++i)
37
72.2M
      {
38
72.2M
      const word y_i = y[i];
39
72.2M
40
72.2M
      word carry = 0;
41
72.2M
42
234M
      for(size_t j = 0; j != x_size_8; j += 8)
43
162M
         carry = word8_madd3(z + i + j, x + j, y_i, carry);
44
72.2M
45
345M
      for(size_t j = x_size_8; j != x_size; ++j)
46
273M
         z[i+j] = word_madd3(x[j], y_i, z[i+j], &carry);
47
72.2M
48
72.2M
      z[x_size+i] = carry;
49
72.2M
      }
50
6.22M
   }
51
52
void basecase_sqr(word z[], size_t z_size,
53
                  const word x[], size_t x_size)
54
1.29M
   {
55
1.29M
   if(z_size < 2*x_size)
56
0
      throw Invalid_Argument("basecase_sqr z_size too small");
57
1.29M
58
1.29M
   const size_t x_size_8 = x_size - (x_size % 8);
59
1.29M
60
1.29M
   clear_mem(z, z_size);
61
1.29M
62
23.3M
   for(size_t i = 0; i != x_size; ++i)
63
22.0M
      {
64
22.0M
      const word x_i = x[i];
65
22.0M
66
22.0M
      word carry = 0;
67
22.0M
68
79.5M
      for(size_t j = 0; j != x_size_8; j += 8)
69
57.4M
         carry = word8_madd3(z + i + j, x + j, x_i, carry);
70
22.0M
71
126M
      for(size_t j = x_size_8; j != x_size; ++j)
72
104M
         z[i+j] = word_madd3(x[j], x_i, z[i+j], &carry);
73
22.0M
74
22.0M
      z[x_size+i] = carry;
75
22.0M
      }
76
1.29M
   }
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.04M
   {
84
7.04M
   if(N < KARATSUBA_MULTIPLY_THRESHOLD || N % 2)
85
5.13M
      {
86
5.13M
      switch(N)
87
5.13M
         {
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
4.79M
         case 16:
95
4.79M
            return bigint_comba_mul16(z, x, y);
96
335k
         case 24:
97
335k
            return bigint_comba_mul24(z, x, y);
98
9.58k
         default:
99
9.58k
            return basecase_mul(z, 2*N, x, N, y, N);
100
1.90M
         }
101
1.90M
      }
102
1.90M
103
1.90M
   const size_t N2 = N / 2;
104
1.90M
105
1.90M
   const word* x0 = x;
106
1.90M
   const word* x1 = x + N2;
107
1.90M
   const word* y0 = y;
108
1.90M
   const word* y1 = y + N2;
109
1.90M
   word* z0 = z;
110
1.90M
   word* z1 = z + N;
111
1.90M
112
1.90M
   word* ws0 = workspace;
113
1.90M
   word* ws1 = workspace + N;
114
1.90M
115
1.90M
   clear_mem(workspace, 2*N);
116
1.90M
117
1.90M
   /*
118
1.90M
   * If either of cmp0 or cmp1 is zero then z0 or z1 resp is zero here,
119
1.90M
   * resulting in a no-op - z0*z1 will be equal to zero so we don't need to do
120
1.90M
   * anything, clear_mem above already set the correct result.
121
1.90M
   *
122
1.90M
   * However we ignore the result of the comparisons and always perform the
123
1.90M
   * subtractions and recursively multiply to avoid the timing channel.
124
1.90M
   */
125
1.90M
126
1.90M
   // First compute (X_lo - X_hi)*(Y_hi - Y_lo)
127
1.90M
   const auto cmp0 = bigint_sub_abs(z0, x0, x1, N2, workspace);
128
1.90M
   const auto cmp1 = bigint_sub_abs(z1, y1, y0, N2, workspace);
129
1.90M
   const auto neg_mask = ~(cmp0 ^ cmp1);
130
1.90M
131
1.90M
   karatsuba_mul(ws0, z0, z1, N2, ws1);
132
1.90M
133
1.90M
   // Compute X_lo * Y_lo
134
1.90M
   karatsuba_mul(z0, x0, y0, N2, ws1);
135
1.90M
136
1.90M
   // Compute X_hi * Y_hi
137
1.90M
   karatsuba_mul(z1, x1, y1, N2, ws1);
138
1.90M
139
1.90M
   const word ws_carry = bigint_add3_nc(ws1, z0, N, z1, N);
140
1.90M
   word z_carry = bigint_add2_nc(z + N2, N, ws1, N);
141
1.90M
142
1.90M
   z_carry += bigint_add2_nc(z + N + N2, N2, &ws_carry, 1);
143
1.90M
   bigint_add2_nc(z + N + N2, N2, &z_carry, 1);
144
1.90M
145
1.90M
   clear_mem(workspace + N, N2);
146
1.90M
147
1.90M
   bigint_cnd_add_or_sub(neg_mask, z + N2, workspace, 2*N-N2);
148
1.90M
   }
149
150
/*
151
* Karatsuba Squaring Operation
152
*/
153
void karatsuba_sqr(word z[], const word x[], size_t N, word workspace[])
154
20.2M
   {
155
20.2M
   if(N < KARATSUBA_SQUARE_THRESHOLD || N % 2)
156
14.8M
      {
157
14.8M
      switch(N)
158
14.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
13.2M
         case 16:
166
13.2M
            return bigint_comba_sqr16(z, x);
167
1.10M
         case 24:
168
1.10M
            return bigint_comba_sqr24(z, x);
169
531k
         default:
170
531k
            return basecase_sqr(z, 2*N, x, N);
171
5.42M
         }
172
5.42M
      }
173
5.42M
174
5.42M
   const size_t N2 = N / 2;
175
5.42M
176
5.42M
   const word* x0 = x;
177
5.42M
   const word* x1 = x + N2;
178
5.42M
   word* z0 = z;
179
5.42M
   word* z1 = z + N;
180
5.42M
181
5.42M
   word* ws0 = workspace;
182
5.42M
   word* ws1 = workspace + N;
183
5.42M
184
5.42M
   clear_mem(workspace, 2*N);
185
5.42M
186
5.42M
   // See comment in karatsuba_mul
187
5.42M
   bigint_sub_abs(z0, x0, x1, N2, workspace);
188
5.42M
   karatsuba_sqr(ws0, z0, N2, ws1);
189
5.42M
190
5.42M
   karatsuba_sqr(z0, x0, N2, ws1);
191
5.42M
   karatsuba_sqr(z1, x1, N2, ws1);
192
5.42M
193
5.42M
   const word ws_carry = bigint_add3_nc(ws1, z0, N, z1, N);
194
5.42M
   word z_carry = bigint_add2_nc(z + N2, N, ws1, N);
195
5.42M
196
5.42M
   z_carry += bigint_add2_nc(z + N + N2, N2, &ws_carry, 1);
197
5.42M
   bigint_add2_nc(z + N + N2, N2, &z_carry, 1);
198
5.42M
199
5.42M
   /*
200
5.42M
   * This is only actually required if cmp (result of bigint_sub_abs) is != 0,
201
5.42M
   * however if cmp==0 then ws0[0:N] == 0 and avoiding the jump hides a
202
5.42M
   * timing channel.
203
5.42M
   */
204
5.42M
   bigint_sub2(z + N2, 2*N-N2, ws0, N);
205
5.42M
   }
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.58M
   {
214
1.58M
   if(x_sw > x_size || x_sw > y_size || y_sw > x_size || y_sw > y_size)
215
391
      return 0;
216
1.58M
217
1.58M
   if(((x_size == x_sw) && (x_size % 2)) ||
218
1.58M
      ((y_size == y_sw) && (y_size % 2)))
219
0
      return 0;
220
1.58M
221
1.58M
   const size_t start = (x_sw > y_sw) ? x_sw : y_sw;
222
1.58M
   const size_t end = (x_size < y_size) ? x_size : y_size;
223
1.58M
224
1.58M
   if(start == end)
225
995k
      {
226
995k
      if(start % 2)
227
0
         return 0;
228
995k
      return start;
229
995k
      }
230
591k
231
845k
   for(size_t j = start; j <= end; ++j)
232
845k
      {
233
845k
      if(j % 2)
234
253k
         continue;
235
591k
236
591k
      if(2*j > z_size)
237
252k
         return 0;
238
339k
239
339k
      if(x_sw <= j && j <= x_size && y_sw <= j && j <= y_size)
240
339k
         {
241
339k
         if(j % 4 == 2 &&
242
339k
            (j+2) <= x_size && (j+2) <= y_size && 2*(j+2) <= z_size)
243
340
            return j+2;
244
338k
         return j;
245
338k
         }
246
339k
      }
247
591k
248
591k
   return 0;
249
591k
   }
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.23M
   {
256
4.23M
   if(x_sw == x_size)
257
13.6k
      {
258
13.6k
      if(x_sw % 2)
259
0
         return 0;
260
13.6k
      return x_sw;
261
13.6k
      }
262
4.21M
263
4.83M
   for(size_t j = x_sw; j <= x_size; ++j)
264
4.83M
      {
265
4.83M
      if(j % 2)
266
618k
         continue;
267
4.21M
268
4.21M
      if(2*j > z_size)
269
212k
         return 0;
270
4.00M
271
4.00M
      if(j % 4 == 2 && (j+2) <= x_size && 2*(j+2) <= z_size)
272
66
         return j+2;
273
4.00M
      return j;
274
4.00M
      }
275
4.21M
276
4.21M
   return 0;
277
4.21M
   }
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
436M
   {
284
436M
   return (x_sw <= SZ && x_size >= SZ &&
285
436M
           y_sw <= SZ && y_size >= SZ &&
286
436M
           z_size >= 2*SZ);
287
436M
   }
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
151M
   {
284
151M
   return (x_sw <= SZ && x_size >= SZ &&
285
151M
           y_sw <= SZ && y_size >= SZ &&
286
151M
           z_size >= 2*SZ);
287
151M
   }
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
111M
   {
284
111M
   return (x_sw <= SZ && x_size >= SZ &&
285
111M
           y_sw <= SZ && y_size >= SZ &&
286
111M
           z_size >= 2*SZ);
287
111M
   }
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
82.7M
   {
284
82.7M
   return (x_sw <= SZ && x_size >= SZ &&
285
82.7M
           y_sw <= SZ && y_size >= SZ &&
286
82.7M
           z_size >= 2*SZ);
287
82.7M
   }
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
75.7M
   {
284
75.7M
   return (x_sw <= SZ && x_size >= SZ &&
285
75.7M
           y_sw <= SZ && y_size >= SZ &&
286
75.7M
           z_size >= 2*SZ);
287
75.7M
   }
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
7.89M
   {
284
7.89M
   return (x_sw <= SZ && x_size >= SZ &&
285
7.89M
           y_sw <= SZ && y_size >= SZ &&
286
7.89M
           z_size >= 2*SZ);
287
7.89M
   }
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
7.70M
   {
284
7.70M
   return (x_sw <= SZ && x_size >= SZ &&
285
7.70M
           y_sw <= SZ && y_size >= SZ &&
286
7.70M
           z_size >= 2*SZ);
287
7.70M
   }
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
396M
   {
293
396M
   return (x_sw <= SZ && x_size >= SZ && z_size >= 2*SZ);
294
396M
   }
mp_karat.cpp:bool Botan::(anonymous namespace)::sized_for_comba_sqr<4ul>(unsigned long, unsigned long, unsigned long)
Line
Count
Source
292
136M
   {
293
136M
   return (x_sw <= SZ && x_size >= SZ && z_size >= 2*SZ);
294
136M
   }
mp_karat.cpp:bool Botan::(anonymous namespace)::sized_for_comba_sqr<6ul>(unsigned long, unsigned long, unsigned long)
Line
Count
Source
292
104M
   {
293
104M
   return (x_sw <= SZ && x_size >= SZ && z_size >= 2*SZ);
294
104M
   }
mp_karat.cpp:bool Botan::(anonymous namespace)::sized_for_comba_sqr<8ul>(unsigned long, unsigned long, unsigned long)
Line
Count
Source
292
75.5M
   {
293
75.5M
   return (x_sw <= SZ && x_size >= SZ && z_size >= 2*SZ);
294
75.5M
   }
mp_karat.cpp:bool Botan::(anonymous namespace)::sized_for_comba_sqr<9ul>(unsigned long, unsigned long, unsigned long)
Line
Count
Source
292
69.7M
   {
293
69.7M
   return (x_sw <= SZ && x_size >= SZ && z_size >= 2*SZ);
294
69.7M
   }
mp_karat.cpp:bool Botan::(anonymous namespace)::sized_for_comba_sqr<16ul>(unsigned long, unsigned long, unsigned long)
Line
Count
Source
292
5.08M
   {
293
5.08M
   return (x_sw <= SZ && x_size >= SZ && z_size >= 2*SZ);
294
5.08M
   }
mp_karat.cpp:bool Botan::(anonymous namespace)::sized_for_comba_sqr<24ul>(unsigned long, unsigned long, unsigned long)
Line
Count
Source
292
4.91M
   {
293
4.91M
   return (x_sw <= SZ && x_size >= SZ && z_size >= 2*SZ);
294
4.91M
   }
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
151M
   {
303
151M
   clear_mem(z, z_size);
304
151M
305
151M
   if(x_sw == 1)
306
222k
      {
307
222k
      bigint_linmul3(z, y, y_sw, x[0]);
308
222k
      }
309
151M
   else if(y_sw == 1)
310
1.27k
      {
311
1.27k
      bigint_linmul3(z, x, x_sw, y[0]);
312
1.27k
      }
313
151M
   else if(sized_for_comba_mul<4>(x_sw, x_size, y_sw, y_size, z_size))
314
40.2M
      {
315
40.2M
      bigint_comba_mul4(z, x, y);
316
40.2M
      }
317
111M
   else if(sized_for_comba_mul<6>(x_sw, x_size, y_sw, y_size, z_size))
318
28.4M
      {
319
28.4M
      bigint_comba_mul6(z, x, y);
320
28.4M
      }
321
82.7M
   else if(sized_for_comba_mul<8>(x_sw, x_size, y_sw, y_size, z_size))
322
7.00M
      {
323
7.00M
      bigint_comba_mul8(z, x, y);
324
7.00M
      }
325
75.7M
   else if(sized_for_comba_mul<9>(x_sw, x_size, y_sw, y_size, z_size))
326
67.8M
      {
327
67.8M
      bigint_comba_mul9(z, x, y);
328
67.8M
      }
329
7.89M
   else if(sized_for_comba_mul<16>(x_sw, x_size, y_sw, y_size, z_size))
330
190k
      {
331
190k
      bigint_comba_mul16(z, x, y);
332
190k
      }
333
7.70M
   else if(sized_for_comba_mul<24>(x_sw, x_size, y_sw, y_size, z_size))
334
156k
      {
335
156k
      bigint_comba_mul24(z, x, y);
336
156k
      }
337
7.54M
   else if(x_sw < KARATSUBA_MULTIPLY_THRESHOLD ||
338
7.54M
           y_sw < KARATSUBA_MULTIPLY_THRESHOLD ||
339
7.54M
           !workspace)
340
5.95M
      {
341
5.95M
      basecase_mul(z, z_size, x, x_sw, y, y_sw);
342
5.95M
      }
343
1.58M
   else
344
1.58M
      {
345
1.58M
      const size_t N = karatsuba_size(z_size, x_size, x_sw, y_size, y_sw);
346
1.58M
347
1.58M
      if(N && z_size >= 2*N && ws_size >= 2*N)
348
1.33M
         karatsuba_mul(z, x, y, N, workspace);
349
254k
      else
350
254k
         basecase_mul(z, z_size, x, x_sw, y, y_sw);
351
1.58M
      }
352
151M
   }
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
136M
   {
361
136M
   clear_mem(z, z_size);
362
136M
363
136M
   BOTAN_ASSERT(z_size/2 >= x_sw, "Output size is sufficient");
364
136M
365
136M
   if(x_sw == 1)
366
216k
      {
367
216k
      bigint_linmul3(z, x, x_sw, x[0]);
368
216k
      }
369
136M
   else if(sized_for_comba_sqr<4>(x_sw, x_size, z_size))
370
32.5M
      {
371
32.5M
      bigint_comba_sqr4(z, x);
372
32.5M
      }
373
104M
   else if(sized_for_comba_sqr<6>(x_sw, x_size, z_size))
374
28.6M
      {
375
28.6M
      bigint_comba_sqr6(z, x);
376
28.6M
      }
377
75.5M
   else if(sized_for_comba_sqr<8>(x_sw, x_size, z_size))
378
5.83M
      {
379
5.83M
      bigint_comba_sqr8(z, x);
380
5.83M
      }
381
69.7M
   else if(sized_for_comba_sqr<9>(x_sw, x_size, z_size))
382
64.6M
      {
383
64.6M
      bigint_comba_sqr9(z, x);
384
64.6M
      }
385
5.08M
   else if(sized_for_comba_sqr<16>(x_sw, x_size, z_size))
386
176k
      {
387
176k
      bigint_comba_sqr16(z, x);
388
176k
      }
389
4.91M
   else if(sized_for_comba_sqr<24>(x_sw, x_size, z_size))
390
130k
      {
391
130k
      bigint_comba_sqr24(z, x);
392
130k
      }
393
4.77M
   else if(x_size < KARATSUBA_SQUARE_THRESHOLD || !workspace)
394
546k
      {
395
546k
      basecase_sqr(z, z_size, x, x_sw);
396
546k
      }
397
4.23M
   else
398
4.23M
      {
399
4.23M
      const size_t N = karatsuba_size(z_size, x_size, x_sw);
400
4.23M
401
4.23M
      if(N && z_size >= 2*N && ws_size >= 2*N)
402
4.02M
         karatsuba_sqr(z, x, N, workspace);
403
212k
      else
404
212k
         basecase_sqr(z, z_size, x, x_sw);
405
4.23M
      }
406
136M
   }
407
408
}