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