crypto: special case ∞+a, a+∞ and a+a in p224.

crypto/p224.cc

 // Copyright (c) 2012 The Chromium Authors. All rights reserved.
 // Use of this source code is governed by a BSD-style license that can be
 // found in the LICENSE file.
 
 // This is an implementation of the P224 elliptic curve group. It's written to
 // be short and simple rather than fast, although it's still constant-time.
 //
 // See http://www.imperialviolet.org/2010/12/04/ecc.html ([1]) for background.
 
 #include "crypto/p224.h"
@@ skipping 14 matching lines @@
 // Field elements are represented by a FieldElement, which is a typedef to an
 // array of 8 uint32's. The value of a FieldElement, a, is:
 //   a[0] + 2**28·a[1] + 2**56·a[1] + ... + 2**196·a[7]
 //
 // Using 28-bit limbs means that there's only 4 bits of headroom, which is less
 // than we would really like. But it has the useful feature that we hit 2**224
 // exactly, making the reflections during a reduce much nicer.
 
 using crypto::p224::FieldElement;
 
+// kP is the P224 prime.
+const FieldElement kP = {
+  1, 0, 0, 268431360,
+  268435455, 268435455, 268435455, 268435455,
+};
+
+void Contract(FieldElement* inout);
+
+// IsZero returns 0xffffffff if a == 0 mod p and 0 otherwise.
+uint32_t IsZero(const FieldElement& a) {
+  FieldElement minimal;
+  memcpy(&minimal, &a, sizeof(minimal));
+  Contract(&minimal);
+
+  uint32 is_zero = 0, is_p = 0;
+  for (unsigned i = 0; i < 8; i++) {
+    is_zero |= minimal[i];
+    is_p |= minimal[i] - kP[i];
+  }
+
+  // If either is_zero or is_p is 0, then we should return 1.
+  is_zero |= is_zero >> 16;
+  is_zero |= is_zero >> 8;
+  is_zero |= is_zero >> 4;
+  is_zero |= is_zero >> 2;
+  is_zero |= is_zero >> 1;
+
+  is_p |= is_p >> 16;
+  is_p |= is_p >> 8;
+  is_p |= is_p >> 4;
+  is_p |= is_p >> 2;
+  is_p |= is_p >> 1;
+
+  // For is_zero and is_p, the LSB is 0 iff all the bits are zero.
+  is_zero &= is_p & 1;
+  is_zero = (~is_zero) << 31;
+  is_zero = static_cast<int32>(is_zero) >> 31;
+  return is_zero;
+}
+
 // Add computes *out = a+b
 //
 // a[i] + b[i] < 2**32
 void Add(FieldElement* out, const FieldElement& a, const FieldElement& b) {
   for (int i = 0; i < 8; i++) {
     (*out)[i] = a[i] + b[i];
   }
 }
 
 static const uint32 kTwo31p3 = (1u<<31) + (1u<<3);
@@ skipping 264 matching lines @@
     out[i+1] -= 1 & mask;
   }
 
   // The value is < 2**224, but maybe greater than p. In order to reduce to a
   // unique, minimal value we see if the value is >= p and, if so, subtract p.
 
   // First we build a mask from the top four limbs, which must all be
   // equal to bottom28Bits if the whole value is >= p. If top4AllOnes
   // ends up with any zero bits in the bottom 28 bits, then this wasn't
   // true.
   uint32 top4AllOnes = 0xffffffffu;

[willchan no longer on Chromium, 2012/07/25 23:35:26]

I know you're going to hate me but...
nit: Top4AllOnes is proper camel casing naming style in Chromium.

[agl, 2012/07/27 17:21:06]

No, that's my screw up. I debugged this code in Go and translated across. I've
changed to top_4_all_ones as its a local variable.

   for (int i = 4; i < 8; i++) {
-    top4AllOnes &= (out[i] & kBottom28Bits) - 1;
+    top4AllOnes &= out[i];
   }
   top4AllOnes |= 0xf0000000;
   // Now we replicate any zero bits to all the bits in top4AllOnes.
   top4AllOnes &= top4AllOnes >> 16;
   top4AllOnes &= top4AllOnes >> 8;
   top4AllOnes &= top4AllOnes >> 4;
   top4AllOnes &= top4AllOnes >> 2;
   top4AllOnes &= top4AllOnes >> 1;
   top4AllOnes =
       static_cast<uint32>(static_cast<int32>(top4AllOnes << 31) >> 31);
 
   // Now we test whether the bottom three limbs are non-zero.
   uint32 bottom3NonZero = out[0] | out[1] | out[2];

