Coverage Report

Created: 2020-05-23 13:54

/src/botan/src/lib/pubkey/curve25519/donna.cpp
Line
Count
Source
1
/*
2
* Based on curve25519-donna-c64.c from github.com/agl/curve25519-donna
3
* revision 80ad9b9930c9baef5829dd2a235b6b7646d32a8e
4
*
5
* Further changes
6
* (C) 2014,2018 Jack Lloyd
7
*
8
* Botan is released under the Simplified BSD License (see license.txt)
9
*/
10
11
/* Copyright 2008, Google Inc.
12
* All rights reserved.
13
*
14
* Code released into the public domain.
15
*
16
* curve25519-donna: Curve25519 elliptic curve, public key function
17
*
18
* https://code.google.com/p/curve25519-donna/
19
*
20
* Adam Langley <agl@imperialviolet.org>
21
*
22
* Derived from public domain C code by Daniel J. Bernstein <djb@cr.yp.to>
23
*
24
* More information about curve25519 can be found here
25
*   https://cr.yp.to/ecdh.html
26
*
27
* djb's sample implementation of curve25519 is written in a special assembly
28
* language called qhasm and uses the floating point registers.
29
*
30
* This is, almost, a clean room reimplementation from the curve25519 paper. It
31
* uses many of the tricks described therein. Only the crecip function is taken
32
* from the sample implementation.
33
*/
34
35
#include <botan/curve25519.h>
36
#include <botan/mul128.h>
37
#include <botan/internal/ct_utils.h>
38
#include <botan/internal/donna128.h>
39
#include <botan/loadstor.h>
40
41
namespace Botan {
42
43
namespace {
44
45
#if !defined(BOTAN_TARGET_HAS_NATIVE_UINT128)
46
typedef donna128 uint128_t;
47
#endif
48
49
/* Sum two numbers: output += in */
50
inline void fsum(uint64_t out[5], const uint64_t in[5])
51
445k
   {
52
445k
   out[0] += in[0];
53
445k
   out[1] += in[1];
54
445k
   out[2] += in[2];
55
445k
   out[3] += in[3];
56
445k
   out[4] += in[4];
57
445k
   }
58
59
/* Find the difference of two numbers: out = in - out
60
* (note the order of the arguments!)
61
*
62
* Assumes that out[i] < 2**52
63
* On return, out[i] < 2**55
64
*/
65
inline void fdifference_backwards(uint64_t out[5], const uint64_t in[5])
66
445k
   {
67
445k
   /* 152 is 19 << 3 */
68
445k
   const uint64_t two54m152 = (static_cast<uint64_t>(1) << 54) - 152;
69
445k
   const uint64_t two54m8   = (static_cast<uint64_t>(1) << 54) - 8;
70
445k
71
445k
   out[0] = in[0] + two54m152 - out[0];
72
445k
   out[1] = in[1] + two54m8 - out[1];
73
445k
   out[2] = in[2] + two54m8 - out[2];
74
445k
   out[3] = in[3] + two54m8 - out[3];
75
445k
   out[4] = in[4] + two54m8 - out[4];
76
445k
   }
77
78
inline void fadd_sub(uint64_t x[5],
79
                     uint64_t y[5])
80
334k
   {
81
334k
   // TODO merge these and avoid the tmp array
82
334k
   uint64_t tmp[5];
83
334k
   copy_mem(tmp, y, 5);
84
334k
   fsum(y, x);
85
334k
   fdifference_backwards(x, tmp);  // does x - z
86
334k
   }
87
88
/* Multiply a number by a scalar: out = in * scalar */
89
inline void fscalar_product(uint64_t out[5], const uint64_t in[5], const uint64_t scalar)
90
111k
   {
91
111k
   uint128_t a = uint128_t(in[0]) * scalar;
92
111k
   out[0] = a & 0x7ffffffffffff;
93
111k
94
111k
   a = uint128_t(in[1]) * scalar + carry_shift(a, 51);
95
111k
   out[1] = a & 0x7ffffffffffff;
96
111k
97
111k
   a = uint128_t(in[2]) * scalar + carry_shift(a, 51);
98
111k
   out[2] = a & 0x7ffffffffffff;
99
111k
100
111k
   a = uint128_t(in[3]) * scalar + carry_shift(a, 51);
101
111k
   out[3] = a & 0x7ffffffffffff;
102
111k
103
111k
   a = uint128_t(in[4]) * scalar + carry_shift(a, 51);
104
111k
   out[4] = a & 0x7ffffffffffff;
105
111k
106
111k
   out[0] += carry_shift(a, 51) * 19;
107
111k
   }
108
109
/* Multiply two numbers: out = in2 * in
110
*
111
* out must be distinct to both inputs. The inputs are reduced coefficient
112
* form, the output is not.
113
*
114
* Assumes that in[i] < 2**55 and likewise for in2.
115
* On return, out[i] < 2**52
116
*/
117
inline void fmul(uint64_t out[5], const uint64_t in[5], const uint64_t in2[5])
118
562k
   {
119
562k
   const uint128_t s0 = in2[0];
120
562k
   const uint128_t s1 = in2[1];
121
562k
   const uint128_t s2 = in2[2];
122
562k
   const uint128_t s3 = in2[3];
123
562k
   const uint128_t s4 = in2[4];
124
562k
125
562k
   uint64_t r0 = in[0];
126
562k
   uint64_t r1 = in[1];
127
562k
   uint64_t r2 = in[2];
128
562k
   uint64_t r3 = in[3];
129
562k
   uint64_t r4 = in[4];
130
562k
131
562k
   uint128_t t0 = r0 * s0;
132
562k
   uint128_t t1 = r0 * s1 + r1 * s0;
133
562k
   uint128_t t2 = r0 * s2 + r2 * s0 + r1 * s1;
134
562k
   uint128_t t3 = r0 * s3 + r3 * s0 + r1 * s2 + r2 * s1;
135
562k
   uint128_t t4 = r0 * s4 + r4 * s0 + r3 * s1 + r1 * s3 + r2 * s2;
136
562k
137
562k
   r4 *= 19;
138
562k
   r1 *= 19;
139
562k
   r2 *= 19;
140
562k
   r3 *= 19;
141
562k
142
562k
   t0 += r4 * s1 + r1 * s4 + r2 * s3 + r3 * s2;
143
562k
   t1 += r4 * s2 + r2 * s4 + r3 * s3;
144
562k
   t2 += r4 * s3 + r3 * s4;
145
562k
   t3 += r4 * s4;
146
562k
147
562k
   r0 = t0 & 0x7ffffffffffff; t1 += carry_shift(t0, 51);
148
562k
   r1 = t1 & 0x7ffffffffffff; t2 += carry_shift(t1, 51);
149
562k
   r2 = t2 & 0x7ffffffffffff; t3 += carry_shift(t2, 51);
150
562k
   r3 = t3 & 0x7ffffffffffff; t4 += carry_shift(t3, 51);
151
562k
   r4 = t4 & 0x7ffffffffffff; uint64_t c = carry_shift(t4, 51);
152
562k
153
562k
   r0 += c * 19; c = r0 >> 51; r0 = r0 & 0x7ffffffffffff;
154
562k
   r1 += c;      c = r1 >> 51; r1 = r1 & 0x7ffffffffffff;
155
562k
   r2 += c;
156
562k
157
562k
   out[0] = r0;
158
562k
   out[1] = r1;
159
562k
   out[2] = r2;
160
562k
   out[3] = r3;
161
562k
   out[4] = r4;
162
562k
   }
163
164
inline void fsquare(uint64_t out[5], const uint64_t in[5], size_t count = 1)
165
450k
   {
166
450k
   uint64_t r0 = in[0];
167
450k
   uint64_t r1 = in[1];
168
450k
   uint64_t r2 = in[2];
169
450k
   uint64_t r3 = in[3];
170
450k
   uint64_t r4 = in[4];
171
450k
172
1.00M
   for(size_t i = 0; i != count; ++i)
173
555k
      {
174
555k
      const uint64_t d0 = r0 * 2;
175
555k
      const uint64_t d1 = r1 * 2;
176
555k
      const uint64_t d2 = r2 * 2 * 19;
177
555k
      const uint64_t d419 = r4 * 19;
178
555k
      const uint64_t d4 = d419 * 2;
179
555k
180
555k
      uint128_t t0 = uint128_t(r0) * r0 + uint128_t(d4) * r1 + uint128_t(d2) * (r3     );
181
555k
      uint128_t t1 = uint128_t(d0) * r1 + uint128_t(d4) * r2 + uint128_t(r3) * (r3 * 19);
182
555k
      uint128_t t2 = uint128_t(d0) * r2 + uint128_t(r1) * r1 + uint128_t(d4) * (r3     );
183
555k
      uint128_t t3 = uint128_t(d0) * r3 + uint128_t(d1) * r2 + uint128_t(r4) * (d419   );
184
555k
      uint128_t t4 = uint128_t(d0) * r4 + uint128_t(d1) * r3 + uint128_t(r2) * (r2     );
185
555k
186
555k
      r0 = t0 & 0x7ffffffffffff; t1 += carry_shift(t0, 51);
187
555k
      r1 = t1 & 0x7ffffffffffff; t2 += carry_shift(t1, 51);
188
555k
      r2 = t2 & 0x7ffffffffffff; t3 += carry_shift(t2, 51);
189
555k
      r3 = t3 & 0x7ffffffffffff; t4 += carry_shift(t3, 51);
190
555k
      r4 = t4 & 0x7ffffffffffff; uint64_t c = carry_shift(t4, 51);
191
555k
192
555k
      r0 += c * 19; c = r0 >> 51; r0 = r0 & 0x7ffffffffffff;
193
555k
      r1 += c;      c = r1 >> 51; r1 = r1 & 0x7ffffffffffff;
194
555k
      r2 += c;
195
555k
      }
196
450k
197
450k
   out[0] = r0;
198
450k
   out[1] = r1;
199
450k
   out[2] = r2;
200
450k
   out[3] = r3;
201
450k
   out[4] = r4;
202
450k
   }
203
204
/* Take a little-endian, 32-byte number and expand it into polynomial form */
205
inline void fexpand(uint64_t *out, const uint8_t *in)
206
435
   {
207
435
   out[0] = load_le<uint64_t>(in, 0) & 0x7ffffffffffff;
208
435
   out[1] = (load_le<uint64_t>(in+6, 0) >> 3) & 0x7ffffffffffff;
209
435
   out[2] = (load_le<uint64_t>(in+12, 0) >> 6) & 0x7ffffffffffff;
210
435
   out[3] = (load_le<uint64_t>(in+19, 0) >> 1) & 0x7ffffffffffff;
211
435
   out[4] = (load_le<uint64_t>(in+24, 0) >> 12) & 0x7ffffffffffff;
212
435
   }
213
214
/* Take a fully reduced polynomial form number and contract it into a
215
* little-endian, 32-byte array
216
*/
217
inline void fcontract(uint8_t *out, const uint64_t input[5])
218
435
   {
219
435
   uint128_t t0 = input[0];
220
435
   uint128_t t1 = input[1];
221
435
   uint128_t t2 = input[2];
222
435
   uint128_t t3 = input[3];
223
435
   uint128_t t4 = input[4];
224
435
225
1.30k
   for(size_t i = 0; i != 2; ++i)
226
870
      {
227
870
      t1 += t0 >> 51;        t0 &= 0x7ffffffffffff;
228
870
      t2 += t1 >> 51;        t1 &= 0x7ffffffffffff;
229
870
      t3 += t2 >> 51;        t2 &= 0x7ffffffffffff;
230
870
      t4 += t3 >> 51;        t3 &= 0x7ffffffffffff;
231
870
      t0 += (t4 >> 51) * 19; t4 &= 0x7ffffffffffff;
232
870
      }
233
435
234
435
   /* now t is between 0 and 2^255-1, properly carried. */
235
435
   /* case 1: between 0 and 2^255-20. case 2: between 2^255-19 and 2^255-1. */
236
435
237
435
   t0 += 19;
238
435
239
435
   t1 += t0 >> 51; t0 &= 0x7ffffffffffff;
240
435
   t2 += t1 >> 51; t1 &= 0x7ffffffffffff;
241
435
   t3 += t2 >> 51; t2 &= 0x7ffffffffffff;
242
435
   t4 += t3 >> 51; t3 &= 0x7ffffffffffff;
243
435
   t0 += (t4 >> 51) * 19; t4 &= 0x7ffffffffffff;
244
435
245
435
   /* now between 19 and 2^255-1 in both cases, and offset by 19. */
246
435
247
435
   t0 += 0x8000000000000 - 19;
248
435
   t1 += 0x8000000000000 - 1;
249
435
   t2 += 0x8000000000000 - 1;
250
435
   t3 += 0x8000000000000 - 1;
251
435
   t4 += 0x8000000000000 - 1;
252
435
253
435
   /* now between 2^255 and 2^256-20, and offset by 2^255. */
254
435
255
435
   t1 += t0 >> 51; t0 &= 0x7ffffffffffff;
256
435
   t2 += t1 >> 51; t1 &= 0x7ffffffffffff;
257
435
   t3 += t2 >> 51; t2 &= 0x7ffffffffffff;
258
435
   t4 += t3 >> 51; t3 &= 0x7ffffffffffff;
259
435
   t4 &= 0x7ffffffffffff;
260
435
261
435
   store_le(out,
262
435
            combine_lower(t0,  0, t1, 51),
263
435
            combine_lower(t1, 13, t2, 38),
264
435
            combine_lower(t2, 26, t3, 25),
265
435
            combine_lower(t3, 39, t4, 12));
266
435
   }
267
268
/* Input: Q, Q', Q-Q'
269
* Out: 2Q, Q+Q'
270
*
271
*   result.two_q (2*Q): long form
272
*   result.q_plus_q_dash (Q + Q): long form
273
*   in_q: short form, destroyed
274
*   in_q_dash: short form, destroyed
275
*   in_q_minus_q_dash: short form, preserved
276
*/
277
void fmonty(uint64_t result_two_q_x[5],
278
            uint64_t result_two_q_z[5],
279
            uint64_t result_q_plus_q_dash_x[5],
280
            uint64_t result_q_plus_q_dash_z[5],
281
            uint64_t in_q_x[5],
282
            uint64_t in_q_z[5],
283
            uint64_t in_q_dash_x[5],
284
            uint64_t in_q_dash_z[5],
285
            const uint64_t q_minus_q_dash[5])
286
111k
   {
287
111k
   uint64_t zzz[5];
288
111k
   uint64_t xx[5];
289
111k
   uint64_t zz[5];
290
111k
   uint64_t xxprime[5];
291
111k
   uint64_t zzprime[5];
292
111k
   uint64_t zzzprime[5];
293
111k
294
111k
   fadd_sub(in_q_z, in_q_x);
295
111k
   fadd_sub(in_q_dash_z, in_q_dash_x);
296
111k
297
111k
   fmul(xxprime, in_q_dash_x, in_q_z);
298
111k
   fmul(zzprime, in_q_dash_z, in_q_x);
299
111k
300
111k
   fadd_sub(zzprime, xxprime);
301
111k
302
111k
   fsquare(result_q_plus_q_dash_x, xxprime);
303
111k
   fsquare(zzzprime, zzprime);
304
111k
   fmul(result_q_plus_q_dash_z, zzzprime, q_minus_q_dash);
305
111k
306
111k
   fsquare(xx, in_q_x);
307
111k
   fsquare(zz, in_q_z);
308
111k
   fmul(result_two_q_x, xx, zz);
309
111k
310
111k
   fdifference_backwards(zz, xx);  // does zz = xx - zz
311
111k
   fscalar_product(zzz, zz, 121665);
312
111k
   fsum(zzz, xx);
313
111k
314
111k
   fmul(result_two_q_z, zz, zzz);
315
111k
   }
316
317
/*
318
* Maybe swap the contents of two uint64_t arrays (@a and @b),
319
* Param @iswap is assumed to be either 0 or 1
320
*
321
* This function performs the swap without leaking any side-channel
322
* information.
323
*/
324
inline void swap_conditional(uint64_t a[5], uint64_t b[5],
325
                             uint64_t c[5], uint64_t d[5],
326
                             uint64_t iswap)
327
125k
   {
328
125k
   const uint64_t swap = 0 - iswap;
329
125k
330
751k
   for(size_t i = 0; i < 5; ++i)
331
626k
      {
332
626k
      const uint64_t x0 = swap & (a[i] ^ b[i]);
333
626k
      const uint64_t x1 = swap & (c[i] ^ d[i]);
334
626k
      a[i] ^= x0;
335
626k
      b[i] ^= x0;
336
626k
      c[i] ^= x1;
337
626k
      d[i] ^= x1;
338
626k
      }
339
125k
   }
340
341
/* Calculates nQ where Q is the x-coordinate of a point on the curve
342
*
343
*   resultx/resultz: the x/z coordinate of the resulting curve point (short form)
344
*   n: a little endian, 32-byte number
345
*   q: a point of the curve (short form)
346
*/
347
void cmult(uint64_t resultx[5], uint64_t resultz[5], const uint8_t n[32], const uint64_t q[5])
348
435
   {
349
435
   uint64_t a[5] = {0}; // nqpqx
350
435
   uint64_t b[5] = {1}; // npqpz
351
435
   uint64_t c[5] = {1}; // nqx
352
435
   uint64_t d[5] = {0}; // nqz
353
435
   uint64_t e[5] = {0}; // npqqx2
354
435
   uint64_t f[5] = {1}; // npqqz2
355
435
   uint64_t g[5] = {0}; // nqx2
356
435
   uint64_t h[5] = {1}; // nqz2
357
435
358
435
   copy_mem(a, q, 5);
359
435
360
14.3k
   for(size_t i = 0; i < 32; ++i)
361
13.9k
      {
362
13.9k
      const uint64_t bit0 = (n[31 - i] >> 7) & 1;
363
13.9k
      const uint64_t bit1 = (n[31 - i] >> 6) & 1;
364
13.9k
      const uint64_t bit2 = (n[31 - i] >> 5) & 1;
365
13.9k
      const uint64_t bit3 = (n[31 - i] >> 4) & 1;
366
13.9k
      const uint64_t bit4 = (n[31 - i] >> 3) & 1;
367
13.9k
      const uint64_t bit5 = (n[31 - i] >> 2) & 1;
368
13.9k
      const uint64_t bit6 = (n[31 - i] >> 1) & 1;
369
13.9k
      const uint64_t bit7 = (n[31 - i] >> 0) & 1;
370
13.9k
371
13.9k
      swap_conditional(c, a, d, b, bit0);
372
13.9k
      fmonty(g, h, e, f, c, d, a, b, q);
373
13.9k
374
13.9k
      swap_conditional(g, e, h, f, bit0 ^ bit1);
375
13.9k
      fmonty(c, d, a, b, g, h, e, f, q);
376
13.9k
377
13.9k
      swap_conditional(c, a, d, b, bit1 ^ bit2);
378
13.9k
      fmonty(g, h, e, f, c, d, a, b, q);
379
13.9k
380
13.9k
      swap_conditional(g, e, h, f, bit2 ^ bit3);
381
13.9k
      fmonty(c, d, a, b, g, h, e, f, q);
382
13.9k
383
13.9k
      swap_conditional(c, a, d, b, bit3 ^ bit4);
384
13.9k
      fmonty(g, h, e, f, c, d, a, b, q);
385
13.9k
386
13.9k
      swap_conditional(g, e, h, f, bit4 ^ bit5);
387
13.9k
      fmonty(c, d, a, b, g, h, e, f, q);
388
13.9k
389
13.9k
      swap_conditional(c, a, d, b, bit5 ^ bit6);
390
13.9k
      fmonty(g, h, e, f, c, d, a, b, q);
391
13.9k
392
13.9k
      swap_conditional(g, e, h, f, bit6 ^ bit7);
393
13.9k
      fmonty(c, d, a, b, g, h, e, f, q);
394
13.9k
395
13.9k
      swap_conditional(c, a, d, b, bit7);
396
13.9k
      }
397
435
398
435
   copy_mem(resultx, c, 5);
399
435
   copy_mem(resultz, d, 5);
400
435
   }
401
402
403
// -----------------------------------------------------------------------------
404
// Shamelessly copied from djb's code, tightened a little
405
// -----------------------------------------------------------------------------
406
void crecip(uint64_t out[5], const uint64_t z[5])
407
435
   {
408
435
   uint64_t a[5];
409
435
   uint64_t b[5];
410
435
   uint64_t c[5];
411
435
   uint64_t t0[5];
412
435
413
435
   fsquare(a, z);       // 2
414
435
   fsquare(t0, a, 2);   // 8
415
435
   fmul(b, t0, z);      // 9
416
435
   fmul(a, b, a);       // 11
417
435
   fsquare(t0, a);      // 22
418
435
   fmul(b, t0, b);      // 2^5 - 2^0 = 31
419
435
   fsquare(t0, b, 5);   // 2^10 - 2^5
420
435
   fmul(b, t0, b);      // 2^10 - 2^0
421
435
   fsquare(t0, b, 10);  // 2^20 - 2^10
422
435
   fmul(c, t0, b);      // 2^20 - 2^0
423
435
   fsquare(t0, c, 20);  // 2^40 - 2^20
424
435
   fmul(t0, t0, c);     // 2^40 - 2^0
425
435
   fsquare(t0, t0, 10); // 2^50 - 2^10
426
435
   fmul(b, t0, b);      // 2^50 - 2^0
427
435
   fsquare(t0, b, 50);  // 2^100 - 2^50
428
435
   fmul(c, t0, b);      // 2^100 - 2^0
429
435
   fsquare(t0, c, 100); // 2^200 - 2^100
430
435
   fmul(t0, t0, c);     // 2^200 - 2^0
431
435
   fsquare(t0, t0, 50); // 2^250 - 2^50
432
435
   fmul(t0, t0, b);     // 2^250 - 2^0
433
435
   fsquare(t0, t0, 5);  // 2^255 - 2^5
434
435
   fmul(out, t0, a);    // 2^255 - 21
435
435
   }
436
437
}
438
439
void
440
curve25519_donna(uint8_t mypublic[32], const uint8_t secret[32], const uint8_t basepoint[32])
441
435
   {
442
435
   CT::poison(secret, 32);
443
435
   CT::poison(basepoint, 32);
444
435
445
435
   uint64_t bp[5], x[5], z[5], zmone[5];
446
435
   uint8_t e[32];
447
435
448
435
   copy_mem(e, secret, 32);
449
435
   e[ 0] &= 248;
450
435
   e[31] &= 127;
451
435
   e[31] |= 64;
452
435
453
435
   fexpand(bp, basepoint);
454
435
   cmult(x, z, e, bp);
455
435
   crecip(zmone, z);
456
435
   fmul(z, x, zmone);
457
435
   fcontract(mypublic, z);
458
435
459
435
   CT::unpoison(secret, 32);
460
435
   CT::unpoison(basepoint, 32);
461
435
   CT::unpoison(mypublic, 32);
462
435
   }
463
464
}