Tor 0.4.9.0-alpha-dev
curve25519-donna.c
1/* Copyright 2008, Google Inc.
2 * All rights reserved.
3 *
4 * Redistribution and use in source and binary forms, with or without
5 * modification, are permitted provided that the following conditions are
6 * met:
7 *
8 * * Redistributions of source code must retain the above copyright
9 * notice, this list of conditions and the following disclaimer.
10 * * Redistributions in binary form must reproduce the above
11 * copyright notice, this list of conditions and the following disclaimer
12 * in the documentation and/or other materials provided with the
13 * distribution.
14 * * Neither the name of Google Inc. nor the names of its
15 * contributors may be used to endorse or promote products derived from
16 * this software without specific prior written permission.
17 *
18 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
19 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
20 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
21 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
22 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
23 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
24 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
25 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
26 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
27 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
28 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29 *
30 * curve25519-donna: Curve25519 elliptic curve, public key function
31 *
32 * http://code.google.com/p/curve25519-donna/
33 *
34 * Adam Langley <agl@imperialviolet.org>
35 *
36 * Derived from public domain C code by Daniel J. Bernstein <djb@cr.yp.to>
37 *
38 * More information about curve25519 can be found here
39 * http://cr.yp.to/ecdh.html
40 *
41 * djb's sample implementation of curve25519 is written in a special assembly
42 * language called qhasm and uses the floating point registers.
43 *
44 * This is, almost, a clean room reimplementation from the curve25519 paper. It
45 * uses many of the tricks described therein. Only the crecip function is taken
46 * from the sample implementation. */
47
48#include "orconfig.h"
49
50#include <string.h>
51#include "lib/cc/torint.h"
52
53typedef uint8_t u8;
54typedef int32_t s32;
55typedef int64_t limb;
56
57/* Field element representation:
58 *
59 * Field elements are written as an array of signed, 64-bit limbs, least
60 * significant first. The value of the field element is:
61 * x[0] + 2^26·x[1] + x^51·x[2] + 2^102·x[3] + ...
62 *
63 * i.e. the limbs are 26, 25, 26, 25, ... bits wide. */
64
65/* Sum two numbers: output += in */
66static void fsum(limb *output, const limb *in) {
67 unsigned i;
68 for (i = 0; i < 10; i += 2) {
69 output[0+i] = output[0+i] + in[0+i];
70 output[1+i] = output[1+i] + in[1+i];
71 }
72}
73
74/* Find the difference of two numbers: output = in - output
75 * (note the order of the arguments!). */
76static void fdifference(limb *output, const limb *in) {
77 unsigned i;
78 for (i = 0; i < 10; ++i) {
79 output[i] = in[i] - output[i];
80 }
81}
82
83/* Multiply a number by a scalar: output = in * scalar */
84static void fscalar_product(limb *output, const limb *in, const limb scalar) {
85 unsigned i;
86 for (i = 0; i < 10; ++i) {
87 output[i] = in[i] * scalar;
88 }
89}
90
91/* Multiply two numbers: output = in2 * in
92 *
93 * output must be distinct to both inputs. The inputs are reduced coefficient
94 * form, the output is not.
95 *
96 * output[x] <= 14 * the largest product of the input limbs. */
97static void fproduct(limb *output, const limb *in2, const limb *in) {
98 output[0] = ((limb) ((s32) in2[0])) * ((s32) in[0]);
99 output[1] = ((limb) ((s32) in2[0])) * ((s32) in[1]) +
100 ((limb) ((s32) in2[1])) * ((s32) in[0]);
101 output[2] = 2 * ((limb) ((s32) in2[1])) * ((s32) in[1]) +
102 ((limb) ((s32) in2[0])) * ((s32) in[2]) +
103 ((limb) ((s32) in2[2])) * ((s32) in[0]);
104 output[3] = ((limb) ((s32) in2[1])) * ((s32) in[2]) +
105 ((limb) ((s32) in2[2])) * ((s32) in[1]) +
106 ((limb) ((s32) in2[0])) * ((s32) in[3]) +
107 ((limb) ((s32) in2[3])) * ((s32) in[0]);
108 output[4] = ((limb) ((s32) in2[2])) * ((s32) in[2]) +
109 2 * (((limb) ((s32) in2[1])) * ((s32) in[3]) +
110 ((limb) ((s32) in2[3])) * ((s32) in[1])) +
111 ((limb) ((s32) in2[0])) * ((s32) in[4]) +
112 ((limb) ((s32) in2[4])) * ((s32) in[0]);
113 output[5] = ((limb) ((s32) in2[2])) * ((s32) in[3]) +
114 ((limb) ((s32) in2[3])) * ((s32) in[2]) +
115 ((limb) ((s32) in2[1])) * ((s32) in[4]) +
116 ((limb) ((s32) in2[4])) * ((s32) in[1]) +
117 ((limb) ((s32) in2[0])) * ((s32) in[5]) +
118 ((limb) ((s32) in2[5])) * ((s32) in[0]);
119 output[6] = 2 * (((limb) ((s32) in2[3])) * ((s32) in[3]) +
120 ((limb) ((s32) in2[1])) * ((s32) in[5]) +
121 ((limb) ((s32) in2[5])) * ((s32) in[1])) +
122 ((limb) ((s32) in2[2])) * ((s32) in[4]) +
123 ((limb) ((s32) in2[4])) * ((s32) in[2]) +
124 ((limb) ((s32) in2[0])) * ((s32) in[6]) +
125 ((limb) ((s32) in2[6])) * ((s32) in[0]);
126 output[7] = ((limb) ((s32) in2[3])) * ((s32) in[4]) +
127 ((limb) ((s32) in2[4])) * ((s32) in[3]) +
128 ((limb) ((s32) in2[2])) * ((s32) in[5]) +
129 ((limb) ((s32) in2[5])) * ((s32) in[2]) +
130 ((limb) ((s32) in2[1])) * ((s32) in[6]) +
131 ((limb) ((s32) in2[6])) * ((s32) in[1]) +
132 ((limb) ((s32) in2[0])) * ((s32) in[7]) +
133 ((limb) ((s32) in2[7])) * ((s32) in[0]);
134 output[8] = ((limb) ((s32) in2[4])) * ((s32) in[4]) +
135 2 * (((limb) ((s32) in2[3])) * ((s32) in[5]) +
136 ((limb) ((s32) in2[5])) * ((s32) in[3]) +
137 ((limb) ((s32) in2[1])) * ((s32) in[7]) +
138 ((limb) ((s32) in2[7])) * ((s32) in[1])) +
139 ((limb) ((s32) in2[2])) * ((s32) in[6]) +
140 ((limb) ((s32) in2[6])) * ((s32) in[2]) +
141 ((limb) ((s32) in2[0])) * ((s32) in[8]) +
142 ((limb) ((s32) in2[8])) * ((s32) in[0]);
143 output[9] = ((limb) ((s32) in2[4])) * ((s32) in[5]) +
144 ((limb) ((s32) in2[5])) * ((s32) in[4]) +
145 ((limb) ((s32) in2[3])) * ((s32) in[6]) +
146 ((limb) ((s32) in2[6])) * ((s32) in[3]) +
147 ((limb) ((s32) in2[2])) * ((s32) in[7]) +
148 ((limb) ((s32) in2[7])) * ((s32) in[2]) +
149 ((limb) ((s32) in2[1])) * ((s32) in[8]) +
150 ((limb) ((s32) in2[8])) * ((s32) in[1]) +
151 ((limb) ((s32) in2[0])) * ((s32) in[9]) +
152 ((limb) ((s32) in2[9])) * ((s32) in[0]);
153 output[10] = 2 * (((limb) ((s32) in2[5])) * ((s32) in[5]) +
154 ((limb) ((s32) in2[3])) * ((s32) in[7]) +
155 ((limb) ((s32) in2[7])) * ((s32) in[3]) +
156 ((limb) ((s32) in2[1])) * ((s32) in[9]) +
157 ((limb) ((s32) in2[9])) * ((s32) in[1])) +
158 ((limb) ((s32) in2[4])) * ((s32) in[6]) +
159 ((limb) ((s32) in2[6])) * ((s32) in[4]) +
160 ((limb) ((s32) in2[2])) * ((s32) in[8]) +
161 ((limb) ((s32) in2[8])) * ((s32) in[2]);
162 output[11] = ((limb) ((s32) in2[5])) * ((s32) in[6]) +
163 ((limb) ((s32) in2[6])) * ((s32) in[5]) +
164 ((limb) ((s32) in2[4])) * ((s32) in[7]) +
165 ((limb) ((s32) in2[7])) * ((s32) in[4]) +
166 ((limb) ((s32) in2[3])) * ((s32) in[8]) +
167 ((limb) ((s32) in2[8])) * ((s32) in[3]) +
168 ((limb) ((s32) in2[2])) * ((s32) in[9]) +
169 ((limb) ((s32) in2[9])) * ((s32) in[2]);
170 output[12] = ((limb) ((s32) in2[6])) * ((s32) in[6]) +
171 2 * (((limb) ((s32) in2[5])) * ((s32) in[7]) +
172 ((limb) ((s32) in2[7])) * ((s32) in[5]) +
173 ((limb) ((s32) in2[3])) * ((s32) in[9]) +
174 ((limb) ((s32) in2[9])) * ((s32) in[3])) +
175 ((limb) ((s32) in2[4])) * ((s32) in[8]) +
176 ((limb) ((s32) in2[8])) * ((s32) in[4]);
177 output[13] = ((limb) ((s32) in2[6])) * ((s32) in[7]) +
178 ((limb) ((s32) in2[7])) * ((s32) in[6]) +
179 ((limb) ((s32) in2[5])) * ((s32) in[8]) +
180 ((limb) ((s32) in2[8])) * ((s32) in[5]) +
181 ((limb) ((s32) in2[4])) * ((s32) in[9]) +
182 ((limb) ((s32) in2[9])) * ((s32) in[4]);
183 output[14] = 2 * (((limb) ((s32) in2[7])) * ((s32) in[7]) +
184 ((limb) ((s32) in2[5])) * ((s32) in[9]) +
185 ((limb) ((s32) in2[9])) * ((s32) in[5])) +
186 ((limb) ((s32) in2[6])) * ((s32) in[8]) +
187 ((limb) ((s32) in2[8])) * ((s32) in[6]);
188 output[15] = ((limb) ((s32) in2[7])) * ((s32) in[8]) +
189 ((limb) ((s32) in2[8])) * ((s32) in[7]) +
190 ((limb) ((s32) in2[6])) * ((s32) in[9]) +
191 ((limb) ((s32) in2[9])) * ((s32) in[6]);
192 output[16] = ((limb) ((s32) in2[8])) * ((s32) in[8]) +
193 2 * (((limb) ((s32) in2[7])) * ((s32) in[9]) +
194 ((limb) ((s32) in2[9])) * ((s32) in[7]));
195 output[17] = ((limb) ((s32) in2[8])) * ((s32) in[9]) +
196 ((limb) ((s32) in2[9])) * ((s32) in[8]);
197 output[18] = 2 * ((limb) ((s32) in2[9])) * ((s32) in[9]);
198}
199
200/* Reduce a long form to a short form by taking the input mod 2^255 - 19.
201 *
202 * On entry: |output[i]| < 14*2^54
203 * On exit: |output[0..8]| < 280*2^54 */
204static void freduce_degree(limb *output) {
205 /* Each of these shifts and adds ends up multiplying the value by 19.
206 *
207 * For output[0..8], the absolute entry value is < 14*2^54 and we add, at
208 * most, 19*14*2^54 thus, on exit, |output[0..8]| < 280*2^54. */
209 output[8] += output[18] << 4;
210 output[8] += output[18] << 1;
211 output[8] += output[18];
212 output[7] += output[17] << 4;
213 output[7] += output[17] << 1;
214 output[7] += output[17];
215 output[6] += output[16] << 4;
216 output[6] += output[16] << 1;
217 output[6] += output[16];
218 output[5] += output[15] << 4;
219 output[5] += output[15] << 1;
220 output[5] += output[15];
221 output[4] += output[14] << 4;
222 output[4] += output[14] << 1;
223 output[4] += output[14];
224 output[3] += output[13] << 4;
225 output[3] += output[13] << 1;
226 output[3] += output[13];
227 output[2] += output[12] << 4;
228 output[2] += output[12] << 1;
229 output[2] += output[12];
230 output[1] += output[11] << 4;
231 output[1] += output[11] << 1;
232 output[1] += output[11];
233 output[0] += output[10] << 4;
234 output[0] += output[10] << 1;
235 output[0] += output[10];
236}
237
238#if (-1 & 3) != 3
239#error "This code only works on a two's complement system"
240#endif
241
242/* return v / 2^26, using only shifts and adds.
243 *
244 * On entry: v can take any value. */
245static inline limb
246div_by_2_26(const limb v)
247{
248 /* High word of v; no shift needed. */
249 const uint32_t highword = (uint32_t) (((uint64_t) v) >> 32);
250 /* Set to all 1s if v was negative; else set to 0s. */
251 const int32_t sign = ((int32_t) highword) >> 31;
252 /* Set to 0x3ffffff if v was negative; else set to 0. */
253 const int32_t roundoff = ((uint32_t) sign) >> 6;
254 /* Should return v / (1<<26) */
255 return (v + roundoff) >> 26;
256}
257
258/* return v / (2^25), using only shifts and adds.
259 *
260 * On entry: v can take any value. */
261static inline limb
262div_by_2_25(const limb v)
263{
264 /* High word of v; no shift needed*/
265 const uint32_t highword = (uint32_t) (((uint64_t) v) >> 32);
266 /* Set to all 1s if v was negative; else set to 0s. */
267 const int32_t sign = ((int32_t) highword) >> 31;
268 /* Set to 0x1ffffff if v was negative; else set to 0. */
269 const int32_t roundoff = ((uint32_t) sign) >> 7;
270 /* Should return v / (1<<25) */
271 return (v + roundoff) >> 25;
272}
273
274#if 0
275/* return v / (2^25), using only shifts and adds.
276 *
277 * On entry: v can take any value. */
278static inline s32
279div_s32_by_2_25(const s32 v)
280{
281 const s32 roundoff = ((uint32_t)(v >> 31)) >> 7;
282 return (v + roundoff) >> 25;
283}
284#endif
285
286/* Reduce all coefficients of the short form input so that |x| < 2^26.
287 *
288 * On entry: |output[i]| < 280*2^54 */
289static void freduce_coefficients(limb *output) {
290 unsigned i;
291
292 output[10] = 0;
293
294 for (i = 0; i < 10; i += 2) {
295 limb over = div_by_2_26(output[i]);
296 /* The entry condition (that |output[i]| < 280*2^54) means that over is, at
297 * most, 280*2^28 in the first iteration of this loop. This is added to the
298 * next limb and we can approximate the resulting bound of that limb by
299 * 281*2^54. */
300 output[i] -= over << 26;
301 output[i+1] += over;
302
303 /* For the first iteration, |output[i+1]| < 281*2^54, thus |over| <
304 * 281*2^29. When this is added to the next limb, the resulting bound can
305 * be approximated as 281*2^54.
306 *
307 * For subsequent iterations of the loop, 281*2^54 remains a conservative
308 * bound and no overflow occurs. */
309 over = div_by_2_25(output[i+1]);
310 output[i+1] -= over << 25;
311 output[i+2] += over;
312 }
313 /* Now |output[10]| < 281*2^29 and all other coefficients are reduced. */
314 output[0] += output[10] << 4;
315 output[0] += output[10] << 1;
316 output[0] += output[10];
317
318 output[10] = 0;
319
320 /* Now output[1..9] are reduced, and |output[0]| < 2^26 + 19*281*2^29
321 * So |over| will be no more than 2^16. */
322 {
323 limb over = div_by_2_26(output[0]);
324 output[0] -= over << 26;
325 output[1] += over;
326 }
327
328 /* Now output[0,2..9] are reduced, and |output[1]| < 2^25 + 2^16 < 2^26. The
329 * bound on |output[1]| is sufficient to meet our needs. */
330}
331
332/* A helpful wrapper around fproduct: output = in * in2.
333 *
334 * On entry: |in[i]| < 2^27 and |in2[i]| < 2^27.
335 *
336 * output must be distinct to both inputs. The output is reduced degree
337 * (indeed, one need only provide storage for 10 limbs) and |output[i]| < 2^26. */
338static void
339fmul(limb *output, const limb *in, const limb *in2) {
340 limb t[19];
341 fproduct(t, in, in2);
342 /* |t[i]| < 14*2^54 */
343 freduce_degree(t);
344 freduce_coefficients(t);
345 /* |t[i]| < 2^26 */
346 memcpy(output, t, sizeof(limb) * 10);
347}
348
349/* Square a number: output = in**2
350 *
351 * output must be distinct from the input. The inputs are reduced coefficient
352 * form, the output is not.
353 *
354 * output[x] <= 14 * the largest product of the input limbs. */
355static void fsquare_inner(limb *output, const limb *in) {
356 output[0] = ((limb) ((s32) in[0])) * ((s32) in[0]);
357 output[1] = 2 * ((limb) ((s32) in[0])) * ((s32) in[1]);
358 output[2] = 2 * (((limb) ((s32) in[1])) * ((s32) in[1]) +
359 ((limb) ((s32) in[0])) * ((s32) in[2]));
360 output[3] = 2 * (((limb) ((s32) in[1])) * ((s32) in[2]) +
361 ((limb) ((s32) in[0])) * ((s32) in[3]));
362 output[4] = ((limb) ((s32) in[2])) * ((s32) in[2]) +
363 4 * ((limb) ((s32) in[1])) * ((s32) in[3]) +
364 2 * ((limb) ((s32) in[0])) * ((s32) in[4]);
365 output[5] = 2 * (((limb) ((s32) in[2])) * ((s32) in[3]) +
366 ((limb) ((s32) in[1])) * ((s32) in[4]) +
367 ((limb) ((s32) in[0])) * ((s32) in[5]));
368 output[6] = 2 * (((limb) ((s32) in[3])) * ((s32) in[3]) +
369 ((limb) ((s32) in[2])) * ((s32) in[4]) +
370 ((limb) ((s32) in[0])) * ((s32) in[6]) +
371 2 * ((limb) ((s32) in[1])) * ((s32) in[5]));
372 output[7] = 2 * (((limb) ((s32) in[3])) * ((s32) in[4]) +
373 ((limb) ((s32) in[2])) * ((s32) in[5]) +
374 ((limb) ((s32) in[1])) * ((s32) in[6]) +
375 ((limb) ((s32) in[0])) * ((s32) in[7]));
376 output[8] = ((limb) ((s32) in[4])) * ((s32) in[4]) +
377 2 * (((limb) ((s32) in[2])) * ((s32) in[6]) +
378 ((limb) ((s32) in[0])) * ((s32) in[8]) +
379 2 * (((limb) ((s32) in[1])) * ((s32) in[7]) +
380 ((limb) ((s32) in[3])) * ((s32) in[5])));
381 output[9] = 2 * (((limb) ((s32) in[4])) * ((s32) in[5]) +
382 ((limb) ((s32) in[3])) * ((s32) in[6]) +
383 ((limb) ((s32) in[2])) * ((s32) in[7]) +
384 ((limb) ((s32) in[1])) * ((s32) in[8]) +
385 ((limb) ((s32) in[0])) * ((s32) in[9]));
386 output[10] = 2 * (((limb) ((s32) in[5])) * ((s32) in[5]) +
387 ((limb) ((s32) in[4])) * ((s32) in[6]) +
388 ((limb) ((s32) in[2])) * ((s32) in[8]) +
389 2 * (((limb) ((s32) in[3])) * ((s32) in[7]) +
390 ((limb) ((s32) in[1])) * ((s32) in[9])));
391 output[11] = 2 * (((limb) ((s32) in[5])) * ((s32) in[6]) +
392 ((limb) ((s32) in[4])) * ((s32) in[7]) +
393 ((limb) ((s32) in[3])) * ((s32) in[8]) +
394 ((limb) ((s32) in[2])) * ((s32) in[9]));
395 output[12] = ((limb) ((s32) in[6])) * ((s32) in[6]) +
396 2 * (((limb) ((s32) in[4])) * ((s32) in[8]) +
397 2 * (((limb) ((s32) in[5])) * ((s32) in[7]) +
398 ((limb) ((s32) in[3])) * ((s32) in[9])));
399 output[13] = 2 * (((limb) ((s32) in[6])) * ((s32) in[7]) +
400 ((limb) ((s32) in[5])) * ((s32) in[8]) +
401 ((limb) ((s32) in[4])) * ((s32) in[9]));
402 output[14] = 2 * (((limb) ((s32) in[7])) * ((s32) in[7]) +
403 ((limb) ((s32) in[6])) * ((s32) in[8]) +
404 2 * ((limb) ((s32) in[5])) * ((s32) in[9]));
405 output[15] = 2 * (((limb) ((s32) in[7])) * ((s32) in[8]) +
406 ((limb) ((s32) in[6])) * ((s32) in[9]));
407 output[16] = ((limb) ((s32) in[8])) * ((s32) in[8]) +
408 4 * ((limb) ((s32) in[7])) * ((s32) in[9]);
409 output[17] = 2 * ((limb) ((s32) in[8])) * ((s32) in[9]);
410 output[18] = 2 * ((limb) ((s32) in[9])) * ((s32) in[9]);
411}
412
413/* fsquare sets output = in^2.
414 *
415 * On entry: The |in| argument is in reduced coefficients form and |in[i]| <
416 * 2^27.
417 *
418 * On exit: The |output| argument is in reduced coefficients form (indeed, one
419 * need only provide storage for 10 limbs) and |out[i]| < 2^26. */
420static void
421fsquare(limb *output, const limb *in) {
422 limb t[19];
423 fsquare_inner(t, in);
424 /* |t[i]| < 14*2^54 because the largest product of two limbs will be <
425 * 2^(27+27) and fsquare_inner adds together, at most, 14 of those
426 * products. */
427 freduce_degree(t);
428 freduce_coefficients(t);
429 /* |t[i]| < 2^26 */
430 memcpy(output, t, sizeof(limb) * 10);
431}
432
433/* Take a little-endian, 32-byte number and expand it into polynomial form */
434static void
435fexpand(limb *output, const u8 *input) {
436#define F(n,start,shift,mask) \
437 output[n] = ((((limb) input[start + 0]) | \
438 ((limb) input[start + 1]) << 8 | \
439 ((limb) input[start + 2]) << 16 | \
440 ((limb) input[start + 3]) << 24) >> shift) & mask;
441 F(0, 0, 0, 0x3ffffff);
442 F(1, 3, 2, 0x1ffffff);
443 F(2, 6, 3, 0x3ffffff);
444 F(3, 9, 5, 0x1ffffff);
445 F(4, 12, 6, 0x3ffffff);
446 F(5, 16, 0, 0x1ffffff);
447 F(6, 19, 1, 0x3ffffff);
448 F(7, 22, 3, 0x1ffffff);
449 F(8, 25, 4, 0x3ffffff);
450 F(9, 28, 6, 0x1ffffff);
451#undef F
452}
453
454#if (-32 >> 1) != -16
455#error "This code only works when >> does sign-extension on negative numbers"
456#endif
457
458/* s32_eq returns 0xffffffff iff a == b and zero otherwise. */
459static s32 s32_eq(s32 a, s32 b) {
460 a = ~(a ^ b);
461 a &= a << 16;
462 a &= a << 8;
463 a &= a << 4;
464 a &= a << 2;
465 a &= a << 1;
466 return a >> 31;
467}
468
469/* s32_gte returns 0xffffffff if a >= b and zero otherwise, where a and b are
470 * both non-negative. */
471static s32 s32_gte(s32 a, s32 b) {
472 a -= b;
473 /* a >= 0 iff a >= b. */
474 return ~(a >> 31);
475}
476
477/* Take a fully reduced polynomial form number and contract it into a
478 * little-endian, 32-byte array.
479 *
480 * On entry: |input_limbs[i]| < 2^26 */
481static void
482fcontract(u8 *output, limb *input_limbs) {
483 int i;
484 int j;
485 s32 input[10];
486
487 /* |input_limbs[i]| < 2^26, so it's valid to convert to an s32. */
488 for (i = 0; i < 10; i++) {
489 input[i] = (s32) input_limbs[i];
490 }
491
492 for (j = 0; j < 2; ++j) {
493 for (i = 0; i < 9; ++i) {
494 if ((i & 1) == 1) {
495 /* This calculation is a time-invariant way to make input[i]
496 * non-negative by borrowing from the next-larger limb. */
497 const s32 mask = input[i] >> 31;
498 const s32 carry = -((input[i] & mask) >> 25);
499 input[i] = input[i] + (carry << 25);
500 input[i+1] = input[i+1] - carry;
501 } else {
502 const s32 mask = input[i] >> 31;
503 const s32 carry = -((input[i] & mask) >> 26);
504 input[i] = input[i] + (carry << 26);
505 input[i+1] = input[i+1] - carry;
506 }
507 }
508
509 /* There's no greater limb for input[9] to borrow from, but we can multiply
510 * by 19 and borrow from input[0], which is valid mod 2^255-19. */
511 {
512 const s32 mask = input[9] >> 31;
513 const s32 carry = -((input[9] & mask) >> 25);
514 input[9] = input[9] + (carry << 25);
515 input[0] = input[0] - (carry * 19);
516 }
517
518 /* After the first iteration, input[1..9] are non-negative and fit within
519 * 25 or 26 bits, depending on position. However, input[0] may be
520 * negative. */
521 }
522
523 /* The first borrow-propagation pass above ended with every limb
524 except (possibly) input[0] non-negative.
525
526 If input[0] was negative after the first pass, then it was because of a
527 carry from input[9]. On entry, input[9] < 2^26 so the carry was, at most,
528 one, since (2**26-1) >> 25 = 1. Thus input[0] >= -19.
529
530 In the second pass, each limb is decreased by at most one. Thus the second
531 borrow-propagation pass could only have wrapped around to decrease
532 input[0] again if the first pass left input[0] negative *and* input[1]
533 through input[9] were all zero. In that case, input[1] is now 2^25 - 1,
534 and this last borrow-propagation step will leave input[1] non-negative. */
535 {
536 const s32 mask = input[0] >> 31;
537 const s32 carry = -((input[0] & mask) >> 26);
538 input[0] = input[0] + (carry << 26);
539 input[1] = input[1] - carry;
540 }
541
542 /* All input[i] are now non-negative. However, there might be values between
543 * 2^25 and 2^26 in a limb which is, nominally, 25 bits wide. */
544 for (j = 0; j < 2; j++) {
545 for (i = 0; i < 9; i++) {
546 if ((i & 1) == 1) {
547 const s32 carry = input[i] >> 25;
548 input[i] &= 0x1ffffff;
549 input[i+1] += carry;
550 } else {
551 const s32 carry = input[i] >> 26;
552 input[i] &= 0x3ffffff;
553 input[i+1] += carry;
554 }
555 }
556
557 {
558 const s32 carry = input[9] >> 25;
559 input[9] &= 0x1ffffff;
560 input[0] += 19*carry;
561 }
562 }
563
564 /* If the first carry-chain pass, just above, ended up with a carry from
565 * input[9], and that caused input[0] to be out-of-bounds, then input[0] was
566 * < 2^26 + 2*19, because the carry was, at most, two.
567 *
568 * If the second pass carried from input[9] again then input[0] is < 2*19 and
569 * the input[9] -> input[0] carry didn't push input[0] out of bounds. */
570
571 /* It still remains the case that input might be between 2^255-19 and 2^255.
572 * In this case, input[1..9] must take their maximum value and input[0] must
573 * be >= (2^255-19) & 0x3ffffff, which is 0x3ffffed. */
574 s32 mask = s32_gte(input[0], 0x3ffffed);
575 for (i = 1; i < 10; i++) {
576 if ((i & 1) == 1) {
577 mask &= s32_eq(input[i], 0x1ffffff);
578 } else {
579 mask &= s32_eq(input[i], 0x3ffffff);
580 }
581 }
582
583 /* mask is either 0xffffffff (if input >= 2^255-19) and zero otherwise. Thus
584 * this conditionally subtracts 2^255-19. */
585 input[0] -= mask & 0x3ffffed;
586
587 for (i = 1; i < 10; i++) {
588 if ((i & 1) == 1) {
589 input[i] -= mask & 0x1ffffff;
590 } else {
591 input[i] -= mask & 0x3ffffff;
592 }
593 }
594
595 input[1] <<= 2;
596 input[2] <<= 3;
597 input[3] <<= 5;
598 input[4] <<= 6;
599 input[6] <<= 1;
600 input[7] <<= 3;
601 input[8] <<= 4;
602 input[9] <<= 6;
603#define F(i, s) \
604 output[s+0] |= input[i] & 0xff; \
605 output[s+1] = (input[i] >> 8) & 0xff; \
606 output[s+2] = (input[i] >> 16) & 0xff; \
607 output[s+3] = (input[i] >> 24) & 0xff;
608 output[0] = 0;
609 output[16] = 0;
610 F(0,0);
611 F(1,3);
612 F(2,6);
613 F(3,9);
614 F(4,12);
615 F(5,16);
616 F(6,19);
617 F(7,22);
618 F(8,25);
619 F(9,28);
620#undef F
621}
622
623/* Input: Q, Q', Q-Q'
624 * Output: 2Q, Q+Q'
625 *
626 * x2 z3: long form
627 * x3 z3: long form
628 * x z: short form, destroyed
629 * xprime zprime: short form, destroyed
630 * qmqp: short form, preserved
631 *
632 * On entry and exit, the absolute value of the limbs of all inputs and outputs
633 * are < 2^26. */
634static void fmonty(limb *x2, limb *z2, /* output 2Q */
635 limb *x3, limb *z3, /* output Q + Q' */
636 limb *x, limb *z, /* input Q */
637 limb *xprime, limb *zprime, /* input Q' */
638 const limb *qmqp /* input Q - Q' */) {
639 limb origx[10], origxprime[10], zzz[19], xx[19], zz[19], xxprime[19],
640 zzprime[19], zzzprime[19], xxxprime[19];
641
642 memcpy(origx, x, 10 * sizeof(limb));
643 fsum(x, z);
644 /* |x[i]| < 2^27 */
645 fdifference(z, origx); /* does x - z */
646 /* |z[i]| < 2^27 */
647
648 memcpy(origxprime, xprime, sizeof(limb) * 10);
649 fsum(xprime, zprime);
650 /* |xprime[i]| < 2^27 */
651 fdifference(zprime, origxprime);
652 /* |zprime[i]| < 2^27 */
653 fproduct(xxprime, xprime, z);
654 /* |xxprime[i]| < 14*2^54: the largest product of two limbs will be <
655 * 2^(27+27) and fproduct adds together, at most, 14 of those products.
656 * (Approximating that to 2^58 doesn't work out.) */
657 fproduct(zzprime, x, zprime);
658 /* |zzprime[i]| < 14*2^54 */
659 freduce_degree(xxprime);
660 freduce_coefficients(xxprime);
661 /* |xxprime[i]| < 2^26 */
662 freduce_degree(zzprime);
663 freduce_coefficients(zzprime);
664 /* |zzprime[i]| < 2^26 */
665 memcpy(origxprime, xxprime, sizeof(limb) * 10);
666 fsum(xxprime, zzprime);
667 /* |xxprime[i]| < 2^27 */
668 fdifference(zzprime, origxprime);
669 /* |zzprime[i]| < 2^27 */
670 fsquare(xxxprime, xxprime);
671 /* |xxxprime[i]| < 2^26 */
672 fsquare(zzzprime, zzprime);
673 /* |zzzprime[i]| < 2^26 */
674 fproduct(zzprime, zzzprime, qmqp);
675 /* |zzprime[i]| < 14*2^52 */
676 freduce_degree(zzprime);
677 freduce_coefficients(zzprime);
678 /* |zzprime[i]| < 2^26 */
679 memcpy(x3, xxxprime, sizeof(limb) * 10);
680 memcpy(z3, zzprime, sizeof(limb) * 10);
681
682 fsquare(xx, x);
683 /* |xx[i]| < 2^26 */
684 fsquare(zz, z);
685 /* |zz[i]| < 2^26 */
686 fproduct(x2, xx, zz);
687 /* |x2[i]| < 14*2^52 */
688 freduce_degree(x2);
689 freduce_coefficients(x2);
690 /* |x2[i]| < 2^26 */
691 fdifference(zz, xx); // does zz = xx - zz
692 /* |zz[i]| < 2^27 */
693 memset(zzz + 10, 0, sizeof(limb) * 9);
694 fscalar_product(zzz, zz, 121665);
695 /* |zzz[i]| < 2^(27+17) */
696 /* No need to call freduce_degree here:
697 fscalar_product doesn't increase the degree of its input. */
698 freduce_coefficients(zzz);
699 /* |zzz[i]| < 2^26 */
700 fsum(zzz, xx);
701 /* |zzz[i]| < 2^27 */
702 fproduct(z2, zz, zzz);
703 /* |z2[i]| < 14*2^(26+27) */
704 freduce_degree(z2);
705 freduce_coefficients(z2);
706 /* |z2|i| < 2^26 */
707}
708
709/* Conditionally swap two reduced-form limb arrays if 'iswap' is 1, but leave
710 * them unchanged if 'iswap' is 0. Runs in data-invariant time to avoid
711 * side-channel attacks.
712 *
713 * NOTE that this function requires that 'iswap' be 1 or 0; other values give
714 * wrong results. Also, the two limb arrays must be in reduced-coefficient,
715 * reduced-degree form: the values in a[10..19] or b[10..19] aren't swapped,
716 * and all all values in a[0..9],b[0..9] must have magnitude less than
717 * INT32_MAX. */
718static void
719swap_conditional(limb a[19], limb b[19], limb iswap) {
720 unsigned i;
721 const s32 swap = (s32) -iswap;
722
723 for (i = 0; i < 10; ++i) {
724 const s32 x = swap & ( ((s32)a[i]) ^ ((s32)b[i]) );
725 a[i] = ((s32)a[i]) ^ x;
726 b[i] = ((s32)b[i]) ^ x;
727 }
728}
729
730/* Calculates nQ where Q is the x-coordinate of a point on the curve
731 *
732 * resultx/resultz: the x coordinate of the resulting curve point (short form)
733 * n: a little endian, 32-byte number
734 * q: a point of the curve (short form) */
735static void
736cmult(limb *resultx, limb *resultz, const u8 *n, const limb *q) {
737 limb a[19] = {0}, b[19] = {1}, c[19] = {1}, d[19] = {0};
738 limb *nqpqx = a, *nqpqz = b, *nqx = c, *nqz = d, *t;
739 limb e[19] = {0}, f[19] = {1}, g[19] = {0}, h[19] = {1};
740 limb *nqpqx2 = e, *nqpqz2 = f, *nqx2 = g, *nqz2 = h;
741
742 unsigned i, j;
743
744 memcpy(nqpqx, q, sizeof(limb) * 10);
745
746 for (i = 0; i < 32; ++i) {
747 u8 byte = n[31 - i];
748 for (j = 0; j < 8; ++j) {
749 const limb bit = byte >> 7;
750
751 swap_conditional(nqx, nqpqx, bit);
752 swap_conditional(nqz, nqpqz, bit);
753 fmonty(nqx2, nqz2,
754 nqpqx2, nqpqz2,
755 nqx, nqz,
756 nqpqx, nqpqz,
757 q);
758 swap_conditional(nqx2, nqpqx2, bit);
759 swap_conditional(nqz2, nqpqz2, bit);
760
761 t = nqx;
762 nqx = nqx2;
763 nqx2 = t;
764 t = nqz;
765 nqz = nqz2;
766 nqz2 = t;
767 t = nqpqx;
768 nqpqx = nqpqx2;
769 nqpqx2 = t;
770 t = nqpqz;
771 nqpqz = nqpqz2;
772 nqpqz2 = t;
773
774 byte <<= 1;
775 }
776 }
777
778 memcpy(resultx, nqx, sizeof(limb) * 10);
779 memcpy(resultz, nqz, sizeof(limb) * 10);
780}
781
782// -----------------------------------------------------------------------------
783// Shamelessly copied from djb's code
784// -----------------------------------------------------------------------------
785static void
786crecip(limb *out, const limb *z) {
787 limb z2[10];
788 limb z9[10];
789 limb z11[10];
790 limb z2_5_0[10];
791 limb z2_10_0[10];
792 limb z2_20_0[10];
793 limb z2_50_0[10];
794 limb z2_100_0[10];
795 limb t0[10];
796 limb t1[10];
797 int i;
798
799 /* 2 */ fsquare(z2,z);
800 /* 4 */ fsquare(t1,z2);
801 /* 8 */ fsquare(t0,t1);
802 /* 9 */ fmul(z9,t0,z);
803 /* 11 */ fmul(z11,z9,z2);
804 /* 22 */ fsquare(t0,z11);
805 /* 2^5 - 2^0 = 31 */ fmul(z2_5_0,t0,z9);
806
807 /* 2^6 - 2^1 */ fsquare(t0,z2_5_0);
808 /* 2^7 - 2^2 */ fsquare(t1,t0);
809 /* 2^8 - 2^3 */ fsquare(t0,t1);
810 /* 2^9 - 2^4 */ fsquare(t1,t0);
811 /* 2^10 - 2^5 */ fsquare(t0,t1);
812 /* 2^10 - 2^0 */ fmul(z2_10_0,t0,z2_5_0);
813
814 /* 2^11 - 2^1 */ fsquare(t0,z2_10_0);
815 /* 2^12 - 2^2 */ fsquare(t1,t0);
816 /* 2^20 - 2^10 */ for (i = 2;i < 10;i += 2) { fsquare(t0,t1); fsquare(t1,t0); }
817 /* 2^20 - 2^0 */ fmul(z2_20_0,t1,z2_10_0);
818
819 /* 2^21 - 2^1 */ fsquare(t0,z2_20_0);
820 /* 2^22 - 2^2 */ fsquare(t1,t0);
821 /* 2^40 - 2^20 */ for (i = 2;i < 20;i += 2) { fsquare(t0,t1); fsquare(t1,t0); }
822 /* 2^40 - 2^0 */ fmul(t0,t1,z2_20_0);
823
824 /* 2^41 - 2^1 */ fsquare(t1,t0);
825 /* 2^42 - 2^2 */ fsquare(t0,t1);
826 /* 2^50 - 2^10 */ for (i = 2;i < 10;i += 2) { fsquare(t1,t0); fsquare(t0,t1); }
827 /* 2^50 - 2^0 */ fmul(z2_50_0,t0,z2_10_0);
828
829 /* 2^51 - 2^1 */ fsquare(t0,z2_50_0);
830 /* 2^52 - 2^2 */ fsquare(t1,t0);
831 /* 2^100 - 2^50 */ for (i = 2;i < 50;i += 2) { fsquare(t0,t1); fsquare(t1,t0); }
832 /* 2^100 - 2^0 */ fmul(z2_100_0,t1,z2_50_0);
833
834 /* 2^101 - 2^1 */ fsquare(t1,z2_100_0);
835 /* 2^102 - 2^2 */ fsquare(t0,t1);
836 /* 2^200 - 2^100 */ for (i = 2;i < 100;i += 2) { fsquare(t1,t0); fsquare(t0,t1); }
837 /* 2^200 - 2^0 */ fmul(t1,t0,z2_100_0);
838
839 /* 2^201 - 2^1 */ fsquare(t0,t1);
840 /* 2^202 - 2^2 */ fsquare(t1,t0);
841 /* 2^250 - 2^50 */ for (i = 2;i < 50;i += 2) { fsquare(t0,t1); fsquare(t1,t0); }
842 /* 2^250 - 2^0 */ fmul(t0,t1,z2_50_0);
843
844 /* 2^251 - 2^1 */ fsquare(t1,t0);
845 /* 2^252 - 2^2 */ fsquare(t0,t1);
846 /* 2^253 - 2^3 */ fsquare(t1,t0);
847 /* 2^254 - 2^4 */ fsquare(t0,t1);
848 /* 2^255 - 2^5 */ fsquare(t1,t0);
849 /* 2^255 - 21 */ fmul(out,t1,z11);
850}
851
852int curve25519_donna(u8 *mypublic, const u8 *secret, const u8 *basepoint);
853
854int
855curve25519_donna(u8 *mypublic, const u8 *secret, const u8 *basepoint) {
856 limb bp[10], x[10], z[11], zmone[10];
857 uint8_t e[32];
858 int i;
859
860 for (i = 0; i < 32; ++i) e[i] = secret[i];
861 e[0] &= 248;
862 e[31] &= 127;
863 e[31] |= 64;
864
865 fexpand(bp, basepoint);
866 cmult(x, z, e, bp);
867 crecip(zmone, z);
868 fmul(z, x, zmone);
869 fcontract(mypublic, z);
870 return 0;
871}
torint.h
Integer definitions used throughout Tor.
Generated by doxygen 1.9.4