a0f40fa0b8c2e0c983383d3c9542a43d87d3a06d
3 * Arithmetic modulo 2^255 - 19
5 * (c) 2017 Straylight/Edgeware
8 /*----- Licensing notice --------------------------------------------------*
10 * This file is part of Catacomb.
12 * Catacomb is free software; you can redistribute it and/or modify
13 * it under the terms of the GNU Library General Public License as
14 * published by the Free Software Foundation; either version 2 of the
15 * License, or (at your option) any later version.
17 * Catacomb is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU Library General Public License for more details.
22 * You should have received a copy of the GNU Library General Public
23 * License along with Catacomb; if not, write to the Free
24 * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
28 /*----- Header files ------------------------------------------------------*/
35 /*----- Basic setup -------------------------------------------------------*/
37 typedef f25519_piece piece
;
40 /* Elements x of GF(2^255 - 19) are represented by ten signed integers x_i: x
41 * = SUM_{0<=i<10} x_i 2^ceil(51i/2), mostly following Bernstein's original
45 typedef int64 dblpiece
;
46 typedef uint32 upiece
; typedef uint64 udblpiece
;
48 #define PIECEWD(i) ((i)%2 ? 25 : 26)
51 #define M26 0x03ffffffu
52 #define M25 0x01ffffffu
53 #define B25 0x02000000u
54 #define B24 0x01000000u
56 #define PIECES(v) v##0, v##1, v##2, v##3, v##4, v##5, v##6, v##7, v##8, v##9
57 #define FETCH(v, w) do { \
58 v##0 = (w)->P[0]; v##1 = (w)->P[1]; \
59 v##2 = (w)->P[2]; v##3 = (w)->P[3]; \
60 v##4 = (w)->P[4]; v##5 = (w)->P[5]; \
61 v##6 = (w)->P[6]; v##7 = (w)->P[7]; \
62 v##8 = (w)->P[8]; v##9 = (w)->P[9]; \
64 #define STASH(w, v) do { \
65 (w)->P[0] = v##0; (w)->P[1] = v##1; \
66 (w)->P[2] = v##2; (w)->P[3] = v##3; \
67 (w)->P[4] = v##4; (w)->P[5] = v##5; \
68 (w)->P[6] = v##6; (w)->P[7] = v##7; \
69 (w)->P[8] = v##8; (w)->P[9] = v##9; \
72 #elif F25519_IMPL == 10
73 /* Elements x of GF(2^255 - 19) are represented by 26 signed integers x_i: x
74 * = SUM_{0<=i<26} x_i 2^ceil(255i/26); i.e., most pieces are 10 bits wide,
75 * except for pieces 5, 10, 15, 20, and 25 which have 9 bits.
78 typedef int32 dblpiece
;
79 typedef uint16 upiece
; typedef uint32 udblpiece
;
82 ((i) == 5 || (i) == 10 || (i) == 15 || (i) == 20 || (i) == 25 ? 9 : 10)
92 /*----- Debugging machinery -----------------------------------------------*/
94 #if defined(F25519_DEBUG) || defined(TEST_RIG)
101 static mp
*get_2p255m91(void)
106 p
= mp_setbit(p
, MP_ZERO
, 255);
107 mp_build(&m19
, &w19
, &w19
+ 1);
108 p
= mp_sub(p
, p
, &m19
);
112 DEF_FDUMP(fdump
, piece
, PIECEWD
, NPIECE
, 32, get_2p255m91())
116 /*----- Loading and storing -----------------------------------------------*/
118 /* --- @f25519_load@ --- *
120 * Arguments: @f25519 *z@ = where to store the result
121 * @const octet xv[32]@ = source to read
125 * Use: Reads an element of %$\gf{2^{255} - 19}$% in external
126 * representation from @xv@ and stores it in @z@.
128 * External representation is little-endian base-256. Some
129 * elements have multiple encodings, which are not produced by
130 * correct software; use of noncanonical encodings is not an
131 * error, and toleration of them is considered a performance
135 void f25519_load(f25519
*z
, const octet xv
[32])
137 #if F25519_IMPL == 26
139 uint32 xw0
= LOAD32_L(xv
+ 0), xw1
= LOAD32_L(xv
+ 4),
140 xw2
= LOAD32_L(xv
+ 8), xw3
= LOAD32_L(xv
+ 12),
141 xw4
= LOAD32_L(xv
+ 16), xw5
= LOAD32_L(xv
+ 20),
142 xw6
= LOAD32_L(xv
+ 24), xw7
= LOAD32_L(xv
+ 28);
143 piece
PIECES(x
), b
, c
;
145 /* First, split the 32-bit words into the irregularly-sized pieces we need
146 * for the field representation. These pieces are all positive. We'll do
147 * the sign correction afterwards.
149 * It may be that the top bit of the input is set. Avoid trouble by
150 * folding that back round into the bottom piece of the representation.
152 * Here, we briefly have 0 <= x_0 < 2^26 + 19, but will resolve this later.
153 * Otherwise, we have 0 <= x_{2i} < 2^26, and 0 <= x_{2i+1} < 2^25.
155 x0
= ((xw0
<< 0)&0x03ffffff) + 19*((xw7
>> 31)&0x00000001);
156 x1
= ((xw1
<< 6)&0x01ffffc0) | ((xw0
>> 26)&0x0000003f);
157 x2
= ((xw2
<< 13)&0x03ffe000) | ((xw1
>> 19)&0x00001fff);
158 x3
= ((xw3
<< 19)&0x01f80000) | ((xw2
>> 13)&0x0007ffff);
159 x4
= ((xw3
>> 6)&0x03ffffff);
160 x5
= (xw4
<< 0)&0x01ffffff;
161 x6
= ((xw5
<< 7)&0x03ffff80) | ((xw4
>> 25)&0x0000007f);
162 x7
= ((xw6
<< 13)&0x01ffe000) | ((xw5
>> 19)&0x00001fff);
163 x8
= ((xw7
<< 20)&0x03f00000) | ((xw6
>> 12)&0x000fffff);
164 x9
= ((xw7
>> 6)&0x01ffffff);
166 /* Next, we convert these pieces into a roughly balanced signed
167 * representation. For each piece, check to see if its top bit is set. If
168 * it is, then lend a bit to the next piece over. For x_9, this needs to
169 * be carried around, which is a little fiddly. In particular, we delay
170 * the carry until after all of the pieces have been balanced. If we don't
171 * do this, then we have to do a more expensive test for nonzeroness to
172 * decide whether to lend a bit leftwards rather than just testing a single
175 * This fixes up the anomalous size of x_0: the loan of a bit becomes an
176 * actual carry if x_0 >= 2^26. By the end, then, we have:
182 * Note that we don't try for a canonical representation here: both upper
183 * and lower bounds are achievable.
185 * All of the x_i at this point are positive, so we don't need to do
186 * anything weird when masking them.
188 b
= x9
&B24
; c
= 19&((b
>> 19) - (b
>> 24)); x9
-= b
<< 1;
189 b
= x8
&B25
; x9
+= b
>> 25; x8
-= b
<< 1;
190 b
= x7
&B24
; x8
+= b
>> 24; x7
-= b
<< 1;
191 b
= x6
&B25
; x7
+= b
>> 25; x6
-= b
<< 1;
192 b
= x5
&B24
; x6
+= b
>> 24; x5
-= b
<< 1;
193 b
= x4
&B25
; x5
+= b
>> 25; x4
-= b
<< 1;
194 b
= x3
&B24
; x4
+= b
>> 24; x3
-= b
<< 1;
195 b
= x2
&B25
; x3
+= b
>> 25; x2
-= b
<< 1;
196 b
= x1
&B24
; x2
+= b
>> 24; x1
-= b
<< 1;
197 b
= x0
&B25
; x1
+= (b
>> 25) + (x0
>> 26); x0
= (x0
&M26
) - (b
<< 1);
200 /* And with that, we're done. */
203 #elif F25519_IMPL == 10
206 unsigned i
, j
, n
, wd
;
210 /* First, just get the content out of the buffer. */
211 for (i
= j
= a
= n
= 0, wd
= 10; j
< NPIECE
; i
++) {
212 a
|= (uint32
)xv
[i
] << n
; n
+= 8;
220 /* There's a little bit left over from the top byte. Carry it into the low
223 x
[0] += 19*(int)(a
&MASK(n
));
225 /* Next, convert the pieces into a roughly balanced signed representation.
226 * If a piece's top bit is set, lend a bit to the next piece over. For
227 * x_25, this needs to be carried around, which is a bit fiddly.
229 b
= x
[NPIECE
- 1]&B8
;
230 c
= 19&((b
>> 3) - (b
>> 8));
231 x
[NPIECE
- 1] -= b
<< 1;
232 for (i
= NPIECE
- 2; i
> 0; i
--) {
239 x
[1] += (b
>> 9) + (x
[0] >> 10);
240 x
[0] = (x
[0]&M10
) - (b
<< 1) + c
;
242 /* And we're done. */
243 for (i
= 0; i
< NPIECE
; i
++) z
->P
[i
] = x
[i
];
248 /* --- @f25519_store@ --- *
250 * Arguments: @octet zv[32]@ = where to write the result
251 * @const f25519 *x@ = the field element to write
255 * Use: Stores a field element in the given octet vector in external
256 * representation. A canonical encoding is always stored, so,
257 * in particular, the top bit of @xv[31]@ is always left clear.
260 void f25519_store(octet zv
[32], const f25519
*x
)
262 #if F25519_IMPL == 26
264 piece
PIECES(x
), PIECES(y
), c
, d
;
265 uint32 zw0
, zw1
, zw2
, zw3
, zw4
, zw5
, zw6
, zw7
;
270 /* First, propagate the carries throughout the pieces. By the end of this,
271 * we'll have all of the pieces canonically sized and positive, and maybe
272 * there'll be (signed) carry out. The carry c is in { -1, 0, +1 }, and
273 * the remaining value will be in the half-open interval [0, 2^255). The
274 * whole represented value is then x + 2^255 c.
276 * It's worth paying careful attention to the bounds. We assume that we
277 * start out with |x_i| <= 2^30. We start by cutting off and reducing the
278 * carry c_9 from the topmost piece, x_9. This leaves 0 <= x_9 < 2^25; and
279 * we'll have |c_9| <= 2^5. We multiply this by 19 and we'll add it onto
280 * x_0 and propagate the carries: but what bounds can we calculate on x
283 * Let o_i = floor(51 i/2). We have X_i = SUM_{0<=j<i} x_j 2^{o_i}, so
284 * x = X_10. We see, inductively, that |X_i| < 2^{31+o_{i-1}}: X_0 = 0;
285 * |x_i| <= 2^30; and |X_{i+1}| = |X_i + x_i 2^{o_i}| <= |X_i| + 2^{30+o_i}
286 * < 2^{31+o_i}. Then x = X_9 + 2^230 x_9, and we have better bounds for
289 * -2^235 < x + 19 c_9 < 2^255 + 2^235
291 * Here, the x_i are signed, so we must be cautious about bithacking them.
293 c
= ASR(piece
, x9
, 25); x9
= (upiece
)x9
&M25
;
294 x0
+= 19*c
; c
= ASR(piece
, x0
, 26); x0
= (upiece
)x0
&M26
;
295 x1
+= c
; c
= ASR(piece
, x1
, 25); x1
= (upiece
)x1
&M25
;
296 x2
+= c
; c
= ASR(piece
, x2
, 26); x2
= (upiece
)x2
&M26
;
297 x3
+= c
; c
= ASR(piece
, x3
, 25); x3
= (upiece
)x3
&M25
;
298 x4
+= c
; c
= ASR(piece
, x4
, 26); x4
= (upiece
)x4
&M26
;
299 x5
+= c
; c
= ASR(piece
, x5
, 25); x5
= (upiece
)x5
&M25
;
300 x6
+= c
; c
= ASR(piece
, x6
, 26); x6
= (upiece
)x6
&M26
;
301 x7
+= c
; c
= ASR(piece
, x7
, 25); x7
= (upiece
)x7
&M25
;
302 x8
+= c
; c
= ASR(piece
, x8
, 26); x8
= (upiece
)x8
&M26
;
303 x9
+= c
; c
= ASR(piece
, x9
, 25); x9
= (upiece
)x9
&M25
;
305 /* Now we have a slightly fiddly job to do. If c = +1, or if c = 0 and
306 * x >= 2^255 - 19, then we should subtract 2^255 - 19 from the whole
307 * value; if c = -1 then we should add 2^255 - 19; and otherwise we should
310 * But conditional behaviour is bad, m'kay. So here's what we do instead.
312 * The first job is to sort out what we wanted to do. If c = -1 then we
313 * want to (a) invert the constant addend and (b) feed in a carry-in;
314 * otherwise, we don't.
319 /* Now do the addition/subtraction. Remember that all of the x_i are
320 * nonnegative, so shifting and masking are safe and easy.
322 d
+= x0
+ (19 ^ (M26
&m
)); y0
= d
&M26
; d
>>= 26;
323 d
+= x1
+ (M25
&m
); y1
= d
&M25
; d
>>= 25;
324 d
+= x2
+ (M26
&m
); y2
= d
&M26
; d
>>= 26;
325 d
+= x3
+ (M25
&m
); y3
= d
&M25
; d
>>= 25;
326 d
+= x4
+ (M26
&m
); y4
= d
&M26
; d
>>= 26;
327 d
+= x5
+ (M25
&m
); y5
= d
&M25
; d
>>= 25;
328 d
+= x6
+ (M26
&m
); y6
= d
&M26
; d
>>= 26;
329 d
+= x7
+ (M25
&m
); y7
= d
&M25
; d
>>= 25;
330 d
+= x8
+ (M26
&m
); y8
= d
&M26
; d
>>= 26;
331 d
+= x9
+ (M25
&m
); y9
= d
&M25
; d
>>= 25;
333 /* The final carry-out is in d; since we only did addition, and the x_i are
334 * nonnegative, then d is in { 0, 1 }. We want to keep y, rather than x,
335 * if (a) c /= 0 (in which case we know that the old value was
336 * unsatisfactory), or (b) if d = 1 (in which case, if c = 0, we know that
337 * the subtraction didn't cause a borrow, so we must be in the case where
338 * 2^255 - 19 <= x < 2^255).
340 m
= NONZEROP(c
) | ~NONZEROP(d
- 1);
341 x0
= (y0
&m
) | (x0
&~m
); x1
= (y1
&m
) | (x1
&~m
);
342 x2
= (y2
&m
) | (x2
&~m
); x3
= (y3
&m
) | (x3
&~m
);
343 x4
= (y4
&m
) | (x4
&~m
); x5
= (y5
&m
) | (x5
&~m
);
344 x6
= (y6
&m
) | (x6
&~m
); x7
= (y7
&m
) | (x7
&~m
);
345 x8
= (y8
&m
) | (x8
&~m
); x9
= (y9
&m
) | (x9
&~m
);
347 /* Extract 32-bit words from the value. */
348 zw0
= ((x0
>> 0)&0x03ffffff) | (((uint32
)x1
<< 26)&0xfc000000);
349 zw1
= ((x1
>> 6)&0x0007ffff) | (((uint32
)x2
<< 19)&0xfff80000);
350 zw2
= ((x2
>> 13)&0x00001fff) | (((uint32
)x3
<< 13)&0xffffe000);
351 zw3
= ((x3
>> 19)&0x0000003f) | (((uint32
)x4
<< 6)&0xffffffc0);
352 zw4
= ((x5
>> 0)&0x01ffffff) | (((uint32
)x6
<< 25)&0xfe000000);
353 zw5
= ((x6
>> 7)&0x0007ffff) | (((uint32
)x7
<< 19)&0xfff80000);
354 zw6
= ((x7
>> 13)&0x00000fff) | (((uint32
)x8
<< 12)&0xfffff000);
355 zw7
= ((x8
>> 20)&0x0000003f) | (((uint32
)x9
<< 6)&0x7fffffc0);
357 /* Store the result as an octet string. */
358 STORE32_L(zv
+ 0, zw0
); STORE32_L(zv
+ 4, zw1
);
359 STORE32_L(zv
+ 8, zw2
); STORE32_L(zv
+ 12, zw3
);
360 STORE32_L(zv
+ 16, zw4
); STORE32_L(zv
+ 20, zw5
);
361 STORE32_L(zv
+ 24, zw6
); STORE32_L(zv
+ 28, zw7
);
363 #elif F25519_IMPL == 10
365 piece y
[NPIECE
], yy
[NPIECE
], c
, d
;
366 unsigned i
, j
, n
, wd
;
369 /* Before we do anything, copy the input so we can hack on it. */
370 for (i
= 0; i
< NPIECE
; i
++) y
[i
] = x
->P
[i
];
372 /* First, propagate the carries throughout the pieces.
374 * It's worth paying careful attention to the bounds. We assume that we
375 * start out with |y_i| <= 2^14. We start by cutting off and reducing the
376 * carry c_25 from the topmost piece, y_25. This leaves 0 <= y_25 < 2^9;
377 * and we'll have |c_25| <= 2^5. We multiply this by 19 and we'll ad it
378 * onto y_0 and propagte the carries: but what bounds can we calculate on
381 * Let o_i = floor(255 i/26). We have Y_i = SUM_{0<=j<i} y_j 2^{o_i}, so
382 * y = Y_26. We see, inductively, that |Y_i| < 2^{31+o_{i-1}}: Y_0 = 0;
383 * |y_i| <= 2^14; and |Y_{i+1}| = |Y_i + y_i 2^{o_i}| <= |Y_i| + 2^{14+o_i}
384 * < 2^{15+o_i}. Then x = Y_25 + 2^246 y_25, and we have better bounds for
387 * -2^251 < y + 19 c_25 < 2^255 + 2^251
389 * Here, the y_i are signed, so we must be cautious about bithacking them.
391 * (Rather closer than the 10-piece case above, but still doable in one
394 c
= 19*ASR(piece
, y
[NPIECE
- 1], 9);
395 y
[NPIECE
- 1] = (upiece
)y
[NPIECE
- 1]&M9
;
396 for (i
= 0; i
< NPIECE
; i
++) {
399 c
= ASR(piece
, y
[i
], wd
);
400 y
[i
] = (upiece
)y
[i
]&MASK(wd
);
403 /* Now the addition or subtraction. */
407 d
+= y
[0] + (19 ^ (M10
&m
));
410 for (i
= 1; i
< NPIECE
; i
++) {
412 d
+= y
[i
] + (MASK(wd
)&m
);
417 /* Choose which value to keep. */
418 m
= NONZEROP(c
) | ~NONZEROP(d
- 1);
419 for (i
= 0; i
< NPIECE
; i
++) y
[i
] = (yy
[i
]&m
) | (y
[i
]&~m
);
421 /* Store the result as an octet string. */
422 for (i
= j
= a
= n
= 0; i
< NPIECE
; i
++) {
423 a
|= (upiece
)y
[i
] << n
; n
+= PIECEWD(i
);
434 /* --- @f25519_set@ --- *
436 * Arguments: @f25519 *z@ = where to write the result
437 * @int a@ = a small-ish constant
441 * Use: Sets @z@ to equal @a@.
444 void f25519_set(f25519
*x
, int a
)
449 for (i
= 1; i
< NPIECE
; i
++) x
->P
[i
] = 0;
452 /*----- Basic arithmetic --------------------------------------------------*/
454 /* --- @f25519_add@ --- *
456 * Arguments: @f25519 *z@ = where to put the result (may alias @x@ or @y@)
457 * @const f25519 *x, *y@ = two operands
461 * Use: Set @z@ to the sum %$x + y$%.
464 void f25519_add(f25519
*z
, const f25519
*x
, const f25519
*y
)
466 #if F25519_IMPL == 26
467 z
->P
[0] = x
->P
[0] + y
->P
[0]; z
->P
[1] = x
->P
[1] + y
->P
[1];
468 z
->P
[2] = x
->P
[2] + y
->P
[2]; z
->P
[3] = x
->P
[3] + y
->P
[3];
469 z
->P
[4] = x
->P
[4] + y
->P
[4]; z
->P
[5] = x
->P
[5] + y
->P
[5];
470 z
->P
[6] = x
->P
[6] + y
->P
[6]; z
->P
[7] = x
->P
[7] + y
->P
[7];
471 z
->P
[8] = x
->P
[8] + y
->P
[8]; z
->P
[9] = x
->P
[9] + y
->P
[9];
472 #elif F25519_IMPL == 10
474 for (i
= 0; i
< NPIECE
; i
++) z
->P
[i
] = x
->P
[i
] + y
->P
[i
];
478 /* --- @f25519_sub@ --- *
480 * Arguments: @f25519 *z@ = where to put the result (may alias @x@ or @y@)
481 * @const f25519 *x, *y@ = two operands
485 * Use: Set @z@ to the difference %$x - y$%.
488 void f25519_sub(f25519
*z
, const f25519
*x
, const f25519
*y
)
490 #if F25519_IMPL == 26
491 z
->P
[0] = x
->P
[0] - y
->P
[0]; z
->P
[1] = x
->P
[1] - y
->P
[1];
492 z
->P
[2] = x
->P
[2] - y
->P
[2]; z
->P
[3] = x
->P
[3] - y
->P
[3];
493 z
->P
[4] = x
->P
[4] - y
->P
[4]; z
->P
[5] = x
->P
[5] - y
->P
[5];
494 z
->P
[6] = x
->P
[6] - y
->P
[6]; z
->P
[7] = x
->P
[7] - y
->P
[7];
495 z
->P
[8] = x
->P
[8] - y
->P
[8]; z
->P
[9] = x
->P
[9] - y
->P
[9];
496 #elif F25519_IMPL == 10
498 for (i
= 0; i
< NPIECE
; i
++) z
->P
[i
] = x
->P
[i
] - y
->P
[i
];
502 /* --- @f25519_neg@ --- *
504 * Arguments: @f25519 *z@ = where to put the result (may alias @x@)
505 * @const f25519 *x@ = an operand
512 void f25519_neg(f25519
*z
, const f25519
*x
)
514 #if F25519_IMPL == 26
515 z
->P
[0] = -x
->P
[0]; z
->P
[1] = -x
->P
[1];
516 z
->P
[2] = -x
->P
[2]; z
->P
[3] = -x
->P
[3];
517 z
->P
[4] = -x
->P
[4]; z
->P
[5] = -x
->P
[5];
518 z
->P
[6] = -x
->P
[6]; z
->P
[7] = -x
->P
[7];
519 z
->P
[8] = -x
->P
[8]; z
->P
[9] = -x
->P
[9];
520 #elif F25519_IMPL == 10
522 for (i
= 0; i
< NPIECE
; i
++) z
->P
[i
] = -x
->P
[i
];
526 /*----- Constant-time utilities -------------------------------------------*/
528 /* --- @f25519_pick2@ --- *
530 * Arguments: @f25519 *z@ = where to put the result (may alias @x@ or @y@)
531 * @const f25519 *x, *y@ = two operands
532 * @uint32 m@ = a mask
536 * Use: If @m@ is zero, set @z = y@; if @m@ is all-bits-set, then set
537 * @z = x@. If @m@ has some other value, then scramble @z@ in
541 void f25519_pick2(f25519
*z
, const f25519
*x
, const f25519
*y
, uint32 m
)
543 mask32 mm
= FIX_MASK32(m
);
545 #if F25519_IMPL == 26
546 z
->P
[0] = PICK2(x
->P
[0], y
->P
[0], mm
);
547 z
->P
[1] = PICK2(x
->P
[1], y
->P
[1], mm
);
548 z
->P
[2] = PICK2(x
->P
[2], y
->P
[2], mm
);
549 z
->P
[3] = PICK2(x
->P
[3], y
->P
[3], mm
);
550 z
->P
[4] = PICK2(x
->P
[4], y
->P
[4], mm
);
551 z
->P
[5] = PICK2(x
->P
[5], y
->P
[5], mm
);
552 z
->P
[6] = PICK2(x
->P
[6], y
->P
[6], mm
);
553 z
->P
[7] = PICK2(x
->P
[7], y
->P
[7], mm
);
554 z
->P
[8] = PICK2(x
->P
[8], y
->P
[8], mm
);
555 z
->P
[9] = PICK2(x
->P
[9], y
->P
[9], mm
);
556 #elif F25519_IMPL == 10
558 for (i
= 0; i
< NPIECE
; i
++) z
->P
[i
] = PICK2(x
->P
[i
], y
->P
[i
], mm
);
562 /* --- @f25519_pickn@ --- *
564 * Arguments: @f25519 *z@ = where to put the result
565 * @const f25519 *v@ = a table of entries
566 * @size_t n@ = the number of entries in @v@
567 * @size_t i@ = an index
571 * Use: If @0 <= i < n < 32@ then set @z = v[i]@. If @n >= 32@ then
572 * do something unhelpful; otherwise, if @i >= n@ then set @z@
576 void f25519_pickn(f25519
*z
, const f25519
*v
, size_t n
, size_t i
)
578 uint32 b
= (uint32
)1 << (31 - i
);
581 #if F25519_IMPL == 26
582 z
->P
[0] = z
->P
[1] = z
->P
[2] = z
->P
[3] = z
->P
[4] =
583 z
->P
[5] = z
->P
[6] = z
->P
[7] = z
->P
[8] = z
->P
[9] = 0;
586 CONDPICK(z
->P
[0], v
->P
[0], m
);
587 CONDPICK(z
->P
[1], v
->P
[1], m
);
588 CONDPICK(z
->P
[2], v
->P
[2], m
);
589 CONDPICK(z
->P
[3], v
->P
[3], m
);
590 CONDPICK(z
->P
[4], v
->P
[4], m
);
591 CONDPICK(z
->P
[5], v
->P
[5], m
);
592 CONDPICK(z
->P
[6], v
->P
[6], m
);
593 CONDPICK(z
->P
[7], v
->P
[7], m
);
594 CONDPICK(z
->P
[8], v
->P
[8], m
);
595 CONDPICK(z
->P
[9], v
->P
[9], m
);
598 #elif F25519_IMPL == 10
601 for (j
= 0; j
< NPIECE
; j
++) z
->P
[j
] = 0;
604 for (j
= 0; j
< NPIECE
; j
++) CONDPICK(z
->P
[j
], v
->P
[j
], m
);
610 /* --- @f25519_condswap@ --- *
612 * Arguments: @f25519 *x, *y@ = two operands
613 * @uint32 m@ = a mask
617 * Use: If @m@ is zero, do nothing; if @m@ is all-bits-set, then
618 * exchange @x@ and @y@. If @m@ has some other value, then
619 * scramble @x@ and @y@ in an unhelpful way.
622 void f25519_condswap(f25519
*x
, f25519
*y
, uint32 m
)
624 mask32 mm
= FIX_MASK32(m
);
626 #if F25519_IMPL == 26
627 CONDSWAP(x
->P
[0], y
->P
[0], mm
);
628 CONDSWAP(x
->P
[1], y
->P
[1], mm
);
629 CONDSWAP(x
->P
[2], y
->P
[2], mm
);
630 CONDSWAP(x
->P
[3], y
->P
[3], mm
);
631 CONDSWAP(x
->P
[4], y
->P
[4], mm
);
632 CONDSWAP(x
->P
[5], y
->P
[5], mm
);
633 CONDSWAP(x
->P
[6], y
->P
[6], mm
);
634 CONDSWAP(x
->P
[7], y
->P
[7], mm
);
635 CONDSWAP(x
->P
[8], y
->P
[8], mm
);
636 CONDSWAP(x
->P
[9], y
->P
[9], mm
);
637 #elif F25519_IMPL == 10
639 for (i
= 0; i
< NPIECE
; i
++) CONDSWAP(x
->P
[i
], y
->P
[i
], mm
);
643 /* --- @f25519_condneg@ --- *
645 * Arguments: @f25519 *z@ = where to put the result (may alias @x@)
646 * @const f25519 *x@ = an operand
647 * @uint32 m@ = a mask
651 * Use: If @m@ is zero, set @z = x@; if @m@ is all-bits-set, then set
652 * @z = -x@. If @m@ has some other value then scramble @z@ in
656 void f25519_condneg(f25519
*z
, const f25519
*x
, uint32 m
)
659 mask32 m_xor
= FIX_MASK32(m
);
661 # define CONDNEG(x) (((x) ^ m_xor) + m_add)
663 int s
= PICK2(-1, +1, m
);
664 # define CONDNEG(x) (s*(x))
667 #if F25519_IMPL == 26
668 z
->P
[0] = CONDNEG(x
->P
[0]);
669 z
->P
[1] = CONDNEG(x
->P
[1]);
670 z
->P
[2] = CONDNEG(x
->P
[2]);
671 z
->P
[3] = CONDNEG(x
->P
[3]);
672 z
->P
[4] = CONDNEG(x
->P
[4]);
673 z
->P
[5] = CONDNEG(x
->P
[5]);
674 z
->P
[6] = CONDNEG(x
->P
[6]);
675 z
->P
[7] = CONDNEG(x
->P
[7]);
676 z
->P
[8] = CONDNEG(x
->P
[8]);
677 z
->P
[9] = CONDNEG(x
->P
[9]);
678 #elif F25519_IMPL == 10
680 for (i
= 0; i
< NPIECE
; i
++) z
->P
[i
] = CONDNEG(x
->P
[i
]);
686 /*----- Multiplication ----------------------------------------------------*/
688 #if F25519_IMPL == 26
690 /* Let B = 2^63 - 1 be the largest value such that +B and -B can be
691 * represented in a double-precision piece. On entry, it must be the case
692 * that |X_i| <= M <= B - 2^25 for some M. If this is the case, then, on
693 * exit, we will have |Z_i| <= 2^25 + 19 M/2^25.
695 #define CARRYSTEP(z, x, m, b, f, xx, n) do { \
696 (z) = (dblpiece)((udblpiece)(x)&(m)) - (b) + \
697 (f)*ASR(dblpiece, (xx), (n)); \
699 #define CARRY_REDUCE(z, x) do { \
700 dblpiece PIECES(_t); \
702 /* Bias the input pieces. This keeps the carries and so on centred \
703 * around zero rather than biased positive. \
705 _t0 = (x##0) + B25; _t1 = (x##1) + B24; \
706 _t2 = (x##2) + B25; _t3 = (x##3) + B24; \
707 _t4 = (x##4) + B25; _t5 = (x##5) + B24; \
708 _t6 = (x##6) + B25; _t7 = (x##7) + B24; \
709 _t8 = (x##8) + B25; _t9 = (x##9) + B24; \
711 /* Calculate the reduced pieces. Careful with the bithacking. */ \
712 CARRYSTEP(z##0, _t0, M26, B25, 19, _t9, 25); \
713 CARRYSTEP(z##1, _t1, M25, B24, 1, _t0, 26); \
714 CARRYSTEP(z##2, _t2, M26, B25, 1, _t1, 25); \
715 CARRYSTEP(z##3, _t3, M25, B24, 1, _t2, 26); \
716 CARRYSTEP(z##4, _t4, M26, B25, 1, _t3, 25); \
717 CARRYSTEP(z##5, _t5, M25, B24, 1, _t4, 26); \
718 CARRYSTEP(z##6, _t6, M26, B25, 1, _t5, 25); \
719 CARRYSTEP(z##7, _t7, M25, B24, 1, _t6, 26); \
720 CARRYSTEP(z##8, _t8, M26, B25, 1, _t7, 25); \
721 CARRYSTEP(z##9, _t9, M25, B24, 1, _t8, 26); \
724 #elif F25519_IMPL == 10
726 /* Perform carry propagation on X. */
727 static void carry_reduce(dblpiece x
[NPIECE
])
729 /* Initial bounds: we assume |x_i| < 2^31 - 2^27. */
734 /* The result is nearly canonical, because we do sequential carry
735 * propagation, because smaller processors are more likely to prefer the
736 * smaller working set than the instruction-level parallelism.
738 * Start at x_23; truncate it to 10 bits, and propagate the carry to x_24.
739 * Truncate x_24 to 10 bits, and add the carry onto x_25. Truncate x_25 to
740 * 9 bits, and add 19 times the carry onto x_0. And so on.
742 * Let c_i be the portion of x_i to be carried onto x_{i+1}. I claim that
743 * |c_i| <= 2^22. Then the carry /into/ any x_i has magnitude at most
744 * 19*2^22 < 2^27 (allowing for the reduction as we carry from x_25 to
745 * x_0), and x_i after carry is bounded above by 2^31. Hence, the carry
746 * out is at most 2^22, as claimed.
748 * Once we reach x_23 for the second time, we start with |x_23| <= 2^9.
749 * The carry into x_23 is at most 2^27 as calculated above; so the carry
750 * out into x_24 has magnitude at most 2^17. In turn, |x_24| <= 2^9 before
751 * the carry, so is now no more than 2^18 in magnitude, and the carry out
752 * into x_25 is at most 2^8. This leaves |x_25| < 2^9 after carry
755 * Be careful with the bit hacking because the quantities involved are
759 /* For each piece, we bias it so that floor division (as done by an
760 * arithmetic right shift) and modulus (as done by bitwise-AND) does the
763 #define CARRY(i, wd, b, m) do { \
765 c = ASR(dblpiece, x[i], (wd)); \
766 x[i] = (dblpiece)((udblpiece)x[i]&(m)) - (b); \
769 { CARRY(23, 10, B9
, M10
); }
770 { x
[24] += c
; CARRY(24, 10, B9
, M10
); }
771 { x
[25] += c
; CARRY(25, 9, B8
, M9
); }
772 { x
[0] += 19*c
; CARRY( 0, 10, B9
, M10
); }
773 for (i
= 1; i
< 21; ) {
774 for (j
= i
+ 4; i
< j
; ) { x
[i
] += c
; CARRY( i
, 10, B9
, M10
); i
++; }
775 { x
[i
] += c
; CARRY( i
, 9, B8
, M9
); i
++; }
777 while (i
< 25) { x
[i
] += c
; CARRY( i
, 10, B9
, M10
); i
++; }
785 /* --- @f25519_mulconst@ --- *
787 * Arguments: @f25519 *z@ = where to put the result (may alias @x@)
788 * @const f25519 *x@ = an operand
789 * @long a@ = a small-ish constant; %$|a| < 2^{20}$%.
793 * Use: Set @z@ to the product %$a x$%.
796 void f25519_mulconst(f25519
*z
, const f25519
*x
, long a
)
798 #if F25519_IMPL == 26
801 dblpiece
PIECES(z
), aa
= a
;
805 /* Suppose that |x_i| <= 2^27, and |a| <= 2^23. Then we'll have
808 z0
= aa
*x0
; z1
= aa
*x1
; z2
= aa
*x2
; z3
= aa
*x3
; z4
= aa
*x4
;
809 z5
= aa
*x5
; z6
= aa
*x6
; z7
= aa
*x7
; z8
= aa
*x8
; z9
= aa
*x9
;
811 /* Following `CARRY_REDUCE', we'll have |z_i| <= 2^26. */
815 #elif F25519_IMPL == 10
820 for (i
= 0; i
< NPIECE
; i
++) y
[i
] = a
*x
->P
[i
];
822 for (i
= 0; i
< NPIECE
; i
++) z
->P
[i
] = y
[i
];
827 /* --- @f25519_mul@ --- *
829 * Arguments: @f25519 *z@ = where to put the result (may alias @x@ or @y@)
830 * @const f25519 *x, *y@ = two operands
834 * Use: Set @z@ to the product %$x y$%.
837 void f25519_mul(f25519
*z
, const f25519
*x
, const f25519
*y
)
839 #if F25519_IMPL == 26
841 piece
PIECES(x
), PIECES(y
);
845 FETCH(x
, x
); FETCH(y
, y
);
847 /* Suppose that |x_i|, |y_i| <= 2^27. Then we'll have
860 * all of which are less than 2^63 - 2^25.
863 #define M(a, b) ((dblpiece)(a)*(b))
865 19*(M(x2
, y8
) + M(x4
, y6
) + M(x6
, y4
) + M(x8
, y2
)) +
866 38*(M(x1
, y9
) + M(x3
, y7
) + M(x5
, y5
) + M(x7
, y3
) + M(x9
, y1
));
867 z1
= M(x0
, y1
) + M(x1
, y0
) +
868 19*(M(x2
, y9
) + M(x3
, y8
) + M(x4
, y7
) + M(x5
, y6
) +
869 M(x6
, y5
) + M(x7
, y4
) + M(x8
, y3
) + M(x9
, y2
));
870 z2
= M(x0
, y2
) + M(x2
, y0
) +
872 19*(M(x4
, y8
) + M(x6
, y6
) + M(x8
, y4
)) +
873 38*(M(x3
, y9
) + M(x5
, y7
) + M(x7
, y5
) + M(x9
, y3
));
874 z3
= M(x0
, y3
) + M(x1
, y2
) + M(x2
, y1
) + M(x3
, y0
) +
875 19*(M(x4
, y9
) + M(x5
, y8
) + M(x6
, y7
) +
876 M(x7
, y6
) + M(x8
, y5
) + M(x9
, y4
));
877 z4
= M(x0
, y4
) + M(x2
, y2
) + M(x4
, y0
) +
878 2*(M(x1
, y3
) + M(x3
, y1
)) +
879 19*(M(x6
, y8
) + M(x8
, y6
)) +
880 38*(M(x5
, y9
) + M(x7
, y7
) + M(x9
, y5
));
881 z5
= M(x0
, y5
) + M(x1
, y4
) + M(x2
, y3
) +
882 M(x3
, y2
) + M(x4
, y1
) + M(x5
, y0
) +
883 19*(M(x6
, y9
) + M(x7
, y8
) + M(x8
, y7
) + M(x9
, y6
));
884 z6
= M(x0
, y6
) + M(x2
, y4
) + M(x4
, y2
) + M(x6
, y0
) +
885 2*(M(x1
, y5
) + M(x3
, y3
) + M(x5
, y1
)) +
887 38*(M(x7
, y9
) + M(x9
, y7
));
888 z7
= M(x0
, y7
) + M(x1
, y6
) + M(x2
, y5
) + M(x3
, y4
) +
889 M(x4
, y3
) + M(x5
, y2
) + M(x6
, y1
) + M(x7
, y0
) +
890 19*(M(x8
, y9
) + M(x9
, y8
));
891 z8
= M(x0
, y8
) + M(x2
, y6
) + M(x4
, y4
) + M(x6
, y2
) + M(x8
, y0
) +
892 2*(M(x1
, y7
) + M(x3
, y5
) + M(x5
, y3
) + M(x7
, y1
)) +
894 z9
= M(x0
, y9
) + M(x1
, y8
) + M(x2
, y7
) + M(x3
, y6
) + M(x4
, y5
) +
895 M(x5
, y4
) + M(x6
, y3
) + M(x7
, y2
) + M(x8
, y1
) + M(x9
, y0
);
898 /* From above, we have |z_i| <= 2^63 - 2^25. A pass of `CARRY_REDUCE' will
899 * leave |z_i| <= 2^38 + 2^25; and a second pass will leave |z_i| <= 2^25 +
900 * 2^13, which is comfortable for an addition prior to the next
903 for (i
= 0; i
< 2; i
++) CARRY_REDUCE(z
, z
);
906 #elif F25519_IMPL == 10
908 dblpiece u
[NPIECE
], t
, tt
, p
;
911 /* This is unpleasant. Honestly, this table seems to be the best way of
914 static const unsigned short off
[NPIECE
] = {
915 0, 10, 20, 30, 40, 50, 59, 69, 79, 89, 99, 108, 118,
916 128, 138, 148, 157, 167, 177, 187, 197, 206, 216, 226, 236, 246
919 /* First pass: things we must multiply by 19 or 38. */
920 for (i
= 0; i
< NPIECE
- 1; i
++) {
922 for (j
= i
+ 1; j
< NPIECE
; j
++) {
923 k
= NPIECE
+ i
- j
; p
= (dblpiece
)x
->P
[j
]*y
->P
[k
];
924 if (off
[i
] < off
[j
] + off
[k
] - 255) tt
+= p
;
927 u
[i
] = 19*(t
+ 2*tt
);
931 /* Second pass: things we must multiply by 1 or 2. */
932 for (i
= 0; i
< NPIECE
; i
++) {
934 for (j
= 0; j
<= i
; j
++) {
935 k
= i
- j
; p
= (dblpiece
)x
->P
[j
]*y
->P
[k
];
936 if (off
[i
] < off
[j
] + off
[k
]) tt
+= p
;
942 /* And we're done. */
944 for (i
= 0; i
< NPIECE
; i
++) z
->P
[i
] = u
[i
];
949 /* --- @f25519_sqr@ --- *
951 * Arguments: @f25519 *z@ = where to put the result (may alias @x@ or @y@)
952 * @const f25519 *x@ = an operand
956 * Use: Set @z@ to the square %$x^2$%.
959 void f25519_sqr(f25519
*z
, const f25519
*x
)
961 #if F25519_IMPL == 26
969 /* See `f25519_mul' for bounds. */
971 #define M(a, b) ((dblpiece)(a)*(b))
973 38*(M(x2
, x8
) + M(x4
, x6
) + M(x5
, x5
)) +
974 76*(M(x1
, x9
) + M(x3
, x7
));
976 38*(M(x2
, x9
) + M(x3
, x8
) + M(x4
, x7
) + M(x5
, x6
));
977 z2
= 2*(M(x0
, x2
) + M(x1
, x1
)) +
980 76*(M(x3
, x9
) + M(x5
, x7
));
981 z3
= 2*(M(x0
, x3
) + M(x1
, x2
)) +
982 38*(M(x4
, x9
) + M(x5
, x8
) + M(x6
, x7
));
986 38*(M(x6
, x8
) + M(x7
, x7
)) +
988 z5
= 2*(M(x0
, x5
) + M(x1
, x4
) + M(x2
, x3
)) +
989 38*(M(x6
, x9
) + M(x7
, x8
));
990 z6
= 2*(M(x0
, x6
) + M(x2
, x4
) + M(x3
, x3
)) +
994 z7
= 2*(M(x0
, x7
) + M(x1
, x6
) + M(x2
, x5
) + M(x3
, x4
)) +
997 2*(M(x0
, x8
) + M(x2
, x6
)) +
998 4*(M(x1
, x7
) + M(x3
, x5
)) +
1000 z9
= 2*(M(x0
, x9
) + M(x1
, x8
) + M(x2
, x7
) + M(x3
, x6
) + M(x4
, x5
));
1003 /* See `f25519_mul' for details. */
1004 for (i
= 0; i
< 2; i
++) CARRY_REDUCE(z
, z
);
1007 #elif F25519_IMPL == 10
1008 f25519_mul(z
, x
, x
);
1012 /*----- More complicated things -------------------------------------------*/
1014 /* --- @f25519_inv@ --- *
1016 * Arguments: @f25519 *z@ = where to put the result (may alias @x@)
1017 * @const f25519 *x@ = an operand
1021 * Use: Stores in @z@ the multiplicative inverse %$x^{-1}$%. If
1022 * %$x = 0$% then @z@ is set to zero. This is considered a
1026 void f25519_inv(f25519
*z
, const f25519
*x
)
1028 f25519 t
, u
, t2
, t11
, t2p10m1
, t2p50m1
;
1031 #define SQRN(z, x, n) do { \
1032 f25519_sqr((z), (x)); \
1033 for (i = 1; i < (n); i++) f25519_sqr((z), (z)); \
1036 /* Calculate x^-1 = x^(p - 2) = x^(2^255 - 21), which also handles x = 0 as
1037 * intended. The addition chain here is from Bernstein's implementation; I
1038 * couldn't find a better one.
1039 */ /* step | value */
1040 f25519_sqr(&t2
, x
); /* 1 | 2 */
1041 SQRN(&u
, &t2
, 2); /* 3 | 8 */
1042 f25519_mul(&t
, &u
, x
); /* 4 | 9 */
1043 f25519_mul(&t11
, &t
, &t2
); /* 5 | 11 = 2^5 - 21 */
1044 f25519_sqr(&u
, &t11
); /* 6 | 22 */
1045 f25519_mul(&t
, &t
, &u
); /* 7 | 31 = 2^5 - 1 */
1046 SQRN(&u
, &t
, 5); /* 12 | 2^10 - 2^5 */
1047 f25519_mul(&t2p10m1
, &t
, &u
); /* 13 | 2^10 - 1 */
1048 SQRN(&u
, &t2p10m1
, 10); /* 23 | 2^20 - 2^10 */
1049 f25519_mul(&t
, &t2p10m1
, &u
); /* 24 | 2^20 - 1 */
1050 SQRN(&u
, &t
, 20); /* 44 | 2^40 - 2^20 */
1051 f25519_mul(&t
, &t
, &u
); /* 45 | 2^40 - 1 */
1052 SQRN(&u
, &t
, 10); /* 55 | 2^50 - 2^10 */
1053 f25519_mul(&t2p50m1
, &t2p10m1
, &u
); /* 56 | 2^50 - 1 */
1054 SQRN(&u
, &t2p50m1
, 50); /* 106 | 2^100 - 2^50 */
1055 f25519_mul(&t
, &t2p50m1
, &u
); /* 107 | 2^100 - 1 */
1056 SQRN(&u
, &t
, 100); /* 207 | 2^200 - 2^100 */
1057 f25519_mul(&t
, &t
, &u
); /* 208 | 2^200 - 1 */
1058 SQRN(&u
, &t
, 50); /* 258 | 2^250 - 2^50 */
1059 f25519_mul(&t
, &t2p50m1
, &u
); /* 259 | 2^250 - 1 */
1060 SQRN(&u
, &t
, 5); /* 264 | 2^255 - 2^5 */
1061 f25519_mul(z
, &u
, &t11
); /* 265 | 2^255 - 21 */
1066 /* --- @f25519_quosqrt@ --- *
1068 * Arguments: @f25519 *z@ = where to put the result (may alias @x@ or @y@)
1069 * @const f25519 *x, *y@ = two operands
1071 * Returns: Zero if successful, @-1@ if %$x/y$% is not a square.
1073 * Use: Stores in @z@ the one of the square roots %$\pm\sqrt{x/y}$%.
1074 * If %$x = y = 0% then the result is zero; if %$y = 0$% but %$x
1075 * \ne 0$% then the operation fails. If you wanted a specific
1076 * square root then you'll have to pick it yourself.
1079 static const piece sqrtm1_pieces
[NPIECE
] = {
1080 #if F25519_IMPL == 26
1081 -32595792, -7943725, 9377950, 3500415, 12389472,
1082 -272473, -25146209, -2005654, 326686, 11406482
1083 #elif F25519_IMPL == 10
1084 176, -88, 161, 157, -485, -196, -231, -220, -416,
1085 -169, -255, 50, 189, -89, -266, -32, 202, -511,
1086 423, 357, 248, -249, 80, 288, 50, 174
1089 #define SQRTM1 ((const f25519 *)sqrtm1_pieces)
1091 int f25519_quosqrt(f25519
*z
, const f25519
*x
, const f25519
*y
)
1093 f25519 t
, u
, v
, w
, t15
;
1094 octet xb
[32], b0
[32], b1
[32];
1099 #define SQRN(z, x, n) do { \
1100 f25519_sqr((z), (x)); \
1101 for (i = 1; i < (n); i++) f25519_sqr((z), (z)); \
1104 /* This is a bit tricky; the algorithm is loosely based on Bernstein, Duif,
1105 * Lange, Schwabe, and Yang, `High-speed high-security signatures',
1106 * 2011-09-26, https://ed25519.cr.yp.to/ed25519-20110926.pdf.
1108 f25519_mul(&v
, x
, y
);
1110 /* Now for an addition chain. */ /* step | value */
1111 f25519_sqr(&u
, &v
); /* 1 | 2 */
1112 f25519_mul(&t
, &u
, &v
); /* 2 | 3 */
1113 SQRN(&u
, &t
, 2); /* 4 | 12 */
1114 f25519_mul(&t15
, &u
, &t
); /* 5 | 15 */
1115 f25519_sqr(&u
, &t15
); /* 6 | 30 */
1116 f25519_mul(&t
, &u
, &v
); /* 7 | 31 = 2^5 - 1 */
1117 SQRN(&u
, &t
, 5); /* 12 | 2^10 - 2^5 */
1118 f25519_mul(&t
, &u
, &t
); /* 13 | 2^10 - 1 */
1119 SQRN(&u
, &t
, 10); /* 23 | 2^20 - 2^10 */
1120 f25519_mul(&u
, &u
, &t
); /* 24 | 2^20 - 1 */
1121 SQRN(&u
, &u
, 10); /* 34 | 2^30 - 2^10 */
1122 f25519_mul(&t
, &u
, &t
); /* 35 | 2^30 - 1 */
1123 f25519_sqr(&u
, &t
); /* 36 | 2^31 - 2 */
1124 f25519_mul(&t
, &u
, &v
); /* 37 | 2^31 - 1 */
1125 SQRN(&u
, &t
, 31); /* 68 | 2^62 - 2^31 */
1126 f25519_mul(&t
, &u
, &t
); /* 69 | 2^62 - 1 */
1127 SQRN(&u
, &t
, 62); /* 131 | 2^124 - 2^62 */
1128 f25519_mul(&t
, &u
, &t
); /* 132 | 2^124 - 1 */
1129 SQRN(&u
, &t
, 124); /* 256 | 2^248 - 2^124 */
1130 f25519_mul(&t
, &u
, &t
); /* 257 | 2^248 - 1 */
1131 f25519_sqr(&u
, &t
); /* 258 | 2^249 - 2 */
1132 f25519_mul(&t
, &u
, &v
); /* 259 | 2^249 - 1 */
1133 SQRN(&t
, &t
, 3); /* 262 | 2^252 - 8 */
1134 f25519_sqr(&u
, &t
); /* 263 | 2^253 - 16 */
1135 f25519_mul(&t
, &u
, &t
); /* 264 | 3*2^252 - 24 */
1136 f25519_mul(&t
, &t
, &t15
); /* 265 | 3*2^252 - 9 */
1137 f25519_mul(&w
, &t
, &v
); /* 266 | 3*2^252 - 8 */
1139 /* Awesome. Now let me explain. Let v be a square in GF(p), and let w =
1140 * v^(3*2^252 - 8). In particular, let's consider
1142 * v^2 w^4 = v^2 v^{3*2^254 - 32} = (v^{2^254 - 10})^3
1144 * But 2^254 - 10 = ((2^255 - 19) - 1)/2 = (p - 1)/2. Since v is a square,
1145 * it has order dividing (p - 1)/2, and therefore v^2 w^4 = 1 and
1149 * That in turn implies that w^2 = ±1/v. Now, recall that v = x y, and let
1150 * w' = w x. Then w'^2 = ±x^2/v = ±x/y. If y w'^2 = x then we set
1151 * z = w', since we have z^2 = x/y; otherwise let z = i w', where i^2 = -1,
1152 * so z^2 = -w^2 = x/y, and we're done.
1154 * The easiest way to compare is to encode. This isn't as wasteful as it
1155 * sounds: the hard part is normalizing the representations, which we have
1158 f25519_mul(&w
, &w
, x
);
1160 f25519_mul(&t
, &t
, y
);
1162 f25519_store(xb
, x
);
1163 f25519_store(b0
, &t
);
1164 f25519_store(b1
, &u
);
1165 f25519_mul(&u
, &w
, SQRTM1
);
1167 m
= -ct_memeq(b0
, xb
, 32);
1168 rc
= PICK2(0, rc
, m
);
1169 f25519_pick2(z
, &w
, &u
, m
);
1170 m
= -ct_memeq(b1
, xb
, 32);
1171 rc
= PICK2(0, rc
, m
);
1173 /* And we're done. */
1177 /*----- Test rig ----------------------------------------------------------*/
1181 #include <mLib/macros.h>
1182 #include <mLib/report.h>
1183 #include <mLib/str.h>
1184 #include <mLib/testrig.h>
1186 static void fixdstr(dstr
*d
)
1189 die(1, "invalid length for f25519");
1190 else if (d
->len
< 32) {
1192 memset(d
->buf
+ d
->len
, 0, 32 - d
->len
);
1197 static void cvt_f25519(const char *buf
, dstr
*d
)
1199 dstr dd
= DSTR_INIT
;
1201 type_hex
.cvt(buf
, &dd
); fixdstr(&dd
);
1202 dstr_ensure(d
, sizeof(f25519
)); d
->len
= sizeof(f25519
);
1203 f25519_load((f25519
*)d
->buf
, (const octet
*)dd
.buf
);
1207 static void dump_f25519(dstr
*d
, FILE *fp
)
1208 { fdump(stderr
, "???", (const piece
*)d
->buf
); }
1210 static void cvt_f25519_ref(const char *buf
, dstr
*d
)
1211 { type_hex
.cvt(buf
, d
); fixdstr(d
); }
1213 static void dump_f25519_ref(dstr
*d
, FILE *fp
)
1217 f25519_load(&x
, (const octet
*)d
->buf
);
1218 fdump(stderr
, "???", x
.P
);
1221 static int eq(const f25519
*x
, dstr
*d
)
1222 { octet b
[32]; f25519_store(b
, x
); return (MEMCMP(b
, ==, d
->buf
, 32)); }
1224 static const test_type
1225 type_f25519
= { cvt_f25519
, dump_f25519
},
1226 type_f25519_ref
= { cvt_f25519_ref
, dump_f25519_ref
};
1228 #define TEST_UNOP(op) \
1229 static int vrf_##op(dstr dv[]) \
1231 f25519 *x = (f25519 *)dv[0].buf; \
1235 f25519_##op(&z, x); \
1236 if (!eq(&z, &dv[1])) { \
1238 fprintf(stderr, "failed!\n"); \
1239 fdump(stderr, "x", x->P); \
1240 fdump(stderr, "calc", z.P); \
1241 f25519_load(&zz, (const octet *)dv[1].buf); \
1242 fdump(stderr, "z", zz.P); \
1252 #define TEST_BINOP(op) \
1253 static int vrf_##op(dstr dv[]) \
1255 f25519 *x = (f25519 *)dv[0].buf, *y = (f25519 *)dv[1].buf; \
1259 f25519_##op(&z, x, y); \
1260 if (!eq(&z, &dv[2])) { \
1262 fprintf(stderr, "failed!\n"); \
1263 fdump(stderr, "x", x->P); \
1264 fdump(stderr, "y", y->P); \
1265 fdump(stderr, "calc", z.P); \
1266 f25519_load(&zz, (const octet *)dv[2].buf); \
1267 fdump(stderr, "z", zz.P); \
1277 static int vrf_mulc(dstr dv
[])
1279 f25519
*x
= (f25519
*)dv
[0].buf
;
1280 long a
= *(const long *)dv
[1].buf
;
1284 f25519_mulconst(&z
, x
, a
);
1285 if (!eq(&z
, &dv
[2])) {
1287 fprintf(stderr
, "failed!\n");
1288 fdump(stderr
, "x", x
->P
);
1289 fprintf(stderr
, "a = %ld\n", a
);
1290 fdump(stderr
, "calc", z
.P
);
1291 f25519_load(&zz
, (const octet
*)dv
[2].buf
);
1292 fdump(stderr
, "z", zz
.P
);
1298 static int vrf_condneg(dstr dv
[])
1300 f25519
*x
= (f25519
*)dv
[0].buf
;
1301 uint32 m
= *(uint32
*)dv
[1].buf
;
1305 f25519_condneg(&z
, x
, m
);
1306 if (!eq(&z
, &dv
[2])) {
1308 fprintf(stderr
, "failed!\n");
1309 fdump(stderr
, "x", x
->P
);
1310 fprintf(stderr
, "m = 0x%08lx\n", (unsigned long)m
);
1311 fdump(stderr
, "calc z", z
.P
);
1312 f25519_load(&z
, (const octet
*)dv
[1].buf
);
1313 fdump(stderr
, "want z", z
.P
);
1319 static int vrf_pick2(dstr dv
[])
1321 f25519
*x
= (f25519
*)dv
[0].buf
, *y
= (f25519
*)dv
[1].buf
;
1322 uint32 m
= *(uint32
*)dv
[2].buf
;
1326 f25519_pick2(&z
, x
, y
, m
);
1327 if (!eq(&z
, &dv
[3])) {
1329 fprintf(stderr
, "failed!\n");
1330 fdump(stderr
, "x", x
->P
);
1331 fdump(stderr
, "y", y
->P
);
1332 fprintf(stderr
, "m = 0x%08lx\n", (unsigned long)m
);
1333 fdump(stderr
, "calc z", z
.P
);
1334 f25519_load(&z
, (const octet
*)dv
[3].buf
);
1335 fdump(stderr
, "want z", z
.P
);
1341 static int vrf_pickn(dstr dv
[])
1345 size_t i
= *(uint32
*)dv
[1].buf
, j
, n
;
1350 for (q
= dv
[0].buf
, n
= 0; (p
= str_qword(&q
, 0)) != 0; n
++)
1351 { cvt_f25519(p
, &d
); v
[n
] = *(f25519
*)d
.buf
; }
1353 f25519_pickn(&z
, v
, n
, i
);
1354 if (!eq(&z
, &dv
[2])) {
1356 fprintf(stderr
, "failed!\n");
1357 for (j
= 0; j
< n
; j
++) {
1358 fprintf(stderr
, "v[%2u]", (unsigned)j
);
1359 fdump(stderr
, "", v
[j
].P
);
1361 fprintf(stderr
, "i = %u\n", (unsigned)i
);
1362 fdump(stderr
, "calc z", z
.P
);
1363 f25519_load(&z
, (const octet
*)dv
[2].buf
);
1364 fdump(stderr
, "want z", z
.P
);
1371 static int vrf_condswap(dstr dv
[])
1373 f25519
*x
= (f25519
*)dv
[0].buf
, *y
= (f25519
*)dv
[1].buf
;
1374 f25519 xx
= *x
, yy
= *y
;
1375 uint32 m
= *(uint32
*)dv
[2].buf
;
1378 f25519_condswap(&xx
, &yy
, m
);
1379 if (!eq(&xx
, &dv
[3]) || !eq(&yy
, &dv
[4])) {
1381 fprintf(stderr
, "failed!\n");
1382 fdump(stderr
, "x", x
->P
);
1383 fdump(stderr
, "y", y
->P
);
1384 fprintf(stderr
, "m = 0x%08lx\n", (unsigned long)m
);
1385 fdump(stderr
, "calc xx", xx
.P
);
1386 fdump(stderr
, "calc yy", yy
.P
);
1387 f25519_load(&xx
, (const octet
*)dv
[3].buf
);
1388 f25519_load(&yy
, (const octet
*)dv
[4].buf
);
1389 fdump(stderr
, "want xx", xx
.P
);
1390 fdump(stderr
, "want yy", yy
.P
);
1396 static int vrf_quosqrt(dstr dv
[])
1398 f25519
*x
= (f25519
*)dv
[0].buf
, *y
= (f25519
*)dv
[1].buf
;
1403 if (dv
[2].len
) { fixdstr(&dv
[2]); fixdstr(&dv
[3]); }
1404 rc
= f25519_quosqrt(&z
, x
, y
);
1405 if (!dv
[2].len ?
!rc
: (rc
|| (!eq(&z
, &dv
[2]) && !eq(&z
, &dv
[3])))) {
1407 fprintf(stderr
, "failed!\n");
1408 fdump(stderr
, "x", x
->P
);
1409 fdump(stderr
, "y", y
->P
);
1410 if (rc
) fprintf(stderr
, "calc: FAIL\n");
1411 else fdump(stderr
, "calc", z
.P
);
1413 fprintf(stderr
, "exp: FAIL\n");
1415 f25519_load(&zz
, (const octet
*)dv
[2].buf
);
1416 fdump(stderr
, "z", zz
.P
);
1417 f25519_load(&zz
, (const octet
*)dv
[3].buf
);
1418 fdump(stderr
, "z'", zz
.P
);
1425 static int vrf_sub_mulc_add_sub_mul(dstr dv
[])
1427 f25519
*u
= (f25519
*)dv
[0].buf
, *v
= (f25519
*)dv
[1].buf
,
1428 *w
= (f25519
*)dv
[3].buf
, *x
= (f25519
*)dv
[4].buf
,
1429 *y
= (f25519
*)dv
[5].buf
;
1430 long a
= *(const long *)dv
[2].buf
;
1431 f25519 umv
, aumv
, wpaumv
, xmy
, z
, zz
;
1434 f25519_sub(&umv
, u
, v
);
1435 f25519_mulconst(&aumv
, &umv
, a
);
1436 f25519_add(&wpaumv
, w
, &aumv
);
1437 f25519_sub(&xmy
, x
, y
);
1438 f25519_mul(&z
, &wpaumv
, &xmy
);
1440 if (!eq(&z
, &dv
[6])) {
1442 fprintf(stderr
, "failed!\n");
1443 fdump(stderr
, "u", u
->P
);
1444 fdump(stderr
, "v", v
->P
);
1445 fdump(stderr
, "u - v", umv
.P
);
1446 fprintf(stderr
, "a = %ld\n", a
);
1447 fdump(stderr
, "a (u - v)", aumv
.P
);
1448 fdump(stderr
, "w + a (u - v)", wpaumv
.P
);
1449 fdump(stderr
, "x", x
->P
);
1450 fdump(stderr
, "y", y
->P
);
1451 fdump(stderr
, "x - y", xmy
.P
);
1452 fdump(stderr
, "(x - y) (w + a (u - v))", z
.P
);
1453 f25519_load(&zz
, (const octet
*)dv
[6].buf
); fdump(stderr
, "z", zz
.P
);
1459 static test_chunk tests
[] = {
1460 { "add", vrf_add
, { &type_f25519
, &type_f25519
, &type_f25519_ref
} },
1461 { "sub", vrf_sub
, { &type_f25519
, &type_f25519
, &type_f25519_ref
} },
1462 { "neg", vrf_neg
, { &type_f25519
, &type_f25519_ref
} },
1463 { "condneg", vrf_condneg
,
1464 { &type_f25519
, &type_uint32
, &type_f25519_ref
} },
1465 { "mul", vrf_mul
, { &type_f25519
, &type_f25519
, &type_f25519_ref
} },
1466 { "mulconst", vrf_mulc
, { &type_f25519
, &type_long
, &type_f25519_ref
} },
1467 { "pick2", vrf_pick2
,
1468 { &type_f25519
, &type_f25519
, &type_uint32
, &type_f25519_ref
} },
1469 { "pickn", vrf_pickn
,
1470 { &type_string
, &type_uint32
, &type_f25519_ref
} },
1471 { "condswap", vrf_condswap
,
1472 { &type_f25519
, &type_f25519
, &type_uint32
,
1473 &type_f25519_ref
, &type_f25519_ref
} },
1474 { "sqr", vrf_sqr
, { &type_f25519
, &type_f25519_ref
} },
1475 { "inv", vrf_inv
, { &type_f25519
, &type_f25519_ref
} },
1476 { "quosqrt", vrf_quosqrt
,
1477 { &type_f25519
, &type_f25519
, &type_hex
, &type_hex
} },
1478 { "sub-mulc-add-sub-mul", vrf_sub_mulc_add_sub_mul
,
1479 { &type_f25519
, &type_f25519
, &type_long
, &type_f25519
,
1480 &type_f25519
, &type_f25519
, &type_f25519_ref
} },
1484 int main(int argc
, char *argv
[])
1486 test_run(argc
, argv
, tests
, SRCDIR
"/t/f25519");
1492 /*----- That's all, folks -------------------------------------------------*/