[willchan no longer on Chromium, 2012/07/25 23:35:26]

Ditto here on CamelCase style.

[agl, 2012/07/27 17:21:06]

Done.

   bottom3NonZero |= bottom3NonZero >> 16;
   bottom3NonZero |= bottom3NonZero >> 8;
   bottom3NonZero |= bottom3NonZero >> 4;
   bottom3NonZero |= bottom3NonZero >> 2;
   bottom3NonZero |= bottom3NonZero >> 1;
   bottom3NonZero =
-      static_cast<uint32>(static_cast<int32>(bottom3NonZero << 31) >> 31);
+      static_cast<uint32>(static_cast<int32>(bottom3NonZero) >> 31);
 
   // Everything depends on the value of out[3].
   //    If it's > 0xffff000 and top4AllOnes != 0 then the whole value is >= p
   //    If it's = 0xffff000 and top4AllOnes != 0 and bottom3NonZero != 0,
   //      then the whole value is >= p
   //    If it's < 0xffff000, then the whole value is < p
   uint32 n = out[3] - 0xffff000;
   uint32 out3Equal = n;
   out3Equal |= out3Equal >> 16;
   out3Equal |= out3Equal >> 8;
@@ skipping 16 matching lines @@
 }
 
 
 // Group element functions.
 //
 // These functions deal with group elements. The group is an elliptic curve
 // group with a = -3 defined in FIPS 186-3, section D.2.2.
 
 using crypto::p224::Point;
 
-// kP is the P224 prime.
-const FieldElement kP = {
-  1, 0, 0, 268431360,
-  268435455, 268435455, 268435455, 268435455,
-};
-
 // kB is parameter of the elliptic curve.
 const FieldElement kB = {
   55967668, 11768882, 265861671, 185302395,
   39211076, 180311059, 84673715, 188764328,
 };
 
+void CopyConditional(Point* out, const Point& a, uint32 mask);

[willchan no longer on Chromium, 2012/07/25 23:35:26]

Google style is to put output parameters last. That said, you seem to have
flouted this Google style rule in the whole file, so it's best that you stay
consistently inconsistent :)

[agl, 2012/07/27 17:21:06]

Well, I'm just being consistent with all the memcpys and every other system
function :)

