Coverage Report

Created: 2021-02-21 07:20

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