Tor  0.4.6.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 
53 typedef uint8_t u8;
54 typedef int32_t s32;
55 typedef 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 */
66 static 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!). */
76 static 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 */
84 static 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. */
97 static 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 */
204 static 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. */
245 static inline limb
246 div_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. */
261 static inline limb
262 div_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. */
278 static inline s32
279 div_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 */
289 static 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. */
338 static void
339 fmul(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. */
355 static 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. */
420 static void
421 fsquare(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 */
434 static void
435 fexpand(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. */
459 static 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. */
471 static 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 */
481 static void
482 fcontract(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. */
634 static 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. */
718 static void
719 swap_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) */
735 static void
736 cmult(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 // -----------------------------------------------------------------------------
785 static void
786 crecip(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 
852 int curve25519_donna(u8 *mypublic, const u8 *secret, const u8 *basepoint);
853 
854 int
855 curve25519_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.