+void DoubleJacobian(Point* out, const Point& a);
+
 // AddJacobian computes *out = a+b where a != b.
 void AddJacobian(Point *out,
                  const Point& a,
                  const Point& b) {
   // See http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#addition-add-2007-bl
   FieldElement z1z1, z2z2, u1, u2, s1, s2, h, i, j, r, v;
 
+  uint32 z1_is_zero = IsZero(a.z);
+  uint32 z2_is_zero = IsZero(b.z);
+
   // Z1Z1 = Z1²
   Square(&z1z1, a.z);
 
   // Z2Z2 = Z2²
   Square(&z2z2, b.z);
 
   // U1 = X1*Z2Z2
   Mul(&u1, a.x, z2z2);
 
   // U2 = X2*Z1Z1
   Mul(&u2, b.x, z1z1);
 
   // S1 = Y1*Z2*Z2Z2
   Mul(&s1, b.z, z2z2);
   Mul(&s1, a.y, s1);
 
   // S2 = Y2*Z1*Z1Z1
   Mul(&s2, a.z, z1z1);
   Mul(&s2, b.y, s2);
 
   // H = U2-U1
   Subtract(&h, u2, u1);
   Reduce(&h);
+  uint32 x_equal = IsZero(h);
 
   // I = (2*H)²
   for (int j = 0; j < 8; j++) {
     i[j] = h[j] << 1;
   }
   Reduce(&i);
   Square(&i, i);
 
   // J = H*I
   Mul(&j, h, i);
   // r = 2*(S2-S1)
   Subtract(&r, s2, s1);
   Reduce(&r);
+  uint32 y_equal = IsZero(r);
+
+  if (x_equal && y_equal && !z1_is_zero && !z2_is_zero) {
+    // The two input points are the same therefore we must use the dedicated
+    // doubling function as the slope of the line is undefined.
+    DoubleJacobian(out, a);
+    return;
+  }
+
   for (int i = 0; i < 8; i++) {
     r[i] <<= 1;
   }
   Reduce(&r);
 
   // V = U1*I
   Mul(&v, u1, i);
 
   // Z3 = ((Z1+Z2)²-Z1Z1-Z2Z2)*H
   Add(&z1z1, z1z1, z2z2);
@@ skipping 17 matching lines @@
   // Y3 = r*(V-X3)-2*S1*J
   for (int i = 0; i < 8; i++) {
     s1[i] <<= 1;
   }
   Mul(&s1, s1, j);
   Subtract(&z1z1, v, out->x);
   Reduce(&z1z1);
   Mul(&z1z1, z1z1, r);
   Subtract(&out->y, z1z1, s1);
   Reduce(&out->y);
+
+  CopyConditional(out, a, z2_is_zero);
+  CopyConditional(out, b, z1_is_zero);
 }
 
 // DoubleJacobian computes *out = a+a.
 void DoubleJacobian(Point* out, const Point& a) {
   // See http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#doubling-dbl-2001-b
   FieldElement delta, gamma, beta, alpha, t;
 
   Square(&delta, a.z);
   Square(&gamma, a.y);
   Mul(&beta, a.x, gamma);
@@ skipping 55 matching lines @@
   }
 }
 
 // ScalarMult calculates *out = a*scalar where scalar is a big-endian number of
 // length scalar_len and != 0.
 void ScalarMult(Point* out, const Point& a,
                 const uint8* scalar, size_t scalar_len) {
   memset(out, 0, sizeof(*out));
   Point tmp;
 
-  uint32 first_bit = 0xffffffff;
   for (size_t i = 0; i < scalar_len; i++) {
     for (unsigned int bit_num = 0; bit_num < 8; bit_num++) {
       DoubleJacobian(out, *out);
       uint32 bit = static_cast<uint32>(static_cast<int32>(
           (((scalar[i] >> (7 - bit_num)) & 1) << 31) >> 31));
       AddJacobian(&tmp, a, *out);
-      CopyConditional(out, a, first_bit & bit);
-      CopyConditional(out, tmp, ~first_bit & bit);
-      first_bit = first_bit & ~bit;
+      CopyConditional(out, tmp, bit);
     }
   }
 }
 
 // Get224Bits reads 7 words from in and scatters their contents in
 // little-endian form into 8 words at out, 28 bits per output word.
 void Get224Bits(uint32* out, const uint32* in) {
   out[0] = NetToHost32(in[6]) & kBottom28Bits;
   out[1] = ((NetToHost32(in[5]) << 4) |
             (NetToHost32(in[6]) >> 28)) & kBottom28Bits;
@@ skipping 56 matching lines @@
   Reduce(&rhs);
 
   ::Add(&rhs, rhs, kB);
   Contract(&rhs);
   return memcmp(&lhs, &rhs, sizeof(lhs)) == 0;
 }
 
 std::string Point::ToString() const {
   FieldElement zinv, zinv_sq, x, y;
 
+  // If this is the point at infinity we return a string of all zeros.
+  if (IsZero(this->z)) {
+    char zeros[56];

[willchan no longer on Chromium, 2012/07/25 23:35:26]

Why isn't this just a static const? It'd POD, so the linker should initialize it
to 0. That said, I guess this is only used by tests, so it's not perf critical.
So whatever.

[agl, 2012/07/27 17:21:06]

Good point. Done.

+    memset(zeros, 0, sizeof(zeros));
+    return std::string(zeros, sizeof(zeros));
+  }
+
   Invert(&zinv, this->z);
   Square(&zinv_sq, zinv);
   Mul(&x, this->x, zinv_sq);
   Mul(&zinv_sq, zinv_sq, zinv);
   Mul(&y, this->y, zinv_sq);
 
   Contract(&x);
   Contract(&y);
 
   uint32 outwords[14];
@@ skipping 37 matching lines @@
   Subtract(&out->y, kP, y);
   Reduce(&out->y);
 
   memset(&out->z, 0, sizeof(out->z));
   out->z[0] = 1;
 }
 
 }  // namespace p224
 
 }  // namespace crypto