github.com/FISCO-BCOS/crypto@v0.0.0-20200202032121-bd8ab0b5d4f1/elliptic/p224.go (about)

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