1 // Copyright (c) 2012 The Chromium Authors. All rights reserved.
2 // Use of this source code is governed by a BSD-style license that can be
3 // found in the LICENSE file.
5 // This is an implementation of the P224 elliptic curve group. It's written to
6 // be short and simple rather than fast, although it's still constant-time.
8 // See http://www.imperialviolet.org/2010/12/04/ecc.html ([1]) for background.
10 #include "crypto/p224.h"
16 #include "base/sys_byteorder.h"
20 using base::HostToNet32;
21 using base::NetToHost32;
23 // Field element functions.
25 // The field that we're dealing with is ℤ/pℤ where p = 2**224 - 2**96 + 1.
27 // Field elements are represented by a FieldElement, which is a typedef to an
28 // array of 8 uint32_t's. The value of a FieldElement, a, is:
29 // a[0] + 2**28·a[1] + 2**56·a[1] + ... + 2**196·a[7]
31 // Using 28-bit limbs means that there's only 4 bits of headroom, which is less
32 // than we would really like. But it has the useful feature that we hit 2**224
33 // exactly, making the reflections during a reduce much nicer.
35 using crypto::p224::FieldElement;
37 // kP is the P224 prime.
38 const FieldElement kP = {
40 268435455, 268435455, 268435455, 268435455,
43 void Contract(FieldElement* inout);
45 // IsZero returns 0xffffffff if a == 0 mod p and 0 otherwise.
46 uint32_t IsZero(const FieldElement& a) {
48 memcpy(&minimal, &a, sizeof(minimal));
51 uint32_t is_zero = 0, is_p = 0;
52 for (unsigned i = 0; i < 8; i++) {
53 is_zero |= minimal[i];
54 is_p |= minimal[i] - kP[i];
57 // If either is_zero or is_p is 0, then we should return 1.
58 is_zero |= is_zero >> 16;
59 is_zero |= is_zero >> 8;
60 is_zero |= is_zero >> 4;
61 is_zero |= is_zero >> 2;
62 is_zero |= is_zero >> 1;
70 // For is_zero and is_p, the LSB is 0 iff all the bits are zero.
72 is_zero = (~is_zero) << 31;
73 is_zero = static_cast<int32_t>(is_zero) >> 31;
77 // Add computes *out = a+b
79 // a[i] + b[i] < 2**32
80 void Add(FieldElement* out, const FieldElement& a, const FieldElement& b) {
81 for (int i = 0; i < 8; i++) {
82 (*out)[i] = a[i] + b[i];
86 static const uint32_t kTwo31p3 = (1u << 31) + (1u << 3);
87 static const uint32_t kTwo31m3 = (1u << 31) - (1u << 3);
88 static const uint32_t kTwo31m15m3 = (1u << 31) - (1u << 15) - (1u << 3);
89 // kZero31ModP is 0 mod p where bit 31 is set in all limbs so that we can
90 // subtract smaller amounts without underflow. See the section "Subtraction" in
92 static const FieldElement kZero31ModP = {
93 kTwo31p3, kTwo31m3, kTwo31m3, kTwo31m15m3,
94 kTwo31m3, kTwo31m3, kTwo31m3, kTwo31m3
97 // Subtract computes *out = a-b
101 void Subtract(FieldElement* out, const FieldElement& a, const FieldElement& b) {
102 for (int i = 0; i < 8; i++) {
103 // See the section on "Subtraction" in [1] for details.
104 (*out)[i] = a[i] + kZero31ModP[i] - b[i];
108 static const uint64_t kTwo63p35 = (1ull << 63) + (1ull << 35);
109 static const uint64_t kTwo63m35 = (1ull << 63) - (1ull << 35);
110 static const uint64_t kTwo63m35m19 = (1ull << 63) - (1ull << 35) - (1ull << 19);
111 // kZero63ModP is 0 mod p where bit 63 is set in all limbs. See the section
112 // "Subtraction" in [1] for why.
113 static const uint64_t kZero63ModP[8] = {
114 kTwo63p35, kTwo63m35, kTwo63m35, kTwo63m35,
115 kTwo63m35m19, kTwo63m35, kTwo63m35, kTwo63m35,
118 static const uint32_t kBottom28Bits = 0xfffffff;
120 // LargeFieldElement also represents an element of the field. The limbs are
121 // still spaced 28-bits apart and in little-endian order. So the limbs are at
122 // 0, 28, 56, ..., 392 bits, each 64-bits wide.
123 typedef uint64_t LargeFieldElement[15];
125 // ReduceLarge converts a LargeFieldElement to a FieldElement.
128 void ReduceLarge(FieldElement* out, LargeFieldElement* inptr) {
129 LargeFieldElement& in(*inptr);
131 for (int i = 0; i < 8; i++) {
132 in[i] += kZero63ModP[i];
135 // Eliminate the coefficients at 2**224 and greater while maintaining the
137 for (int i = 14; i >= 8; i--) {
138 in[i-8] -= in[i]; // reflection off the "+1" term of p.
139 in[i-5] += (in[i] & 0xffff) << 12; // part of the "-2**96" reflection.
140 in[i-4] += in[i] >> 16; // the rest of the "-2**96" reflection.
145 // As the values become small enough, we start to store them in |out| and use
146 // 32-bit operations.
147 for (int i = 1; i < 8; i++) {
148 in[i+1] += in[i] >> 28;
149 (*out)[i] = static_cast<uint32_t>(in[i] & kBottom28Bits);
151 // Eliminate the term at 2*224 that we introduced while keeping the same
153 in[0] -= in[8]; // reflection off the "+1" term of p.
154 (*out)[3] += static_cast<uint32_t>(in[8] & 0xffff) << 12; // "-2**96" term
155 (*out)[4] += static_cast<uint32_t>(in[8] >> 16); // rest of "-2**96" term
159 // out[1,2,5..7] < 2**28
161 (*out)[0] = static_cast<uint32_t>(in[0] & kBottom28Bits);
162 (*out)[1] += static_cast<uint32_t>((in[0] >> 28) & kBottom28Bits);
163 (*out)[2] += static_cast<uint32_t>(in[0] >> 56);
169 // Mul computes *out = a*b
171 // a[i] < 2**29, b[i] < 2**30 (or vice versa)
173 void Mul(FieldElement* out, const FieldElement& a, const FieldElement& b) {
174 LargeFieldElement tmp;
175 memset(&tmp, 0, sizeof(tmp));
177 for (int i = 0; i < 8; i++) {
178 for (int j = 0; j < 8; j++) {
179 tmp[i + j] += static_cast<uint64_t>(a[i]) * static_cast<uint64_t>(b[j]);
183 ReduceLarge(out, &tmp);
186 // Square computes *out = a*a
190 void Square(FieldElement* out, const FieldElement& a) {
191 LargeFieldElement tmp;
192 memset(&tmp, 0, sizeof(tmp));
194 for (int i = 0; i < 8; i++) {
195 for (int j = 0; j <= i; j++) {
196 uint64_t r = static_cast<uint64_t>(a[i]) * static_cast<uint64_t>(a[j]);
205 ReduceLarge(out, &tmp);
208 // Reduce reduces the coefficients of in_out to smaller bounds.
210 // On entry: a[i] < 2**31 + 2**30
211 // On exit: a[i] < 2**29
212 void Reduce(FieldElement* in_out) {
213 FieldElement& a = *in_out;
215 for (int i = 0; i < 7; i++) {
216 a[i+1] += a[i] >> 28;
217 a[i] &= kBottom28Bits;
219 uint32_t top = a[7] >> 28;
220 a[7] &= kBottom28Bits;
223 // Constant-time: mask = (top != 0) ? 0xffffffff : 0
228 mask = static_cast<uint32_t>(static_cast<int32_t>(mask) >> 31);
230 // Eliminate top while maintaining the same value mod p.
234 // We may have just made a[0] negative but, if we did, then we must
235 // have added something to a[3], thus it's > 2**12. Therefore we can
236 // carry down to a[0].
238 a[2] += mask & ((1<<28) - 1);
239 a[1] += mask & ((1<<28) - 1);
240 a[0] += mask & (1<<28);
243 // Invert calcuates *out = in**-1 by computing in**(2**224 - 2**96 - 1), i.e.
244 // Fermat's little theorem.
245 void Invert(FieldElement* out, const FieldElement& in) {
246 FieldElement f1, f2, f3, f4;
248 Square(&f1, in); // 2
249 Mul(&f1, f1, in); // 2**2 - 1
250 Square(&f1, f1); // 2**3 - 2
251 Mul(&f1, f1, in); // 2**3 - 1
252 Square(&f2, f1); // 2**4 - 2
253 Square(&f2, f2); // 2**5 - 4
254 Square(&f2, f2); // 2**6 - 8
255 Mul(&f1, f1, f2); // 2**6 - 1
256 Square(&f2, f1); // 2**7 - 2
257 for (int i = 0; i < 5; i++) { // 2**12 - 2**6
260 Mul(&f2, f2, f1); // 2**12 - 1
261 Square(&f3, f2); // 2**13 - 2
262 for (int i = 0; i < 11; i++) { // 2**24 - 2**12
265 Mul(&f2, f3, f2); // 2**24 - 1
266 Square(&f3, f2); // 2**25 - 2
267 for (int i = 0; i < 23; i++) { // 2**48 - 2**24
270 Mul(&f3, f3, f2); // 2**48 - 1
271 Square(&f4, f3); // 2**49 - 2
272 for (int i = 0; i < 47; i++) { // 2**96 - 2**48
275 Mul(&f3, f3, f4); // 2**96 - 1
276 Square(&f4, f3); // 2**97 - 2
277 for (int i = 0; i < 23; i++) { // 2**120 - 2**24
280 Mul(&f2, f4, f2); // 2**120 - 1
281 for (int i = 0; i < 6; i++) { // 2**126 - 2**6
284 Mul(&f1, f1, f2); // 2**126 - 1
285 Square(&f1, f1); // 2**127 - 2
286 Mul(&f1, f1, in); // 2**127 - 1
287 for (int i = 0; i < 97; i++) { // 2**224 - 2**97
290 Mul(out, f1, f3); // 2**224 - 2**96 - 1
293 // Contract converts a FieldElement to its minimal, distinguished form.
295 // On entry, in[i] < 2**29
296 // On exit, in[i] < 2**28
297 void Contract(FieldElement* inout) {
298 FieldElement& out = *inout;
300 // Reduce the coefficients to < 2**28.
301 for (int i = 0; i < 7; i++) {
302 out[i+1] += out[i] >> 28;
303 out[i] &= kBottom28Bits;
305 uint32_t top = out[7] >> 28;
306 out[7] &= kBottom28Bits;
308 // Eliminate top while maintaining the same value mod p.
312 // We may just have made out[0] negative. So we carry down. If we made
313 // out[0] negative then we know that out[3] is sufficiently positive
314 // because we just added to it.
315 for (int i = 0; i < 3; i++) {
316 uint32_t mask = static_cast<uint32_t>(static_cast<int32_t>(out[i]) >> 31);
317 out[i] += (1 << 28) & mask;
318 out[i+1] -= 1 & mask;
321 // We might have pushed out[3] over 2**28 so we perform another, partial
323 for (int i = 3; i < 7; i++) {
324 out[i+1] += out[i] >> 28;
325 out[i] &= kBottom28Bits;
328 out[7] &= kBottom28Bits;
330 // Eliminate top while maintaining the same value mod p.
334 // There are two cases to consider for out[3]:
335 // 1) The first time that we eliminated top, we didn't push out[3] over
336 // 2**28. In this case, the partial carry chain didn't change any values
338 // 2) We did push out[3] over 2**28 the first time that we eliminated top.
339 // The first value of top was in [0..16), therefore, prior to eliminating
340 // the first top, 0xfff1000 <= out[3] <= 0xfffffff. Therefore, after
341 // overflowing and being reduced by the second carry chain, out[3] <=
342 // 0xf000. Thus it cannot have overflowed when we eliminated top for the
345 // Again, we may just have made out[0] negative, so do the same carry down.
346 // As before, if we made out[0] negative then we know that out[3] is
347 // sufficiently positive.
348 for (int i = 0; i < 3; i++) {
349 uint32_t mask = static_cast<uint32_t>(static_cast<int32_t>(out[i]) >> 31);
350 out[i] += (1 << 28) & mask;
351 out[i+1] -= 1 & mask;
354 // The value is < 2**224, but maybe greater than p. In order to reduce to a
355 // unique, minimal value we see if the value is >= p and, if so, subtract p.
357 // First we build a mask from the top four limbs, which must all be
358 // equal to bottom28Bits if the whole value is >= p. If top_4_all_ones
359 // ends up with any zero bits in the bottom 28 bits, then this wasn't
361 uint32_t top_4_all_ones = 0xffffffffu;
362 for (int i = 4; i < 8; i++) {
363 top_4_all_ones &= out[i];
365 top_4_all_ones |= 0xf0000000;
366 // Now we replicate any zero bits to all the bits in top_4_all_ones.
367 top_4_all_ones &= top_4_all_ones >> 16;
368 top_4_all_ones &= top_4_all_ones >> 8;
369 top_4_all_ones &= top_4_all_ones >> 4;
370 top_4_all_ones &= top_4_all_ones >> 2;
371 top_4_all_ones &= top_4_all_ones >> 1;
373 static_cast<uint32_t>(static_cast<int32_t>(top_4_all_ones << 31) >> 31);
375 // Now we test whether the bottom three limbs are non-zero.
376 uint32_t bottom_3_non_zero = out[0] | out[1] | out[2];
377 bottom_3_non_zero |= bottom_3_non_zero >> 16;
378 bottom_3_non_zero |= bottom_3_non_zero >> 8;
379 bottom_3_non_zero |= bottom_3_non_zero >> 4;
380 bottom_3_non_zero |= bottom_3_non_zero >> 2;
381 bottom_3_non_zero |= bottom_3_non_zero >> 1;
383 static_cast<uint32_t>(static_cast<int32_t>(bottom_3_non_zero) >> 31);
385 // Everything depends on the value of out[3].
386 // If it's > 0xffff000 and top_4_all_ones != 0 then the whole value is >= p
387 // If it's = 0xffff000 and top_4_all_ones != 0 and bottom_3_non_zero != 0,
388 // then the whole value is >= p
389 // If it's < 0xffff000, then the whole value is < p
390 uint32_t n = out[3] - 0xffff000;
391 uint32_t out_3_equal = n;
392 out_3_equal |= out_3_equal >> 16;
393 out_3_equal |= out_3_equal >> 8;
394 out_3_equal |= out_3_equal >> 4;
395 out_3_equal |= out_3_equal >> 2;
396 out_3_equal |= out_3_equal >> 1;
398 ~static_cast<uint32_t>(static_cast<int32_t>(out_3_equal << 31) >> 31);
400 // If out[3] > 0xffff000 then n's MSB will be zero.
402 ~static_cast<uint32_t>(static_cast<int32_t>(n << 31) >> 31);
405 top_4_all_ones & ((out_3_equal & bottom_3_non_zero) | out_3_gt);
407 out[3] -= 0xffff000 & mask;
408 out[4] -= 0xfffffff & mask;
409 out[5] -= 0xfffffff & mask;
410 out[6] -= 0xfffffff & mask;
411 out[7] -= 0xfffffff & mask;
415 // Group element functions.
417 // These functions deal with group elements. The group is an elliptic curve
418 // group with a = -3 defined in FIPS 186-3, section D.2.2.
420 using crypto::p224::Point;
422 // kB is parameter of the elliptic curve.
423 const FieldElement kB = {
424 55967668, 11768882, 265861671, 185302395,
425 39211076, 180311059, 84673715, 188764328,
428 void CopyConditional(Point* out, const Point& a, uint32_t mask);
429 void DoubleJacobian(Point* out, const Point& a);
431 // AddJacobian computes *out = a+b where a != b.
432 void AddJacobian(Point *out,
435 // See http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#addition-add-2007-bl
436 FieldElement z1z1, z2z2, u1, u2, s1, s2, h, i, j, r, v;
438 uint32_t z1_is_zero = IsZero(a.z);
439 uint32_t z2_is_zero = IsZero(b.z);
462 Subtract(&h, u2, u1);
464 uint32_t x_equal = IsZero(h);
467 for (int k = 0; k < 8; k++) {
476 Subtract(&r, s2, s1);
478 uint32_t y_equal = IsZero(r);
480 if (x_equal && y_equal && !z1_is_zero && !z2_is_zero) {
481 // The two input points are the same therefore we must use the dedicated
482 // doubling function as the slope of the line is undefined.
483 DoubleJacobian(out, a);
487 for (int k = 0; k < 8; k++) {
495 // Z3 = ((Z1+Z2)²-Z1Z1-Z2Z2)*H
496 Add(&z1z1, z1z1, z2z2);
497 Add(&z2z2, a.z, b.z);
500 Subtract(&out->z, z2z2, z1z1);
502 Mul(&out->z, out->z, h);
505 for (int k = 0; k < 8; k++) {
511 Subtract(&out->x, out->x, z1z1);
514 // Y3 = r*(V-X3)-2*S1*J
515 for (int k = 0; k < 8; k++) {
519 Subtract(&z1z1, v, out->x);
522 Subtract(&out->y, z1z1, s1);
525 CopyConditional(out, a, z2_is_zero);
526 CopyConditional(out, b, z1_is_zero);
529 // DoubleJacobian computes *out = a+a.
530 void DoubleJacobian(Point* out, const Point& a) {
531 // See http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#doubling-dbl-2001-b
532 FieldElement delta, gamma, beta, alpha, t;
536 Mul(&beta, a.x, gamma);
538 // alpha = 3*(X1-delta)*(X1+delta)
540 for (int i = 0; i < 8; i++) {
544 Subtract(&alpha, a.x, delta);
546 Mul(&alpha, alpha, t);
548 // Z3 = (Y1+Z1)²-gamma-delta
549 Add(&out->z, a.y, a.z);
551 Square(&out->z, out->z);
552 Subtract(&out->z, out->z, gamma);
554 Subtract(&out->z, out->z, delta);
557 // X3 = alpha²-8*beta
558 for (int i = 0; i < 8; i++) {
559 delta[i] = beta[i] << 3;
562 Square(&out->x, alpha);
563 Subtract(&out->x, out->x, delta);
566 // Y3 = alpha*(4*beta-X3)-8*gamma²
567 for (int i = 0; i < 8; i++) {
571 Subtract(&beta, beta, out->x);
573 Square(&gamma, gamma);
574 for (int i = 0; i < 8; i++) {
578 Mul(&out->y, alpha, beta);
579 Subtract(&out->y, out->y, gamma);
583 // CopyConditional sets *out=a if mask is 0xffffffff. mask must be either 0 of
585 void CopyConditional(Point* out, const Point& a, uint32_t mask) {
586 for (int i = 0; i < 8; i++) {
587 out->x[i] ^= mask & (a.x[i] ^ out->x[i]);
588 out->y[i] ^= mask & (a.y[i] ^ out->y[i]);
589 out->z[i] ^= mask & (a.z[i] ^ out->z[i]);
593 // ScalarMult calculates *out = a*scalar where scalar is a big-endian number of
594 // length scalar_len and != 0.
595 void ScalarMult(Point* out,
597 const uint8_t* scalar,
599 memset(out, 0, sizeof(*out));
602 for (size_t i = 0; i < scalar_len; i++) {
603 for (unsigned int bit_num = 0; bit_num < 8; bit_num++) {
604 DoubleJacobian(out, *out);
605 uint32_t bit = static_cast<uint32_t>(static_cast<int32_t>(
606 (((scalar[i] >> (7 - bit_num)) & 1) << 31) >> 31));
607 AddJacobian(&tmp, a, *out);
608 CopyConditional(out, tmp, bit);
613 // Get224Bits reads 7 words from in and scatters their contents in
614 // little-endian form into 8 words at out, 28 bits per output word.
615 void Get224Bits(uint32_t* out, const uint32_t* in) {
616 out[0] = NetToHost32(in[6]) & kBottom28Bits;
617 out[1] = ((NetToHost32(in[5]) << 4) |
618 (NetToHost32(in[6]) >> 28)) & kBottom28Bits;
619 out[2] = ((NetToHost32(in[4]) << 8) |
620 (NetToHost32(in[5]) >> 24)) & kBottom28Bits;
621 out[3] = ((NetToHost32(in[3]) << 12) |
622 (NetToHost32(in[4]) >> 20)) & kBottom28Bits;
623 out[4] = ((NetToHost32(in[2]) << 16) |
624 (NetToHost32(in[3]) >> 16)) & kBottom28Bits;
625 out[5] = ((NetToHost32(in[1]) << 20) |
626 (NetToHost32(in[2]) >> 12)) & kBottom28Bits;
627 out[6] = ((NetToHost32(in[0]) << 24) |
628 (NetToHost32(in[1]) >> 8)) & kBottom28Bits;
629 out[7] = (NetToHost32(in[0]) >> 4) & kBottom28Bits;
632 // Put224Bits performs the inverse operation to Get224Bits: taking 28 bits from
633 // each of 8 input words and writing them in big-endian order to 7 words at
635 void Put224Bits(uint32_t* out, const uint32_t* in) {
636 out[6] = HostToNet32((in[0] >> 0) | (in[1] << 28));
637 out[5] = HostToNet32((in[1] >> 4) | (in[2] << 24));
638 out[4] = HostToNet32((in[2] >> 8) | (in[3] << 20));
639 out[3] = HostToNet32((in[3] >> 12) | (in[4] << 16));
640 out[2] = HostToNet32((in[4] >> 16) | (in[5] << 12));
641 out[1] = HostToNet32((in[5] >> 20) | (in[6] << 8));
642 out[0] = HostToNet32((in[6] >> 24) | (in[7] << 4));
645 } // anonymous namespace
651 bool Point::SetFromString(base::StringPiece in) {
652 if (in.size() != 2*28)
654 const uint32_t* inwords = reinterpret_cast<const uint32_t*>(in.data());
655 Get224Bits(x, inwords);
656 Get224Bits(y, inwords + 7);
657 memset(&z, 0, sizeof(z));
660 // Check that the point is on the curve, i.e. that y² = x³ - 3x + b.
669 FieldElement three_x;
670 for (int i = 0; i < 8; i++) {
671 three_x[i] = x[i] * 3;
674 Subtract(&rhs, rhs, three_x);
677 ::Add(&rhs, rhs, kB);
679 return memcmp(&lhs, &rhs, sizeof(lhs)) == 0;
682 std::string Point::ToString() const {
683 FieldElement zinv, zinv_sq, xx, yy;
685 // If this is the point at infinity we return a string of all zeros.
686 if (IsZero(this->z)) {
687 static const char zeros[56] = {0};
688 return std::string(zeros, sizeof(zeros));
691 Invert(&zinv, this->z);
692 Square(&zinv_sq, zinv);
693 Mul(&xx, x, zinv_sq);
694 Mul(&zinv_sq, zinv_sq, zinv);
695 Mul(&yy, y, zinv_sq);
700 uint32_t outwords[14];
701 Put224Bits(outwords, xx);
702 Put224Bits(outwords + 7, yy);
703 return std::string(reinterpret_cast<const char*>(outwords), sizeof(outwords));
706 void ScalarMult(const Point& in, const uint8_t* scalar, Point* out) {
707 ::ScalarMult(out, in, scalar, 28);
710 // kBasePoint is the base point (generator) of the elliptic curve group.
711 static const Point kBasePoint = {
712 {22813985, 52956513, 34677300, 203240812,
713 12143107, 133374265, 225162431, 191946955},
714 {83918388, 223877528, 122119236, 123340192,
715 266784067, 263504429, 146143011, 198407736},
716 {1, 0, 0, 0, 0, 0, 0, 0},
719 void ScalarBaseMult(const uint8_t* scalar, Point* out) {
720 ::ScalarMult(out, kBasePoint, scalar, 28);
723 void Add(const Point& a, const Point& b, Point* out) {
724 AddJacobian(out, a, b);
727 void Negate(const Point& in, Point* out) {
728 // Guide to elliptic curve cryptography, page 89 suggests that (X : X+Y : Z)
729 // is the negative in Jacobian coordinates, but it doesn't actually appear to
730 // be true in testing so this performs the negation in affine coordinates.
731 FieldElement zinv, zinv_sq, y;
733 Square(&zinv_sq, zinv);
734 Mul(&out->x, in.x, zinv_sq);
735 Mul(&zinv_sq, zinv_sq, zinv);
736 Mul(&y, in.y, zinv_sq);
738 Subtract(&out->y, kP, y);
741 memset(&out->z, 0, sizeof(out->z));
747 } // namespace crypto