libgo/go/crypto/elliptic/p224.go

   1 // Copyright 2012 The Go Authors.  All rights reserved.
   2 // Use of this source code is governed by a BSD-style
   3 // license that can be found in the LICENSE file.
   4
   5 package elliptic
   6
   7 // This is a constant-time, 32-bit implementation of P224. See FIPS 186-3,
   8 // section D.2.2.
   9 //
  10 // See http://www.imperialviolet.org/2010/12/04/ecc.html ([1]) for background.
  11
  12 import (
  13         "math/big"
  14 )
  15
  16 var p224 p224Curve
  17
  18 type p224Curve struct {
  19         *CurveParams
  20         gx, gy, b p224FieldElement
  21 }
  22
  23 func initP224() {
  24         // See FIPS 186-3, section D.2.2
  25         p224.CurveParams = new(CurveParams)
  26         p224.P, _ = new(big.Int).SetString("26959946667150639794667015087019630673557916260026308143510066298881", 10)
  27         p224.N, _ = new(big.Int).SetString("26959946667150639794667015087019625940457807714424391721682722368061", 10)
  28         p224.B, _ = new(big.Int).SetString("b4050a850c04b3abf54132565044b0b7d7bfd8ba270b39432355ffb4", 16)
  29         p224.Gx, _ = new(big.Int).SetString("b70e0cbd6bb4bf7f321390b94a03c1d356c21122343280d6115c1d21", 16)
  30         p224.Gy, _ = new(big.Int).SetString("bd376388b5f723fb4c22dfe6cd4375a05a07476444d5819985007e34", 16)
  31         p224.BitSize = 224
  32
  33         p224FromBig(&p224.gx, p224.Gx)
  34         p224FromBig(&p224.gy, p224.Gy)
  35         p224FromBig(&p224.b, p224.B)
  36 }
  37
  38 // P224 returns a Curve which implements P-224 (see FIPS 186-3, section D.2.2)
  39 func P224() Curve {
  40         initonce.Do(initAll)
  41         return p224
  42 }
  43
  44 func (curve p224Curve) Params() *CurveParams {
  45         return curve.CurveParams
  46 }
  47
  48 func (curve p224Curve) IsOnCurve(bigX, bigY *big.Int) bool {
  49         var x, y p224FieldElement
  50         p224FromBig(&x, bigX)
  51         p224FromBig(&y, bigY)
  52
  53         // y² = x³ - 3x + b
  54         var tmp p224LargeFieldElement
  55         var x3 p224FieldElement
  56         p224Square(&x3, &x, &tmp)
  57         p224Mul(&x3, &x3, &x, &tmp)
  58
  59         for i := 0; i < 8; i++ {
  60                 x[i] *= 3
  61         }
  62         p224Sub(&x3, &x3, &x)
  63         p224Reduce(&x3)
  64         p224Add(&x3, &x3, &curve.b)
  65         p224Contract(&x3, &x3)
  66
  67         p224Square(&y, &y, &tmp)
  68         p224Contract(&y, &y)
  69
  70         for i := 0; i < 8; i++ {
  71                 if y[i] != x3[i] {
  72                         return false
  73                 }
  74         }
  75         return true
  76 }
  77
  78 func (p224Curve) Add(bigX1, bigY1, bigX2, bigY2 *big.Int) (x, y *big.Int) {
  79         var x1, y1, z1, x2, y2, z2, x3, y3, z3 p224FieldElement
  80
  81         p224FromBig(&x1, bigX1)
  82         p224FromBig(&y1, bigY1)
  83         if bigX1.Sign() != 0 || bigY1.Sign() != 0 {
  84                 z1[0] = 1
  85         }
  86         p224FromBig(&x2, bigX2)
  87         p224FromBig(&y2, bigY2)
  88         if bigX2.Sign() != 0 || bigY2.Sign() != 0 {
  89                 z2[0] = 1
  90         }
  91
  92         p224AddJacobian(&x3, &y3, &z3, &x1, &y1, &z1, &x2, &y2, &z2)
  93         return p224ToAffine(&x3, &y3, &z3)
  94 }
  95
  96 func (p224Curve) Double(bigX1, bigY1 *big.Int) (x, y *big.Int) {
  97         var x1, y1, z1, x2, y2, z2 p224FieldElement
  98
  99         p224FromBig(&x1, bigX1)
 100         p224FromBig(&y1, bigY1)
 101         z1[0] = 1
 102
 103         p224DoubleJacobian(&x2, &y2, &z2, &x1, &y1, &z1)
 104         return p224ToAffine(&x2, &y2, &z2)
 105 }
 106
 107 func (p224Curve) ScalarMult(bigX1, bigY1 *big.Int, scalar []byte) (x, y *big.Int) {
 108         var x1, y1, z1, x2, y2, z2 p224FieldElement
 109
 110         p224FromBig(&x1, bigX1)
 111         p224FromBig(&y1, bigY1)
 112         z1[0] = 1
 113
 114         p224ScalarMult(&x2, &y2, &z2, &x1, &y1, &z1, scalar)
 115         return p224ToAffine(&x2, &y2, &z2)
 116 }
 117
 118 func (curve p224Curve) ScalarBaseMult(scalar []byte) (x, y *big.Int) {
 119         var z1, x2, y2, z2 p224FieldElement
 120
 121         z1[0] = 1
 122         p224ScalarMult(&x2, &y2, &z2, &curve.gx, &curve.gy, &z1, scalar)
 123         return p224ToAffine(&x2, &y2, &z2)
 124 }
 125
 126 // Field element functions.
 127 //
 128 // The field that we're dealing with is ℤ/pℤ where p = 2**224 - 2**96 + 1.
 129 //
 130 // Field elements are represented by a FieldElement, which is a typedef to an
 131 // array of 8 uint32's. The value of a FieldElement, a, is:
 132 //   a[0] + 2**28·a[1] + 2**56·a[1] + ... + 2**196·a[7]
 133 //
 134 // Using 28-bit limbs means that there's only 4 bits of headroom, which is less
 135 // than we would really like. But it has the useful feature that we hit 2**224
 136 // exactly, making the reflections during a reduce much nicer.
 137 type p224FieldElement [8]uint32
 138
 139 // p224P is the order of the field, represented as a p224FieldElement.
 140 var p224P = [8]uint32{1, 0, 0, 0xffff000, 0xfffffff, 0xfffffff, 0xfffffff, 0xfffffff}
 141
 142 // p224IsZero returns 1 if a == 0 mod p and 0 otherwise.
 143 //
 144 // a[i] < 2**29
 145 func p224IsZero(a *p224FieldElement) uint32 {
 146         // Since a p224FieldElement contains 224 bits there are two possible
 147         // representations of 0: 0 and p.
 148         var minimal p224FieldElement
 149         p224Contract(&minimal, a)
 150
 151         var isZero, isP uint32
 152         for i, v := range minimal {
 153                 isZero |= v
 154                 isP |= v - p224P[i]
 155         }
 156
 157         // If either isZero or isP is 0, then we should return 1.
 158         isZero |= isZero >> 16
 159         isZero |= isZero >> 8
 160         isZero |= isZero >> 4
 161         isZero |= isZero >> 2
 162         isZero |= isZero >> 1
 163
 164         isP |= isP >> 16
 165         isP |= isP >> 8
 166         isP |= isP >> 4
 167         isP |= isP >> 2
 168         isP |= isP >> 1
 169
 170         // For isZero and isP, the LSB is 0 iff all the bits are zero.
 171         result := isZero & isP
 172         result = (^result) & 1
 173
 174         return result
 175 }
 176
 177 // p224Add computes *out = a+b
 178 //
 179 // a[i] + b[i] < 2**32
 180 func p224Add(out, a, b *p224FieldElement) {
 181         for i := 0; i < 8; i++ {
 182                 out[i] = a[i] + b[i]
 183         }
 184 }
 185
 186 const two31p3 = 1<<31 + 1<<3
 187 const two31m3 = 1<<31 - 1<<3
 188 const two31m15m3 = 1<<31 - 1<<15 - 1<<3
 189
 190 // p224ZeroModP31 is 0 mod p where bit 31 is set in all limbs so that we can
 191 // subtract smaller amounts without underflow. See the section "Subtraction" in
 192 // [1] for reasoning.
 193 var p224ZeroModP31 = []uint32{two31p3, two31m3, two31m3, two31m15m3, two31m3, two31m3, two31m3, two31m3}
 194
 195 // p224Sub computes *out = a-b
 196 //
 197 // a[i], b[i] < 2**30
 198 // out[i] < 2**32
 199 func p224Sub(out, a, b *p224FieldElement) {
 200         for i := 0; i < 8; i++ {
 201                 out[i] = a[i] + p224ZeroModP31[i] - b[i]
 202         }
 203 }
 204
 205 // LargeFieldElement also represents an element of the field. The limbs are
 206 // still spaced 28-bits apart and in little-endian order. So the limbs are at
 207 // 0, 28, 56, ..., 392 bits, each 64-bits wide.
 208 type p224LargeFieldElement [15]uint64
 209
 210 const two63p35 = 1<<63 + 1<<35
 211 const two63m35 = 1<<63 - 1<<35
 212 const two63m35m19 = 1<<63 - 1<<35 - 1<<19
 213
 214 // p224ZeroModP63 is 0 mod p where bit 63 is set in all limbs. See the section
 215 // "Subtraction" in [1] for why.
 216 var p224ZeroModP63 = [8]uint64{two63p35, two63m35, two63m35, two63m35, two63m35m19, two63m35, two63m35, two63m35}
 217
 218 const bottom12Bits = 0xfff
 219 const bottom28Bits = 0xfffffff
 220
 221 // p224Mul computes *out = a*b
 222 //
 223 // a[i] < 2**29, b[i] < 2**30 (or vice versa)
 224 // out[i] < 2**29
 225 func p224Mul(out, a, b *p224FieldElement, tmp *p224LargeFieldElement) {
 226         for i := 0; i < 15; i++ {
 227                 tmp[i] = 0
 228         }
 229
 230         for i := 0; i < 8; i++ {
 231                 for j := 0; j < 8; j++ {
 232                         tmp[i+j] += uint64(a[i]) * uint64(b[j])
 233                 }
 234         }
 235
 236         p224ReduceLarge(out, tmp)
 237 }
 238
 239 // Square computes *out = a*a
 240 //
 241 // a[i] < 2**29
 242 // out[i] < 2**29
 243 func p224Square(out, a *p224FieldElement, tmp *p224LargeFieldElement) {
 244         for i := 0; i < 15; i++ {
 245                 tmp[i] = 0
 246         }
 247
 248         for i := 0; i < 8; i++ {
 249                 for j := 0; j <= i; j++ {
 250                         r := uint64(a[i]) * uint64(a[j])
 251                         if i == j {
 252                                 tmp[i+j] += r
 253                         } else {
 254                                 tmp[i+j] += r << 1
 255                         }
 256                 }
 257         }
 258
 259         p224ReduceLarge(out, tmp)
 260 }
 261
 262 // ReduceLarge converts a p224LargeFieldElement to a p224FieldElement.
 263 //
 264 // in[i] < 2**62
 265 func p224ReduceLarge(out *p224FieldElement, in *p224LargeFieldElement) {
 266         for i := 0; i < 8; i++ {
 267                 in[i] += p224ZeroModP63[i]
 268         }
 269
 270         // Eliminate the coefficients at 2**224 and greater.
 271         for i := 14; i >= 8; i-- {
 272                 in[i-8] -= in[i]
 273                 in[i-5] += (in[i] & 0xffff) << 12
 274                 in[i-4] += in[i] >> 16
 275         }
 276         in[8] = 0
 277         // in[0..8] < 2**64
 278
 279         // As the values become small enough, we start to store them in |out|
 280         // and use 32-bit operations.
 281         for i := 1; i < 8; i++ {
 282                 in[i+1] += in[i] >> 28
 283                 out[i] = uint32(in[i] & bottom28Bits)
 284         }
 285         in[0] -= in[8]
 286         out[3] += uint32(in[8]&0xffff) << 12
 287         out[4] += uint32(in[8] >> 16)
 288         // in[0] < 2**64
 289         // out[3] < 2**29
 290         // out[4] < 2**29
 291         // out[1,2,5..7] < 2**28
 292
 293         out[0] = uint32(in[0] & bottom28Bits)
 294         out[1] += uint32((in[0] >> 28) & bottom28Bits)
 295         out[2] += uint32(in[0] >> 56)
 296         // out[0] < 2**28
 297         // out[1..4] < 2**29
 298         // out[5..7] < 2**28
 299 }
 300
 301 // Reduce reduces the coefficients of a to smaller bounds.
 302 //
 303 // On entry: a[i] < 2**31 + 2**30
 304 // On exit: a[i] < 2**29
 305 func p224Reduce(a *p224FieldElement) {
 306         for i := 0; i < 7; i++ {
 307                 a[i+1] += a[i] >> 28
 308                 a[i] &= bottom28Bits
 309         }
 310         top := a[7] >> 28
 311         a[7] &= bottom28Bits
 312
 313         // top < 2**4
 314         mask := top
 315         mask |= mask >> 2
 316         mask |= mask >> 1
 317         mask <<= 31
 318         mask = uint32(int32(mask) >> 31)
 319         // Mask is all ones if top != 0, all zero otherwise
 320
 321         a[0] -= top
 322         a[3] += top << 12
 323
 324         // We may have just made a[0] negative but, if we did, then we must
 325         // have added something to a[3], this it's > 2**12. Therefore we can
 326         // carry down to a[0].
 327         a[3] -= 1 & mask
 328         a[2] += mask & (1<<28 - 1)
 329         a[1] += mask & (1<<28 - 1)
 330         a[0] += mask & (1 << 28)
 331 }
 332
 333 // p224Invert calculates *out = in**-1 by computing in**(2**224 - 2**96 - 1),
 334 // i.e. Fermat's little theorem.
 335 func p224Invert(out, in *p224FieldElement) {
 336         var f1, f2, f3, f4 p224FieldElement
 337         var c p224LargeFieldElement
 338
 339         p224Square(&f1, in, &c)    // 2
 340         p224Mul(&f1, &f1, in, &c)  // 2**2 - 1
 341         p224Square(&f1, &f1, &c)   // 2**3 - 2
 342         p224Mul(&f1, &f1, in, &c)  // 2**3 - 1
 343         p224Square(&f2, &f1, &c)   // 2**4 - 2
 344         p224Square(&f2, &f2, &c)   // 2**5 - 4
 345         p224Square(&f2, &f2, &c)   // 2**6 - 8
 346         p224Mul(&f1, &f1, &f2, &c) // 2**6 - 1
 347         p224Square(&f2, &f1, &c)   // 2**7 - 2
 348         for i := 0; i < 5; i++ {   // 2**12 - 2**6
 349                 p224Square(&f2, &f2, &c)
 350         }
 351         p224Mul(&f2, &f2, &f1, &c) // 2**12 - 1
 352         p224Square(&f3, &f2, &c)   // 2**13 - 2
 353         for i := 0; i < 11; i++ {  // 2**24 - 2**12
 354                 p224Square(&f3, &f3, &c)
 355         }
 356         p224Mul(&f2, &f3, &f2, &c) // 2**24 - 1
 357         p224Square(&f3, &f2, &c)   // 2**25 - 2
 358         for i := 0; i < 23; i++ {  // 2**48 - 2**24
 359                 p224Square(&f3, &f3, &c)
 360         }
 361         p224Mul(&f3, &f3, &f2, &c) // 2**48 - 1
 362         p224Square(&f4, &f3, &c)   // 2**49 - 2
 363         for i := 0; i < 47; i++ {  // 2**96 - 2**48
 364                 p224Square(&f4, &f4, &c)
 365         }
 366         p224Mul(&f3, &f3, &f4, &c) // 2**96 - 1
 367         p224Square(&f4, &f3, &c)   // 2**97 - 2
 368         for i := 0; i < 23; i++ {  // 2**120 - 2**24
 369                 p224Square(&f4, &f4, &c)
 370         }
 371         p224Mul(&f2, &f4, &f2, &c) // 2**120 - 1
 372         for i := 0; i < 6; i++ {   // 2**126 - 2**6
 373                 p224Square(&f2, &f2, &c)
 374         }
 375         p224Mul(&f1, &f1, &f2, &c) // 2**126 - 1
 376         p224Square(&f1, &f1, &c)   // 2**127 - 2
 377         p224Mul(&f1, &f1, in, &c)  // 2**127 - 1
 378         for i := 0; i < 97; i++ {  // 2**224 - 2**97
 379                 p224Square(&f1, &f1, &c)
 380         }
 381         p224Mul(out, &f1, &f3, &c) // 2**224 - 2**96 - 1
 382 }
 383
 384 // p224Contract converts a FieldElement to its unique, minimal form.
 385 //
 386 // On entry, in[i] < 2**29
 387 // On exit, in[i] < 2**28
 388 func p224Contract(out, in *p224FieldElement) {
 389         copy(out[:], in[:])
 390
 391         for i := 0; i < 7; i++ {
 392                 out[i+1] += out[i] >> 28
 393                 out[i] &= bottom28Bits
 394         }
 395         top := out[7] >> 28
 396         out[7] &= bottom28Bits
 397
 398         out[0] -= top
 399         out[3] += top << 12
 400
 401         // We may just have made out[i] negative. So we carry down. If we made
 402         // out[0] negative then we know that out[3] is sufficiently positive
 403         // because we just added to it.
 404         for i := 0; i < 3; i++ {
 405                 mask := uint32(int32(out[i]) >> 31)
 406                 out[i] += (1 << 28) & mask
 407                 out[i+1] -= 1 & mask
 408         }
 409
 410         // We might have pushed out[3] over 2**28 so we perform another, partial,
 411         // carry chain.
 412         for i := 3; i < 7; i++ {
 413                 out[i+1] += out[i] >> 28
 414                 out[i] &= bottom28Bits
 415         }
 416         top = out[7] >> 28
 417         out[7] &= bottom28Bits
 418
 419         // Eliminate top while maintaining the same value mod p.
 420         out[0] -= top
 421         out[3] += top << 12
 422
 423         // There are two cases to consider for out[3]:
 424         //   1) The first time that we eliminated top, we didn't push out[3] over
 425         //      2**28. In this case, the partial carry chain didn't change any values
 426         //      and top is zero.
 427         //   2) We did push out[3] over 2**28 the first time that we eliminated top.
 428         //      The first value of top was in [0..16), therefore, prior to eliminating
 429         //      the first top, 0xfff1000 <= out[3] <= 0xfffffff. Therefore, after
 430         //      overflowing and being reduced by the second carry chain, out[3] <=
 431         //      0xf000. Thus it cannot have overflowed when we eliminated top for the
 432         //      second time.
 433
 434         // Again, we may just have made out[0] negative, so do the same carry down.
 435         // As before, if we made out[0] negative then we know that out[3] is
 436         // sufficiently positive.
 437         for i := 0; i < 3; i++ {
 438                 mask := uint32(int32(out[i]) >> 31)
 439                 out[i] += (1 << 28) & mask
 440                 out[i+1] -= 1 & mask
 441         }
 442
 443         // Now we see if the value is >= p and, if so, subtract p.
 444
 445         // First we build a mask from the top four limbs, which must all be
 446         // equal to bottom28Bits if the whole value is >= p. If top4AllOnes
 447         // ends up with any zero bits in the bottom 28 bits, then this wasn't
 448         // true.
 449         top4AllOnes := uint32(0xffffffff)
 450         for i := 4; i < 8; i++ {
 451                 top4AllOnes &= out[i]
 452         }
 453         top4AllOnes |= 0xf0000000
 454         // Now we replicate any zero bits to all the bits in top4AllOnes.
 455         top4AllOnes &= top4AllOnes >> 16
 456         top4AllOnes &= top4AllOnes >> 8
 457         top4AllOnes &= top4AllOnes >> 4
 458         top4AllOnes &= top4AllOnes >> 2
 459         top4AllOnes &= top4AllOnes >> 1
 460         top4AllOnes = uint32(int32(top4AllOnes<<31) >> 31)
 461
 462         // Now we test whether the bottom three limbs are non-zero.
 463         bottom3NonZero := out[0] | out[1] | out[2]
 464         bottom3NonZero |= bottom3NonZero >> 16
 465         bottom3NonZero |= bottom3NonZero >> 8
 466         bottom3NonZero |= bottom3NonZero >> 4
 467         bottom3NonZero |= bottom3NonZero >> 2
 468         bottom3NonZero |= bottom3NonZero >> 1
 469         bottom3NonZero = uint32(int32(bottom3NonZero<<31) >> 31)
 470
 471         // Everything depends on the value of out[3].
 472         //    If it's > 0xffff000 and top4AllOnes != 0 then the whole value is >= p
 473         //    If it's = 0xffff000 and top4AllOnes != 0 and bottom3NonZero != 0,
 474         //      then the whole value is >= p
 475         //    If it's < 0xffff000, then the whole value is < p
 476         n := out[3] - 0xffff000
 477         out3Equal := n
 478         out3Equal |= out3Equal >> 16
 479         out3Equal |= out3Equal >> 8
 480         out3Equal |= out3Equal >> 4
 481         out3Equal |= out3Equal >> 2
 482         out3Equal |= out3Equal >> 1
 483         out3Equal = ^uint32(int32(out3Equal<<31) >> 31)
 484
 485         // If out[3] > 0xffff000 then n's MSB will be zero.
 486         out3GT := ^uint32(int32(n) >> 31)
 487
 488         mask := top4AllOnes & ((out3Equal & bottom3NonZero) | out3GT)
 489         out[0] -= 1 & mask
 490         out[3] -= 0xffff000 & mask
 491         out[4] -= 0xfffffff & mask
 492         out[5] -= 0xfffffff & mask
 493         out[6] -= 0xfffffff & mask
 494         out[7] -= 0xfffffff & mask
 495 }
 496
 497 // Group element functions.
 498 //
 499 // These functions deal with group elements. The group is an elliptic curve
 500 // group with a = -3 defined in FIPS 186-3, section D.2.2.
 501
 502 // p224AddJacobian computes *out = a+b where a != b.
 503 func p224AddJacobian(x3, y3, z3, x1, y1, z1, x2, y2, z2 *p224FieldElement) {
 504         // See http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#addition-p224Add-2007-bl
 505         var z1z1, z2z2, u1, u2, s1, s2, h, i, j, r, v p224FieldElement
 506         var c p224LargeFieldElement
 507
 508         z1IsZero := p224IsZero(z1)
 509         z2IsZero := p224IsZero(z2)
 510
 511         // Z1Z1 = Z1²
 512         p224Square(&z1z1, z1, &c)
 513         // Z2Z2 = Z2²
 514         p224Square(&z2z2, z2, &c)
 515         // U1 = X1*Z2Z2
 516         p224Mul(&u1, x1, &z2z2, &c)
 517         // U2 = X2*Z1Z1
 518         p224Mul(&u2, x2, &z1z1, &c)
 519         // S1 = Y1*Z2*Z2Z2
 520         p224Mul(&s1, z2, &z2z2, &c)
 521         p224Mul(&s1, y1, &s1, &c)
 522         // S2 = Y2*Z1*Z1Z1
 523         p224Mul(&s2, z1, &z1z1, &c)
 524         p224Mul(&s2, y2, &s2, &c)
 525         // H = U2-U1
 526         p224Sub(&h, &u2, &u1)
 527         p224Reduce(&h)
 528         xEqual := p224IsZero(&h)
 529         // I = (2*H)²
 530         for j := 0; j < 8; j++ {
 531                 i[j] = h[j] << 1
 532         }
 533         p224Reduce(&i)
 534         p224Square(&i, &i, &c)
 535         // J = H*I
 536         p224Mul(&j, &h, &i, &c)
 537         // r = 2*(S2-S1)
 538         p224Sub(&r, &s2, &s1)
 539         p224Reduce(&r)
 540         yEqual := p224IsZero(&r)
 541         if xEqual == 1 && yEqual == 1 && z1IsZero == 0 && z2IsZero == 0 {
 542                 p224DoubleJacobian(x3, y3, z3, x1, y1, z1)
 543                 return
 544         }
 545         for i := 0; i < 8; i++ {
 546                 r[i] <<= 1
 547         }
 548         p224Reduce(&r)
 549         // V = U1*I
 550         p224Mul(&v, &u1, &i, &c)
 551         // Z3 = ((Z1+Z2)²-Z1Z1-Z2Z2)*H
 552         p224Add(&z1z1, &z1z1, &z2z2)
 553         p224Add(&z2z2, z1, z2)
 554         p224Reduce(&z2z2)
 555         p224Square(&z2z2, &z2z2, &c)
 556         p224Sub(z3, &z2z2, &z1z1)
 557         p224Reduce(z3)
 558         p224Mul(z3, z3, &h, &c)
 559         // X3 = r²-J-2*V
 560         for i := 0; i < 8; i++ {
 561                 z1z1[i] = v[i] << 1
 562         }
 563         p224Add(&z1z1, &j, &z1z1)
 564         p224Reduce(&z1z1)
 565         p224Square(x3, &r, &c)
 566         p224Sub(x3, x3, &z1z1)
 567         p224Reduce(x3)
 568         // Y3 = r*(V-X3)-2*S1*J
 569         for i := 0; i < 8; i++ {
 570                 s1[i] <<= 1
 571         }
 572         p224Mul(&s1, &s1, &j, &c)
 573         p224Sub(&z1z1, &v, x3)
 574         p224Reduce(&z1z1)
 575         p224Mul(&z1z1, &z1z1, &r, &c)
 576         p224Sub(y3, &z1z1, &s1)
 577         p224Reduce(y3)
 578
 579         p224CopyConditional(x3, x2, z1IsZero)
 580         p224CopyConditional(x3, x1, z2IsZero)
 581         p224CopyConditional(y3, y2, z1IsZero)
 582         p224CopyConditional(y3, y1, z2IsZero)
 583         p224CopyConditional(z3, z2, z1IsZero)
 584         p224CopyConditional(z3, z1, z2IsZero)
 585 }
 586
 587 // p224DoubleJacobian computes *out = a+a.
 588 func p224DoubleJacobian(x3, y3, z3, x1, y1, z1 *p224FieldElement) {
 589         var delta, gamma, beta, alpha, t p224FieldElement
 590         var c p224LargeFieldElement
 591
 592         p224Square(&delta, z1, &c)
 593         p224Square(&gamma, y1, &c)
 594         p224Mul(&beta, x1, &gamma, &c)
 595
 596         // alpha = 3*(X1-delta)*(X1+delta)
 597         p224Add(&t, x1, &delta)
 598         for i := 0; i < 8; i++ {
 599                 t[i] += t[i] << 1
 600         }
 601         p224Reduce(&t)
 602         p224Sub(&alpha, x1, &delta)
 603         p224Reduce(&alpha)
 604         p224Mul(&alpha, &alpha, &t, &c)
 605
 606         // Z3 = (Y1+Z1)²-gamma-delta
 607         p224Add(z3, y1, z1)
 608         p224Reduce(z3)
 609         p224Square(z3, z3, &c)
 610         p224Sub(z3, z3, &gamma)
 611         p224Reduce(z3)
 612         p224Sub(z3, z3, &delta)
 613         p224Reduce(z3)
 614
 615         // X3 = alpha²-8*beta
 616         for i := 0; i < 8; i++ {
 617                 delta[i] = beta[i] << 3
 618         }
 619         p224Reduce(&delta)
 620         p224Square(x3, &alpha, &c)
 621         p224Sub(x3, x3, &delta)
 622         p224Reduce(x3)
 623
 624         // Y3 = alpha*(4*beta-X3)-8*gamma²
 625         for i := 0; i < 8; i++ {
 626                 beta[i] <<= 2
 627         }
 628         p224Sub(&beta, &beta, x3)
 629         p224Reduce(&beta)
 630         p224Square(&gamma, &gamma, &c)
 631         for i := 0; i < 8; i++ {
 632                 gamma[i] <<= 3
 633         }
 634         p224Reduce(&gamma)
 635         p224Mul(y3, &alpha, &beta, &c)
 636         p224Sub(y3, y3, &gamma)
 637         p224Reduce(y3)
 638 }
 639
 640 // p224CopyConditional sets *out = *in iff the least-significant-bit of control
 641 // is true, and it runs in constant time.
 642 func p224CopyConditional(out, in *p224FieldElement, control uint32) {
 643         control <<= 31
 644         control = uint32(int32(control) >> 31)
 645
 646         for i := 0; i < 8; i++ {
 647                 out[i] ^= (out[i] ^ in[i]) & control
 648         }
 649 }
 650
 651 func p224ScalarMult(outX, outY, outZ, inX, inY, inZ *p224FieldElement, scalar []byte) {
 652         var xx, yy, zz p224FieldElement
 653         for i := 0; i < 8; i++ {
 654                 outX[i] = 0
 655                 outY[i] = 0
 656                 outZ[i] = 0
 657         }
 658
 659         for _, byte := range scalar {
 660                 for bitNum := uint(0); bitNum < 8; bitNum++ {
 661                         p224DoubleJacobian(outX, outY, outZ, outX, outY, outZ)
 662                         bit := uint32((byte >> (7 - bitNum)) & 1)
 663                         p224AddJacobian(&xx, &yy, &zz, inX, inY, inZ, outX, outY, outZ)
 664                         p224CopyConditional(outX, &xx, bit)
 665                         p224CopyConditional(outY, &yy, bit)
 666                         p224CopyConditional(outZ, &zz, bit)
 667                 }
 668         }
 669 }
 670
 671 // p224ToAffine converts from Jacobian to affine form.
 672 func p224ToAffine(x, y, z *p224FieldElement) (*big.Int, *big.Int) {
 673         var zinv, zinvsq, outx, outy p224FieldElement
 674         var tmp p224LargeFieldElement
 675
 676         if isPointAtInfinity := p224IsZero(z); isPointAtInfinity == 1 {
 677                 return new(big.Int), new(big.Int)
 678         }
 679
 680         p224Invert(&zinv, z)
 681         p224Square(&zinvsq, &zinv, &tmp)
 682         p224Mul(x, x, &zinvsq, &tmp)
 683         p224Mul(&zinvsq, &zinvsq, &zinv, &tmp)
 684         p224Mul(y, y, &zinvsq, &tmp)
 685
 686         p224Contract(&outx, x)
 687         p224Contract(&outy, y)
 688         return p224ToBig(&outx), p224ToBig(&outy)
 689 }
 690
 691 // get28BitsFromEnd returns the least-significant 28 bits from buf>>shift,
 692 // where buf is interpreted as a big-endian number.
 693 func get28BitsFromEnd(buf []byte, shift uint) (uint32, []byte) {
 694         var ret uint32
 695
 696         for i := uint(0); i < 4; i++ {
 697                 var b byte
 698                 if l := len(buf); l > 0 {
 699                         b = buf[l-1]
 700                         // We don't remove the byte if we're about to return and we're not
 701                         // reading all of it.
 702                         if i != 3 || shift == 4 {
 703                                 buf = buf[:l-1]
 704                         }
 705                 }
 706                 ret |= uint32(b) << (8 * i) >> shift
 707         }
 708         ret &= bottom28Bits
 709         return ret, buf
 710 }
 711
 712 // p224FromBig sets *out = *in.
 713 func p224FromBig(out *p224FieldElement, in *big.Int) {
 714         bytes := in.Bytes()
 715         out[0], bytes = get28BitsFromEnd(bytes, 0)
 716         out[1], bytes = get28BitsFromEnd(bytes, 4)
 717         out[2], bytes = get28BitsFromEnd(bytes, 0)
 718         out[3], bytes = get28BitsFromEnd(bytes, 4)
 719         out[4], bytes = get28BitsFromEnd(bytes, 0)
 720         out[5], bytes = get28BitsFromEnd(bytes, 4)
 721         out[6], bytes = get28BitsFromEnd(bytes, 0)
 722         out[7], bytes = get28BitsFromEnd(bytes, 4)
 723 }
 724
 725 // p224ToBig returns in as a big.Int.
 726 func p224ToBig(in *p224FieldElement) *big.Int {
 727         var buf [28]byte
 728         buf[27] = byte(in[0])
 729         buf[26] = byte(in[0] >> 8)
 730         buf[25] = byte(in[0] >> 16)
 731         buf[24] = byte(((in[0] >> 24) & 0x0f) | (in[1]<<4)&0xf0)
 732
 733         buf[23] = byte(in[1] >> 4)
 734         buf[22] = byte(in[1] >> 12)
 735         buf[21] = byte(in[1] >> 20)
 736
 737         buf[20] = byte(in[2])
 738         buf[19] = byte(in[2] >> 8)
 739         buf[18] = byte(in[2] >> 16)
 740         buf[17] = byte(((in[2] >> 24) & 0x0f) | (in[3]<<4)&0xf0)
 741
 742         buf[16] = byte(in[3] >> 4)
 743         buf[15] = byte(in[3] >> 12)
 744         buf[14] = byte(in[3] >> 20)
 745
 746         buf[13] = byte(in[4])
 747         buf[12] = byte(in[4] >> 8)
 748         buf[11] = byte(in[4] >> 16)
 749         buf[10] = byte(((in[4] >> 24) & 0x0f) | (in[5]<<4)&0xf0)
 750
 751         buf[9] = byte(in[5] >> 4)
 752         buf[8] = byte(in[5] >> 12)
 753         buf[7] = byte(in[5] >> 20)
 754
 755         buf[6] = byte(in[6])
 756         buf[5] = byte(in[6] >> 8)
 757         buf[4] = byte(in[6] >> 16)
 758         buf[3] = byte(((in[6] >> 24) & 0x0f) | (in[7]<<4)&0xf0)
 759
 760         buf[2] = byte(in[7] >> 4)
 761         buf[1] = byte(in[7] >> 12)
 762         buf[0] = byte(in[7] >> 20)
 763
 764         return new(big.Int).SetBytes(buf[:])
 765 }