...
Run Format

Source file src/crypto/elliptic/p256.go

Documentation: crypto/elliptic

     1  // Copyright 2013 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  // +build !amd64,!arm64
     6  
     7  package elliptic
     8  
     9  // This file contains a constant-time, 32-bit implementation of P256.
    10  
    11  import (
    12  	"math/big"
    13  )
    14  
    15  type p256Curve struct {
    16  	*CurveParams
    17  }
    18  
    19  var (
    20  	p256Params *CurveParams
    21  
    22  	// RInverse contains 1/R mod p - the inverse of the Montgomery constant
    23  	// (2**257).
    24  	p256RInverse *big.Int
    25  )
    26  
    27  func initP256() {
    28  	// See FIPS 186-3, section D.2.3
    29  	p256Params = &CurveParams{Name: "P-256"}
    30  	p256Params.P, _ = new(big.Int).SetString("115792089210356248762697446949407573530086143415290314195533631308867097853951", 10)
    31  	p256Params.N, _ = new(big.Int).SetString("115792089210356248762697446949407573529996955224135760342422259061068512044369", 10)
    32  	p256Params.B, _ = new(big.Int).SetString("5ac635d8aa3a93e7b3ebbd55769886bc651d06b0cc53b0f63bce3c3e27d2604b", 16)
    33  	p256Params.Gx, _ = new(big.Int).SetString("6b17d1f2e12c4247f8bce6e563a440f277037d812deb33a0f4a13945d898c296", 16)
    34  	p256Params.Gy, _ = new(big.Int).SetString("4fe342e2fe1a7f9b8ee7eb4a7c0f9e162bce33576b315ececbb6406837bf51f5", 16)
    35  	p256Params.BitSize = 256
    36  
    37  	p256RInverse, _ = new(big.Int).SetString("7fffffff00000001fffffffe8000000100000000ffffffff0000000180000000", 16)
    38  
    39  	// Arch-specific initialization, i.e. let a platform dynamically pick a P256 implementation
    40  	initP256Arch()
    41  }
    42  
    43  func (curve p256Curve) Params() *CurveParams {
    44  	return curve.CurveParams
    45  }
    46  
    47  // p256GetScalar endian-swaps the big-endian scalar value from in and writes it
    48  // to out. If the scalar is equal or greater than the order of the group, it's
    49  // reduced modulo that order.
    50  func p256GetScalar(out *[32]byte, in []byte) {
    51  	n := new(big.Int).SetBytes(in)
    52  	var scalarBytes []byte
    53  
    54  	if n.Cmp(p256Params.N) >= 0 {
    55  		n.Mod(n, p256Params.N)
    56  		scalarBytes = n.Bytes()
    57  	} else {
    58  		scalarBytes = in
    59  	}
    60  
    61  	for i, v := range scalarBytes {
    62  		out[len(scalarBytes)-(1+i)] = v
    63  	}
    64  }
    65  
    66  func (p256Curve) ScalarBaseMult(scalar []byte) (x, y *big.Int) {
    67  	var scalarReversed [32]byte
    68  	p256GetScalar(&scalarReversed, scalar)
    69  
    70  	var x1, y1, z1 [p256Limbs]uint32
    71  	p256ScalarBaseMult(&x1, &y1, &z1, &scalarReversed)
    72  	return p256ToAffine(&x1, &y1, &z1)
    73  }
    74  
    75  func (p256Curve) ScalarMult(bigX, bigY *big.Int, scalar []byte) (x, y *big.Int) {
    76  	var scalarReversed [32]byte
    77  	p256GetScalar(&scalarReversed, scalar)
    78  
    79  	var px, py, x1, y1, z1 [p256Limbs]uint32
    80  	p256FromBig(&px, bigX)
    81  	p256FromBig(&py, bigY)
    82  	p256ScalarMult(&x1, &y1, &z1, &px, &py, &scalarReversed)
    83  	return p256ToAffine(&x1, &y1, &z1)
    84  }
    85  
    86  // Field elements are represented as nine, unsigned 32-bit words.
    87  //
    88  // The value of an field element is:
    89  //   x[0] + (x[1] * 2**29) + (x[2] * 2**57) + ... + (x[8] * 2**228)
    90  //
    91  // That is, each limb is alternately 29 or 28-bits wide in little-endian
    92  // order.
    93  //
    94  // This means that a field element hits 2**257, rather than 2**256 as we would
    95  // like. A 28, 29, ... pattern would cause us to hit 2**256, but that causes
    96  // problems when multiplying as terms end up one bit short of a limb which
    97  // would require much bit-shifting to correct.
    98  //
    99  // Finally, the values stored in a field element are in Montgomery form. So the
   100  // value |y| is stored as (y*R) mod p, where p is the P-256 prime and R is
   101  // 2**257.
   102  
   103  const (
   104  	p256Limbs    = 9
   105  	bottom29Bits = 0x1fffffff
   106  )
   107  
   108  var (
   109  	// p256One is the number 1 as a field element.
   110  	p256One  = [p256Limbs]uint32{2, 0, 0, 0xffff800, 0x1fffffff, 0xfffffff, 0x1fbfffff, 0x1ffffff, 0}
   111  	p256Zero = [p256Limbs]uint32{0, 0, 0, 0, 0, 0, 0, 0, 0}
   112  	// p256P is the prime modulus as a field element.
   113  	p256P = [p256Limbs]uint32{0x1fffffff, 0xfffffff, 0x1fffffff, 0x3ff, 0, 0, 0x200000, 0xf000000, 0xfffffff}
   114  	// p2562P is the twice prime modulus as a field element.
   115  	p2562P = [p256Limbs]uint32{0x1ffffffe, 0xfffffff, 0x1fffffff, 0x7ff, 0, 0, 0x400000, 0xe000000, 0x1fffffff}
   116  )
   117  
   118  // p256Precomputed contains precomputed values to aid the calculation of scalar
   119  // multiples of the base point, G. It's actually two, equal length, tables
   120  // concatenated.
   121  //
   122  // The first table contains (x,y) field element pairs for 16 multiples of the
   123  // base point, G.
   124  //
   125  //   Index  |  Index (binary) | Value
   126  //       0  |           0000  | 0G (all zeros, omitted)
   127  //       1  |           0001  | G
   128  //       2  |           0010  | 2**64G
   129  //       3  |           0011  | 2**64G + G
   130  //       4  |           0100  | 2**128G
   131  //       5  |           0101  | 2**128G + G
   132  //       6  |           0110  | 2**128G + 2**64G
   133  //       7  |           0111  | 2**128G + 2**64G + G
   134  //       8  |           1000  | 2**192G
   135  //       9  |           1001  | 2**192G + G
   136  //      10  |           1010  | 2**192G + 2**64G
   137  //      11  |           1011  | 2**192G + 2**64G + G
   138  //      12  |           1100  | 2**192G + 2**128G
   139  //      13  |           1101  | 2**192G + 2**128G + G
   140  //      14  |           1110  | 2**192G + 2**128G + 2**64G
   141  //      15  |           1111  | 2**192G + 2**128G + 2**64G + G
   142  //
   143  // The second table follows the same style, but the terms are 2**32G,
   144  // 2**96G, 2**160G, 2**224G.
   145  //
   146  // This is ~2KB of data.
   147  var p256Precomputed = [p256Limbs * 2 * 15 * 2]uint32{
   148  	0x11522878, 0xe730d41, 0xdb60179, 0x4afe2ff, 0x12883add, 0xcaddd88, 0x119e7edc, 0xd4a6eab, 0x3120bee,
   149  	0x1d2aac15, 0xf25357c, 0x19e45cdd, 0x5c721d0, 0x1992c5a5, 0xa237487, 0x154ba21, 0x14b10bb, 0xae3fe3,
   150  	0xd41a576, 0x922fc51, 0x234994f, 0x60b60d3, 0x164586ae, 0xce95f18, 0x1fe49073, 0x3fa36cc, 0x5ebcd2c,
   151  	0xb402f2f, 0x15c70bf, 0x1561925c, 0x5a26704, 0xda91e90, 0xcdc1c7f, 0x1ea12446, 0xe1ade1e, 0xec91f22,
   152  	0x26f7778, 0x566847e, 0xa0bec9e, 0x234f453, 0x1a31f21a, 0xd85e75c, 0x56c7109, 0xa267a00, 0xb57c050,
   153  	0x98fb57, 0xaa837cc, 0x60c0792, 0xcfa5e19, 0x61bab9e, 0x589e39b, 0xa324c5, 0x7d6dee7, 0x2976e4b,
   154  	0x1fc4124a, 0xa8c244b, 0x1ce86762, 0xcd61c7e, 0x1831c8e0, 0x75774e1, 0x1d96a5a9, 0x843a649, 0xc3ab0fa,
   155  	0x6e2e7d5, 0x7673a2a, 0x178b65e8, 0x4003e9b, 0x1a1f11c2, 0x7816ea, 0xf643e11, 0x58c43df, 0xf423fc2,
   156  	0x19633ffa, 0x891f2b2, 0x123c231c, 0x46add8c, 0x54700dd, 0x59e2b17, 0x172db40f, 0x83e277d, 0xb0dd609,
   157  	0xfd1da12, 0x35c6e52, 0x19ede20c, 0xd19e0c0, 0x97d0f40, 0xb015b19, 0x449e3f5, 0xe10c9e, 0x33ab581,
   158  	0x56a67ab, 0x577734d, 0x1dddc062, 0xc57b10d, 0x149b39d, 0x26a9e7b, 0xc35df9f, 0x48764cd, 0x76dbcca,
   159  	0xca4b366, 0xe9303ab, 0x1a7480e7, 0x57e9e81, 0x1e13eb50, 0xf466cf3, 0x6f16b20, 0x4ba3173, 0xc168c33,
   160  	0x15cb5439, 0x6a38e11, 0x73658bd, 0xb29564f, 0x3f6dc5b, 0x53b97e, 0x1322c4c0, 0x65dd7ff, 0x3a1e4f6,
   161  	0x14e614aa, 0x9246317, 0x1bc83aca, 0xad97eed, 0xd38ce4a, 0xf82b006, 0x341f077, 0xa6add89, 0x4894acd,
   162  	0x9f162d5, 0xf8410ef, 0x1b266a56, 0xd7f223, 0x3e0cb92, 0xe39b672, 0x6a2901a, 0x69a8556, 0x7e7c0,
   163  	0x9b7d8d3, 0x309a80, 0x1ad05f7f, 0xc2fb5dd, 0xcbfd41d, 0x9ceb638, 0x1051825c, 0xda0cf5b, 0x812e881,
   164  	0x6f35669, 0x6a56f2c, 0x1df8d184, 0x345820, 0x1477d477, 0x1645db1, 0xbe80c51, 0xc22be3e, 0xe35e65a,
   165  	0x1aeb7aa0, 0xc375315, 0xf67bc99, 0x7fdd7b9, 0x191fc1be, 0x61235d, 0x2c184e9, 0x1c5a839, 0x47a1e26,
   166  	0xb7cb456, 0x93e225d, 0x14f3c6ed, 0xccc1ac9, 0x17fe37f3, 0x4988989, 0x1a90c502, 0x2f32042, 0xa17769b,
   167  	0xafd8c7c, 0x8191c6e, 0x1dcdb237, 0x16200c0, 0x107b32a1, 0x66c08db, 0x10d06a02, 0x3fc93, 0x5620023,
   168  	0x16722b27, 0x68b5c59, 0x270fcfc, 0xfad0ecc, 0xe5de1c2, 0xeab466b, 0x2fc513c, 0x407f75c, 0xbaab133,
   169  	0x9705fe9, 0xb88b8e7, 0x734c993, 0x1e1ff8f, 0x19156970, 0xabd0f00, 0x10469ea7, 0x3293ac0, 0xcdc98aa,
   170  	0x1d843fd, 0xe14bfe8, 0x15be825f, 0x8b5212, 0xeb3fb67, 0x81cbd29, 0xbc62f16, 0x2b6fcc7, 0xf5a4e29,
   171  	0x13560b66, 0xc0b6ac2, 0x51ae690, 0xd41e271, 0xf3e9bd4, 0x1d70aab, 0x1029f72, 0x73e1c35, 0xee70fbc,
   172  	0xad81baf, 0x9ecc49a, 0x86c741e, 0xfe6be30, 0x176752e7, 0x23d416, 0x1f83de85, 0x27de188, 0x66f70b8,
   173  	0x181cd51f, 0x96b6e4c, 0x188f2335, 0xa5df759, 0x17a77eb6, 0xfeb0e73, 0x154ae914, 0x2f3ec51, 0x3826b59,
   174  	0xb91f17d, 0x1c72949, 0x1362bf0a, 0xe23fddf, 0xa5614b0, 0xf7d8f, 0x79061, 0x823d9d2, 0x8213f39,
   175  	0x1128ae0b, 0xd095d05, 0xb85c0c2, 0x1ecb2ef, 0x24ddc84, 0xe35e901, 0x18411a4a, 0xf5ddc3d, 0x3786689,
   176  	0x52260e8, 0x5ae3564, 0x542b10d, 0x8d93a45, 0x19952aa4, 0x996cc41, 0x1051a729, 0x4be3499, 0x52b23aa,
   177  	0x109f307e, 0x6f5b6bb, 0x1f84e1e7, 0x77a0cfa, 0x10c4df3f, 0x25a02ea, 0xb048035, 0xe31de66, 0xc6ecaa3,
   178  	0x28ea335, 0x2886024, 0x1372f020, 0xf55d35, 0x15e4684c, 0xf2a9e17, 0x1a4a7529, 0xcb7beb1, 0xb2a78a1,
   179  	0x1ab21f1f, 0x6361ccf, 0x6c9179d, 0xb135627, 0x1267b974, 0x4408bad, 0x1cbff658, 0xe3d6511, 0xc7d76f,
   180  	0x1cc7a69, 0xe7ee31b, 0x54fab4f, 0x2b914f, 0x1ad27a30, 0xcd3579e, 0xc50124c, 0x50daa90, 0xb13f72,
   181  	0xb06aa75, 0x70f5cc6, 0x1649e5aa, 0x84a5312, 0x329043c, 0x41c4011, 0x13d32411, 0xb04a838, 0xd760d2d,
   182  	0x1713b532, 0xbaa0c03, 0x84022ab, 0x6bcf5c1, 0x2f45379, 0x18ae070, 0x18c9e11e, 0x20bca9a, 0x66f496b,
   183  	0x3eef294, 0x67500d2, 0xd7f613c, 0x2dbbeb, 0xb741038, 0xe04133f, 0x1582968d, 0xbe985f7, 0x1acbc1a,
   184  	0x1a6a939f, 0x33e50f6, 0xd665ed4, 0xb4b7bd6, 0x1e5a3799, 0x6b33847, 0x17fa56ff, 0x65ef930, 0x21dc4a,
   185  	0x2b37659, 0x450fe17, 0xb357b65, 0xdf5efac, 0x15397bef, 0x9d35a7f, 0x112ac15f, 0x624e62e, 0xa90ae2f,
   186  	0x107eecd2, 0x1f69bbe, 0x77d6bce, 0x5741394, 0x13c684fc, 0x950c910, 0x725522b, 0xdc78583, 0x40eeabb,
   187  	0x1fde328a, 0xbd61d96, 0xd28c387, 0x9e77d89, 0x12550c40, 0x759cb7d, 0x367ef34, 0xae2a960, 0x91b8bdc,
   188  	0x93462a9, 0xf469ef, 0xb2e9aef, 0xd2ca771, 0x54e1f42, 0x7aaa49, 0x6316abb, 0x2413c8e, 0x5425bf9,
   189  	0x1bed3e3a, 0xf272274, 0x1f5e7326, 0x6416517, 0xea27072, 0x9cedea7, 0x6e7633, 0x7c91952, 0xd806dce,
   190  	0x8e2a7e1, 0xe421e1a, 0x418c9e1, 0x1dbc890, 0x1b395c36, 0xa1dc175, 0x1dc4ef73, 0x8956f34, 0xe4b5cf2,
   191  	0x1b0d3a18, 0x3194a36, 0x6c2641f, 0xe44124c, 0xa2f4eaa, 0xa8c25ba, 0xf927ed7, 0x627b614, 0x7371cca,
   192  	0xba16694, 0x417bc03, 0x7c0a7e3, 0x9c35c19, 0x1168a205, 0x8b6b00d, 0x10e3edc9, 0x9c19bf2, 0x5882229,
   193  	0x1b2b4162, 0xa5cef1a, 0x1543622b, 0x9bd433e, 0x364e04d, 0x7480792, 0x5c9b5b3, 0xe85ff25, 0x408ef57,
   194  	0x1814cfa4, 0x121b41b, 0xd248a0f, 0x3b05222, 0x39bb16a, 0xc75966d, 0xa038113, 0xa4a1769, 0x11fbc6c,
   195  	0x917e50e, 0xeec3da8, 0x169d6eac, 0x10c1699, 0xa416153, 0xf724912, 0x15cd60b7, 0x4acbad9, 0x5efc5fa,
   196  	0xf150ed7, 0x122b51, 0x1104b40a, 0xcb7f442, 0xfbb28ff, 0x6ac53ca, 0x196142cc, 0x7bf0fa9, 0x957651,
   197  	0x4e0f215, 0xed439f8, 0x3f46bd5, 0x5ace82f, 0x110916b6, 0x6db078, 0xffd7d57, 0xf2ecaac, 0xca86dec,
   198  	0x15d6b2da, 0x965ecc9, 0x1c92b4c2, 0x1f3811, 0x1cb080f5, 0x2d8b804, 0x19d1c12d, 0xf20bd46, 0x1951fa7,
   199  	0xa3656c3, 0x523a425, 0xfcd0692, 0xd44ddc8, 0x131f0f5b, 0xaf80e4a, 0xcd9fc74, 0x99bb618, 0x2db944c,
   200  	0xa673090, 0x1c210e1, 0x178c8d23, 0x1474383, 0x10b8743d, 0x985a55b, 0x2e74779, 0x576138, 0x9587927,
   201  	0x133130fa, 0xbe05516, 0x9f4d619, 0xbb62570, 0x99ec591, 0xd9468fe, 0x1d07782d, 0xfc72e0b, 0x701b298,
   202  	0x1863863b, 0x85954b8, 0x121a0c36, 0x9e7fedf, 0xf64b429, 0x9b9d71e, 0x14e2f5d8, 0xf858d3a, 0x942eea8,
   203  	0xda5b765, 0x6edafff, 0xa9d18cc, 0xc65e4ba, 0x1c747e86, 0xe4ea915, 0x1981d7a1, 0x8395659, 0x52ed4e2,
   204  	0x87d43b7, 0x37ab11b, 0x19d292ce, 0xf8d4692, 0x18c3053f, 0x8863e13, 0x4c146c0, 0x6bdf55a, 0x4e4457d,
   205  	0x16152289, 0xac78ec2, 0x1a59c5a2, 0x2028b97, 0x71c2d01, 0x295851f, 0x404747b, 0x878558d, 0x7d29aa4,
   206  	0x13d8341f, 0x8daefd7, 0x139c972d, 0x6b7ea75, 0xd4a9dde, 0xff163d8, 0x81d55d7, 0xa5bef68, 0xb7b30d8,
   207  	0xbe73d6f, 0xaa88141, 0xd976c81, 0x7e7a9cc, 0x18beb771, 0xd773cbd, 0x13f51951, 0x9d0c177, 0x1c49a78,
   208  }
   209  
   210  // Field element operations:
   211  
   212  // nonZeroToAllOnes returns:
   213  //   0xffffffff for 0 < x <= 2**31
   214  //   0 for x == 0 or x > 2**31.
   215  func nonZeroToAllOnes(x uint32) uint32 {
   216  	return ((x - 1) >> 31) - 1
   217  }
   218  
   219  // p256ReduceCarry adds a multiple of p in order to cancel |carry|,
   220  // which is a term at 2**257.
   221  //
   222  // On entry: carry < 2**3, inout[0,2,...] < 2**29, inout[1,3,...] < 2**28.
   223  // On exit: inout[0,2,..] < 2**30, inout[1,3,...] < 2**29.
   224  func p256ReduceCarry(inout *[p256Limbs]uint32, carry uint32) {
   225  	carry_mask := nonZeroToAllOnes(carry)
   226  
   227  	inout[0] += carry << 1
   228  	inout[3] += 0x10000000 & carry_mask
   229  	// carry < 2**3 thus (carry << 11) < 2**14 and we added 2**28 in the
   230  	// previous line therefore this doesn't underflow.
   231  	inout[3] -= carry << 11
   232  	inout[4] += (0x20000000 - 1) & carry_mask
   233  	inout[5] += (0x10000000 - 1) & carry_mask
   234  	inout[6] += (0x20000000 - 1) & carry_mask
   235  	inout[6] -= carry << 22
   236  	// This may underflow if carry is non-zero but, if so, we'll fix it in the
   237  	// next line.
   238  	inout[7] -= 1 & carry_mask
   239  	inout[7] += carry << 25
   240  }
   241  
   242  // p256Sum sets out = in+in2.
   243  //
   244  // On entry, in[i]+in2[i] must not overflow a 32-bit word.
   245  // On exit: out[0,2,...] < 2**30, out[1,3,...] < 2**29
   246  func p256Sum(out, in, in2 *[p256Limbs]uint32) {
   247  	carry := uint32(0)
   248  	for i := 0; ; i++ {
   249  		out[i] = in[i] + in2[i]
   250  		out[i] += carry
   251  		carry = out[i] >> 29
   252  		out[i] &= bottom29Bits
   253  
   254  		i++
   255  		if i == p256Limbs {
   256  			break
   257  		}
   258  
   259  		out[i] = in[i] + in2[i]
   260  		out[i] += carry
   261  		carry = out[i] >> 28
   262  		out[i] &= bottom28Bits
   263  	}
   264  
   265  	p256ReduceCarry(out, carry)
   266  }
   267  
   268  const (
   269  	two30m2    = 1<<30 - 1<<2
   270  	two30p13m2 = 1<<30 + 1<<13 - 1<<2
   271  	two31m2    = 1<<31 - 1<<2
   272  	two31p24m2 = 1<<31 + 1<<24 - 1<<2
   273  	two30m27m2 = 1<<30 - 1<<27 - 1<<2
   274  )
   275  
   276  // p256Zero31 is 0 mod p.
   277  var p256Zero31 = [p256Limbs]uint32{two31m3, two30m2, two31m2, two30p13m2, two31m2, two30m2, two31p24m2, two30m27m2, two31m2}
   278  
   279  // p256Diff sets out = in-in2.
   280  //
   281  // On entry: in[0,2,...] < 2**30, in[1,3,...] < 2**29 and
   282  //           in2[0,2,...] < 2**30, in2[1,3,...] < 2**29.
   283  // On exit: out[0,2,...] < 2**30, out[1,3,...] < 2**29.
   284  func p256Diff(out, in, in2 *[p256Limbs]uint32) {
   285  	var carry uint32
   286  
   287  	for i := 0; ; i++ {
   288  		out[i] = in[i] - in2[i]
   289  		out[i] += p256Zero31[i]
   290  		out[i] += carry
   291  		carry = out[i] >> 29
   292  		out[i] &= bottom29Bits
   293  
   294  		i++
   295  		if i == p256Limbs {
   296  			break
   297  		}
   298  
   299  		out[i] = in[i] - in2[i]
   300  		out[i] += p256Zero31[i]
   301  		out[i] += carry
   302  		carry = out[i] >> 28
   303  		out[i] &= bottom28Bits
   304  	}
   305  
   306  	p256ReduceCarry(out, carry)
   307  }
   308  
   309  // p256ReduceDegree sets out = tmp/R mod p where tmp contains 64-bit words with
   310  // the same 29,28,... bit positions as an field element.
   311  //
   312  // The values in field elements are in Montgomery form: x*R mod p where R =
   313  // 2**257. Since we just multiplied two Montgomery values together, the result
   314  // is x*y*R*R mod p. We wish to divide by R in order for the result also to be
   315  // in Montgomery form.
   316  //
   317  // On entry: tmp[i] < 2**64
   318  // On exit: out[0,2,...] < 2**30, out[1,3,...] < 2**29
   319  func p256ReduceDegree(out *[p256Limbs]uint32, tmp [17]uint64) {
   320  	// The following table may be helpful when reading this code:
   321  	//
   322  	// Limb number:   0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10...
   323  	// Width (bits):  29| 28| 29| 28| 29| 28| 29| 28| 29| 28| 29
   324  	// Start bit:     0 | 29| 57| 86|114|143|171|200|228|257|285
   325  	//   (odd phase): 0 | 28| 57| 85|114|142|171|199|228|256|285
   326  	var tmp2 [18]uint32
   327  	var carry, x, xMask uint32
   328  
   329  	// tmp contains 64-bit words with the same 29,28,29-bit positions as an
   330  	// field element. So the top of an element of tmp might overlap with
   331  	// another element two positions down. The following loop eliminates
   332  	// this overlap.
   333  	tmp2[0] = uint32(tmp[0]) & bottom29Bits
   334  
   335  	tmp2[1] = uint32(tmp[0]) >> 29
   336  	tmp2[1] |= (uint32(tmp[0]>>32) << 3) & bottom28Bits
   337  	tmp2[1] += uint32(tmp[1]) & bottom28Bits
   338  	carry = tmp2[1] >> 28
   339  	tmp2[1] &= bottom28Bits
   340  
   341  	for i := 2; i < 17; i++ {
   342  		tmp2[i] = (uint32(tmp[i-2] >> 32)) >> 25
   343  		tmp2[i] += (uint32(tmp[i-1])) >> 28
   344  		tmp2[i] += (uint32(tmp[i-1]>>32) << 4) & bottom29Bits
   345  		tmp2[i] += uint32(tmp[i]) & bottom29Bits
   346  		tmp2[i] += carry
   347  		carry = tmp2[i] >> 29
   348  		tmp2[i] &= bottom29Bits
   349  
   350  		i++
   351  		if i == 17 {
   352  			break
   353  		}
   354  		tmp2[i] = uint32(tmp[i-2]>>32) >> 25
   355  		tmp2[i] += uint32(tmp[i-1]) >> 29
   356  		tmp2[i] += ((uint32(tmp[i-1] >> 32)) << 3) & bottom28Bits
   357  		tmp2[i] += uint32(tmp[i]) & bottom28Bits
   358  		tmp2[i] += carry
   359  		carry = tmp2[i] >> 28
   360  		tmp2[i] &= bottom28Bits
   361  	}
   362  
   363  	tmp2[17] = uint32(tmp[15]>>32) >> 25
   364  	tmp2[17] += uint32(tmp[16]) >> 29
   365  	tmp2[17] += uint32(tmp[16]>>32) << 3
   366  	tmp2[17] += carry
   367  
   368  	// Montgomery elimination of terms:
   369  	//
   370  	// Since R is 2**257, we can divide by R with a bitwise shift if we can
   371  	// ensure that the right-most 257 bits are all zero. We can make that true
   372  	// by adding multiplies of p without affecting the value.
   373  	//
   374  	// So we eliminate limbs from right to left. Since the bottom 29 bits of p
   375  	// are all ones, then by adding tmp2[0]*p to tmp2 we'll make tmp2[0] == 0.
   376  	// We can do that for 8 further limbs and then right shift to eliminate the
   377  	// extra factor of R.
   378  	for i := 0; ; i += 2 {
   379  		tmp2[i+1] += tmp2[i] >> 29
   380  		x = tmp2[i] & bottom29Bits
   381  		xMask = nonZeroToAllOnes(x)
   382  		tmp2[i] = 0
   383  
   384  		// The bounds calculations for this loop are tricky. Each iteration of
   385  		// the loop eliminates two words by adding values to words to their
   386  		// right.
   387  		//
   388  		// The following table contains the amounts added to each word (as an
   389  		// offset from the value of i at the top of the loop). The amounts are
   390  		// accounted for from the first and second half of the loop separately
   391  		// and are written as, for example, 28 to mean a value <2**28.
   392  		//
   393  		// Word:                   3   4   5   6   7   8   9   10
   394  		// Added in top half:     28  11      29  21  29  28
   395  		//                                        28  29
   396  		//                                            29
   397  		// Added in bottom half:      29  10      28  21  28   28
   398  		//                                            29
   399  		//
   400  		// The value that is currently offset 7 will be offset 5 for the next
   401  		// iteration and then offset 3 for the iteration after that. Therefore
   402  		// the total value added will be the values added at 7, 5 and 3.
   403  		//
   404  		// The following table accumulates these values. The sums at the bottom
   405  		// are written as, for example, 29+28, to mean a value < 2**29+2**28.
   406  		//
   407  		// Word:                   3   4   5   6   7   8   9  10  11  12  13
   408  		//                        28  11  10  29  21  29  28  28  28  28  28
   409  		//                            29  28  11  28  29  28  29  28  29  28
   410  		//                                    29  28  21  21  29  21  29  21
   411  		//                                        10  29  28  21  28  21  28
   412  		//                                        28  29  28  29  28  29  28
   413  		//                                            11  10  29  10  29  10
   414  		//                                            29  28  11  28  11
   415  		//                                                    29      29
   416  		//                        --------------------------------------------
   417  		//                                                30+ 31+ 30+ 31+ 30+
   418  		//                                                28+ 29+ 28+ 29+ 21+
   419  		//                                                21+ 28+ 21+ 28+ 10
   420  		//                                                10  21+ 10  21+
   421  		//                                                    11      11
   422  		//
   423  		// So the greatest amount is added to tmp2[10] and tmp2[12]. If
   424  		// tmp2[10/12] has an initial value of <2**29, then the maximum value
   425  		// will be < 2**31 + 2**30 + 2**28 + 2**21 + 2**11, which is < 2**32,
   426  		// as required.
   427  		tmp2[i+3] += (x << 10) & bottom28Bits
   428  		tmp2[i+4] += (x >> 18)
   429  
   430  		tmp2[i+6] += (x << 21) & bottom29Bits
   431  		tmp2[i+7] += x >> 8
   432  
   433  		// At position 200, which is the starting bit position for word 7, we
   434  		// have a factor of 0xf000000 = 2**28 - 2**24.
   435  		tmp2[i+7] += 0x10000000 & xMask
   436  		tmp2[i+8] += (x - 1) & xMask
   437  		tmp2[i+7] -= (x << 24) & bottom28Bits
   438  		tmp2[i+8] -= x >> 4
   439  
   440  		tmp2[i+8] += 0x20000000 & xMask
   441  		tmp2[i+8] -= x
   442  		tmp2[i+8] += (x << 28) & bottom29Bits
   443  		tmp2[i+9] += ((x >> 1) - 1) & xMask
   444  
   445  		if i+1 == p256Limbs {
   446  			break
   447  		}
   448  		tmp2[i+2] += tmp2[i+1] >> 28
   449  		x = tmp2[i+1] & bottom28Bits
   450  		xMask = nonZeroToAllOnes(x)
   451  		tmp2[i+1] = 0
   452  
   453  		tmp2[i+4] += (x << 11) & bottom29Bits
   454  		tmp2[i+5] += (x >> 18)
   455  
   456  		tmp2[i+7] += (x << 21) & bottom28Bits
   457  		tmp2[i+8] += x >> 7
   458  
   459  		// At position 199, which is the starting bit of the 8th word when
   460  		// dealing with a context starting on an odd word, we have a factor of
   461  		// 0x1e000000 = 2**29 - 2**25. Since we have not updated i, the 8th
   462  		// word from i+1 is i+8.
   463  		tmp2[i+8] += 0x20000000 & xMask
   464  		tmp2[i+9] += (x - 1) & xMask
   465  		tmp2[i+8] -= (x << 25) & bottom29Bits
   466  		tmp2[i+9] -= x >> 4
   467  
   468  		tmp2[i+9] += 0x10000000 & xMask
   469  		tmp2[i+9] -= x
   470  		tmp2[i+10] += (x - 1) & xMask
   471  	}
   472  
   473  	// We merge the right shift with a carry chain. The words above 2**257 have
   474  	// widths of 28,29,... which we need to correct when copying them down.
   475  	carry = 0
   476  	for i := 0; i < 8; i++ {
   477  		// The maximum value of tmp2[i + 9] occurs on the first iteration and
   478  		// is < 2**30+2**29+2**28. Adding 2**29 (from tmp2[i + 10]) is
   479  		// therefore safe.
   480  		out[i] = tmp2[i+9]
   481  		out[i] += carry
   482  		out[i] += (tmp2[i+10] << 28) & bottom29Bits
   483  		carry = out[i] >> 29
   484  		out[i] &= bottom29Bits
   485  
   486  		i++
   487  		out[i] = tmp2[i+9] >> 1
   488  		out[i] += carry
   489  		carry = out[i] >> 28
   490  		out[i] &= bottom28Bits
   491  	}
   492  
   493  	out[8] = tmp2[17]
   494  	out[8] += carry
   495  	carry = out[8] >> 29
   496  	out[8] &= bottom29Bits
   497  
   498  	p256ReduceCarry(out, carry)
   499  }
   500  
   501  // p256Square sets out=in*in.
   502  //
   503  // On entry: in[0,2,...] < 2**30, in[1,3,...] < 2**29.
   504  // On exit: out[0,2,...] < 2**30, out[1,3,...] < 2**29.
   505  func p256Square(out, in *[p256Limbs]uint32) {
   506  	var tmp [17]uint64
   507  
   508  	tmp[0] = uint64(in[0]) * uint64(in[0])
   509  	tmp[1] = uint64(in[0]) * (uint64(in[1]) << 1)
   510  	tmp[2] = uint64(in[0])*(uint64(in[2])<<1) +
   511  		uint64(in[1])*(uint64(in[1])<<1)
   512  	tmp[3] = uint64(in[0])*(uint64(in[3])<<1) +
   513  		uint64(in[1])*(uint64(in[2])<<1)
   514  	tmp[4] = uint64(in[0])*(uint64(in[4])<<1) +
   515  		uint64(in[1])*(uint64(in[3])<<2) +
   516  		uint64(in[2])*uint64(in[2])
   517  	tmp[5] = uint64(in[0])*(uint64(in[5])<<1) +
   518  		uint64(in[1])*(uint64(in[4])<<1) +
   519  		uint64(in[2])*(uint64(in[3])<<1)
   520  	tmp[6] = uint64(in[0])*(uint64(in[6])<<1) +
   521  		uint64(in[1])*(uint64(in[5])<<2) +
   522  		uint64(in[2])*(uint64(in[4])<<1) +
   523  		uint64(in[3])*(uint64(in[3])<<1)
   524  	tmp[7] = uint64(in[0])*(uint64(in[7])<<1) +
   525  		uint64(in[1])*(uint64(in[6])<<1) +
   526  		uint64(in[2])*(uint64(in[5])<<1) +
   527  		uint64(in[3])*(uint64(in[4])<<1)
   528  	// tmp[8] has the greatest value of 2**61 + 2**60 + 2**61 + 2**60 + 2**60,
   529  	// which is < 2**64 as required.
   530  	tmp[8] = uint64(in[0])*(uint64(in[8])<<1) +
   531  		uint64(in[1])*(uint64(in[7])<<2) +
   532  		uint64(in[2])*(uint64(in[6])<<1) +
   533  		uint64(in[3])*(uint64(in[5])<<2) +
   534  		uint64(in[4])*uint64(in[4])
   535  	tmp[9] = uint64(in[1])*(uint64(in[8])<<1) +
   536  		uint64(in[2])*(uint64(in[7])<<1) +
   537  		uint64(in[3])*(uint64(in[6])<<1) +
   538  		uint64(in[4])*(uint64(in[5])<<1)
   539  	tmp[10] = uint64(in[2])*(uint64(in[8])<<1) +
   540  		uint64(in[3])*(uint64(in[7])<<2) +
   541  		uint64(in[4])*(uint64(in[6])<<1) +
   542  		uint64(in[5])*(uint64(in[5])<<1)
   543  	tmp[11] = uint64(in[3])*(uint64(in[8])<<1) +
   544  		uint64(in[4])*(uint64(in[7])<<1) +
   545  		uint64(in[5])*(uint64(in[6])<<1)
   546  	tmp[12] = uint64(in[4])*(uint64(in[8])<<1) +
   547  		uint64(in[5])*(uint64(in[7])<<2) +
   548  		uint64(in[6])*uint64(in[6])
   549  	tmp[13] = uint64(in[5])*(uint64(in[8])<<1) +
   550  		uint64(in[6])*(uint64(in[7])<<1)
   551  	tmp[14] = uint64(in[6])*(uint64(in[8])<<1) +
   552  		uint64(in[7])*(uint64(in[7])<<1)
   553  	tmp[15] = uint64(in[7]) * (uint64(in[8]) << 1)
   554  	tmp[16] = uint64(in[8]) * uint64(in[8])
   555  
   556  	p256ReduceDegree(out, tmp)
   557  }
   558  
   559  // p256Mul sets out=in*in2.
   560  //
   561  // On entry: in[0,2,...] < 2**30, in[1,3,...] < 2**29 and
   562  //           in2[0,2,...] < 2**30, in2[1,3,...] < 2**29.
   563  // On exit: out[0,2,...] < 2**30, out[1,3,...] < 2**29.
   564  func p256Mul(out, in, in2 *[p256Limbs]uint32) {
   565  	var tmp [17]uint64
   566  
   567  	tmp[0] = uint64(in[0]) * uint64(in2[0])
   568  	tmp[1] = uint64(in[0])*(uint64(in2[1])<<0) +
   569  		uint64(in[1])*(uint64(in2[0])<<0)
   570  	tmp[2] = uint64(in[0])*(uint64(in2[2])<<0) +
   571  		uint64(in[1])*(uint64(in2[1])<<1) +
   572  		uint64(in[2])*(uint64(in2[0])<<0)
   573  	tmp[3] = uint64(in[0])*(uint64(in2[3])<<0) +
   574  		uint64(in[1])*(uint64(in2[2])<<0) +
   575  		uint64(in[2])*(uint64(in2[1])<<0) +
   576  		uint64(in[3])*(uint64(in2[0])<<0)
   577  	tmp[4] = uint64(in[0])*(uint64(in2[4])<<0) +
   578  		uint64(in[1])*(uint64(in2[3])<<1) +
   579  		uint64(in[2])*(uint64(in2[2])<<0) +
   580  		uint64(in[3])*(uint64(in2[1])<<1) +
   581  		uint64(in[4])*(uint64(in2[0])<<0)
   582  	tmp[5] = uint64(in[0])*(uint64(in2[5])<<0) +
   583  		uint64(in[1])*(uint64(in2[4])<<0) +
   584  		uint64(in[2])*(uint64(in2[3])<<0) +
   585  		uint64(in[3])*(uint64(in2[2])<<0) +
   586  		uint64(in[4])*(uint64(in2[1])<<0) +
   587  		uint64(in[5])*(uint64(in2[0])<<0)
   588  	tmp[6] = uint64(in[0])*(uint64(in2[6])<<0) +
   589  		uint64(in[1])*(uint64(in2[5])<<1) +
   590  		uint64(in[2])*(uint64(in2[4])<<0) +
   591  		uint64(in[3])*(uint64(in2[3])<<1) +
   592  		uint64(in[4])*(uint64(in2[2])<<0) +
   593  		uint64(in[5])*(uint64(in2[1])<<1) +
   594  		uint64(in[6])*(uint64(in2[0])<<0)
   595  	tmp[7] = uint64(in[0])*(uint64(in2[7])<<0) +
   596  		uint64(in[1])*(uint64(in2[6])<<0) +
   597  		uint64(in[2])*(uint64(in2[5])<<0) +
   598  		uint64(in[3])*(uint64(in2[4])<<0) +
   599  		uint64(in[4])*(uint64(in2[3])<<0) +
   600  		uint64(in[5])*(uint64(in2[2])<<0) +
   601  		uint64(in[6])*(uint64(in2[1])<<0) +
   602  		uint64(in[7])*(uint64(in2[0])<<0)
   603  	// tmp[8] has the greatest value but doesn't overflow. See logic in
   604  	// p256Square.
   605  	tmp[8] = uint64(in[0])*(uint64(in2[8])<<0) +
   606  		uint64(in[1])*(uint64(in2[7])<<1) +
   607  		uint64(in[2])*(uint64(in2[6])<<0) +
   608  		uint64(in[3])*(uint64(in2[5])<<1) +
   609  		uint64(in[4])*(uint64(in2[4])<<0) +
   610  		uint64(in[5])*(uint64(in2[3])<<1) +
   611  		uint64(in[6])*(uint64(in2[2])<<0) +
   612  		uint64(in[7])*(uint64(in2[1])<<1) +
   613  		uint64(in[8])*(uint64(in2[0])<<0)
   614  	tmp[9] = uint64(in[1])*(uint64(in2[8])<<0) +
   615  		uint64(in[2])*(uint64(in2[7])<<0) +
   616  		uint64(in[3])*(uint64(in2[6])<<0) +
   617  		uint64(in[4])*(uint64(in2[5])<<0) +
   618  		uint64(in[5])*(uint64(in2[4])<<0) +
   619  		uint64(in[6])*(uint64(in2[3])<<0) +
   620  		uint64(in[7])*(uint64(in2[2])<<0) +
   621  		uint64(in[8])*(uint64(in2[1])<<0)
   622  	tmp[10] = uint64(in[2])*(uint64(in2[8])<<0) +
   623  		uint64(in[3])*(uint64(in2[7])<<1) +
   624  		uint64(in[4])*(uint64(in2[6])<<0) +
   625  		uint64(in[5])*(uint64(in2[5])<<1) +
   626  		uint64(in[6])*(uint64(in2[4])<<0) +
   627  		uint64(in[7])*(uint64(in2[3])<<1) +
   628  		uint64(in[8])*(uint64(in2[2])<<0)
   629  	tmp[11] = uint64(in[3])*(uint64(in2[8])<<0) +
   630  		uint64(in[4])*(uint64(in2[7])<<0) +
   631  		uint64(in[5])*(uint64(in2[6])<<0) +
   632  		uint64(in[6])*(uint64(in2[5])<<0) +
   633  		uint64(in[7])*(uint64(in2[4])<<0) +
   634  		uint64(in[8])*(uint64(in2[3])<<0)
   635  	tmp[12] = uint64(in[4])*(uint64(in2[8])<<0) +
   636  		uint64(in[5])*(uint64(in2[7])<<1) +
   637  		uint64(in[6])*(uint64(in2[6])<<0) +
   638  		uint64(in[7])*(uint64(in2[5])<<1) +
   639  		uint64(in[8])*(uint64(in2[4])<<0)
   640  	tmp[13] = uint64(in[5])*(uint64(in2[8])<<0) +
   641  		uint64(in[6])*(uint64(in2[7])<<0) +
   642  		uint64(in[7])*(uint64(in2[6])<<0) +
   643  		uint64(in[8])*(uint64(in2[5])<<0)
   644  	tmp[14] = uint64(in[6])*(uint64(in2[8])<<0) +
   645  		uint64(in[7])*(uint64(in2[7])<<1) +
   646  		uint64(in[8])*(uint64(in2[6])<<0)
   647  	tmp[15] = uint64(in[7])*(uint64(in2[8])<<0) +
   648  		uint64(in[8])*(uint64(in2[7])<<0)
   649  	tmp[16] = uint64(in[8]) * (uint64(in2[8]) << 0)
   650  
   651  	p256ReduceDegree(out, tmp)
   652  }
   653  
   654  func p256Assign(out, in *[p256Limbs]uint32) {
   655  	*out = *in
   656  }
   657  
   658  // p256Invert calculates |out| = |in|^{-1}
   659  //
   660  // Based on Fermat's Little Theorem:
   661  //   a^p = a (mod p)
   662  //   a^{p-1} = 1 (mod p)
   663  //   a^{p-2} = a^{-1} (mod p)
   664  func p256Invert(out, in *[p256Limbs]uint32) {
   665  	var ftmp, ftmp2 [p256Limbs]uint32
   666  
   667  	// each e_I will hold |in|^{2^I - 1}
   668  	var e2, e4, e8, e16, e32, e64 [p256Limbs]uint32
   669  
   670  	p256Square(&ftmp, in)     // 2^1
   671  	p256Mul(&ftmp, in, &ftmp) // 2^2 - 2^0
   672  	p256Assign(&e2, &ftmp)
   673  	p256Square(&ftmp, &ftmp)   // 2^3 - 2^1
   674  	p256Square(&ftmp, &ftmp)   // 2^4 - 2^2
   675  	p256Mul(&ftmp, &ftmp, &e2) // 2^4 - 2^0
   676  	p256Assign(&e4, &ftmp)
   677  	p256Square(&ftmp, &ftmp)   // 2^5 - 2^1
   678  	p256Square(&ftmp, &ftmp)   // 2^6 - 2^2
   679  	p256Square(&ftmp, &ftmp)   // 2^7 - 2^3
   680  	p256Square(&ftmp, &ftmp)   // 2^8 - 2^4
   681  	p256Mul(&ftmp, &ftmp, &e4) // 2^8 - 2^0
   682  	p256Assign(&e8, &ftmp)
   683  	for i := 0; i < 8; i++ {
   684  		p256Square(&ftmp, &ftmp)
   685  	} // 2^16 - 2^8
   686  	p256Mul(&ftmp, &ftmp, &e8) // 2^16 - 2^0
   687  	p256Assign(&e16, &ftmp)
   688  	for i := 0; i < 16; i++ {
   689  		p256Square(&ftmp, &ftmp)
   690  	} // 2^32 - 2^16
   691  	p256Mul(&ftmp, &ftmp, &e16) // 2^32 - 2^0
   692  	p256Assign(&e32, &ftmp)
   693  	for i := 0; i < 32; i++ {
   694  		p256Square(&ftmp, &ftmp)
   695  	} // 2^64 - 2^32
   696  	p256Assign(&e64, &ftmp)
   697  	p256Mul(&ftmp, &ftmp, in) // 2^64 - 2^32 + 2^0
   698  	for i := 0; i < 192; i++ {
   699  		p256Square(&ftmp, &ftmp)
   700  	} // 2^256 - 2^224 + 2^192
   701  
   702  	p256Mul(&ftmp2, &e64, &e32) // 2^64 - 2^0
   703  	for i := 0; i < 16; i++ {
   704  		p256Square(&ftmp2, &ftmp2)
   705  	} // 2^80 - 2^16
   706  	p256Mul(&ftmp2, &ftmp2, &e16) // 2^80 - 2^0
   707  	for i := 0; i < 8; i++ {
   708  		p256Square(&ftmp2, &ftmp2)
   709  	} // 2^88 - 2^8
   710  	p256Mul(&ftmp2, &ftmp2, &e8) // 2^88 - 2^0
   711  	for i := 0; i < 4; i++ {
   712  		p256Square(&ftmp2, &ftmp2)
   713  	} // 2^92 - 2^4
   714  	p256Mul(&ftmp2, &ftmp2, &e4) // 2^92 - 2^0
   715  	p256Square(&ftmp2, &ftmp2)   // 2^93 - 2^1
   716  	p256Square(&ftmp2, &ftmp2)   // 2^94 - 2^2
   717  	p256Mul(&ftmp2, &ftmp2, &e2) // 2^94 - 2^0
   718  	p256Square(&ftmp2, &ftmp2)   // 2^95 - 2^1
   719  	p256Square(&ftmp2, &ftmp2)   // 2^96 - 2^2
   720  	p256Mul(&ftmp2, &ftmp2, in)  // 2^96 - 3
   721  
   722  	p256Mul(out, &ftmp2, &ftmp) // 2^256 - 2^224 + 2^192 + 2^96 - 3
   723  }
   724  
   725  // p256Scalar3 sets out=3*out.
   726  //
   727  // On entry: out[0,2,...] < 2**30, out[1,3,...] < 2**29.
   728  // On exit: out[0,2,...] < 2**30, out[1,3,...] < 2**29.
   729  func p256Scalar3(out *[p256Limbs]uint32) {
   730  	var carry uint32
   731  
   732  	for i := 0; ; i++ {
   733  		out[i] *= 3
   734  		out[i] += carry
   735  		carry = out[i] >> 29
   736  		out[i] &= bottom29Bits
   737  
   738  		i++
   739  		if i == p256Limbs {
   740  			break
   741  		}
   742  
   743  		out[i] *= 3
   744  		out[i] += carry
   745  		carry = out[i] >> 28
   746  		out[i] &= bottom28Bits
   747  	}
   748  
   749  	p256ReduceCarry(out, carry)
   750  }
   751  
   752  // p256Scalar4 sets out=4*out.
   753  //
   754  // On entry: out[0,2,...] < 2**30, out[1,3,...] < 2**29.
   755  // On exit: out[0,2,...] < 2**30, out[1,3,...] < 2**29.
   756  func p256Scalar4(out *[p256Limbs]uint32) {
   757  	var carry, nextCarry uint32
   758  
   759  	for i := 0; ; i++ {
   760  		nextCarry = out[i] >> 27
   761  		out[i] <<= 2
   762  		out[i] &= bottom29Bits
   763  		out[i] += carry
   764  		carry = nextCarry + (out[i] >> 29)
   765  		out[i] &= bottom29Bits
   766  
   767  		i++
   768  		if i == p256Limbs {
   769  			break
   770  		}
   771  		nextCarry = out[i] >> 26
   772  		out[i] <<= 2
   773  		out[i] &= bottom28Bits
   774  		out[i] += carry
   775  		carry = nextCarry + (out[i] >> 28)
   776  		out[i] &= bottom28Bits
   777  	}
   778  
   779  	p256ReduceCarry(out, carry)
   780  }
   781  
   782  // p256Scalar8 sets out=8*out.
   783  //
   784  // On entry: out[0,2,...] < 2**30, out[1,3,...] < 2**29.
   785  // On exit: out[0,2,...] < 2**30, out[1,3,...] < 2**29.
   786  func p256Scalar8(out *[p256Limbs]uint32) {
   787  	var carry, nextCarry uint32
   788  
   789  	for i := 0; ; i++ {
   790  		nextCarry = out[i] >> 26
   791  		out[i] <<= 3
   792  		out[i] &= bottom29Bits
   793  		out[i] += carry
   794  		carry = nextCarry + (out[i] >> 29)
   795  		out[i] &= bottom29Bits
   796  
   797  		i++
   798  		if i == p256Limbs {
   799  			break
   800  		}
   801  		nextCarry = out[i] >> 25
   802  		out[i] <<= 3
   803  		out[i] &= bottom28Bits
   804  		out[i] += carry
   805  		carry = nextCarry + (out[i] >> 28)
   806  		out[i] &= bottom28Bits
   807  	}
   808  
   809  	p256ReduceCarry(out, carry)
   810  }
   811  
   812  // Group operations:
   813  //
   814  // Elements of the elliptic curve group are represented in Jacobian
   815  // coordinates: (x, y, z). An affine point (x', y') is x'=x/z**2, y'=y/z**3 in
   816  // Jacobian form.
   817  
   818  // p256PointDouble sets {xOut,yOut,zOut} = 2*{x,y,z}.
   819  //
   820  // See https://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#doubling-dbl-2009-l
   821  func p256PointDouble(xOut, yOut, zOut, x, y, z *[p256Limbs]uint32) {
   822  	var delta, gamma, alpha, beta, tmp, tmp2 [p256Limbs]uint32
   823  
   824  	p256Square(&delta, z)
   825  	p256Square(&gamma, y)
   826  	p256Mul(&beta, x, &gamma)
   827  
   828  	p256Sum(&tmp, x, &delta)
   829  	p256Diff(&tmp2, x, &delta)
   830  	p256Mul(&alpha, &tmp, &tmp2)
   831  	p256Scalar3(&alpha)
   832  
   833  	p256Sum(&tmp, y, z)
   834  	p256Square(&tmp, &tmp)
   835  	p256Diff(&tmp, &tmp, &gamma)
   836  	p256Diff(zOut, &tmp, &delta)
   837  
   838  	p256Scalar4(&beta)
   839  	p256Square(xOut, &alpha)
   840  	p256Diff(xOut, xOut, &beta)
   841  	p256Diff(xOut, xOut, &beta)
   842  
   843  	p256Diff(&tmp, &beta, xOut)
   844  	p256Mul(&tmp, &alpha, &tmp)
   845  	p256Square(&tmp2, &gamma)
   846  	p256Scalar8(&tmp2)
   847  	p256Diff(yOut, &tmp, &tmp2)
   848  }
   849  
   850  // p256PointAddMixed sets {xOut,yOut,zOut} = {x1,y1,z1} + {x2,y2,1}.
   851  // (i.e. the second point is affine.)
   852  //
   853  // See https://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#addition-add-2007-bl
   854  //
   855  // Note that this function does not handle P+P, infinity+P nor P+infinity
   856  // correctly.
   857  func p256PointAddMixed(xOut, yOut, zOut, x1, y1, z1, x2, y2 *[p256Limbs]uint32) {
   858  	var z1z1, z1z1z1, s2, u2, h, i, j, r, rr, v, tmp [p256Limbs]uint32
   859  
   860  	p256Square(&z1z1, z1)
   861  	p256Sum(&tmp, z1, z1)
   862  
   863  	p256Mul(&u2, x2, &z1z1)
   864  	p256Mul(&z1z1z1, z1, &z1z1)
   865  	p256Mul(&s2, y2, &z1z1z1)
   866  	p256Diff(&h, &u2, x1)
   867  	p256Sum(&i, &h, &h)
   868  	p256Square(&i, &i)
   869  	p256Mul(&j, &h, &i)
   870  	p256Diff(&r, &s2, y1)
   871  	p256Sum(&r, &r, &r)
   872  	p256Mul(&v, x1, &i)
   873  
   874  	p256Mul(zOut, &tmp, &h)
   875  	p256Square(&rr, &r)
   876  	p256Diff(xOut, &rr, &j)
   877  	p256Diff(xOut, xOut, &v)
   878  	p256Diff(xOut, xOut, &v)
   879  
   880  	p256Diff(&tmp, &v, xOut)
   881  	p256Mul(yOut, &tmp, &r)
   882  	p256Mul(&tmp, y1, &j)
   883  	p256Diff(yOut, yOut, &tmp)
   884  	p256Diff(yOut, yOut, &tmp)
   885  }
   886  
   887  // p256PointAdd sets {xOut,yOut,zOut} = {x1,y1,z1} + {x2,y2,z2}.
   888  //
   889  // See https://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#addition-add-2007-bl
   890  //
   891  // Note that this function does not handle P+P, infinity+P nor P+infinity
   892  // correctly.
   893  func p256PointAdd(xOut, yOut, zOut, x1, y1, z1, x2, y2, z2 *[p256Limbs]uint32) {
   894  	var z1z1, z1z1z1, z2z2, z2z2z2, s1, s2, u1, u2, h, i, j, r, rr, v, tmp [p256Limbs]uint32
   895  
   896  	p256Square(&z1z1, z1)
   897  	p256Square(&z2z2, z2)
   898  	p256Mul(&u1, x1, &z2z2)
   899  
   900  	p256Sum(&tmp, z1, z2)
   901  	p256Square(&tmp, &tmp)
   902  	p256Diff(&tmp, &tmp, &z1z1)
   903  	p256Diff(&tmp, &tmp, &z2z2)
   904  
   905  	p256Mul(&z2z2z2, z2, &z2z2)
   906  	p256Mul(&s1, y1, &z2z2z2)
   907  
   908  	p256Mul(&u2, x2, &z1z1)
   909  	p256Mul(&z1z1z1, z1, &z1z1)
   910  	p256Mul(&s2, y2, &z1z1z1)
   911  	p256Diff(&h, &u2, &u1)
   912  	p256Sum(&i, &h, &h)
   913  	p256Square(&i, &i)
   914  	p256Mul(&j, &h, &i)
   915  	p256Diff(&r, &s2, &s1)
   916  	p256Sum(&r, &r, &r)
   917  	p256Mul(&v, &u1, &i)
   918  
   919  	p256Mul(zOut, &tmp, &h)
   920  	p256Square(&rr, &r)
   921  	p256Diff(xOut, &rr, &j)
   922  	p256Diff(xOut, xOut, &v)
   923  	p256Diff(xOut, xOut, &v)
   924  
   925  	p256Diff(&tmp, &v, xOut)
   926  	p256Mul(yOut, &tmp, &r)
   927  	p256Mul(&tmp, &s1, &j)
   928  	p256Diff(yOut, yOut, &tmp)
   929  	p256Diff(yOut, yOut, &tmp)
   930  }
   931  
   932  // p256CopyConditional sets out=in if mask = 0xffffffff in constant time.
   933  //
   934  // On entry: mask is either 0 or 0xffffffff.
   935  func p256CopyConditional(out, in *[p256Limbs]uint32, mask uint32) {
   936  	for i := 0; i < p256Limbs; i++ {
   937  		tmp := mask & (in[i] ^ out[i])
   938  		out[i] ^= tmp
   939  	}
   940  }
   941  
   942  // p256SelectAffinePoint sets {out_x,out_y} to the index'th entry of table.
   943  // On entry: index < 16, table[0] must be zero.
   944  func p256SelectAffinePoint(xOut, yOut *[p256Limbs]uint32, table []uint32, index uint32) {
   945  	for i := range xOut {
   946  		xOut[i] = 0
   947  	}
   948  	for i := range yOut {
   949  		yOut[i] = 0
   950  	}
   951  
   952  	for i := uint32(1); i < 16; i++ {
   953  		mask := i ^ index
   954  		mask |= mask >> 2
   955  		mask |= mask >> 1
   956  		mask &= 1
   957  		mask--
   958  		for j := range xOut {
   959  			xOut[j] |= table[0] & mask
   960  			table = table[1:]
   961  		}
   962  		for j := range yOut {
   963  			yOut[j] |= table[0] & mask
   964  			table = table[1:]
   965  		}
   966  	}
   967  }
   968  
   969  // p256SelectJacobianPoint sets {out_x,out_y,out_z} to the index'th entry of
   970  // table.
   971  // On entry: index < 16, table[0] must be zero.
   972  func p256SelectJacobianPoint(xOut, yOut, zOut *[p256Limbs]uint32, table *[16][3][p256Limbs]uint32, index uint32) {
   973  	for i := range xOut {
   974  		xOut[i] = 0
   975  	}
   976  	for i := range yOut {
   977  		yOut[i] = 0
   978  	}
   979  	for i := range zOut {
   980  		zOut[i] = 0
   981  	}
   982  
   983  	// The implicit value at index 0 is all zero. We don't need to perform that
   984  	// iteration of the loop because we already set out_* to zero.
   985  	for i := uint32(1); i < 16; i++ {
   986  		mask := i ^ index
   987  		mask |= mask >> 2
   988  		mask |= mask >> 1
   989  		mask &= 1
   990  		mask--
   991  		for j := range xOut {
   992  			xOut[j] |= table[i][0][j] & mask
   993  		}
   994  		for j := range yOut {
   995  			yOut[j] |= table[i][1][j] & mask
   996  		}
   997  		for j := range zOut {
   998  			zOut[j] |= table[i][2][j] & mask
   999  		}
  1000  	}
  1001  }
  1002  
  1003  // p256GetBit returns the bit'th bit of scalar.
  1004  func p256GetBit(scalar *[32]uint8, bit uint) uint32 {
  1005  	return uint32(((scalar[bit>>3]) >> (bit & 7)) & 1)
  1006  }
  1007  
  1008  // p256ScalarBaseMult sets {xOut,yOut,zOut} = scalar*G where scalar is a
  1009  // little-endian number. Note that the value of scalar must be less than the
  1010  // order of the group.
  1011  func p256ScalarBaseMult(xOut, yOut, zOut *[p256Limbs]uint32, scalar *[32]uint8) {
  1012  	nIsInfinityMask := ^uint32(0)
  1013  	var pIsNoninfiniteMask, mask, tableOffset uint32
  1014  	var px, py, tx, ty, tz [p256Limbs]uint32
  1015  
  1016  	for i := range xOut {
  1017  		xOut[i] = 0
  1018  	}
  1019  	for i := range yOut {
  1020  		yOut[i] = 0
  1021  	}
  1022  	for i := range zOut {
  1023  		zOut[i] = 0
  1024  	}
  1025  
  1026  	// The loop adds bits at positions 0, 64, 128 and 192, followed by
  1027  	// positions 32,96,160 and 224 and does this 32 times.
  1028  	for i := uint(0); i < 32; i++ {
  1029  		if i != 0 {
  1030  			p256PointDouble(xOut, yOut, zOut, xOut, yOut, zOut)
  1031  		}
  1032  		tableOffset = 0
  1033  		for j := uint(0); j <= 32; j += 32 {
  1034  			bit0 := p256GetBit(scalar, 31-i+j)
  1035  			bit1 := p256GetBit(scalar, 95-i+j)
  1036  			bit2 := p256GetBit(scalar, 159-i+j)
  1037  			bit3 := p256GetBit(scalar, 223-i+j)
  1038  			index := bit0 | (bit1 << 1) | (bit2 << 2) | (bit3 << 3)
  1039  
  1040  			p256SelectAffinePoint(&px, &py, p256Precomputed[tableOffset:], index)
  1041  			tableOffset += 30 * p256Limbs
  1042  
  1043  			// Since scalar is less than the order of the group, we know that
  1044  			// {xOut,yOut,zOut} != {px,py,1}, unless both are zero, which we handle
  1045  			// below.
  1046  			p256PointAddMixed(&tx, &ty, &tz, xOut, yOut, zOut, &px, &py)
  1047  			// The result of pointAddMixed is incorrect if {xOut,yOut,zOut} is zero
  1048  			// (a.k.a.  the point at infinity). We handle that situation by
  1049  			// copying the point from the table.
  1050  			p256CopyConditional(xOut, &px, nIsInfinityMask)
  1051  			p256CopyConditional(yOut, &py, nIsInfinityMask)
  1052  			p256CopyConditional(zOut, &p256One, nIsInfinityMask)
  1053  
  1054  			// Equally, the result is also wrong if the point from the table is
  1055  			// zero, which happens when the index is zero. We handle that by
  1056  			// only copying from {tx,ty,tz} to {xOut,yOut,zOut} if index != 0.
  1057  			pIsNoninfiniteMask = nonZeroToAllOnes(index)
  1058  			mask = pIsNoninfiniteMask & ^nIsInfinityMask
  1059  			p256CopyConditional(xOut, &tx, mask)
  1060  			p256CopyConditional(yOut, &ty, mask)
  1061  			p256CopyConditional(zOut, &tz, mask)
  1062  			// If p was not zero, then n is now non-zero.
  1063  			nIsInfinityMask &^= pIsNoninfiniteMask
  1064  		}
  1065  	}
  1066  }
  1067  
  1068  // p256PointToAffine converts a Jacobian point to an affine point. If the input
  1069  // is the point at infinity then it returns (0, 0) in constant time.
  1070  func p256PointToAffine(xOut, yOut, x, y, z *[p256Limbs]uint32) {
  1071  	var zInv, zInvSq [p256Limbs]uint32
  1072  
  1073  	p256Invert(&zInv, z)
  1074  	p256Square(&zInvSq, &zInv)
  1075  	p256Mul(xOut, x, &zInvSq)
  1076  	p256Mul(&zInv, &zInv, &zInvSq)
  1077  	p256Mul(yOut, y, &zInv)
  1078  }
  1079  
  1080  // p256ToAffine returns a pair of *big.Int containing the affine representation
  1081  // of {x,y,z}.
  1082  func p256ToAffine(x, y, z *[p256Limbs]uint32) (xOut, yOut *big.Int) {
  1083  	var xx, yy [p256Limbs]uint32
  1084  	p256PointToAffine(&xx, &yy, x, y, z)
  1085  	return p256ToBig(&xx), p256ToBig(&yy)
  1086  }
  1087  
  1088  // p256ScalarMult sets {xOut,yOut,zOut} = scalar*{x,y}.
  1089  func p256ScalarMult(xOut, yOut, zOut, x, y *[p256Limbs]uint32, scalar *[32]uint8) {
  1090  	var px, py, pz, tx, ty, tz [p256Limbs]uint32
  1091  	var precomp [16][3][p256Limbs]uint32
  1092  	var nIsInfinityMask, index, pIsNoninfiniteMask, mask uint32
  1093  
  1094  	// We precompute 0,1,2,... times {x,y}.
  1095  	precomp[1][0] = *x
  1096  	precomp[1][1] = *y
  1097  	precomp[1][2] = p256One
  1098  
  1099  	for i := 2; i < 16; i += 2 {
  1100  		p256PointDouble(&precomp[i][0], &precomp[i][1], &precomp[i][2], &precomp[i/2][0], &precomp[i/2][1], &precomp[i/2][2])
  1101  		p256PointAddMixed(&precomp[i+1][0], &precomp[i+1][1], &precomp[i+1][2], &precomp[i][0], &precomp[i][1], &precomp[i][2], x, y)
  1102  	}
  1103  
  1104  	for i := range xOut {
  1105  		xOut[i] = 0
  1106  	}
  1107  	for i := range yOut {
  1108  		yOut[i] = 0
  1109  	}
  1110  	for i := range zOut {
  1111  		zOut[i] = 0
  1112  	}
  1113  	nIsInfinityMask = ^uint32(0)
  1114  
  1115  	// We add in a window of four bits each iteration and do this 64 times.
  1116  	for i := 0; i < 64; i++ {
  1117  		if i != 0 {
  1118  			p256PointDouble(xOut, yOut, zOut, xOut, yOut, zOut)
  1119  			p256PointDouble(xOut, yOut, zOut, xOut, yOut, zOut)
  1120  			p256PointDouble(xOut, yOut, zOut, xOut, yOut, zOut)
  1121  			p256PointDouble(xOut, yOut, zOut, xOut, yOut, zOut)
  1122  		}
  1123  
  1124  		index = uint32(scalar[31-i/2])
  1125  		if (i & 1) == 1 {
  1126  			index &= 15
  1127  		} else {
  1128  			index >>= 4
  1129  		}
  1130  
  1131  		// See the comments in scalarBaseMult about handling infinities.
  1132  		p256SelectJacobianPoint(&px, &py, &pz, &precomp, index)
  1133  		p256PointAdd(&tx, &ty, &tz, xOut, yOut, zOut, &px, &py, &pz)
  1134  		p256CopyConditional(xOut, &px, nIsInfinityMask)
  1135  		p256CopyConditional(yOut, &py, nIsInfinityMask)
  1136  		p256CopyConditional(zOut, &pz, nIsInfinityMask)
  1137  
  1138  		pIsNoninfiniteMask = nonZeroToAllOnes(index)
  1139  		mask = pIsNoninfiniteMask & ^nIsInfinityMask
  1140  		p256CopyConditional(xOut, &tx, mask)
  1141  		p256CopyConditional(yOut, &ty, mask)
  1142  		p256CopyConditional(zOut, &tz, mask)
  1143  		nIsInfinityMask &^= pIsNoninfiniteMask
  1144  	}
  1145  }
  1146  
  1147  // p256FromBig sets out = R*in.
  1148  func p256FromBig(out *[p256Limbs]uint32, in *big.Int) {
  1149  	tmp := new(big.Int).Lsh(in, 257)
  1150  	tmp.Mod(tmp, p256Params.P)
  1151  
  1152  	for i := 0; i < p256Limbs; i++ {
  1153  		if bits := tmp.Bits(); len(bits) > 0 {
  1154  			out[i] = uint32(bits[0]) & bottom29Bits
  1155  		} else {
  1156  			out[i] = 0
  1157  		}
  1158  		tmp.Rsh(tmp, 29)
  1159  
  1160  		i++
  1161  		if i == p256Limbs {
  1162  			break
  1163  		}
  1164  
  1165  		if bits := tmp.Bits(); len(bits) > 0 {
  1166  			out[i] = uint32(bits[0]) & bottom28Bits
  1167  		} else {
  1168  			out[i] = 0
  1169  		}
  1170  		tmp.Rsh(tmp, 28)
  1171  	}
  1172  }
  1173  
  1174  // p256ToBig returns a *big.Int containing the value of in.
  1175  func p256ToBig(in *[p256Limbs]uint32) *big.Int {
  1176  	result, tmp := new(big.Int), new(big.Int)
  1177  
  1178  	result.SetInt64(int64(in[p256Limbs-1]))
  1179  	for i := p256Limbs - 2; i >= 0; i-- {
  1180  		if (i & 1) == 0 {
  1181  			result.Lsh(result, 29)
  1182  		} else {
  1183  			result.Lsh(result, 28)
  1184  		}
  1185  		tmp.SetInt64(int64(in[i]))
  1186  		result.Add(result, tmp)
  1187  	}
  1188  
  1189  	result.Mul(result, p256RInverse)
  1190  	result.Mod(result, p256Params.P)
  1191  	return result
  1192  }
  1193  

View as plain text