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
16namespace {
17
18using base::HostToNet32;
19using 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
33using crypto::p224::FieldElement;
34
35// kP is the P224 prime.
36const FieldElement kP = {
37  1, 0, 0, 268431360,
38  268435455, 268435455, 268435455, 268435455,
39};
40
41void Contract(FieldElement* inout);
42
43// IsZero returns 0xffffffff if a == 0 mod p and 0 otherwise.
44uint32 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
78void 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
84static const uint32 kTwo31p3 = (1u<<31) + (1u<<3);
85static const uint32 kTwo31m3 = (1u<<31) - (1u<<3);
86static 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.
90static 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
99void 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
106static const uint64 kTwo63p35 = (1ull<<63) + (1ull<<35);
107static const uint64 kTwo63m35 = (1ull<<63) - (1ull<<35);
108static 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.
111static const uint64 kZero63ModP[8] = {
112  kTwo63p35, kTwo63m35, kTwo63m35, kTwo63m35,
113  kTwo63m35m19, kTwo63m35, kTwo63m35, kTwo63m35,
114};
115
116static 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.
121typedef uint64 LargeFieldElement[15];
122
123// ReduceLarge converts a LargeFieldElement to a FieldElement.
124//
125// in[i] < 2**62
126void 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
171void 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
188void 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
210void 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.
243void 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
295void 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
416using crypto::p224::Point;
417
418// kB is parameter of the elliptic curve.
419const FieldElement kB = {
420  55967668, 11768882, 265861671, 185302395,
421  39211076, 180311059, 84673715, 188764328,
422};
423
424void CopyConditional(Point* out, const Point& a, uint32 mask);
425void DoubleJacobian(Point* out, const Point& a);
426
427// AddJacobian computes *out = a+b where a != b.
428void 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.
526void 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.
581void 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.
593void 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.
611void 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.
631void 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
643namespace crypto {
644
645namespace p224 {
646
647bool 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
678std::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
702void 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.
707static 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
715void ScalarBaseMult(const uint8* scalar, Point* out) {
716  ::ScalarMult(out, kBasePoint, scalar, 28);
717}
718
719void Add(const Point& a, const Point& b, Point* out) {
720  AddJacobian(out, a, b);
721}
722
723void 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
744