src/crypto/p224.cc

   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.
   4
   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.
   7 //
   8 // See http://www.imperialviolet.org/2010/12/04/ecc.html ([1]) for background.
   9
  10 #include "crypto/p224.h"
  11
  12 #include <string.h>
  13
  14 #include "base/sys_byteorder.h"
  15
  16 namespace {
  17
  18 using base::HostToNet32;
  19 using base::NetToHost32;
  20
  21 // Field element functions.
  22 //
  23 // The field that we're dealing with is ℤ/pℤ where p = 2**224 - 2**96 + 1.
  24 //
  25 // Field elements are represented by a FieldElement, which is a typedef to an
  26 // array of 8 uint32's. The value of a FieldElement, a, is:
  27 //   a[0] + 2**28·a[1] + 2**56·a[1] + ... + 2**196·a[7]
  28 //
  29 // Using 28-bit limbs means that there's only 4 bits of headroom, which is less
  30 // than we would really like. But it has the useful feature that we hit 2**224
  31 // exactly, making the reflections during a reduce much nicer.
  32
  33 using crypto::p224::FieldElement;
  34
  35 // kP is the P224 prime.
  36 const FieldElement kP = {
  37   1, 0, 0, 268431360,
  38   268435455, 268435455, 268435455, 268435455,
  39 };
  40
  41 void Contract(FieldElement* inout);
  42
  43 // IsZero returns 0xffffffff if a == 0 mod p and 0 otherwise.
  44 uint32 IsZero(const FieldElement& a) {
  45   FieldElement minimal;
  46   memcpy(&minimal, &a, sizeof(minimal));
  47   Contract(&minimal);
  48
  49   uint32 is_zero = 0, is_p = 0;
  50   for (unsigned i = 0; i < 8; i++) {
  51     is_zero |= minimal[i];
  52     is_p |= minimal[i] - kP[i];
  53   }
  54
  55   // If either is_zero or is_p is 0, then we should return 1.
  56   is_zero |= is_zero >> 16;
  57   is_zero |= is_zero >> 8;
  58   is_zero |= is_zero >> 4;
  59   is_zero |= is_zero >> 2;
  60   is_zero |= is_zero >> 1;
  61
  62   is_p |= is_p >> 16;
  63   is_p |= is_p >> 8;
  64   is_p |= is_p >> 4;
  65   is_p |= is_p >> 2;
  66   is_p |= is_p >> 1;
  67
  68   // For is_zero and is_p, the LSB is 0 iff all the bits are zero.
  69   is_zero &= is_p & 1;
  70   is_zero = (~is_zero) << 31;
  71   is_zero = static_cast<int32>(is_zero) >> 31;
  72   return is_zero;
  73 }
  74
  75 // Add computes *out = a+b
  76 //
  77 // a[i] + b[i] < 2**32
  78 void Add(FieldElement* out, const FieldElement& a, const FieldElement& b) {
  79   for (int i = 0; i < 8; i++) {
  80     (*out)[i] = a[i] + b[i];
  81   }
  82 }
  83
  84 static const uint32 kTwo31p3 = (1u<<31) + (1u<<3);
  85 static const uint32 kTwo31m3 = (1u<<31) - (1u<<3);
  86 static const uint32 kTwo31m15m3 = (1u<<31) - (1u<<15) - (1u<<3);
  87 // kZero31ModP is 0 mod p where bit 31 is set in all limbs so that we can
  88 // subtract smaller amounts without underflow. See the section "Subtraction" in
  89 // [1] for why.
  90 static const FieldElement kZero31ModP = {
  91   kTwo31p3, kTwo31m3, kTwo31m3, kTwo31m15m3,
  92   kTwo31m3, kTwo31m3, kTwo31m3, kTwo31m3
  93 };
  94
  95 // Subtract computes *out = a-b
  96 //
  97 // a[i], b[i] < 2**30
  98 // out[i] < 2**32
  99 void Subtract(FieldElement* out, const FieldElement& a, const FieldElement& b) {
 100   for (int i = 0; i < 8; i++) {
 101     // See the section on "Subtraction" in [1] for details.
 102     (*out)[i] = a[i] + kZero31ModP[i] - b[i];
 103   }
 104 }
 105
 106 static const uint64 kTwo63p35 = (1ull<<63) + (1ull<<35);
 107 static const uint64 kTwo63m35 = (1ull<<63) - (1ull<<35);
 108 static const uint64 kTwo63m35m19 = (1ull<<63) - (1ull<<35) - (1ull<<19);
 109 // kZero63ModP is 0 mod p where bit 63 is set in all limbs. See the section
 110 // "Subtraction" in [1] for why.
 111 static const uint64 kZero63ModP[8] = {
 112   kTwo63p35, kTwo63m35, kTwo63m35, kTwo63m35,
 113   kTwo63m35m19, kTwo63m35, kTwo63m35, kTwo63m35,
 114 };
 115
 116 static const uint32 kBottom28Bits = 0xfffffff;
 117
 118 // LargeFieldElement also represents an element of the field. The limbs are
 119 // still spaced 28-bits apart and in little-endian order. So the limbs are at
 120 // 0, 28, 56, ..., 392 bits, each 64-bits wide.
 121 typedef uint64 LargeFieldElement[15];
 122
 123 // ReduceLarge converts a LargeFieldElement to a FieldElement.
 124 //
 125 // in[i] < 2**62
 126 void ReduceLarge(FieldElement* out, LargeFieldElement* inptr) {
 127   LargeFieldElement& in(*inptr);
 128
 129   for (int i = 0; i < 8; i++) {
 130     in[i] += kZero63ModP[i];
 131   }
 132
 133   // Eliminate the coefficients at 2**224 and greater while maintaining the
 134   // same value mod p.
 135   for (int i = 14; i >= 8; i--) {
 136     in[i-8] -= in[i];  // reflection off the "+1" term of p.
 137     in[i-5] += (in[i] & 0xffff) << 12;  // part of the "-2**96" reflection.
 138     in[i-4] += in[i] >> 16;  // the rest of the "-2**96" reflection.
 139   }
 140   in[8] = 0;
 141   // in[0..8] < 2**64
 142
 143   // As the values become small enough, we start to store them in |out| and use
 144   // 32-bit operations.
 145   for (int i = 1; i < 8; i++) {
 146     in[i+1] += in[i] >> 28;
 147     (*out)[i] = static_cast<uint32>(in[i] & kBottom28Bits);
 148   }
 149   // Eliminate the term at 2*224 that we introduced while keeping the same
 150   // value mod p.
 151   in[0] -= in[8];  // reflection off the "+1" term of p.
 152   (*out)[3] += static_cast<uint32>(in[8] & 0xffff) << 12; // "-2**96" term
 153   (*out)[4] += static_cast<uint32>(in[8] >> 16);  // rest of "-2**96" term
 154   // in[0] < 2**64
 155   // out[3] < 2**29
 156   // out[4] < 2**29
 157   // out[1,2,5..7] < 2**28
 158
 159   (*out)[0] = static_cast<uint32>(in[0] & kBottom28Bits);
 160   (*out)[1] += static_cast<uint32>((in[0] >> 28) & kBottom28Bits);
 161   (*out)[2] += static_cast<uint32>(in[0] >> 56);
 162   // out[0] < 2**28
 163   // out[1..4] < 2**29
 164   // out[5..7] < 2**28
 165 }
 166
 167 // Mul computes *out = a*b
 168 //
 169 // a[i] < 2**29, b[i] < 2**30 (or vice versa)
 170 // out[i] < 2**29
 171 void Mul(FieldElement* out, const FieldElement& a, const FieldElement& b) {
 172   LargeFieldElement tmp;
 173   memset(&tmp, 0, sizeof(tmp));
 174
 175   for (int i = 0; i < 8; i++) {
 176     for (int j = 0; j < 8; j++) {
 177       tmp[i+j] += static_cast<uint64>(a[i]) * static_cast<uint64>(b[j]);
 178     }
 179   }
 180
 181   ReduceLarge(out, &tmp);
 182 }
 183
 184 // Square computes *out = a*a
 185 //
 186 // a[i] < 2**29
 187 // out[i] < 2**29
 188 void Square(FieldElement* out, const FieldElement& a) {
 189   LargeFieldElement tmp;
 190   memset(&tmp, 0, sizeof(tmp));
 191
 192   for (int i = 0; i < 8; i++) {
 193     for (int j = 0; j <= i; j++) {
 194       uint64 r = static_cast<uint64>(a[i]) * static_cast<uint64>(a[j]);
 195       if (i == j) {
 196         tmp[i+j] += r;
 197       } else {
 198         tmp[i+j] += r << 1;
 199       }
 200     }
 201   }
 202
 203   ReduceLarge(out, &tmp);
 204 }
 205
 206 // Reduce reduces the coefficients of in_out to smaller bounds.
 207 //
 208 // On entry: a[i] < 2**31 + 2**30
 209 // On exit: a[i] < 2**29
 210 void Reduce(FieldElement* in_out) {
 211   FieldElement& a = *in_out;
 212
 213   for (int i = 0; i < 7; i++) {
 214     a[i+1] += a[i] >> 28;
 215     a[i] &= kBottom28Bits;
 216   }
 217   uint32 top = a[7] >> 28;
 218   a[7] &= kBottom28Bits;
 219
 220   // top < 2**4
 221   // Constant-time: mask = (top != 0) ? 0xffffffff : 0
 222   uint32 mask = top;
 223   mask |= mask >> 2;
 224   mask |= mask >> 1;
 225   mask <<= 31;
 226   mask = static_cast<uint32>(static_cast<int32>(mask) >> 31);
 227
 228   // Eliminate top while maintaining the same value mod p.
 229   a[0] -= top;
 230   a[3] += top << 12;
 231
 232   // We may have just made a[0] negative but, if we did, then we must
 233   // have added something to a[3], thus it's > 2**12. Therefore we can
 234   // carry down to a[0].
 235   a[3] -= 1 & mask;
 236   a[2] += mask & ((1<<28) - 1);
 237   a[1] += mask & ((1<<28) - 1);
 238   a[0] += mask & (1<<28);
 239 }
 240
 241 // Invert calcuates *out = in**-1 by computing in**(2**224 - 2**96 - 1), i.e.
 242 // Fermat's little theorem.
 243 void Invert(FieldElement* out, const FieldElement& in) {
 244   FieldElement f1, f2, f3, f4;
 245
 246   Square(&f1, in);                        // 2
 247   Mul(&f1, f1, in);                       // 2**2 - 1
 248   Square(&f1, f1);                        // 2**3 - 2
 249   Mul(&f1, f1, in);                       // 2**3 - 1
 250   Square(&f2, f1);                        // 2**4 - 2
 251   Square(&f2, f2);                        // 2**5 - 4
 252   Square(&f2, f2);                        // 2**6 - 8
 253   Mul(&f1, f1, f2);                       // 2**6 - 1
 254   Square(&f2, f1);                        // 2**7 - 2
 255   for (int i = 0; i < 5; i++) {           // 2**12 - 2**6
 256     Square(&f2, f2);
 257   }
 258   Mul(&f2, f2, f1);                       // 2**12 - 1
 259   Square(&f3, f2);                        // 2**13 - 2
 260   for (int i = 0; i < 11; i++) {          // 2**24 - 2**12
 261     Square(&f3, f3);
 262   }
 263   Mul(&f2, f3, f2);                       // 2**24 - 1
 264   Square(&f3, f2);                        // 2**25 - 2
 265   for (int i = 0; i < 23; i++) {          // 2**48 - 2**24
 266     Square(&f3, f3);
 267   }
 268   Mul(&f3, f3, f2);                       // 2**48 - 1
 269   Square(&f4, f3);                        // 2**49 - 2
 270   for (int i = 0; i < 47; i++) {          // 2**96 - 2**48
 271     Square(&f4, f4);
 272   }
 273   Mul(&f3, f3, f4);                       // 2**96 - 1
 274   Square(&f4, f3);                        // 2**97 - 2
 275   for (int i = 0; i < 23; i++) {          // 2**120 - 2**24
 276     Square(&f4, f4);
 277   }
 278   Mul(&f2, f4, f2);                       // 2**120 - 1
 279   for (int i = 0; i < 6; i++) {           // 2**126 - 2**6
 280     Square(&f2, f2);
 281   }
 282   Mul(&f1, f1, f2);                       // 2**126 - 1
 283   Square(&f1, f1);                        // 2**127 - 2
 284   Mul(&f1, f1, in);                       // 2**127 - 1
 285   for (int i = 0; i < 97; i++) {          // 2**224 - 2**97
 286     Square(&f1, f1);
 287   }
 288   Mul(out, f1, f3);                       // 2**224 - 2**96 - 1
 289 }
 290
 291 // Contract converts a FieldElement to its minimal, distinguished form.
 292 //
 293 // On entry, in[i] < 2**29
 294 // On exit, in[i] < 2**28
 295 void Contract(FieldElement* inout) {
 296   FieldElement& out = *inout;
 297
 298   // Reduce the coefficients to < 2**28.
 299   for (int i = 0; i < 7; i++) {
 300     out[i+1] += out[i] >> 28;
 301     out[i] &= kBottom28Bits;
 302   }
 303   uint32 top = out[7] >> 28;
 304   out[7] &= kBottom28Bits;
 305
 306   // Eliminate top while maintaining the same value mod p.
 307   out[0] -= top;
 308   out[3] += top << 12;
 309
 310   // We may just have made out[0] negative. So we carry down. If we made
 311   // out[0] negative then we know that out[3] is sufficiently positive
 312   // because we just added to it.
 313   for (int i = 0; i < 3; i++) {
 314     uint32 mask = static_cast<uint32>(static_cast<int32>(out[i]) >> 31);
 315     out[i] += (1 << 28) & mask;
 316     out[i+1] -= 1 & mask;
 317   }
 318
 319   // We might have pushed out[3] over 2**28 so we perform another, partial
 320   // carry chain.
 321   for (int i = 3; i < 7; i++) {
 322     out[i+1] += out[i] >> 28;
 323     out[i] &= kBottom28Bits;
 324   }
 325   top = out[7] >> 28;
 326   out[7] &= kBottom28Bits;
 327
 328   // Eliminate top while maintaining the same value mod p.
 329   out[0] -= top;
 330   out[3] += top << 12;
 331
 332   // There are two cases to consider for out[3]:
 333   //   1) The first time that we eliminated top, we didn't push out[3] over
 334   //      2**28. In this case, the partial carry chain didn't change any values
 335   //      and top is zero.
 336   //   2) We did push out[3] over 2**28 the first time that we eliminated top.
 337   //      The first value of top was in [0..16), therefore, prior to eliminating
 338   //      the first top, 0xfff1000 <= out[3] <= 0xfffffff. Therefore, after
 339   //      overflowing and being reduced by the second carry chain, out[3] <=
 340   //      0xf000. Thus it cannot have overflowed when we eliminated top for the
 341   //      second time.
 342
 343   // Again, we may just have made out[0] negative, so do the same carry down.
 344   // As before, if we made out[0] negative then we know that out[3] is
 345   // sufficiently positive.
 346   for (int i = 0; i < 3; i++) {
 347     uint32 mask = static_cast<uint32>(static_cast<int32>(out[i]) >> 31);
 348     out[i] += (1 << 28) & mask;
 349     out[i+1] -= 1 & mask;
 350   }
 351
 352   // The value is < 2**224, but maybe greater than p. In order to reduce to a
 353   // unique, minimal value we see if the value is >= p and, if so, subtract p.
 354
 355   // First we build a mask from the top four limbs, which must all be
 356   // equal to bottom28Bits if the whole value is >= p. If top_4_all_ones
 357   // ends up with any zero bits in the bottom 28 bits, then this wasn't
 358   // true.
 359   uint32 top_4_all_ones = 0xffffffffu;
 360   for (int i = 4; i < 8; i++) {
 361     top_4_all_ones &= out[i];
 362   }
 363   top_4_all_ones |= 0xf0000000;
 364   // Now we replicate any zero bits to all the bits in top_4_all_ones.
 365   top_4_all_ones &= top_4_all_ones >> 16;
 366   top_4_all_ones &= top_4_all_ones >> 8;
 367   top_4_all_ones &= top_4_all_ones >> 4;
 368   top_4_all_ones &= top_4_all_ones >> 2;
 369   top_4_all_ones &= top_4_all_ones >> 1;
 370   top_4_all_ones =
 371       static_cast<uint32>(static_cast<int32>(top_4_all_ones << 31) >> 31);
 372
 373   // Now we test whether the bottom three limbs are non-zero.
 374   uint32 bottom_3_non_zero = out[0] | out[1] | out[2];
 375   bottom_3_non_zero |= bottom_3_non_zero >> 16;
 376   bottom_3_non_zero |= bottom_3_non_zero >> 8;
 377   bottom_3_non_zero |= bottom_3_non_zero >> 4;
 378   bottom_3_non_zero |= bottom_3_non_zero >> 2;
 379   bottom_3_non_zero |= bottom_3_non_zero >> 1;
 380   bottom_3_non_zero =
 381       static_cast<uint32>(static_cast<int32>(bottom_3_non_zero) >> 31);
 382
 383   // Everything depends on the value of out[3].
 384   //    If it's > 0xffff000 and top_4_all_ones != 0 then the whole value is >= p
 385   //    If it's = 0xffff000 and top_4_all_ones != 0 and bottom_3_non_zero != 0,
 386   //      then the whole value is >= p
 387   //    If it's < 0xffff000, then the whole value is < p
 388   uint32 n = out[3] - 0xffff000;
 389   uint32 out_3_equal = n;
 390   out_3_equal |= out_3_equal >> 16;
 391   out_3_equal |= out_3_equal >> 8;
 392   out_3_equal |= out_3_equal >> 4;
 393   out_3_equal |= out_3_equal >> 2;
 394   out_3_equal |= out_3_equal >> 1;
 395   out_3_equal =
 396       ~static_cast<uint32>(static_cast<int32>(out_3_equal << 31) >> 31);
 397
 398   // If out[3] > 0xffff000 then n's MSB will be zero.
 399   uint32 out_3_gt = ~static_cast<uint32>(static_cast<int32>(n << 31) >> 31);
 400
 401   uint32 mask = top_4_all_ones & ((out_3_equal & bottom_3_non_zero) | out_3_gt);
 402   out[0] -= 1 & mask;
 403   out[3] -= 0xffff000 & mask;
 404   out[4] -= 0xfffffff & mask;
 405   out[5] -= 0xfffffff & mask;
 406   out[6] -= 0xfffffff & mask;
 407   out[7] -= 0xfffffff & mask;
 408 }
 409
 410
 411 // Group element functions.
 412 //
 413 // These functions deal with group elements. The group is an elliptic curve
 414 // group with a = -3 defined in FIPS 186-3, section D.2.2.
 415
 416 using crypto::p224::Point;
 417
 418 // kB is parameter of the elliptic curve.
 419 const FieldElement kB = {
 420   55967668, 11768882, 265861671, 185302395,
 421   39211076, 180311059, 84673715, 188764328,
 422 };
 423
 424 void CopyConditional(Point* out, const Point& a, uint32 mask);
 425 void DoubleJacobian(Point* out, const Point& a);
 426
 427 // AddJacobian computes *out = a+b where a != b.
 428 void AddJacobian(Point *out,
 429                  const Point& a,
 430                  const Point& b) {
 431   // See http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#addition-add-2007-bl
 432   FieldElement z1z1, z2z2, u1, u2, s1, s2, h, i, j, r, v;
 433
 434   uint32 z1_is_zero = IsZero(a.z);
 435   uint32 z2_is_zero = IsZero(b.z);
 436
 437   // Z1Z1 = Z1²
 438   Square(&z1z1, a.z);
 439
 440   // Z2Z2 = Z2²
 441   Square(&z2z2, b.z);
 442
 443   // U1 = X1*Z2Z2
 444   Mul(&u1, a.x, z2z2);
 445
 446   // U2 = X2*Z1Z1
 447   Mul(&u2, b.x, z1z1);
 448
 449   // S1 = Y1*Z2*Z2Z2
 450   Mul(&s1, b.z, z2z2);
 451   Mul(&s1, a.y, s1);
 452
 453   // S2 = Y2*Z1*Z1Z1
 454   Mul(&s2, a.z, z1z1);
 455   Mul(&s2, b.y, s2);
 456
 457   // H = U2-U1
 458   Subtract(&h, u2, u1);
 459   Reduce(&h);
 460   uint32 x_equal = IsZero(h);
 461
 462   // I = (2*H)²
 463   for (int j = 0; j < 8; j++) {
 464     i[j] = h[j] << 1;
 465   }
 466   Reduce(&i);
 467   Square(&i, i);
 468
 469   // J = H*I
 470   Mul(&j, h, i);
 471   // r = 2*(S2-S1)
 472   Subtract(&r, s2, s1);
 473   Reduce(&r);
 474   uint32 y_equal = IsZero(r);
 475
 476   if (x_equal && y_equal && !z1_is_zero && !z2_is_zero) {
 477     // The two input points are the same therefore we must use the dedicated
 478     // doubling function as the slope of the line is undefined.
 479     DoubleJacobian(out, a);
 480     return;
 481   }
 482
 483   for (int i = 0; i < 8; i++) {
 484     r[i] <<= 1;
 485   }
 486   Reduce(&r);
 487
 488   // V = U1*I
 489   Mul(&v, u1, i);
 490
 491   // Z3 = ((Z1+Z2)²-Z1Z1-Z2Z2)*H
 492   Add(&z1z1, z1z1, z2z2);
 493   Add(&z2z2, a.z, b.z);
 494   Reduce(&z2z2);
 495   Square(&z2z2, z2z2);
 496   Subtract(&out->z, z2z2, z1z1);
 497   Reduce(&out->z);
 498   Mul(&out->z, out->z, h);
 499
 500   // X3 = r²-J-2*V
 501   for (int i = 0; i < 8; i++) {
 502     z1z1[i] = v[i] << 1;
 503   }
 504   Add(&z1z1, j, z1z1);
 505   Reduce(&z1z1);
 506   Square(&out->x, r);
 507   Subtract(&out->x, out->x, z1z1);
 508   Reduce(&out->x);
 509
 510   // Y3 = r*(V-X3)-2*S1*J
 511   for (int i = 0; i < 8; i++) {
 512     s1[i] <<= 1;
 513   }
 514   Mul(&s1, s1, j);
 515   Subtract(&z1z1, v, out->x);
 516   Reduce(&z1z1);
 517   Mul(&z1z1, z1z1, r);
 518   Subtract(&out->y, z1z1, s1);
 519   Reduce(&out->y);
 520
 521   CopyConditional(out, a, z2_is_zero);
 522   CopyConditional(out, b, z1_is_zero);
 523 }
 524
 525 // DoubleJacobian computes *out = a+a.
 526 void DoubleJacobian(Point* out, const Point& a) {
 527   // See http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#doubling-dbl-2001-b
 528   FieldElement delta, gamma, beta, alpha, t;
 529
 530   Square(&delta, a.z);
 531   Square(&gamma, a.y);
 532   Mul(&beta, a.x, gamma);
 533
 534   // alpha = 3*(X1-delta)*(X1+delta)
 535   Add(&t, a.x, delta);
 536   for (int i = 0; i < 8; i++) {
 537           t[i] += t[i] << 1;
 538   }
 539   Reduce(&t);
 540   Subtract(&alpha, a.x, delta);
 541   Reduce(&alpha);
 542   Mul(&alpha, alpha, t);
 543
 544   // Z3 = (Y1+Z1)²-gamma-delta
 545   Add(&out->z, a.y, a.z);
 546   Reduce(&out->z);
 547   Square(&out->z, out->z);
 548   Subtract(&out->z, out->z, gamma);
 549   Reduce(&out->z);
 550   Subtract(&out->z, out->z, delta);
 551   Reduce(&out->z);
 552
 553   // X3 = alpha²-8*beta
 554   for (int i = 0; i < 8; i++) {
 555           delta[i] = beta[i] << 3;
 556   }
 557   Reduce(&delta);
 558   Square(&out->x, alpha);
 559   Subtract(&out->x, out->x, delta);
 560   Reduce(&out->x);
 561
 562   // Y3 = alpha*(4*beta-X3)-8*gamma²
 563   for (int i = 0; i < 8; i++) {
 564           beta[i] <<= 2;
 565   }
 566   Reduce(&beta);
 567   Subtract(&beta, beta, out->x);
 568   Reduce(&beta);
 569   Square(&gamma, gamma);
 570   for (int i = 0; i < 8; i++) {
 571           gamma[i] <<= 3;
 572   }
 573   Reduce(&gamma);
 574   Mul(&out->y, alpha, beta);
 575   Subtract(&out->y, out->y, gamma);
 576   Reduce(&out->y);
 577 }
 578
 579 // CopyConditional sets *out=a if mask is 0xffffffff. mask must be either 0 of
 580 // 0xffffffff.
 581 void CopyConditional(Point* out,
 582                      const Point& a,
 583                      uint32 mask) {
 584   for (int i = 0; i < 8; i++) {
 585     out->x[i] ^= mask & (a.x[i] ^ out->x[i]);
 586     out->y[i] ^= mask & (a.y[i] ^ out->y[i]);
 587     out->z[i] ^= mask & (a.z[i] ^ out->z[i]);
 588   }
 589 }
 590
 591 // ScalarMult calculates *out = a*scalar where scalar is a big-endian number of
 592 // length scalar_len and != 0.
 593 void ScalarMult(Point* out, const Point& a,
 594                 const uint8* scalar, size_t scalar_len) {
 595   memset(out, 0, sizeof(*out));
 596   Point tmp;
 597
 598   for (size_t i = 0; i < scalar_len; i++) {
 599     for (unsigned int bit_num = 0; bit_num < 8; bit_num++) {
 600       DoubleJacobian(out, *out);
 601       uint32 bit = static_cast<uint32>(static_cast<int32>(
 602           (((scalar[i] >> (7 - bit_num)) & 1) << 31) >> 31));
 603       AddJacobian(&tmp, a, *out);
 604       CopyConditional(out, tmp, bit);
 605     }
 606   }
 607 }
 608
 609 // Get224Bits reads 7 words from in and scatters their contents in
 610 // little-endian form into 8 words at out, 28 bits per output word.
 611 void Get224Bits(uint32* out, const uint32* in) {
 612   out[0] = NetToHost32(in[6]) & kBottom28Bits;
 613   out[1] = ((NetToHost32(in[5]) << 4) |
 614             (NetToHost32(in[6]) >> 28)) & kBottom28Bits;
 615   out[2] = ((NetToHost32(in[4]) << 8) |
 616             (NetToHost32(in[5]) >> 24)) & kBottom28Bits;
 617   out[3] = ((NetToHost32(in[3]) << 12) |
 618             (NetToHost32(in[4]) >> 20)) & kBottom28Bits;
 619   out[4] = ((NetToHost32(in[2]) << 16) |
 620             (NetToHost32(in[3]) >> 16)) & kBottom28Bits;
 621   out[5] = ((NetToHost32(in[1]) << 20) |
 622             (NetToHost32(in[2]) >> 12)) & kBottom28Bits;
 623   out[6] = ((NetToHost32(in[0]) << 24) |
 624             (NetToHost32(in[1]) >> 8)) & kBottom28Bits;
 625   out[7] = (NetToHost32(in[0]) >> 4) & kBottom28Bits;
 626 }
 627
 628 // Put224Bits performs the inverse operation to Get224Bits: taking 28 bits from
 629 // each of 8 input words and writing them in big-endian order to 7 words at
 630 // out.
 631 void Put224Bits(uint32* out, const uint32* in) {
 632   out[6] = HostToNet32((in[0] >> 0) | (in[1] << 28));
 633   out[5] = HostToNet32((in[1] >> 4) | (in[2] << 24));
 634   out[4] = HostToNet32((in[2] >> 8) | (in[3] << 20));
 635   out[3] = HostToNet32((in[3] >> 12) | (in[4] << 16));
 636   out[2] = HostToNet32((in[4] >> 16) | (in[5] << 12));
 637   out[1] = HostToNet32((in[5] >> 20) | (in[6] << 8));
 638   out[0] = HostToNet32((in[6] >> 24) | (in[7] << 4));
 639 }
 640
 641 } // anonymous namespace
 642
 643 namespace crypto {
 644
 645 namespace p224 {
 646
 647 bool Point::SetFromString(const base::StringPiece& in) {
 648   if (in.size() != 2*28)
 649     return false;
 650   const uint32* inwords = reinterpret_cast<const uint32*>(in.data());
 651   Get224Bits(x, inwords);
 652   Get224Bits(y, inwords + 7);
 653   memset(&z, 0, sizeof(z));
 654   z[0] = 1;
 655
 656   // Check that the point is on the curve, i.e. that y² = x³ - 3x + b.
 657   FieldElement lhs;
 658   Square(&lhs, y);
 659   Contract(&lhs);
 660
 661   FieldElement rhs;
 662   Square(&rhs, x);
 663   Mul(&rhs, x, rhs);
 664
 665   FieldElement three_x;
 666   for (int i = 0; i < 8; i++) {
 667     three_x[i] = x[i] * 3;
 668   }
 669   Reduce(&three_x);
 670   Subtract(&rhs, rhs, three_x);
 671   Reduce(&rhs);
 672
 673   ::Add(&rhs, rhs, kB);
 674   Contract(&rhs);
 675   return memcmp(&lhs, &rhs, sizeof(lhs)) == 0;
 676 }
 677
 678 std::string Point::ToString() const {
 679   FieldElement zinv, zinv_sq, x, y;
 680
 681   // If this is the point at infinity we return a string of all zeros.
 682   if (IsZero(this->z)) {
 683     static const char zeros[56] = {0};
 684     return std::string(zeros, sizeof(zeros));
 685   }
 686
 687   Invert(&zinv, this->z);
 688   Square(&zinv_sq, zinv);
 689   Mul(&x, this->x, zinv_sq);
 690   Mul(&zinv_sq, zinv_sq, zinv);
 691   Mul(&y, this->y, zinv_sq);
 692
 693   Contract(&x);
 694   Contract(&y);
 695
 696   uint32 outwords[14];
 697   Put224Bits(outwords, x);
 698   Put224Bits(outwords + 7, y);
 699   return std::string(reinterpret_cast<const char*>(outwords), sizeof(outwords));
 700 }
 701
 702 void ScalarMult(const Point& in, const uint8* scalar, Point* out) {
 703   ::ScalarMult(out, in, scalar, 28);
 704 }
 705
 706 // kBasePoint is the base point (generator) of the elliptic curve group.
 707 static const Point kBasePoint = {
 708   {22813985, 52956513, 34677300, 203240812,
 709    12143107, 133374265, 225162431, 191946955},
 710   {83918388, 223877528, 122119236, 123340192,
 711    266784067, 263504429, 146143011, 198407736},
 712   {1, 0, 0, 0, 0, 0, 0, 0},
 713 };
 714
 715 void ScalarBaseMult(const uint8* scalar, Point* out) {
 716   ::ScalarMult(out, kBasePoint, scalar, 28);
 717 }
 718
 719 void Add(const Point& a, const Point& b, Point* out) {
 720   AddJacobian(out, a, b);
 721 }
 722
 723 void Negate(const Point& in, Point* out) {
 724   // Guide to elliptic curve cryptography, page 89 suggests that (X : X+Y : Z)
 725   // is the negative in Jacobian coordinates, but it doesn't actually appear to
 726   // be true in testing so this performs the negation in affine coordinates.
 727   FieldElement zinv, zinv_sq, y;
 728   Invert(&zinv, in.z);
 729   Square(&zinv_sq, zinv);
 730   Mul(&out->x, in.x, zinv_sq);
 731   Mul(&zinv_sq, zinv_sq, zinv);
 732   Mul(&y, in.y, zinv_sq);
 733
 734   Subtract(&out->y, kP, y);
 735   Reduce(&out->y);
 736
 737   memset(&out->z, 0, sizeof(out->z));
 738   out->z[0] = 1;
 739 }
 740
 741 }  // namespace p224
 742
 743 }  // namespace crypto