2 * f25519.c: arithmetic modulo 2^255 - 19
5 * This file is Free Software. It has been modified to as part of its
6 * incorporation into secnet.
8 * Copyright 2017 Mark Wooding
10 * You may redistribute this file and/or modify it under the terms of
11 * the permissive licence shown below.
13 * You may redistribute secnet as a whole and/or modify it under the
14 * terms of the GNU General Public License as published by the Free
15 * Software Foundation; either version 3, or (at your option) any
18 * This program is distributed in the hope that it will be useful,
19 * but WITHOUT ANY WARRANTY; without even the implied warranty of
20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 * GNU General Public License for more details.
23 * You should have received a copy of the GNU General Public License
24 * along with this program; if not, see
25 * https://www.gnu.org/licenses/gpl.html.
28 * Imported from Catacomb, and modified for Secnet (2017-04-30):
30 * * Use `fake-mLib-bits.h' in place of the real <mLib/bits.h>.
32 * * Remove the 16/32-bit implementation, since C99 always has 64-bit
35 * * Remove the test rig code: a replacement is in a separate source file.
37 * * Disable some of the operations which aren't needed for X25519.
38 * (They're used for Ed25519, which we don't need.)
40 * The file's original comment headers are preserved below.
44 * Arithmetic modulo 2^255 - 19
46 * (c) 2017 Straylight/Edgeware
49 /*----- Licensing notice --------------------------------------------------*
51 * This file is part of Catacomb.
53 * Catacomb is free software; you can redistribute it and/or modify
54 * it under the terms of the GNU Library General Public License as
55 * published by the Free Software Foundation; either version 2 of the
56 * License, or (at your option) any later version.
58 * Catacomb is distributed in the hope that it will be useful,
59 * but WITHOUT ANY WARRANTY; without even the implied warranty of
60 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
61 * GNU Library General Public License for more details.
63 * You should have received a copy of the GNU Library General Public
64 * License along with Catacomb; if not, write to the Free
65 * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
69 /*----- Header files ------------------------------------------------------*/
73 /*----- Basic setup -------------------------------------------------------*/
75 /* Elements x of GF(2^255 - 19) are represented by ten signed integers x_i: x
76 * = SUM_{0<=i<10} x_i 2^ceil(51i/2), mostly following Bernstein's original
80 typedef int32 piece
; typedef int64 dblpiece
;
81 typedef uint32 upiece
; typedef uint64 udblpiece
;
83 #define PIECEWD(i) ((i)%2 ? 25 : 26)
86 #define M26 0x03ffffffu
87 #define M25 0x01ffffffu
88 #define B26 0x04000000u
89 #define B25 0x02000000u
90 #define B24 0x01000000u
92 #define PIECES(v) v##0, v##1, v##2, v##3, v##4, v##5, v##6, v##7, v##8, v##9
93 #define FETCH(v, w) do { \
94 v##0 = (w)->P[0]; v##1 = (w)->P[1]; \
95 v##2 = (w)->P[2]; v##3 = (w)->P[3]; \
96 v##4 = (w)->P[4]; v##5 = (w)->P[5]; \
97 v##6 = (w)->P[6]; v##7 = (w)->P[7]; \
98 v##8 = (w)->P[8]; v##9 = (w)->P[9]; \
100 #define STASH(w, v) do { \
101 (w)->P[0] = v##0; (w)->P[1] = v##1; \
102 (w)->P[2] = v##2; (w)->P[3] = v##3; \
103 (w)->P[4] = v##4; (w)->P[5] = v##5; \
104 (w)->P[6] = v##6; (w)->P[7] = v##7; \
105 (w)->P[8] = v##8; (w)->P[9] = v##9; \
108 /*----- Debugging machinery -----------------------------------------------*/
110 #if defined(F25519_DEBUG)
117 static mp
*get_2p255m91(void)
122 p
= mp_setbit(p
, MP_ZERO
, 255);
123 mp_build(&m19
, &w19
, &w19
+ 1);
124 p
= mp_sub(p
, p
, &m19
);
128 DEF_FDUMP(fdump
, piece
, PIECEWD
, NPIECE
, 32, get_2p255m91())
132 /*----- Loading and storing -----------------------------------------------*/
134 /* --- @f25519_load@ --- *
136 * Arguments: @f25519 *z@ = where to store the result
137 * @const octet xv[32]@ = source to read
141 * Use: Reads an element of %$\gf{2^{255} - 19}$% in external
142 * representation from @xv@ and stores it in @z@.
144 * External representation is little-endian base-256. Some
145 * elements have multiple encodings, which are not produced by
146 * correct software; use of noncanonical encodings is not an
147 * error, and toleration of them is considered a performance
151 void f25519_load(f25519
*z
, const octet xv
[32])
153 uint32 xw0
= LOAD32_L(xv
+ 0), xw1
= LOAD32_L(xv
+ 4),
154 xw2
= LOAD32_L(xv
+ 8), xw3
= LOAD32_L(xv
+ 12),
155 xw4
= LOAD32_L(xv
+ 16), xw5
= LOAD32_L(xv
+ 20),
156 xw6
= LOAD32_L(xv
+ 24), xw7
= LOAD32_L(xv
+ 28);
157 piece
PIECES(x
), b
, c
;
159 /* First, split the 32-bit words into the irregularly-sized pieces we need
160 * for the field representation. These pieces are all positive. We'll do
161 * the sign correction afterwards.
163 * It may be that the top bit of the input is set. Avoid trouble by
164 * folding that back round into the bottom piece of the representation.
166 * Here, we briefly have 0 <= x_0 < 2^26 + 19, but will resolve this later.
167 * Otherwise, we have 0 <= x_{2i} < 2^26, and 0 <= x_{2i+1} < 2^25.
169 x0
= ((xw0
<< 0)&0x03ffffff) + 19*((xw7
>> 31)&0x00000001);
170 x1
= ((xw1
<< 6)&0x01ffffc0) | ((xw0
>> 26)&0x0000003f);
171 x2
= ((xw2
<< 13)&0x03ffe000) | ((xw1
>> 19)&0x00001fff);
172 x3
= ((xw3
<< 19)&0x01f80000) | ((xw2
>> 13)&0x0007ffff);
173 x4
= ((xw3
>> 6)&0x03ffffff);
174 x5
= (xw4
<< 0)&0x01ffffff;
175 x6
= ((xw5
<< 7)&0x03ffff80) | ((xw4
>> 25)&0x0000007f);
176 x7
= ((xw6
<< 13)&0x01ffe000) | ((xw5
>> 19)&0x00001fff);
177 x8
= ((xw7
<< 20)&0x03f00000) | ((xw6
>> 12)&0x000fffff);
178 x9
= ((xw7
>> 6)&0x01ffffff);
180 /* Next, we convert these pieces into a roughly balanced signed
181 * representation. For each piece, check to see if its top bit is set. If
182 * it is, then lend a bit to the next piece over. For x_9, this needs to
183 * be carried around, which is a little fiddly. In particular, we delay
184 * the carry until after all of the pieces have been balanced. If we don't
185 * do this, then we have to do a more expensive test for nonzeroness to
186 * decide whether to lend a bit leftwards rather than just testing a single
189 * This fixes up the anomalous size of x_0: the loan of a bit becomes an
190 * actual carry if x_0 >= 2^26. By the end, then, we have:
196 * Note that we don't try for a canonical representation here: both upper
197 * and lower bounds are achievable.
199 * All of the x_i at this point are positive, so we don't need to do
200 * anything wierd when masking them.
202 b
= x9
&B24
; c
= 19&((b
>> 19) - (b
>> 24)); x9
-= b
<< 1;
203 b
= x8
&B25
; x9
+= b
>> 25; x8
-= b
<< 1;
204 b
= x7
&B24
; x8
+= b
>> 24; x7
-= b
<< 1;
205 b
= x6
&B25
; x7
+= b
>> 25; x6
-= b
<< 1;
206 b
= x5
&B24
; x6
+= b
>> 24; x5
-= b
<< 1;
207 b
= x4
&B25
; x5
+= b
>> 25; x4
-= b
<< 1;
208 b
= x3
&B24
; x4
+= b
>> 24; x3
-= b
<< 1;
209 b
= x2
&B25
; x3
+= b
>> 25; x2
-= b
<< 1;
210 b
= x1
&B24
; x2
+= b
>> 24; x1
-= b
<< 1;
211 b
= x0
&B25
; x1
+= (b
>> 25) + (x0
>> 26); x0
= (x0
&M26
) - (b
<< 1);
214 /* And with that, we're done. */
218 /* --- @f25519_store@ --- *
220 * Arguments: @octet zv[32]@ = where to write the result
221 * @const f25519 *x@ = the field element to write
225 * Use: Stores a field element in the given octet vector in external
226 * representation. A canonical encoding is always stored, so,
227 * in particular, the top bit of @xv[31]@ is always left clear.
230 void f25519_store(octet zv
[32], const f25519
*x
)
232 piece
PIECES(x
), PIECES(y
), c
, d
;
233 uint32 zw0
, zw1
, zw2
, zw3
, zw4
, zw5
, zw6
, zw7
;
238 /* First, propagate the carries throughout the pieces. By the end of this,
239 * we'll have all of the pieces canonically sized and positive, and maybe
240 * there'll be (signed) carry out. The carry c is in { -1, 0, +1 }, and
241 * the remaining value will be in the half-open interval [0, 2^255). The
242 * whole represented value is then x + 2^255 c.
244 * It's worth paying careful attention to the bounds. We assume that we
245 * start out with |x_i| <= 2^30. We start by cutting off and reducing the
246 * carry c_9 from the topmost piece, x_9. This leaves 0 <= x_9 < 2^25; and
247 * we'll have |c_9| <= 2^5. We multiply this by 19 and we'll add it onto
248 * x_0 and propagate the carries: but what bounds can we calculate on x
251 * Let o_i = floor(51 i/2). We have X_i = SUM_{0<=j<i} x_j 2^{o_i}, so
252 * x = X_10. We see, inductively, that |X_i| < 2^{31+o_{i-1}}: X_0 = 0;
253 * |x_i| <= 2^30; and |X_{i+1}| = |X_i + x_i 2^{o_i}| <= |X_i| + 2^{30+o_i}
254 * < 2^{31+o_i}. Then x = X_9 + 2^230 x_9, and we have better bounds for
257 * -2^235 < x + 19 c_9 < 2^255 + 2^235
259 * Here, the x_i are signed, so we must be cautious about bithacking them.
261 c
= ASR(piece
, x9
, 25); x9
= (upiece
)x9
&M25
;
262 x0
+= 19*c
; c
= ASR(piece
, x0
, 26); x0
= (upiece
)x0
&M26
;
263 x1
+= c
; c
= ASR(piece
, x1
, 25); x1
= (upiece
)x1
&M25
;
264 x2
+= c
; c
= ASR(piece
, x2
, 26); x2
= (upiece
)x2
&M26
;
265 x3
+= c
; c
= ASR(piece
, x3
, 25); x3
= (upiece
)x3
&M25
;
266 x4
+= c
; c
= ASR(piece
, x4
, 26); x4
= (upiece
)x4
&M26
;
267 x5
+= c
; c
= ASR(piece
, x5
, 25); x5
= (upiece
)x5
&M25
;
268 x6
+= c
; c
= ASR(piece
, x6
, 26); x6
= (upiece
)x6
&M26
;
269 x7
+= c
; c
= ASR(piece
, x7
, 25); x7
= (upiece
)x7
&M25
;
270 x8
+= c
; c
= ASR(piece
, x8
, 26); x8
= (upiece
)x8
&M26
;
271 x9
+= c
; c
= ASR(piece
, x9
, 25); x9
= (upiece
)x9
&M25
;
273 /* Now we have a slightly fiddly job to do. If c = +1, or if c = 0 and
274 * x >= 2^255 - 19, then we should subtract 2^255 - 19 from the whole
275 * value; if c = -1 then we should add 2^255 - 19; and otherwise we should
278 * But conditional behaviour is bad, m'kay. So here's what we do instead.
280 * The first job is to sort out what we wanted to do. If c = -1 then we
281 * want to (a) invert the constant addend and (b) feed in a carry-in;
282 * otherwise, we don't.
287 /* Now do the addition/subtraction. Remember that all of the x_i are
288 * nonnegative, so shifting and masking are safe and easy.
290 d
+= x0
+ (19 ^ (M26
&m
)); y0
= d
&M26
; d
>>= 26;
291 d
+= x1
+ (M25
&m
); y1
= d
&M25
; d
>>= 25;
292 d
+= x2
+ (M26
&m
); y2
= d
&M26
; d
>>= 26;
293 d
+= x3
+ (M25
&m
); y3
= d
&M25
; d
>>= 25;
294 d
+= x4
+ (M26
&m
); y4
= d
&M26
; d
>>= 26;
295 d
+= x5
+ (M25
&m
); y5
= d
&M25
; d
>>= 25;
296 d
+= x6
+ (M26
&m
); y6
= d
&M26
; d
>>= 26;
297 d
+= x7
+ (M25
&m
); y7
= d
&M25
; d
>>= 25;
298 d
+= x8
+ (M26
&m
); y8
= d
&M26
; d
>>= 26;
299 d
+= x9
+ (M25
&m
); y9
= d
&M25
; d
>>= 25;
301 /* The final carry-out is in d; since we only did addition, and the x_i are
302 * nonnegative, then d is in { 0, 1 }. We want to keep y, rather than x,
303 * if (a) c /= 0 (in which case we know that the old value was
304 * unsatisfactory), or (b) if d = 1 (in which case, if c = 0, we know that
305 * the subtraction didn't cause a borrow, so we must be in the case where
306 * 2^255 - 19 <= x < 2^255).
308 m
= NONZEROP(c
) | ~NONZEROP(d
- 1);
309 x0
= (y0
&m
) | (x0
&~m
); x1
= (y1
&m
) | (x1
&~m
);
310 x2
= (y2
&m
) | (x2
&~m
); x3
= (y3
&m
) | (x3
&~m
);
311 x4
= (y4
&m
) | (x4
&~m
); x5
= (y5
&m
) | (x5
&~m
);
312 x6
= (y6
&m
) | (x6
&~m
); x7
= (y7
&m
) | (x7
&~m
);
313 x8
= (y8
&m
) | (x8
&~m
); x9
= (y9
&m
) | (x9
&~m
);
315 /* Extract 32-bit words from the value. */
316 zw0
= ((x0
>> 0)&0x03ffffff) | (((uint32
)x1
<< 26)&0xfc000000);
317 zw1
= ((x1
>> 6)&0x0007ffff) | (((uint32
)x2
<< 19)&0xfff80000);
318 zw2
= ((x2
>> 13)&0x00001fff) | (((uint32
)x3
<< 13)&0xffffe000);
319 zw3
= ((x3
>> 19)&0x0000003f) | (((uint32
)x4
<< 6)&0xffffffc0);
320 zw4
= ((x5
>> 0)&0x01ffffff) | (((uint32
)x6
<< 25)&0xfe000000);
321 zw5
= ((x6
>> 7)&0x0007ffff) | (((uint32
)x7
<< 19)&0xfff80000);
322 zw6
= ((x7
>> 13)&0x00000fff) | (((uint32
)x8
<< 12)&0xfffff000);
323 zw7
= ((x8
>> 20)&0x0000003f) | (((uint32
)x9
<< 6)&0x7fffffc0);
325 /* Store the result as an octet string. */
326 STORE32_L(zv
+ 0, zw0
); STORE32_L(zv
+ 4, zw1
);
327 STORE32_L(zv
+ 8, zw2
); STORE32_L(zv
+ 12, zw3
);
328 STORE32_L(zv
+ 16, zw4
); STORE32_L(zv
+ 20, zw5
);
329 STORE32_L(zv
+ 24, zw6
); STORE32_L(zv
+ 28, zw7
);
332 /* --- @f25519_set@ --- *
334 * Arguments: @f25519 *z@ = where to write the result
335 * @int a@ = a small-ish constant
339 * Use: Sets @z@ to equal @a@.
342 void f25519_set(f25519
*x
, int a
)
347 for (i
= 1; i
< NPIECE
; i
++) x
->P
[i
] = 0;
350 /*----- Basic arithmetic --------------------------------------------------*/
352 /* --- @f25519_add@ --- *
354 * Arguments: @f25519 *z@ = where to put the result (may alias @x@ or @y@)
355 * @const f25519 *x, *y@ = two operands
359 * Use: Set @z@ to the sum %$x + y$%.
362 void f25519_add(f25519
*z
, const f25519
*x
, const f25519
*y
)
364 z
->P
[0] = x
->P
[0] + y
->P
[0]; z
->P
[1] = x
->P
[1] + y
->P
[1];
365 z
->P
[2] = x
->P
[2] + y
->P
[2]; z
->P
[3] = x
->P
[3] + y
->P
[3];
366 z
->P
[4] = x
->P
[4] + y
->P
[4]; z
->P
[5] = x
->P
[5] + y
->P
[5];
367 z
->P
[6] = x
->P
[6] + y
->P
[6]; z
->P
[7] = x
->P
[7] + y
->P
[7];
368 z
->P
[8] = x
->P
[8] + y
->P
[8]; z
->P
[9] = x
->P
[9] + y
->P
[9];
371 /* --- @f25519_sub@ --- *
373 * Arguments: @f25519 *z@ = where to put the result (may alias @x@ or @y@)
374 * @const f25519 *x, *y@ = two operands
378 * Use: Set @z@ to the difference %$x - y$%.
381 void f25519_sub(f25519
*z
, const f25519
*x
, const f25519
*y
)
383 z
->P
[0] = x
->P
[0] - y
->P
[0]; z
->P
[1] = x
->P
[1] - y
->P
[1];
384 z
->P
[2] = x
->P
[2] - y
->P
[2]; z
->P
[3] = x
->P
[3] - y
->P
[3];
385 z
->P
[4] = x
->P
[4] - y
->P
[4]; z
->P
[5] = x
->P
[5] - y
->P
[5];
386 z
->P
[6] = x
->P
[6] - y
->P
[6]; z
->P
[7] = x
->P
[7] - y
->P
[7];
387 z
->P
[8] = x
->P
[8] - y
->P
[8]; z
->P
[9] = x
->P
[9] - y
->P
[9];
390 #ifndef F25519_TRIM_X25519
392 /* --- @f25519_neg@ --- *
394 * Arguments: @f25519 *z@ = where to put the result (may alias @x@)
395 * @const f25519 *x@ = an operand
402 void f25519_neg(f25519
*z
, const f25519
*x
)
404 z
->P
[0] = -x
->P
[0]; z
->P
[1] = -x
->P
[1];
405 z
->P
[2] = -x
->P
[2]; z
->P
[3] = -x
->P
[3];
406 z
->P
[4] = -x
->P
[4]; z
->P
[5] = -x
->P
[5];
407 z
->P
[6] = -x
->P
[6]; z
->P
[7] = -x
->P
[7];
408 z
->P
[8] = -x
->P
[8]; z
->P
[9] = -x
->P
[9];
413 /*----- Constant-time utilities -------------------------------------------*/
415 #ifndef F25519_TRIM_X25519
417 /* --- @f25519_pick2@ --- *
419 * Arguments: @f25519 *z@ = where to put the result (may alias @x@ or @y@)
420 * @const f25519 *x, *y@ = two operands
421 * @uint32 m@ = a mask
425 * Use: If @m@ is zero, set @z = y@; if @m@ is all-bits-set, then set
426 * @z = x@. If @m@ has some other value, then scramble @z@ in
430 void f25519_pick2(f25519
*z
, const f25519
*x
, const f25519
*y
, uint32 m
)
432 mask32 mm
= FIX_MASK32(m
);
434 z
->P
[0] = PICK2(x
->P
[0], y
->P
[0], mm
);
435 z
->P
[1] = PICK2(x
->P
[1], y
->P
[1], mm
);
436 z
->P
[2] = PICK2(x
->P
[2], y
->P
[2], mm
);
437 z
->P
[3] = PICK2(x
->P
[3], y
->P
[3], mm
);
438 z
->P
[4] = PICK2(x
->P
[4], y
->P
[4], mm
);
439 z
->P
[5] = PICK2(x
->P
[5], y
->P
[5], mm
);
440 z
->P
[6] = PICK2(x
->P
[6], y
->P
[6], mm
);
441 z
->P
[7] = PICK2(x
->P
[7], y
->P
[7], mm
);
442 z
->P
[8] = PICK2(x
->P
[8], y
->P
[8], mm
);
443 z
->P
[9] = PICK2(x
->P
[9], y
->P
[9], mm
);
446 /* --- @f25519_pickn@ --- *
448 * Arguments: @f25519 *z@ = where to put the result
449 * @const f25519 *v@ = a table of entries
450 * @size_t n@ = the number of entries in @v@
451 * @size_t i@ = an index
455 * Use: If @0 <= i < n < 32@ then set @z = v[i]@. If @n >= 32@ then
456 * do something unhelpful; otherwise, if @i >= n@ then set @z@
460 void f25519_pickn(f25519
*z
, const f25519
*v
, size_t n
, size_t i
)
462 uint32 b
= (uint32
)1 << (31 - i
);
465 z
->P
[0] = z
->P
[1] = z
->P
[2] = z
->P
[3] = z
->P
[4] =
466 z
->P
[5] = z
->P
[6] = z
->P
[7] = z
->P
[8] = z
->P
[9] = 0;
469 CONDPICK(z
->P
[0], v
->P
[0], m
);
470 CONDPICK(z
->P
[1], v
->P
[1], m
);
471 CONDPICK(z
->P
[2], v
->P
[2], m
);
472 CONDPICK(z
->P
[3], v
->P
[3], m
);
473 CONDPICK(z
->P
[4], v
->P
[4], m
);
474 CONDPICK(z
->P
[5], v
->P
[5], m
);
475 CONDPICK(z
->P
[6], v
->P
[6], m
);
476 CONDPICK(z
->P
[7], v
->P
[7], m
);
477 CONDPICK(z
->P
[8], v
->P
[8], m
);
478 CONDPICK(z
->P
[9], v
->P
[9], m
);
485 /* --- @f25519_condswap@ --- *
487 * Arguments: @f25519 *x, *y@ = two operands
488 * @uint32 m@ = a mask
492 * Use: If @m@ is zero, do nothing; if @m@ is all-bits-set, then
493 * exchange @x@ and @y@. If @m@ has some other value, then
494 * scramble @x@ and @y@ in an unhelpful way.
497 void f25519_condswap(f25519
*x
, f25519
*y
, uint32 m
)
499 mask32 mm
= FIX_MASK32(m
);
501 CONDSWAP(x
->P
[0], y
->P
[0], mm
);
502 CONDSWAP(x
->P
[1], y
->P
[1], mm
);
503 CONDSWAP(x
->P
[2], y
->P
[2], mm
);
504 CONDSWAP(x
->P
[3], y
->P
[3], mm
);
505 CONDSWAP(x
->P
[4], y
->P
[4], mm
);
506 CONDSWAP(x
->P
[5], y
->P
[5], mm
);
507 CONDSWAP(x
->P
[6], y
->P
[6], mm
);
508 CONDSWAP(x
->P
[7], y
->P
[7], mm
);
509 CONDSWAP(x
->P
[8], y
->P
[8], mm
);
510 CONDSWAP(x
->P
[9], y
->P
[9], mm
);
513 #ifndef F25519_TRIM_X25519
515 /* --- @f25519_condneg@ --- *
517 * Arguments: @f25519 *z@ = where to put the result (may alias @x@)
518 * @const f25519 *x@ = an operand
519 * @uint32 m@ = a mask
523 * Use: If @m@ is zero, set @z = x@; if @m@ is all-bits-set, then set
524 * @z = -x@. If @m@ has some other value then scramble @z@ in
528 void f25519_condneg(f25519
*z
, const f25519
*x
, uint32 m
)
530 mask32 m_xor
= FIX_MASK32(m
);
532 # define CONDNEG(x) (((x) ^ m_xor) + m_add)
534 z
->P
[0] = CONDNEG(x
->P
[0]);
535 z
->P
[1] = CONDNEG(x
->P
[1]);
536 z
->P
[2] = CONDNEG(x
->P
[2]);
537 z
->P
[3] = CONDNEG(x
->P
[3]);
538 z
->P
[4] = CONDNEG(x
->P
[4]);
539 z
->P
[5] = CONDNEG(x
->P
[5]);
540 z
->P
[6] = CONDNEG(x
->P
[6]);
541 z
->P
[7] = CONDNEG(x
->P
[7]);
542 z
->P
[8] = CONDNEG(x
->P
[8]);
543 z
->P
[9] = CONDNEG(x
->P
[9]);
550 /*----- Multiplication ----------------------------------------------------*/
552 /* Let B = 2^63 - 1 be the largest value such that +B and -B can be
553 * represented in a double-precision piece. On entry, it must be the case
554 * that |X_i| <= M <= B - 2^25 for some M. If this is the case, then, on
555 * exit, we will have |Z_i| <= 2^25 + 19 M/2^25.
557 #define CARRYSTEP(z, x, m, b, f, xx, n) do { \
558 (z) = (dblpiece)((udblpiece)(x)&(m)) - (b) + \
559 (f)*ASR(dblpiece, (xx), (n)); \
561 #define CARRY_REDUCE(z, x) do { \
562 dblpiece PIECES(_t); \
564 /* Bias the input pieces. This keeps the carries and so on centred \
565 * around zero rather than biased positive. \
567 _t0 = (x##0) + B25; _t1 = (x##1) + B24; \
568 _t2 = (x##2) + B25; _t3 = (x##3) + B24; \
569 _t4 = (x##4) + B25; _t5 = (x##5) + B24; \
570 _t6 = (x##6) + B25; _t7 = (x##7) + B24; \
571 _t8 = (x##8) + B25; _t9 = (x##9) + B24; \
573 /* Calculate the reduced pieces. Careful with the bithacking. */ \
574 CARRYSTEP(z##0, _t0, M26, B25, 19, _t9, 25); \
575 CARRYSTEP(z##1, _t1, M25, B24, 1, _t0, 26); \
576 CARRYSTEP(z##2, _t2, M26, B25, 1, _t1, 25); \
577 CARRYSTEP(z##3, _t3, M25, B24, 1, _t2, 26); \
578 CARRYSTEP(z##4, _t4, M26, B25, 1, _t3, 25); \
579 CARRYSTEP(z##5, _t5, M25, B24, 1, _t4, 26); \
580 CARRYSTEP(z##6, _t6, M26, B25, 1, _t5, 25); \
581 CARRYSTEP(z##7, _t7, M25, B24, 1, _t6, 26); \
582 CARRYSTEP(z##8, _t8, M26, B25, 1, _t7, 25); \
583 CARRYSTEP(z##9, _t9, M25, B24, 1, _t8, 26); \
586 /* --- @f25519_mulconst@ --- *
588 * Arguments: @f25519 *z@ = where to put the result (may alias @x@)
589 * @const f25519 *x@ = an operand
590 * @long a@ = a small-ish constant; %$|a| < 2^{20}$%.
594 * Use: Set @z@ to the product %$a x$%.
597 void f25519_mulconst(f25519
*z
, const f25519
*x
, long a
)
600 dblpiece
PIECES(z
), aa
= a
;
604 /* Suppose that |x_i| <= 2^27, and |a| <= 2^23. Then we'll have
607 z0
= aa
*x0
; z1
= aa
*x1
; z2
= aa
*x2
; z3
= aa
*x3
; z4
= aa
*x4
;
608 z5
= aa
*x5
; z6
= aa
*x6
; z7
= aa
*x7
; z8
= aa
*x8
; z9
= aa
*x9
;
610 /* Following `CARRY_REDUCE', we'll have |z_i| <= 2^26. */
615 /* --- @f25519_mul@ --- *
617 * Arguments: @f25519 *z@ = where to put the result (may alias @x@ or @y@)
618 * @const f25519 *x, *y@ = two operands
622 * Use: Set @z@ to the product %$x y$%.
625 void f25519_mul(f25519
*z
, const f25519
*x
, const f25519
*y
)
627 piece
PIECES(x
), PIECES(y
);
631 FETCH(x
, x
); FETCH(y
, y
);
633 /* Suppose that |x_i|, |y_i| <= 2^27. Then we'll have
646 * all of which are less than 2^63 - 2^25.
649 #define M(a, b) ((dblpiece)(a)*(b))
651 19*(M(x2
, y8
) + M(x4
, y6
) + M(x6
, y4
) + M(x8
, y2
)) +
652 38*(M(x1
, y9
) + M(x3
, y7
) + M(x5
, y5
) + M(x7
, y3
) + M(x9
, y1
));
653 z1
= M(x0
, y1
) + M(x1
, y0
) +
654 19*(M(x2
, y9
) + M(x3
, y8
) + M(x4
, y7
) + M(x5
, y6
) +
655 M(x6
, y5
) + M(x7
, y4
) + M(x8
, y3
) + M(x9
, y2
));
656 z2
= M(x0
, y2
) + M(x2
, y0
) +
658 19*(M(x4
, y8
) + M(x6
, y6
) + M(x8
, y4
)) +
659 38*(M(x3
, y9
) + M(x5
, y7
) + M(x7
, y5
) + M(x9
, y3
));
660 z3
= M(x0
, y3
) + M(x1
, y2
) + M(x2
, y1
) + M(x3
, y0
) +
661 19*(M(x4
, y9
) + M(x5
, y8
) + M(x6
, y7
) +
662 M(x7
, y6
) + M(x8
, y5
) + M(x9
, y4
));
663 z4
= M(x0
, y4
) + M(x2
, y2
) + M(x4
, y0
) +
664 2*(M(x1
, y3
) + M(x3
, y1
)) +
665 19*(M(x6
, y8
) + M(x8
, y6
)) +
666 38*(M(x5
, y9
) + M(x7
, y7
) + M(x9
, y5
));
667 z5
= M(x0
, y5
) + M(x1
, y4
) + M(x2
, y3
) +
668 M(x3
, y2
) + M(x4
, y1
) + M(x5
, y0
) +
669 19*(M(x6
, y9
) + M(x7
, y8
) + M(x8
, y7
) + M(x9
, y6
));
670 z6
= M(x0
, y6
) + M(x2
, y4
) + M(x4
, y2
) + M(x6
, y0
) +
671 2*(M(x1
, y5
) + M(x3
, y3
) + M(x5
, y1
)) +
673 38*(M(x7
, y9
) + M(x9
, y7
));
674 z7
= M(x0
, y7
) + M(x1
, y6
) + M(x2
, y5
) + M(x3
, y4
) +
675 M(x4
, y3
) + M(x5
, y2
) + M(x6
, y1
) + M(x7
, y0
) +
676 19*(M(x8
, y9
) + M(x9
, y8
));
677 z8
= M(x0
, y8
) + M(x2
, y6
) + M(x4
, y4
) + M(x6
, y2
) + M(x8
, y0
) +
678 2*(M(x1
, y7
) + M(x3
, y5
) + M(x5
, y3
) + M(x7
, y1
)) +
680 z9
= M(x0
, y9
) + M(x1
, y8
) + M(x2
, y7
) + M(x3
, y6
) + M(x4
, y5
) +
681 M(x5
, y4
) + M(x6
, y3
) + M(x7
, y2
) + M(x8
, y1
) + M(x9
, y0
);
684 /* From above, we have |z_i| <= 2^63 - 2^25. A pass of `CARRY_REDUCE' will
685 * leave |z_i| <= 2^38 + 2^25; and a second pass will leave |z_i| <= 2^25 +
686 * 2^13, which is comfortable for an addition prior to the next
689 for (i
= 0; i
< 2; i
++) CARRY_REDUCE(z
, z
);
693 /* --- @f25519_sqr@ --- *
695 * Arguments: @f25519 *z@ = where to put the result (may alias @x@ or @y@)
696 * @const f25519 *x@ = an operand
700 * Use: Set @z@ to the square %$x^2$%.
703 void f25519_sqr(f25519
*z
, const f25519
*x
)
711 /* See `f25519_mul' for bounds. */
713 #define M(a, b) ((dblpiece)(a)*(b))
715 38*(M(x2
, x8
) + M(x4
, x6
) + M(x5
, x5
)) +
716 76*(M(x1
, x9
) + M(x3
, x7
));
718 38*(M(x2
, x9
) + M(x3
, x8
) + M(x4
, x7
) + M(x5
, x6
));
719 z2
= 2*(M(x0
, x2
) + M(x1
, x1
)) +
722 76*(M(x3
, x9
) + M(x5
, x7
));
723 z3
= 2*(M(x0
, x3
) + M(x1
, x2
)) +
724 38*(M(x4
, x9
) + M(x5
, x8
) + M(x6
, x7
));
728 38*(M(x6
, x8
) + M(x7
, x7
)) +
730 z5
= 2*(M(x0
, x5
) + M(x1
, x4
) + M(x2
, x3
)) +
731 38*(M(x6
, x9
) + M(x7
, x8
));
732 z6
= 2*(M(x0
, x6
) + M(x2
, x4
) + M(x3
, x3
)) +
736 z7
= 2*(M(x0
, x7
) + M(x1
, x6
) + M(x2
, x5
) + M(x3
, x4
)) +
739 2*(M(x0
, x8
) + M(x2
, x6
)) +
740 4*(M(x1
, x7
) + M(x3
, x5
)) +
742 z9
= 2*(M(x0
, x9
) + M(x1
, x8
) + M(x2
, x7
) + M(x3
, x6
) + M(x4
, x5
));
745 /* See `f25519_mul' for details. */
746 for (i
= 0; i
< 2; i
++) CARRY_REDUCE(z
, z
);
750 /*----- More complicated things -------------------------------------------*/
752 /* --- @f25519_inv@ --- *
754 * Arguments: @f25519 *z@ = where to put the result (may alias @x@)
755 * @const f25519 *x@ = an operand
759 * Use: Stores in @z@ the multiplicative inverse %$x^{-1}$%. If
760 * %$x = 0$% then @z@ is set to zero. This is considered a
764 void f25519_inv(f25519
*z
, const f25519
*x
)
766 f25519 t
, u
, t2
, t11
, t2p10m1
, t2p50m1
;
769 #define SQRN(z, x, n) do { \
770 f25519_sqr((z), (x)); \
771 for (i = 1; i < (n); i++) f25519_sqr((z), (z)); \
774 /* Calculate x^-1 = x^(p - 2) = x^(2^255 - 21), which also handles x = 0 as
775 * intended. The addition chain here is from Bernstein's implementation; I
776 * couldn't find a better one.
777 */ /* step | value */
778 f25519_sqr(&t2
, x
); /* 1 | 2 */
779 SQRN(&u
, &t2
, 2); /* 3 | 8 */
780 f25519_mul(&t
, &u
, x
); /* 4 | 9 */
781 f25519_mul(&t11
, &t
, &t2
); /* 5 | 11 = 2^5 - 21 */
782 f25519_sqr(&u
, &t11
); /* 6 | 22 */
783 f25519_mul(&t
, &t
, &u
); /* 7 | 31 = 2^5 - 1 */
784 SQRN(&u
, &t
, 5); /* 12 | 2^10 - 2^5 */
785 f25519_mul(&t2p10m1
, &t
, &u
); /* 13 | 2^10 - 1 */
786 SQRN(&u
, &t2p10m1
, 10); /* 23 | 2^20 - 2^10 */
787 f25519_mul(&t
, &t2p10m1
, &u
); /* 24 | 2^20 - 1 */
788 SQRN(&u
, &t
, 20); /* 44 | 2^40 - 2^20 */
789 f25519_mul(&t
, &t
, &u
); /* 45 | 2^40 - 1 */
790 SQRN(&u
, &t
, 10); /* 55 | 2^50 - 2^10 */
791 f25519_mul(&t2p50m1
, &t2p10m1
, &u
); /* 56 | 2^50 - 1 */
792 SQRN(&u
, &t2p50m1
, 50); /* 106 | 2^100 - 2^50 */
793 f25519_mul(&t
, &t2p50m1
, &u
); /* 107 | 2^100 - 1 */
794 SQRN(&u
, &t
, 100); /* 207 | 2^200 - 2^100 */
795 f25519_mul(&t
, &t
, &u
); /* 208 | 2^200 - 1 */
796 SQRN(&u
, &t
, 50); /* 258 | 2^250 - 2^50 */
797 f25519_mul(&t
, &t2p50m1
, &u
); /* 259 | 2^250 - 1 */
798 SQRN(&u
, &t
, 5); /* 264 | 2^255 - 2^5 */
799 f25519_mul(z
, &u
, &t11
); /* 265 | 2^255 - 21 */
804 #ifndef F25519_TRIM_X25519
806 /* --- @f25519_quosqrt@ --- *
808 * Arguments: @f25519 *z@ = where to put the result (may alias @x@ or @y@)
809 * @const f25519 *x, *y@ = two operands
811 * Returns: Zero if successful, @-1@ if %$x/y$% is not a square.
813 * Use: Stores in @z@ the one of the square roots %$\pm\sqrt{x/y}$%.
814 * If %$x = y = 0% then the result is zero; if %$y = 0$% but %$x
815 * \ne 0$% then the operation fails. If you wanted a specific
816 * square root then you'll have to pick it yourself.
819 static const piece sqrtm1_pieces
[NPIECE
] = {
820 #if F25519_IMPL == 26
821 -32595792, -7943725, 9377950, 3500415, 12389472,
822 -272473, -25146209, -2005654, 326686, 11406482
823 #elif F25519_IMPL == 10
824 176, -88, 161, 157, -485, -196, -231, -220, -416,
825 -169, -255, 50, 189, -89, -266, -32, 202, -511,
826 423, 357, 248, -249, 80, 288, 50, 174
829 #define SQRTM1 ((const f25519 *)sqrtm1_pieces)
831 int f25519_quosqrt(f25519
*z
, const f25519
*x
, const f25519
*y
)
833 f25519 t
, u
, w
, beta
, xy3
, t2p50m1
;
834 octet xb
[32], b0
[32], b1
[32];
839 #define SQRN(z, x, n) do { \
840 f25519_sqr((z), (x)); \
841 for (i = 1; i < (n); i++) f25519_sqr((z), (z)); \
844 /* This is a bit tricky; the algorithm is from Bernstein, Duif, Lange,
845 * Schwabe, and Yang, `High-speed high-security signatures', 2011-09-26,
846 * https://ed25519.cr.yp.to/ed25519-20110926.pdf.
848 * First of all, a complicated exponentation. The addition chain here is
849 * mine. We start with some preliminary values.
850 */ /* step | value */
851 SQRN(&u
, y
, 1); /* 1 | 0, 2 */
852 f25519_mul(&t
, &u
, y
); /* 2 | 0, 3 */
853 f25519_mul(&xy3
, &t
, x
); /* 3 | 1, 3 */
854 SQRN(&u
, &u
, 1); /* 4 | 0, 4 */
855 f25519_mul(&w
, &u
, &xy3
); /* 5 | 1, 7 */
857 /* And now we calculate w^((p - 5)/8) = w^(252 - 3). */
858 SQRN(&u
, &w
, 1); /* 6 | 2 */
859 f25519_mul(&t
, &w
, &u
); /* 7 | 3 */
860 SQRN(&u
, &t
, 1); /* 8 | 6 */
861 f25519_mul(&t
, &u
, &w
); /* 9 | 7 */
862 SQRN(&u
, &t
, 3); /* 12 | 56 */
863 f25519_mul(&t
, &t
, &u
); /* 13 | 63 = 2^6 - 1 */
864 SQRN(&u
, &t
, 6); /* 19 | 2^12 - 2^6 */
865 f25519_mul(&t
, &t
, &u
); /* 20 | 2^12 - 1 */
866 SQRN(&u
, &t
, 12); /* 32 | 2^24 - 2^12 */
867 f25519_mul(&t
, &t
, &u
); /* 33 | 2^24 - 1 */
868 SQRN(&u
, &t
, 1); /* 34 | 2^25 - 2 */
869 f25519_mul(&t
, &u
, &w
); /* 35 | 2^25 - 1 */
870 SQRN(&u
, &t
, 25); /* 60 | 2^50 - 2^25 */
871 f25519_mul(&t2p50m1
, &t
, &u
); /* 61 | 2^50 - 1 */
872 SQRN(&u
, &t2p50m1
, 50); /* 111 | 2^100 - 2^50 */
873 f25519_mul(&t
, &t2p50m1
, &u
); /* 112 | 2^100 - 1 */
874 SQRN(&u
, &t
, 100); /* 212 | 2^200 - 2^100 */
875 f25519_mul(&t
, &t
, &u
); /* 213 | 2^200 - 1 */
876 SQRN(&u
, &t
, 50); /* 263 | 2^250 - 2^50 */
877 f25519_mul(&t
, &t2p50m1
, &u
); /* 264 | 2^250 - 1 */
878 SQRN(&u
, &t
, 2); /* 266 | 2^252 - 4 */
879 f25519_mul(&t
, &u
, &w
); /* 267 | 2^252 - 3 */
882 f25519_mul(&beta
, &t
, &xy3
); /* 268 | ... */
884 /* Now we have beta = (x y^3) (x y^7)^((p - 5)/8) = (x/y)^((p + 3)/8), and
885 * we're ready to finish the computation. Suppose that alpha^2 = u/w.
886 * Then beta^4 = (x/y)^((p + 3)/2) = alpha^(p + 3) = alpha^4 = (x/y)^2, so
887 * we have beta^2 = ±x/y. If y beta^2 = x then beta is the one we wanted;
888 * if -y beta^2 = x, then we want beta sqrt(-1), which we already know. Of
889 * course, it might not match either, in which case we fail.
891 * The easiest way to compare is to encode. This isn't as wasteful as it
892 * sounds: the hard part is normalizing the representations, which we have
895 f25519_sqr(&t
, &beta
);
896 f25519_mul(&t
, &t
, y
);
899 f25519_store(b0
, &t
);
900 f25519_store(b1
, &u
);
901 f25519_mul(&u
, &beta
, SQRTM1
);
903 m
= -ct_memeq(b0
, xb
, 32);
904 rc
= PICK2(0, rc
, m
);
905 f25519_pick2(z
, &beta
, &u
, m
);
906 m
= -ct_memeq(b1
, xb
, 32);
907 rc
= PICK2(0, rc
, m
);
909 /* And we're done. */
915 /*----- That's all, folks -------------------------------------------------*/