Botan 2.1.0 Crypto and TLS for C++11

- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

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 * http://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 * http://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)
63 typedef donna128 uint128_t;
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 }
Botan::felem
limb felem[5]
Definition: donna.cpp:40
Botan::donna128
Definition: donna128.h:15
Botan::CT::poison
void poison(const T *p, size_t n)
Definition: ct_utils.h:46
Botan::carry_shift
uint64_t carry_shift(const donna128 &a, size_t shift)
Definition: donna128.h:108
Botan::combine_lower
uint64_t combine_lower(const donna128 &a, size_t s1, const donna128 &b, size_t s2)
Definition: donna128.h:113
Botan::u8
uint8_t u8
Definition: donna.cpp:38
Botan::load_le< uint64_t >
uint64_t load_le< uint64_t >(const uint8_t in[], size_t off)
Definition: loadstor.h:243
Botan::copy_mem
void copy_mem(T *out, const T *in, size_t n)
Definition: mem_ops.h:68
Botan
Definition: alg_id.cpp:13
Botan::limb
uint64_t limb
Definition: donna.cpp:39
Botan::CT::unpoison
void unpoison(const T *p, size_t n)
Definition: ct_utils.h:57
Botan::uint128_t
donna128 uint128_t
Definition: donna.cpp:63
Botan::curve25519_donna
void BOTAN_DLL curve25519_donna(uint8_t mypublic[32], const uint8_t secret[32], const uint8_t basepoint[32])
Botan::store_le
void store_le(uint16_t in, uint8_t out[2])
Definition: loadstor.h:457