2 * Copyright (c) 2018 Thomas Pornin <pornin@bolet.org>
4 * Permission is hereby granted, free of charge, to any person obtaining
5 * a copy of this software and associated documentation files (the
6 * "Software"), to deal in the Software without restriction, including
7 * without limitation the rights to use, copy, modify, merge, publish,
8 * distribute, sublicense, and/or sell copies of the Software, and to
9 * permit persons to whom the Software is furnished to do so, subject to
10 * the following conditions:
12 * The above copyright notice and this permission notice shall be
13 * included in all copies or substantial portions of the Software.
15 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
16 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
17 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
18 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
19 * BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
20 * ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
21 * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
27 #if BR_INT128 || BR_UMUL128
33 static const unsigned char P256_G
[] = {
34 0x04, 0x6B, 0x17, 0xD1, 0xF2, 0xE1, 0x2C, 0x42, 0x47, 0xF8,
35 0xBC, 0xE6, 0xE5, 0x63, 0xA4, 0x40, 0xF2, 0x77, 0x03, 0x7D,
36 0x81, 0x2D, 0xEB, 0x33, 0xA0, 0xF4, 0xA1, 0x39, 0x45, 0xD8,
37 0x98, 0xC2, 0x96, 0x4F, 0xE3, 0x42, 0xE2, 0xFE, 0x1A, 0x7F,
38 0x9B, 0x8E, 0xE7, 0xEB, 0x4A, 0x7C, 0x0F, 0x9E, 0x16, 0x2B,
39 0xCE, 0x33, 0x57, 0x6B, 0x31, 0x5E, 0xCE, 0xCB, 0xB6, 0x40,
40 0x68, 0x37, 0xBF, 0x51, 0xF5
43 static const unsigned char P256_N
[] = {
44 0xFF, 0xFF, 0xFF, 0xFF, 0x00, 0x00, 0x00, 0x00, 0xFF, 0xFF,
45 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xBC, 0xE6, 0xFA, 0xAD,
46 0xA7, 0x17, 0x9E, 0x84, 0xF3, 0xB9, 0xCA, 0xC2, 0xFC, 0x63,
50 static const unsigned char *
51 api_generator(int curve
, size_t *len
)
58 static const unsigned char *
59 api_order(int curve
, size_t *len
)
67 api_xoff(int curve
, size_t *len
)
75 * A field element is encoded as four 64-bit integers, in basis 2^64.
76 * Values may reach up to 2^256-1. Montgomery multiplication is used.
80 static const uint64_t F256_R
[] = {
81 0x0000000000000001, 0xFFFFFFFF00000000,
82 0xFFFFFFFFFFFFFFFF, 0x00000000FFFFFFFE
85 /* Curve equation is y^2 = x^3 - 3*x + B. This constant is B*R mod p
86 (Montgomery representation of B). */
87 static const uint64_t P256_B_MONTY
[] = {
88 0xD89CDF6229C4BDDF, 0xACF005CD78843090,
89 0xE5A220ABF7212ED6, 0xDC30061D04874834
93 * Addition in the field.
96 f256_add(uint64_t *d
, const uint64_t *a
, const uint64_t *b
)
103 * Do the addition, with an extra carry in t.
105 w
= (unsigned __int128
)a
[0] + b
[0];
107 w
= (unsigned __int128
)a
[1] + b
[1] + (w
>> 64);
109 w
= (unsigned __int128
)a
[2] + b
[2] + (w
>> 64);
111 w
= (unsigned __int128
)a
[3] + b
[3] + (w
>> 64);
113 t
= (uint64_t)(w
>> 64);
116 * Fold carry t, using: 2^256 = 2^224 - 2^192 - 2^96 + 1 mod p.
118 w
= (unsigned __int128
)d
[0] + t
;
120 w
= (unsigned __int128
)d
[1] + (w
>> 64) - (t
<< 32);
122 /* Here, carry "w >> 64" can only be 0 or -1 */
123 w
= (unsigned __int128
)d
[2] - ((w
>> 64) & 1);
125 /* Again, carry is 0 or -1. But there can be carry only if t = 1,
126 in which case the addition of (t << 32) - t is positive. */
127 w
= (unsigned __int128
)d
[3] - ((w
>> 64) & 1) + (t
<< 32) - t
;
129 t
= (uint64_t)(w
>> 64);
132 * There can be an extra carry here, which we must fold again.
134 w
= (unsigned __int128
)d
[0] + t
;
136 w
= (unsigned __int128
)d
[1] + (w
>> 64) - (t
<< 32);
138 w
= (unsigned __int128
)d
[2] - ((w
>> 64) & 1);
140 d
[3] += (t
<< 32) - t
- (uint64_t)((w
>> 64) & 1);
147 cc
= _addcarry_u64(0, a
[0], b
[0], &d
[0]);
148 cc
= _addcarry_u64(cc
, a
[1], b
[1], &d
[1]);
149 cc
= _addcarry_u64(cc
, a
[2], b
[2], &d
[2]);
150 cc
= _addcarry_u64(cc
, a
[3], b
[3], &d
[3]);
153 * If there is a carry, then we want to subtract p, which we
154 * do by adding 2^256 - p.
157 cc
= _addcarry_u64(cc
, d
[0], 0, &d
[0]);
158 cc
= _addcarry_u64(cc
, d
[1], -(t
<< 32), &d
[1]);
159 cc
= _addcarry_u64(cc
, d
[2], -t
, &d
[2]);
160 cc
= _addcarry_u64(cc
, d
[3], (t
<< 32) - (t
<< 1), &d
[3]);
163 * We have to do it again if there still is a carry.
166 cc
= _addcarry_u64(cc
, d
[0], 0, &d
[0]);
167 cc
= _addcarry_u64(cc
, d
[1], -(t
<< 32), &d
[1]);
168 cc
= _addcarry_u64(cc
, d
[2], -t
, &d
[2]);
169 (void)_addcarry_u64(cc
, d
[3], (t
<< 32) - (t
<< 1), &d
[3]);
175 * Subtraction in the field.
178 f256_sub(uint64_t *d
, const uint64_t *a
, const uint64_t *b
)
185 w
= (unsigned __int128
)a
[0] - b
[0];
187 w
= (unsigned __int128
)a
[1] - b
[1] - ((w
>> 64) & 1);
189 w
= (unsigned __int128
)a
[2] - b
[2] - ((w
>> 64) & 1);
191 w
= (unsigned __int128
)a
[3] - b
[3] - ((w
>> 64) & 1);
193 t
= (uint64_t)(w
>> 64) & 1;
196 * If there is a borrow (t = 1), then we must add the modulus
197 * p = 2^256 - 2^224 + 2^192 + 2^96 - 1.
199 w
= (unsigned __int128
)d
[0] - t
;
201 w
= (unsigned __int128
)d
[1] + (t
<< 32) - ((w
>> 64) & 1);
203 /* Here, carry "w >> 64" can only be 0 or +1 */
204 w
= (unsigned __int128
)d
[2] + (w
>> 64);
206 /* Again, carry is 0 or +1 */
207 w
= (unsigned __int128
)d
[3] + (w
>> 64) - (t
<< 32) + t
;
209 t
= (uint64_t)(w
>> 64) & 1;
212 * There may be again a borrow, in which case we must add the
215 w
= (unsigned __int128
)d
[0] - t
;
217 w
= (unsigned __int128
)d
[1] + (t
<< 32) - ((w
>> 64) & 1);
219 w
= (unsigned __int128
)d
[2] + (w
>> 64);
221 d
[3] += (uint64_t)(w
>> 64) - (t
<< 32) + t
;
228 cc
= _subborrow_u64(0, a
[0], b
[0], &d
[0]);
229 cc
= _subborrow_u64(cc
, a
[1], b
[1], &d
[1]);
230 cc
= _subborrow_u64(cc
, a
[2], b
[2], &d
[2]);
231 cc
= _subborrow_u64(cc
, a
[3], b
[3], &d
[3]);
234 * If there is a borrow, then we need to add p. We (virtually)
235 * add 2^256, then subtract 2^256 - p.
238 cc
= _subborrow_u64(0, d
[0], t
, &d
[0]);
239 cc
= _subborrow_u64(cc
, d
[1], -(t
<< 32), &d
[1]);
240 cc
= _subborrow_u64(cc
, d
[2], -t
, &d
[2]);
241 cc
= _subborrow_u64(cc
, d
[3], (t
<< 32) - (t
<< 1), &d
[3]);
244 * If there still is a borrow, then we need to add p again.
247 cc
= _subborrow_u64(0, d
[0], t
, &d
[0]);
248 cc
= _subborrow_u64(cc
, d
[1], -(t
<< 32), &d
[1]);
249 cc
= _subborrow_u64(cc
, d
[2], -t
, &d
[2]);
250 (void)_subborrow_u64(cc
, d
[3], (t
<< 32) - (t
<< 1), &d
[3]);
256 * Montgomery multiplication in the field.
259 f256_montymul(uint64_t *d
, const uint64_t *a
, const uint64_t *b
)
263 uint64_t x
, f
, t0
, t1
, t2
, t3
, t4
;
264 unsigned __int128 z
, ff
;
268 * When computing d <- d + a[u]*b, we also add f*p such
269 * that d + a[u]*b + f*p is a multiple of 2^64. Since
270 * p = -1 mod 2^64, we can compute f = d[0] + a[u]*b[0] mod 2^64.
274 * Step 1: t <- (a[0]*b + f*p) / 2^64
275 * We have f = a[0]*b[0] mod 2^64. Since p = -1 mod 2^64, this
276 * ensures that (a[0]*b + f*p) is a multiple of 2^64.
278 * We also have: f*p = f*2^256 - f*2^224 + f*2^192 + f*2^96 - f.
281 z
= (unsigned __int128
)b
[0] * x
;
283 z
= (unsigned __int128
)b
[1] * x
+ (z
>> 64) + (uint64_t)(f
<< 32);
285 z
= (unsigned __int128
)b
[2] * x
+ (z
>> 64) + (uint64_t)(f
>> 32);
287 z
= (unsigned __int128
)b
[3] * x
+ (z
>> 64) + f
;
289 t3
= (uint64_t)(z
>> 64);
290 ff
= ((unsigned __int128
)f
<< 64) - ((unsigned __int128
)f
<< 32);
291 z
= (unsigned __int128
)t2
+ (uint64_t)ff
;
293 z
= (unsigned __int128
)t3
+ (z
>> 64) + (ff
>> 64);
295 t4
= (uint64_t)(z
>> 64);
298 * Steps 2 to 4: t <- (t + a[i]*b + f*p) / 2^64
300 for (i
= 1; i
< 4; i
++) {
303 /* t <- (t + x*b - f) / 2^64 */
304 z
= (unsigned __int128
)b
[0] * x
+ t0
;
306 z
= (unsigned __int128
)b
[1] * x
+ t1
+ (z
>> 64);
308 z
= (unsigned __int128
)b
[2] * x
+ t2
+ (z
>> 64);
310 z
= (unsigned __int128
)b
[3] * x
+ t3
+ (z
>> 64);
314 t4
= (uint64_t)(z
>> 64);
316 /* t <- t + f*2^32, carry in the upper half of z */
317 z
= (unsigned __int128
)t0
+ (uint64_t)(f
<< 32);
319 z
= (z
>> 64) + (unsigned __int128
)t1
+ (uint64_t)(f
>> 32);
322 /* t <- t + f*2^192 - f*2^160 + f*2^128 */
323 ff
= ((unsigned __int128
)f
<< 64)
324 - ((unsigned __int128
)f
<< 32) + f
;
325 z
= (z
>> 64) + (unsigned __int128
)t2
+ (uint64_t)ff
;
327 z
= (unsigned __int128
)t3
+ (z
>> 64) + (ff
>> 64);
329 t4
+= (uint64_t)(z
>> 64);
333 * At that point, we have computed t = (a*b + F*p) / 2^256, where
334 * F is a 256-bit integer whose limbs are the "f" coefficients
335 * in the steps above. We have:
340 * a*b + F*p <= (2^256-1)*(2^256-1) + p*(2^256-1)
341 * a*b + F*p <= 2^256*(2^256 - 2 + p) + 1 - p
344 * Since p < 2^256, it follows that:
345 * t4 can be only 0 or 1
347 * We can therefore subtract p from t, conditionally on t4, to
348 * get a nonnegative result that fits on 256 bits.
350 z
= (unsigned __int128
)t0
+ t4
;
352 z
= (unsigned __int128
)t1
- (t4
<< 32) + (z
>> 64);
354 z
= (unsigned __int128
)t2
- (z
>> 127);
356 t3
= t3
- (uint64_t)(z
>> 127) - t4
+ (t4
<< 32);
365 uint64_t x
, f
, t0
, t1
, t2
, t3
, t4
;
366 uint64_t zl
, zh
, ffl
, ffh
;
371 * When computing d <- d + a[u]*b, we also add f*p such
372 * that d + a[u]*b + f*p is a multiple of 2^64. Since
373 * p = -1 mod 2^64, we can compute f = d[0] + a[u]*b[0] mod 2^64.
377 * Step 1: t <- (a[0]*b + f*p) / 2^64
378 * We have f = a[0]*b[0] mod 2^64. Since p = -1 mod 2^64, this
379 * ensures that (a[0]*b + f*p) is a multiple of 2^64.
381 * We also have: f*p = f*2^256 - f*2^224 + f*2^192 + f*2^96 - f.
385 zl
= _umul128(b
[0], x
, &zh
);
389 zl
= _umul128(b
[1], x
, &zh
);
390 k
= _addcarry_u64(0, zl
, t0
, &zl
);
391 (void)_addcarry_u64(k
, zh
, 0, &zh
);
392 k
= _addcarry_u64(0, zl
, f
<< 32, &zl
);
393 (void)_addcarry_u64(k
, zh
, 0, &zh
);
397 zl
= _umul128(b
[2], x
, &zh
);
398 k
= _addcarry_u64(0, zl
, t1
, &zl
);
399 (void)_addcarry_u64(k
, zh
, 0, &zh
);
400 k
= _addcarry_u64(0, zl
, f
>> 32, &zl
);
401 (void)_addcarry_u64(k
, zh
, 0, &zh
);
405 zl
= _umul128(b
[3], x
, &zh
);
406 k
= _addcarry_u64(0, zl
, t2
, &zl
);
407 (void)_addcarry_u64(k
, zh
, 0, &zh
);
408 k
= _addcarry_u64(0, zl
, f
, &zl
);
409 (void)_addcarry_u64(k
, zh
, 0, &zh
);
413 t4
= _addcarry_u64(0, t3
, f
, &t3
);
414 k
= _subborrow_u64(0, t2
, f
<< 32, &t2
);
415 k
= _subborrow_u64(k
, t3
, f
>> 32, &t3
);
416 (void)_subborrow_u64(k
, t4
, 0, &t4
);
419 * Steps 2 to 4: t <- (t + a[i]*b + f*p) / 2^64
421 for (i
= 1; i
< 4; i
++) {
423 /* f = t0 + x * b[0]; -- computed below */
425 /* t <- (t + x*b - f) / 2^64 */
426 zl
= _umul128(b
[0], x
, &zh
);
427 k
= _addcarry_u64(0, zl
, t0
, &f
);
428 (void)_addcarry_u64(k
, zh
, 0, &t0
);
430 zl
= _umul128(b
[1], x
, &zh
);
431 k
= _addcarry_u64(0, zl
, t0
, &zl
);
432 (void)_addcarry_u64(k
, zh
, 0, &zh
);
433 k
= _addcarry_u64(0, zl
, t1
, &t0
);
434 (void)_addcarry_u64(k
, zh
, 0, &t1
);
436 zl
= _umul128(b
[2], x
, &zh
);
437 k
= _addcarry_u64(0, zl
, t1
, &zl
);
438 (void)_addcarry_u64(k
, zh
, 0, &zh
);
439 k
= _addcarry_u64(0, zl
, t2
, &t1
);
440 (void)_addcarry_u64(k
, zh
, 0, &t2
);
442 zl
= _umul128(b
[3], x
, &zh
);
443 k
= _addcarry_u64(0, zl
, t2
, &zl
);
444 (void)_addcarry_u64(k
, zh
, 0, &zh
);
445 k
= _addcarry_u64(0, zl
, t3
, &t2
);
446 (void)_addcarry_u64(k
, zh
, 0, &t3
);
448 t4
= _addcarry_u64(0, t3
, t4
, &t3
);
450 /* t <- t + f*2^32, carry in k */
451 k
= _addcarry_u64(0, t0
, f
<< 32, &t0
);
452 k
= _addcarry_u64(k
, t1
, f
>> 32, &t1
);
454 /* t <- t + f*2^192 - f*2^160 + f*2^128 */
455 m
= _subborrow_u64(0, f
, f
<< 32, &ffl
);
456 (void)_subborrow_u64(m
, f
, f
>> 32, &ffh
);
457 k
= _addcarry_u64(k
, t2
, ffl
, &t2
);
458 k
= _addcarry_u64(k
, t3
, ffh
, &t3
);
459 (void)_addcarry_u64(k
, t4
, 0, &t4
);
463 * At that point, we have computed t = (a*b + F*p) / 2^256, where
464 * F is a 256-bit integer whose limbs are the "f" coefficients
465 * in the steps above. We have:
470 * a*b + F*p <= (2^256-1)*(2^256-1) + p*(2^256-1)
471 * a*b + F*p <= 2^256*(2^256 - 2 + p) + 1 - p
474 * Since p < 2^256, it follows that:
475 * t4 can be only 0 or 1
477 * We can therefore subtract p from t, conditionally on t4, to
478 * get a nonnegative result that fits on 256 bits.
480 k
= _addcarry_u64(0, t0
, t4
, &t0
);
481 k
= _addcarry_u64(k
, t1
, -(t4
<< 32), &t1
);
482 k
= _addcarry_u64(k
, t2
, -t4
, &t2
);
483 (void)_addcarry_u64(k
, t3
, (t4
<< 32) - (t4
<< 1), &t3
);
494 * Montgomery squaring in the field; currently a basic wrapper around
495 * multiplication (inline, should be optimized away).
496 * TODO: see if some extra speed can be gained here.
499 f256_montysquare(uint64_t *d
, const uint64_t *a
)
501 f256_montymul(d
, a
, a
);
505 * Convert to Montgomery representation.
508 f256_tomonty(uint64_t *d
, const uint64_t *a
)
512 * If R = 2^256 mod p, then R2 = R^2 mod p; and the Montgomery
513 * multiplication of a by R2 is: a*R2/R = a*R mod p, i.e. the
514 * conversion to Montgomery representation.
516 static const uint64_t R2
[] = {
523 f256_montymul(d
, a
, R2
);
527 * Convert from Montgomery representation.
530 f256_frommonty(uint64_t *d
, const uint64_t *a
)
533 * Montgomery multiplication by 1 is division by 2^256 modulo p.
535 static const uint64_t one
[] = { 1, 0, 0, 0 };
537 f256_montymul(d
, a
, one
);
541 * Inversion in the field. If the source value is 0 modulo p, then this
542 * returns 0 or p. This function uses Montgomery representation.
545 f256_invert(uint64_t *d
, const uint64_t *a
)
548 * We compute a^(p-2) mod p. The exponent pattern (from high to
550 * - 32 bits of value 1
551 * - 31 bits of value 0
553 * - 96 bits of value 0
554 * - 94 bits of value 1
557 * To speed up the square-and-multiply algorithm, we precompute
564 memcpy(t
, a
, sizeof t
);
565 for (i
= 0; i
< 30; i
++) {
566 f256_montysquare(t
, t
);
567 f256_montymul(t
, t
, a
);
570 memcpy(r
, t
, sizeof t
);
571 for (i
= 224; i
>= 0; i
--) {
572 f256_montysquare(r
, r
);
578 f256_montymul(r
, r
, a
);
583 f256_montymul(r
, r
, t
);
587 memcpy(d
, r
, sizeof r
);
591 * Finalize reduction.
592 * Input value fits on 256 bits. This function subtracts p if and only
593 * if the input is greater than or equal to p.
596 f256_final_reduce(uint64_t *a
)
600 uint64_t t0
, t1
, t2
, t3
, cc
;
604 * We add 2^224 - 2^192 - 2^96 + 1 to a. If there is no carry,
605 * then a < p; otherwise, the addition result we computed is
606 * the value we must return.
608 z
= (unsigned __int128
)a
[0] + 1;
610 z
= (unsigned __int128
)a
[1] + (z
>> 64) - ((uint64_t)1 << 32);
612 z
= (unsigned __int128
)a
[2] - (z
>> 127);
614 z
= (unsigned __int128
)a
[3] - (z
>> 127) + 0xFFFFFFFF;
616 cc
= -(uint64_t)(z
>> 64);
618 a
[0] ^= cc
& (a
[0] ^ t0
);
619 a
[1] ^= cc
& (a
[1] ^ t1
);
620 a
[2] ^= cc
& (a
[2] ^ t2
);
621 a
[3] ^= cc
& (a
[3] ^ t3
);
625 uint64_t t0
, t1
, t2
, t3
, m
;
628 k
= _addcarry_u64(0, a
[0], (uint64_t)1, &t0
);
629 k
= _addcarry_u64(k
, a
[1], -((uint64_t)1 << 32), &t1
);
630 k
= _addcarry_u64(k
, a
[2], -(uint64_t)1, &t2
);
631 k
= _addcarry_u64(k
, a
[3], ((uint64_t)1 << 32) - 2, &t3
);
634 a
[0] ^= m
& (a
[0] ^ t0
);
635 a
[1] ^= m
& (a
[1] ^ t1
);
636 a
[2] ^= m
& (a
[2] ^ t2
);
637 a
[3] ^= m
& (a
[3] ^ t3
);
643 * Points in affine and Jacobian coordinates.
645 * - In affine coordinates, the point-at-infinity cannot be encoded.
646 * - Jacobian coordinates (X,Y,Z) correspond to affine (X/Z^2,Y/Z^3);
647 * if Z = 0 then this is the point-at-infinity.
661 * Decode a point. The returned point is in Jacobian coordinates, but
662 * with z = 1. If the encoding is invalid, or encodes a point which is
663 * not on the curve, or encodes the point at infinity, then this function
664 * returns 0. Otherwise, 1 is returned.
666 * The buffer is assumed to have length exactly 65 bytes.
669 point_decode(p256_jacobian
*P
, const unsigned char *buf
)
671 uint64_t x
[4], y
[4], t
[4], x3
[4], tt
;
675 * Header byte shall be 0x04.
677 r
= EQ(buf
[0], 0x04);
680 * Decode X and Y coordinates, and convert them into
681 * Montgomery representation.
683 x
[3] = br_dec64be(buf
+ 1);
684 x
[2] = br_dec64be(buf
+ 9);
685 x
[1] = br_dec64be(buf
+ 17);
686 x
[0] = br_dec64be(buf
+ 25);
687 y
[3] = br_dec64be(buf
+ 33);
688 y
[2] = br_dec64be(buf
+ 41);
689 y
[1] = br_dec64be(buf
+ 49);
690 y
[0] = br_dec64be(buf
+ 57);
695 * Verify y^2 = x^3 + A*x + B. In curve P-256, A = -3.
696 * Note that the Montgomery representation of 0 is 0. We must
697 * take care to apply the final reduction to make sure we have
700 f256_montysquare(t
, y
);
701 f256_montysquare(x3
, x
);
702 f256_montymul(x3
, x3
, x
);
707 f256_sub(t
, t
, P256_B_MONTY
);
708 f256_final_reduce(t
);
709 tt
= t
[0] | t
[1] | t
[2] | t
[3];
710 r
&= EQ((uint32_t)(tt
| (tt
>> 32)), 0);
713 * Return the point in Jacobian coordinates (and Montgomery
716 memcpy(P
->x
, x
, sizeof x
);
717 memcpy(P
->y
, y
, sizeof y
);
718 memcpy(P
->z
, F256_R
, sizeof F256_R
);
723 * Final conversion for a point:
724 * - The point is converted back to affine coordinates.
725 * - Final reduction is performed.
726 * - The point is encoded into the provided buffer.
728 * If the point is the point-at-infinity, all operations are performed,
729 * but the buffer contents are indeterminate, and 0 is returned. Otherwise,
730 * the encoded point is written in the buffer, and 1 is returned.
733 point_encode(unsigned char *buf
, const p256_jacobian
*P
)
735 uint64_t t1
[4], t2
[4], z
;
737 /* Set t1 = 1/z^2 and t2 = 1/z^3. */
738 f256_invert(t2
, P
->z
);
739 f256_montysquare(t1
, t2
);
740 f256_montymul(t2
, t2
, t1
);
742 /* Compute affine coordinates x (in t1) and y (in t2). */
743 f256_montymul(t1
, P
->x
, t1
);
744 f256_montymul(t2
, P
->y
, t2
);
746 /* Convert back from Montgomery representation, and finalize
748 f256_frommonty(t1
, t1
);
749 f256_frommonty(t2
, t2
);
750 f256_final_reduce(t1
);
751 f256_final_reduce(t2
);
755 br_enc64be(buf
+ 1, t1
[3]);
756 br_enc64be(buf
+ 9, t1
[2]);
757 br_enc64be(buf
+ 17, t1
[1]);
758 br_enc64be(buf
+ 25, t1
[0]);
759 br_enc64be(buf
+ 33, t2
[3]);
760 br_enc64be(buf
+ 41, t2
[2]);
761 br_enc64be(buf
+ 49, t2
[1]);
762 br_enc64be(buf
+ 57, t2
[0]);
764 /* Return success if and only if P->z != 0. */
765 z
= P
->z
[0] | P
->z
[1] | P
->z
[2] | P
->z
[3];
766 return NEQ((uint32_t)(z
| z
>> 32), 0);
770 * Point doubling in Jacobian coordinates: point P is doubled.
771 * Note: if the source point is the point-at-infinity, then the result is
772 * still the point-at-infinity, which is correct. Moreover, if the three
773 * coordinates were zero, then they still are zero in the returned value.
775 * (Note: this is true even without the final reduction: if the three
776 * coordinates are encoded as four words of value zero each, then the
777 * result will also have all-zero coordinate encodings, not the alternate
778 * encoding as the integer p.)
781 p256_double(p256_jacobian
*P
)
784 * Doubling formulas are:
787 * m = 3*(x + z^2)*(x - z^2)
789 * y' = m*(s - x') - 8*y^4
792 * These formulas work for all points, including points of order 2
793 * and points at infinity:
794 * - If y = 0 then z' = 0. But there is no such point in P-256
796 * - If z = 0 then z' = 0.
798 uint64_t t1
[4], t2
[4], t3
[4], t4
[4];
803 f256_montysquare(t1
, P
->z
);
806 * Compute x-z^2 in t2 and x+z^2 in t1.
808 f256_add(t2
, P
->x
, t1
);
809 f256_sub(t1
, P
->x
, t1
);
812 * Compute 3*(x+z^2)*(x-z^2) in t1.
814 f256_montymul(t3
, t1
, t2
);
815 f256_add(t1
, t3
, t3
);
816 f256_add(t1
, t3
, t1
);
819 * Compute 4*x*y^2 (in t2) and 2*y^2 (in t3).
821 f256_montysquare(t3
, P
->y
);
822 f256_add(t3
, t3
, t3
);
823 f256_montymul(t2
, P
->x
, t3
);
824 f256_add(t2
, t2
, t2
);
827 * Compute x' = m^2 - 2*s.
829 f256_montysquare(P
->x
, t1
);
830 f256_sub(P
->x
, P
->x
, t2
);
831 f256_sub(P
->x
, P
->x
, t2
);
834 * Compute z' = 2*y*z.
836 f256_montymul(t4
, P
->y
, P
->z
);
837 f256_add(P
->z
, t4
, t4
);
840 * Compute y' = m*(s - x') - 8*y^4. Note that we already have
843 f256_sub(t2
, t2
, P
->x
);
844 f256_montymul(P
->y
, t1
, t2
);
845 f256_montysquare(t4
, t3
);
846 f256_add(t4
, t4
, t4
);
847 f256_sub(P
->y
, P
->y
, t4
);
851 * Point addition (Jacobian coordinates): P1 is replaced with P1+P2.
852 * This function computes the wrong result in the following cases:
854 * - If P1 == 0 but P2 != 0
855 * - If P1 != 0 but P2 == 0
858 * In all three cases, P1 is set to the point at infinity.
860 * Returned value is 0 if one of the following occurs:
862 * - P1 and P2 have the same Y coordinate.
863 * - P1 == 0 and P2 == 0.
864 * - The Y coordinate of one of the points is 0 and the other point is
865 * the point at infinity.
867 * The third case cannot actually happen with valid points, since a point
868 * with Y == 0 is a point of order 2, and there is no point of order 2 on
871 * Therefore, assuming that P1 != 0 and P2 != 0 on input, then the caller
872 * can apply the following:
874 * - If the result is not the point at infinity, then it is correct.
875 * - Otherwise, if the returned value is 1, then this is a case of
876 * P1+P2 == 0, so the result is indeed the point at infinity.
877 * - Otherwise, P1 == P2, so a "double" operation should have been
880 * Note that you can get a returned value of 0 with a correct result,
881 * e.g. if P1 and P2 have the same Y coordinate, but distinct X coordinates.
884 p256_add(p256_jacobian
*P1
, const p256_jacobian
*P2
)
887 * Addtions formulas are:
895 * x3 = r^2 - h^3 - 2 * u1 * h^2
896 * y3 = r * (u1 * h^2 - x3) - s1 * h^3
899 uint64_t t1
[4], t2
[4], t3
[4], t4
[4], t5
[4], t6
[4], t7
[4], tt
;
903 * Compute u1 = x1*z2^2 (in t1) and s1 = y1*z2^3 (in t3).
905 f256_montysquare(t3
, P2
->z
);
906 f256_montymul(t1
, P1
->x
, t3
);
907 f256_montymul(t4
, P2
->z
, t3
);
908 f256_montymul(t3
, P1
->y
, t4
);
911 * Compute u2 = x2*z1^2 (in t2) and s2 = y2*z1^3 (in t4).
913 f256_montysquare(t4
, P1
->z
);
914 f256_montymul(t2
, P2
->x
, t4
);
915 f256_montymul(t5
, P1
->z
, t4
);
916 f256_montymul(t4
, P2
->y
, t5
);
919 * Compute h = h2 - u1 (in t2) and r = s2 - s1 (in t4).
920 * We need to test whether r is zero, so we will do some extra
923 f256_sub(t2
, t2
, t1
);
924 f256_sub(t4
, t4
, t3
);
925 f256_final_reduce(t4
);
926 tt
= t4
[0] | t4
[1] | t4
[2] | t4
[3];
927 ret
= (uint32_t)(tt
| (tt
>> 32));
928 ret
= (ret
| -ret
) >> 31;
931 * Compute u1*h^2 (in t6) and h^3 (in t5);
933 f256_montysquare(t7
, t2
);
934 f256_montymul(t6
, t1
, t7
);
935 f256_montymul(t5
, t7
, t2
);
938 * Compute x3 = r^2 - h^3 - 2*u1*h^2.
940 f256_montysquare(P1
->x
, t4
);
941 f256_sub(P1
->x
, P1
->x
, t5
);
942 f256_sub(P1
->x
, P1
->x
, t6
);
943 f256_sub(P1
->x
, P1
->x
, t6
);
946 * Compute y3 = r*(u1*h^2 - x3) - s1*h^3.
948 f256_sub(t6
, t6
, P1
->x
);
949 f256_montymul(P1
->y
, t4
, t6
);
950 f256_montymul(t1
, t5
, t3
);
951 f256_sub(P1
->y
, P1
->y
, t1
);
954 * Compute z3 = h*z1*z2.
956 f256_montymul(t1
, P1
->z
, P2
->z
);
957 f256_montymul(P1
->z
, t1
, t2
);
963 * Point addition (mixed coordinates): P1 is replaced with P1+P2.
964 * This is a specialised function for the case when P2 is a non-zero point
965 * in affine coordinates.
967 * This function computes the wrong result in the following cases:
972 * In both cases, P1 is set to the point at infinity.
974 * Returned value is 0 if one of the following occurs:
976 * - P1 and P2 have the same Y (affine) coordinate.
977 * - The Y coordinate of P2 is 0 and P1 is the point at infinity.
979 * The second case cannot actually happen with valid points, since a point
980 * with Y == 0 is a point of order 2, and there is no point of order 2 on
983 * Therefore, assuming that P1 != 0 on input, then the caller
984 * can apply the following:
986 * - If the result is not the point at infinity, then it is correct.
987 * - Otherwise, if the returned value is 1, then this is a case of
988 * P1+P2 == 0, so the result is indeed the point at infinity.
989 * - Otherwise, P1 == P2, so a "double" operation should have been
992 * Again, a value of 0 may be returned in some cases where the addition
996 p256_add_mixed(p256_jacobian
*P1
, const p256_affine
*P2
)
999 * Addtions formulas are:
1007 * x3 = r^2 - h^3 - 2 * u1 * h^2
1008 * y3 = r * (u1 * h^2 - x3) - s1 * h^3
1011 uint64_t t1
[4], t2
[4], t3
[4], t4
[4], t5
[4], t6
[4], t7
[4], tt
;
1015 * Compute u1 = x1 (in t1) and s1 = y1 (in t3).
1017 memcpy(t1
, P1
->x
, sizeof t1
);
1018 memcpy(t3
, P1
->y
, sizeof t3
);
1021 * Compute u2 = x2*z1^2 (in t2) and s2 = y2*z1^3 (in t4).
1023 f256_montysquare(t4
, P1
->z
);
1024 f256_montymul(t2
, P2
->x
, t4
);
1025 f256_montymul(t5
, P1
->z
, t4
);
1026 f256_montymul(t4
, P2
->y
, t5
);
1029 * Compute h = h2 - u1 (in t2) and r = s2 - s1 (in t4).
1030 * We need to test whether r is zero, so we will do some extra
1033 f256_sub(t2
, t2
, t1
);
1034 f256_sub(t4
, t4
, t3
);
1035 f256_final_reduce(t4
);
1036 tt
= t4
[0] | t4
[1] | t4
[2] | t4
[3];
1037 ret
= (uint32_t)(tt
| (tt
>> 32));
1038 ret
= (ret
| -ret
) >> 31;
1041 * Compute u1*h^2 (in t6) and h^3 (in t5);
1043 f256_montysquare(t7
, t2
);
1044 f256_montymul(t6
, t1
, t7
);
1045 f256_montymul(t5
, t7
, t2
);
1048 * Compute x3 = r^2 - h^3 - 2*u1*h^2.
1050 f256_montysquare(P1
->x
, t4
);
1051 f256_sub(P1
->x
, P1
->x
, t5
);
1052 f256_sub(P1
->x
, P1
->x
, t6
);
1053 f256_sub(P1
->x
, P1
->x
, t6
);
1056 * Compute y3 = r*(u1*h^2 - x3) - s1*h^3.
1058 f256_sub(t6
, t6
, P1
->x
);
1059 f256_montymul(P1
->y
, t4
, t6
);
1060 f256_montymul(t1
, t5
, t3
);
1061 f256_sub(P1
->y
, P1
->y
, t1
);
1064 * Compute z3 = h*z1*z2.
1066 f256_montymul(P1
->z
, P1
->z
, t2
);
1074 * Point addition (mixed coordinates, complete): P1 is replaced with P1+P2.
1075 * This is a specialised function for the case when P2 is a non-zero point
1076 * in affine coordinates.
1078 * This function returns the correct result in all cases.
1081 p256_add_complete_mixed(p256_jacobian
*P1
, const p256_affine
*P2
)
1084 * Addtions formulas, in the general case, are:
1092 * x3 = r^2 - h^3 - 2 * u1 * h^2
1093 * y3 = r * (u1 * h^2 - x3) - s1 * h^3
1096 * These formulas mishandle the two following cases:
1098 * - If P1 is the point-at-infinity (z1 = 0), then z3 is
1099 * incorrectly set to 0.
1101 * - If P1 = P2, then u1 = u2 and s1 = s2, and x3, y3 and z3
1104 * However, if P1 + P2 = 0, then u1 = u2 but s1 != s2, and then
1105 * we correctly get z3 = 0 (the point-at-infinity).
1107 * To fix the case P1 = 0, we perform at the end a copy of P2
1108 * over P1, conditional to z1 = 0.
1110 * For P1 = P2: in that case, both h and r are set to 0, and
1111 * we get x3, y3 and z3 equal to 0. We can test for that
1112 * occurrence to make a mask which will be all-one if P1 = P2,
1113 * or all-zero otherwise; then we can compute the double of P2
1114 * and add it, combined with the mask, to (x3,y3,z3).
1116 * Using the doubling formulas in p256_double() on (x2,y2),
1117 * simplifying since P2 is affine (i.e. z2 = 1, implicitly),
1120 * m = 3*(x2 + 1)*(x2 - 1)
1122 * y' = m*(s - x') - 8*y2^4
1124 * which requires only 6 multiplications. Added to the 11
1125 * multiplications of the normal mixed addition in Jacobian
1126 * coordinates, we get a cost of 17 multiplications in total.
1128 uint64_t t1
[4], t2
[4], t3
[4], t4
[4], t5
[4], t6
[4], t7
[4], tt
, zz
;
1132 * Set zz to -1 if P1 is the point at infinity, 0 otherwise.
1134 zz
= P1
->z
[0] | P1
->z
[1] | P1
->z
[2] | P1
->z
[3];
1135 zz
= ((zz
| -zz
) >> 63) - (uint64_t)1;
1138 * Compute u1 = x1 (in t1) and s1 = y1 (in t3).
1140 memcpy(t1
, P1
->x
, sizeof t1
);
1141 memcpy(t3
, P1
->y
, sizeof t3
);
1144 * Compute u2 = x2*z1^2 (in t2) and s2 = y2*z1^3 (in t4).
1146 f256_montysquare(t4
, P1
->z
);
1147 f256_montymul(t2
, P2
->x
, t4
);
1148 f256_montymul(t5
, P1
->z
, t4
);
1149 f256_montymul(t4
, P2
->y
, t5
);
1152 * Compute h = h2 - u1 (in t2) and r = s2 - s1 (in t4).
1155 f256_sub(t2
, t2
, t1
);
1156 f256_sub(t4
, t4
, t3
);
1159 * If both h = 0 and r = 0, then P1 = P2, and we want to set
1160 * the mask tt to -1; otherwise, the mask will be 0.
1162 f256_final_reduce(t2
);
1163 f256_final_reduce(t4
);
1164 tt
= t2
[0] | t2
[1] | t2
[2] | t2
[3] | t4
[0] | t4
[1] | t4
[2] | t4
[3];
1165 tt
= ((tt
| -tt
) >> 63) - (uint64_t)1;
1168 * Compute u1*h^2 (in t6) and h^3 (in t5);
1170 f256_montysquare(t7
, t2
);
1171 f256_montymul(t6
, t1
, t7
);
1172 f256_montymul(t5
, t7
, t2
);
1175 * Compute x3 = r^2 - h^3 - 2*u1*h^2.
1177 f256_montysquare(P1
->x
, t4
);
1178 f256_sub(P1
->x
, P1
->x
, t5
);
1179 f256_sub(P1
->x
, P1
->x
, t6
);
1180 f256_sub(P1
->x
, P1
->x
, t6
);
1183 * Compute y3 = r*(u1*h^2 - x3) - s1*h^3.
1185 f256_sub(t6
, t6
, P1
->x
);
1186 f256_montymul(P1
->y
, t4
, t6
);
1187 f256_montymul(t1
, t5
, t3
);
1188 f256_sub(P1
->y
, P1
->y
, t1
);
1191 * Compute z3 = h*z1.
1193 f256_montymul(P1
->z
, P1
->z
, t2
);
1196 * The "double" result, in case P1 = P2.
1200 * Compute z' = 2*y2 (in t1).
1202 f256_add(t1
, P2
->y
, P2
->y
);
1205 * Compute 2*(y2^2) (in t2) and s = 4*x2*(y2^2) (in t3).
1207 f256_montysquare(t2
, P2
->y
);
1208 f256_add(t2
, t2
, t2
);
1209 f256_add(t3
, t2
, t2
);
1210 f256_montymul(t3
, P2
->x
, t3
);
1213 * Compute m = 3*(x2^2 - 1) (in t4).
1215 f256_montysquare(t4
, P2
->x
);
1216 f256_sub(t4
, t4
, F256_R
);
1217 f256_add(t5
, t4
, t4
);
1218 f256_add(t4
, t4
, t5
);
1221 * Compute x' = m^2 - 2*s (in t5).
1223 f256_montysquare(t5
, t4
);
1228 * Compute y' = m*(s - x') - 8*y2^4 (in t6).
1230 f256_sub(t6
, t3
, t5
);
1231 f256_montymul(t6
, t6
, t4
);
1232 f256_montysquare(t7
, t2
);
1233 f256_sub(t6
, t6
, t7
);
1234 f256_sub(t6
, t6
, t7
);
1237 * We now have the alternate (doubling) coordinates in (t5,t6,t1).
1238 * We combine them with (x3,y3,z3).
1240 for (i
= 0; i
< 4; i
++) {
1241 P1
->x
[i
] |= tt
& t5
[i
];
1242 P1
->y
[i
] |= tt
& t6
[i
];
1243 P1
->z
[i
] |= tt
& t1
[i
];
1247 * If P1 = 0, then we get z3 = 0 (which is invalid); if z1 is 0,
1248 * then we want to replace the result with a copy of P2. The
1249 * test on z1 was done at the start, in the zz mask.
1251 for (i
= 0; i
< 4; i
++) {
1252 P1
->x
[i
] ^= zz
& (P1
->x
[i
] ^ P2
->x
[i
]);
1253 P1
->y
[i
] ^= zz
& (P1
->y
[i
] ^ P2
->y
[i
]);
1254 P1
->z
[i
] ^= zz
& (P1
->z
[i
] ^ F256_R
[i
]);
1260 * Inner function for computing a point multiplication. A window is
1261 * provided, with points 1*P to 15*P in affine coordinates.
1264 * - All provided points are valid points on the curve.
1265 * - Multiplier is non-zero, and smaller than the curve order.
1266 * - Everything is in Montgomery representation.
1269 point_mul_inner(p256_jacobian
*R
, const p256_affine
*W
,
1270 const unsigned char *k
, size_t klen
)
1275 memset(&Q
, 0, sizeof Q
);
1277 while (klen
-- > 0) {
1282 for (i
= 0; i
< 2; i
++) {
1295 bits
= (bk
>> 4) & 0x0F;
1299 * Lookup point in window. If the bits are 0,
1300 * we get something invalid, which is not a
1301 * problem because we will use it only if the
1302 * bits are non-zero.
1304 memset(&T
, 0, sizeof T
);
1305 for (n
= 0; n
< 15; n
++) {
1306 m
= -(uint64_t)EQ(bits
, n
+ 1);
1307 T
.x
[0] |= m
& W
[n
].x
[0];
1308 T
.x
[1] |= m
& W
[n
].x
[1];
1309 T
.x
[2] |= m
& W
[n
].x
[2];
1310 T
.x
[3] |= m
& W
[n
].x
[3];
1311 T
.y
[0] |= m
& W
[n
].y
[0];
1312 T
.y
[1] |= m
& W
[n
].y
[1];
1313 T
.y
[2] |= m
& W
[n
].y
[2];
1314 T
.y
[3] |= m
& W
[n
].y
[3];
1318 p256_add_mixed(&U
, &T
);
1321 * If qz is still 1, then Q was all-zeros, and this
1322 * is conserved through p256_double().
1324 m
= -(uint64_t)(bnz
& qz
);
1325 for (j
= 0; j
< 4; j
++) {
1326 Q
.x
[j
] |= m
& T
.x
[j
];
1327 Q
.y
[j
] |= m
& T
.y
[j
];
1328 Q
.z
[j
] |= m
& F256_R
[j
];
1330 CCOPY(bnz
& ~qz
, &Q
, &U
, sizeof Q
);
1339 * Convert a window from Jacobian to affine coordinates. A single
1340 * field inversion is used. This function works for windows up to
1343 * The destination array (aff[]) and the source array (jac[]) may
1344 * overlap, provided that the start of aff[] is not after the start of
1345 * jac[]. Even if the arrays do _not_ overlap, the source array is
1349 window_to_affine(p256_affine
*aff
, p256_jacobian
*jac
, int num
)
1352 * Convert the window points to affine coordinates. We use the
1353 * following trick to mutualize the inversion computation: if
1354 * we have z1, z2, z3, and z4, and want to inverse all of them,
1355 * we compute u = 1/(z1*z2*z3*z4), and then we have:
1361 * The partial products are computed recursively:
1363 * - on input (z_1,z_2), return (z_2,z_1) and z_1*z_2
1364 * - on input (z_1,z_2,... z_n):
1365 * recurse on (z_1,z_2,... z_(n/2)) -> r1 and m1
1366 * recurse on (z_(n/2+1),z_(n/2+2)... z_n) -> r2 and m2
1367 * multiply elements of r1 by m2 -> s1
1368 * multiply elements of r2 by m1 -> s2
1369 * return r1||r2 and m1*m2
1371 * In the example below, we suppose that we have 14 elements.
1372 * Let z1, z2,... zE be the 14 values to invert (index noted in
1373 * hexadecimal, starting at 1).
1376 * swap(z1, z2); z12 = z1*z2
1377 * swap(z3, z4); z34 = z3*z4
1378 * swap(z5, z6); z56 = z5*z6
1379 * swap(z7, z8); z78 = z7*z8
1380 * swap(z9, zA); z9A = z9*zA
1381 * swap(zB, zC); zBC = zB*zC
1382 * swap(zD, zE); zDE = zD*zE
1385 * z1 <- z1*z34, z2 <- z2*z34, z3 <- z3*z12, z4 <- z4*z12
1387 * z5 <- z5*z78, z6 <- z6*z78, z7 <- z7*z56, z8 <- z8*z56
1389 * z9 <- z9*zBC, zA <- zA*zBC, zB <- zB*z9A, zC <- zC*z9A
1393 * z1 <- z1*z5678, z2 <- z2*z5678, z3 <- z3*z5678, z4 <- z4*z5678
1394 * z5 <- z5*z1234, z6 <- z6*z1234, z7 <- z7*z1234, z8 <- z8*z1234
1395 * z12345678 = z1234*z5678
1396 * z9 <- z9*zDE, zA <- zA*zDE, zB <- zB*zDE, zC <- zC*zDE
1397 * zD <- zD*z9ABC, zE*z9ABC
1398 * z9ABCDE = z9ABC*zDE
1401 * multiply z1..z8 by z9ABCDE
1402 * multiply z9..zE by z12345678
1403 * final z = z12345678*z9ABCDE
1413 * First recursion step (pairwise swapping and multiplication).
1414 * If there is an odd number of elements, then we "invent" an
1415 * extra one with coordinate Z = 1 (in Montgomery representation).
1417 for (i
= 0; (i
+ 1) < num
; i
+= 2) {
1418 memcpy(zt
, jac
[i
].z
, sizeof zt
);
1419 memcpy(jac
[i
].z
, jac
[i
+ 1].z
, sizeof zt
);
1420 memcpy(jac
[i
+ 1].z
, zt
, sizeof zt
);
1421 f256_montymul(z
[i
>> 1], jac
[i
].z
, jac
[i
+ 1].z
);
1423 if ((num
& 1) != 0) {
1424 memcpy(z
[num
>> 1], jac
[num
- 1].z
, sizeof zt
);
1425 memcpy(jac
[num
- 1].z
, F256_R
, sizeof F256_R
);
1429 * Perform further recursion steps. At the entry of each step,
1430 * the process has been done for groups of 's' points. The
1431 * integer k is the log2 of s.
1433 for (k
= 1, s
= 2; s
< num
; k
++, s
<<= 1) {
1436 for (i
= 0; i
< num
; i
++) {
1437 f256_montymul(jac
[i
].z
, jac
[i
].z
, z
[(i
>> k
) ^ 1]);
1439 n
= (num
+ s
- 1) >> k
;
1440 for (i
= 0; i
< (n
>> 1); i
++) {
1441 f256_montymul(z
[i
], z
[i
<< 1], z
[(i
<< 1) + 1]);
1444 memmove(z
[n
>> 1], z
[n
], sizeof zt
);
1449 * Invert the final result, and convert all points.
1451 f256_invert(zt
, z
[0]);
1452 for (i
= 0; i
< num
; i
++) {
1453 f256_montymul(zv
, jac
[i
].z
, zt
);
1454 f256_montysquare(zu
, zv
);
1455 f256_montymul(zv
, zv
, zu
);
1456 f256_montymul(aff
[i
].x
, jac
[i
].x
, zu
);
1457 f256_montymul(aff
[i
].y
, jac
[i
].y
, zv
);
1462 * Multiply the provided point by an integer.
1464 * - Source point is a valid curve point.
1465 * - Source point is not the point-at-infinity.
1466 * - Integer is not 0, and is lower than the curve order.
1467 * If these conditions are not met, then the result is indeterminate
1468 * (but the process is still constant-time).
1471 p256_mul(p256_jacobian
*P
, const unsigned char *k
, size_t klen
)
1474 p256_affine aff
[15];
1475 p256_jacobian jac
[15];
1480 * Compute window, in Jacobian coordinates.
1483 for (i
= 2; i
< 16; i
++) {
1484 window
.jac
[i
- 1] = window
.jac
[(i
>> 1) - 1];
1486 p256_double(&window
.jac
[i
- 1]);
1488 p256_add(&window
.jac
[i
- 1], &window
.jac
[i
>> 1]);
1493 * Convert the window points to affine coordinates. Point
1494 * window[0] is the source point, already in affine coordinates.
1496 window_to_affine(window
.aff
, window
.jac
, 15);
1499 * Perform point multiplication.
1501 point_mul_inner(P
, window
.aff
, k
, klen
);
1505 * Precomputed window for the conventional generator: P256_Gwin[n]
1506 * contains (n+1)*G (affine coordinates, in Montgomery representation).
1508 static const p256_affine P256_Gwin
[] = {
1510 { 0x79E730D418A9143C, 0x75BA95FC5FEDB601,
1511 0x79FB732B77622510, 0x18905F76A53755C6 },
1512 { 0xDDF25357CE95560A, 0x8B4AB8E4BA19E45C,
1513 0xD2E88688DD21F325, 0x8571FF1825885D85 }
1516 { 0x850046D410DDD64D, 0xAA6AE3C1A433827D,
1517 0x732205038D1490D9, 0xF6BB32E43DCF3A3B },
1518 { 0x2F3648D361BEE1A5, 0x152CD7CBEB236FF8,
1519 0x19A8FB0E92042DBE, 0x78C577510A5B8A3B }
1522 { 0xFFAC3F904EEBC127, 0xB027F84A087D81FB,
1523 0x66AD77DD87CBBC98, 0x26936A3FB6FF747E },
1524 { 0xB04C5C1FC983A7EB, 0x583E47AD0861FE1A,
1525 0x788208311A2EE98E, 0xD5F06A29E587CC07 }
1528 { 0x74B0B50D46918DCC, 0x4650A6EDC623C173,
1529 0x0CDAACACE8100AF2, 0x577362F541B0176B },
1530 { 0x2D96F24CE4CBABA6, 0x17628471FAD6F447,
1531 0x6B6C36DEE5DDD22E, 0x84B14C394C5AB863 }
1534 { 0xBE1B8AAEC45C61F5, 0x90EC649A94B9537D,
1535 0x941CB5AAD076C20C, 0xC9079605890523C8 },
1536 { 0xEB309B4AE7BA4F10, 0x73C568EFE5EB882B,
1537 0x3540A9877E7A1F68, 0x73A076BB2DD1E916 }
1540 { 0x403947373E77664A, 0x55AE744F346CEE3E,
1541 0xD50A961A5B17A3AD, 0x13074B5954213673 },
1542 { 0x93D36220D377E44B, 0x299C2B53ADFF14B5,
1543 0xF424D44CEF639F11, 0xA4C9916D4A07F75F }
1546 { 0x0746354EA0173B4F, 0x2BD20213D23C00F7,
1547 0xF43EAAB50C23BB08, 0x13BA5119C3123E03 },
1548 { 0x2847D0303F5B9D4D, 0x6742F2F25DA67BDD,
1549 0xEF933BDC77C94195, 0xEAEDD9156E240867 }
1552 { 0x27F14CD19499A78F, 0x462AB5C56F9B3455,
1553 0x8F90F02AF02CFC6B, 0xB763891EB265230D },
1554 { 0xF59DA3A9532D4977, 0x21E3327DCF9EBA15,
1555 0x123C7B84BE60BBF0, 0x56EC12F27706DF76 }
1558 { 0x75C96E8F264E20E8, 0xABE6BFED59A7A841,
1559 0x2CC09C0444C8EB00, 0xE05B3080F0C4E16B },
1560 { 0x1EB7777AA45F3314, 0x56AF7BEDCE5D45E3,
1561 0x2B6E019A88B12F1A, 0x086659CDFD835F9B }
1564 { 0x2C18DBD19DC21EC8, 0x98F9868A0FCF8139,
1565 0x737D2CD648250B49, 0xCC61C94724B3428F },
1566 { 0x0C2B407880DD9E76, 0xC43A8991383FBE08,
1567 0x5F7D2D65779BE5D2, 0x78719A54EB3B4AB5 }
1570 { 0xEA7D260A6245E404, 0x9DE407956E7FDFE0,
1571 0x1FF3A4158DAC1AB5, 0x3E7090F1649C9073 },
1572 { 0x1A7685612B944E88, 0x250F939EE57F61C8,
1573 0x0C0DAA891EAD643D, 0x68930023E125B88E }
1576 { 0x04B71AA7D2697768, 0xABDEDEF5CA345A33,
1577 0x2409D29DEE37385E, 0x4EE1DF77CB83E156 },
1578 { 0x0CAC12D91CBB5B43, 0x170ED2F6CA895637,
1579 0x28228CFA8ADE6D66, 0x7FF57C9553238ACA }
1582 { 0xCCC425634B2ED709, 0x0E356769856FD30D,
1583 0xBCBCD43F559E9811, 0x738477AC5395B759 },
1584 { 0x35752B90C00EE17F, 0x68748390742ED2E3,
1585 0x7CD06422BD1F5BC1, 0xFBC08769C9E7B797 }
1588 { 0xA242A35BB0CF664A, 0x126E48F77F9707E3,
1589 0x1717BF54C6832660, 0xFAAE7332FD12C72E },
1590 { 0x27B52DB7995D586B, 0xBE29569E832237C2,
1591 0xE8E4193E2A65E7DB, 0x152706DC2EAA1BBB }
1594 { 0x72BCD8B7BC60055B, 0x03CC23EE56E27E4B,
1595 0xEE337424E4819370, 0xE2AA0E430AD3DA09 },
1596 { 0x40B8524F6383C45D, 0xD766355442A41B25,
1597 0x64EFA6DE778A4797, 0x2042170A7079ADF4 }
1602 * Multiply the conventional generator of the curve by the provided
1603 * integer. Return is written in *P.
1606 * - Integer is not 0, and is lower than the curve order.
1607 * If this conditions is not met, then the result is indeterminate
1608 * (but the process is still constant-time).
1611 p256_mulgen(p256_jacobian
*P
, const unsigned char *k
, size_t klen
)
1613 point_mul_inner(P
, P256_Gwin
, k
, klen
);
1617 * Return 1 if all of the following hold:
1620 * - k is lower than the curve order
1621 * Otherwise, return 0.
1623 * Constant-time behaviour: only klen may be observable.
1626 check_scalar(const unsigned char *k
, size_t klen
)
1636 for (u
= 0; u
< klen
; u
++) {
1641 for (u
= 0; u
< klen
; u
++) {
1642 c
|= -(int32_t)EQ0(c
) & CMP(k
[u
], P256_N
[u
]);
1647 return NEQ(z
, 0) & LT0(c
);
1651 api_mul(unsigned char *G
, size_t Glen
,
1652 const unsigned char *k
, size_t klen
, int curve
)
1661 r
= check_scalar(k
, klen
);
1662 r
&= point_decode(&P
, G
);
1663 p256_mul(&P
, k
, klen
);
1664 r
&= point_encode(G
, &P
);
1669 api_mulgen(unsigned char *R
,
1670 const unsigned char *k
, size_t klen
, int curve
)
1675 p256_mulgen(&P
, k
, klen
);
1676 point_encode(R
, &P
);
1681 api_muladd(unsigned char *A
, const unsigned char *B
, size_t len
,
1682 const unsigned char *x
, size_t xlen
,
1683 const unsigned char *y
, size_t ylen
, int curve
)
1686 * We might want to use Shamir's trick here: make a composite
1687 * window of u*P+v*Q points, to merge the two doubling-ladders
1688 * into one. This, however, has some complications:
1690 * - During the computation, we may hit the point-at-infinity.
1691 * Thus, we would need p256_add_complete_mixed() (complete
1692 * formulas for point addition), with a higher cost (17 muls
1695 * - A 4-bit window would be too large, since it would involve
1696 * 16*16-1 = 255 points. For the same window size as in the
1697 * p256_mul() case, we would need to reduce the window size
1698 * to 2 bits, and thus perform twice as many non-doubling
1701 * - The window may itself contain the point-at-infinity, and
1702 * thus cannot be in all generality be made of affine points.
1703 * Instead, we would need to make it a window of points in
1704 * Jacobian coordinates. Even p256_add_complete_mixed() would
1707 * For these reasons, the code below performs two separate
1708 * point multiplications, then computes the final point addition
1709 * (which is both a "normal" addition, and a doubling, to handle
1721 r
= point_decode(&P
, A
);
1722 p256_mul(&P
, x
, xlen
);
1724 p256_mulgen(&Q
, y
, ylen
);
1726 r
&= point_decode(&Q
, B
);
1727 p256_mul(&Q
, y
, ylen
);
1731 * The final addition may fail in case both points are equal.
1733 t
= p256_add(&P
, &Q
);
1734 f256_final_reduce(P
.z
);
1735 z
= P
.z
[0] | P
.z
[1] | P
.z
[2] | P
.z
[3];
1736 s
= EQ((uint32_t)(z
| (z
>> 32)), 0);
1740 * If s is 1 then either P+Q = 0 (t = 1) or P = Q (t = 0). So we
1741 * have the following:
1743 * s = 0, t = 0 return P (normal addition)
1744 * s = 0, t = 1 return P (normal addition)
1745 * s = 1, t = 0 return Q (a 'double' case)
1746 * s = 1, t = 1 report an error (P+Q = 0)
1748 CCOPY(s
& ~t
, &P
, &Q
, sizeof Q
);
1749 point_encode(A
, &P
);
1754 /* see bearssl_ec.h */
1755 const br_ec_impl br_ec_p256_m64
= {
1756 (uint32_t)0x00800000,
1765 /* see bearssl_ec.h */
1767 br_ec_p256_m64_get(void)
1769 return &br_ec_p256_m64
;
1774 /* see bearssl_ec.h */
1776 br_ec_p256_m64_get(void)