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