field.mx raw

   1  package gnarl
   2  
   3  // Montgomery field arithmetic over Z_P.
   4  // 4x uint64 little-endian limbs in Montgomery form.
   5  // P is 216-bit (27 bytes). R = 2^256.
   6  
   7  import "math/bits"
   8  
   9  type fe [4]uint64
  10  
  11  // These are initialized in init().
  12  var pLimbs fe
  13  var feOne fe
  14  var feR2 fe
  15  var feZero fe
  16  var tsOddQ fe
  17  var tsQm1h fe
  18  var thirdPNorm fe
  19  var halfPNorm fe
  20  var pMinus2 fe
  21  
  22  const pPrime uint64 = 0xf07d39ef3ea058a3
  23  const tsS = 2
  24  const fieldLen = 27
  25  
  26  func initFieldConstants() {
  27  	pLimbs = fe{0x31f9791e0f9ee4f5, 0x121e79ccd682cc99, 0xb6d81bb9b02e5e5d, 0x00000000009563a6}
  28  	feOne = fe{0x761ae156c8e8e77c, 0x023540a7764c229f, 0x315e327188dada36, 0x00000000001012bc}
  29  	feR2 = fe{0xf56927577991be00, 0xf6eb880a87ed3681, 0x18d8b12f949488c0, 0x0000000000684a62}
  30  	tsOddQ = fe{0x4c7e5e4783e7b93d, 0x44879e7335a0b326, 0xadb606ee6c0b9797, 0x00000000002558e9}
  31  	tsQm1h = fe{0x263f2f23c1f3dc9e, 0xa243cf399ad05993, 0xd6db03773605cbcb, 0x000000000012ac74}
  32  	thirdPNorm = fe{0xbb53285f5a8a4c51, 0xb0b4d3444780eedd, 0x3cf2b3e8900f74c9, 0x000000000031cbe2}
  33  	halfPNorm = fe{0x98fcbc8f07cf727a, 0x890f3ce66b41664c, 0x5b6c0ddcd8172f2e, 0x00000000004ab1d3}
  34  	pMinus2 = fe{0x31f9791e0f9ee4f3, 0x121e79ccd682cc99, 0xb6d81bb9b02e5e5d, 0x00000000009563a6}
  35  }
  36  
  37  func feSet(dst, src *fe) {
  38  	*dst = *src
  39  }
  40  
  41  func feIsZero(a *fe) (result int32) {
  42  	z := a[0] | a[1] | a[2] | a[3]
  43  	z = (z | (0 - z)) >> 63
  44  	return 1 - int32(z)
  45  }
  46  
  47  func feEqual(a, b *fe) (result int32) {
  48  	z := (a[0] ^ b[0]) | (a[1] ^ b[1]) | (a[2] ^ b[2]) | (a[3] ^ b[3])
  49  	z = (z | (0 - z)) >> 63
  50  	return 1 - int32(z)
  51  }
  52  
  53  func feReduce(a *fe) {
  54  	var d0, d1, d2, d3, borrow uint64
  55  	d0, borrow = bits.Sub64(a[0], pLimbs[0], 0)
  56  	d1, borrow = bits.Sub64(a[1], pLimbs[1], borrow)
  57  	d2, borrow = bits.Sub64(a[2], pLimbs[2], borrow)
  58  	d3, borrow = bits.Sub64(a[3], pLimbs[3], borrow)
  59  
  60  	mask := uint64(0) - borrow
  61  	a[0] = (a[0] & mask) | (d0 & ^mask)
  62  	a[1] = (a[1] & mask) | (d1 & ^mask)
  63  	a[2] = (a[2] & mask) | (d2 & ^mask)
  64  	a[3] = (a[3] & mask) | (d3 & ^mask)
  65  }
  66  
  67  func feAdd(r, a, b *fe) {
  68  	var carry, borrow, underflow uint64
  69  	var d0, d1, d2, d3 uint64
  70  	r[0], carry = bits.Add64(a[0], b[0], 0)
  71  	r[1], carry = bits.Add64(a[1], b[1], carry)
  72  	r[2], carry = bits.Add64(a[2], b[2], carry)
  73  	r[3], carry = bits.Add64(a[3], b[3], carry)
  74  
  75  	d0, borrow = bits.Sub64(r[0], pLimbs[0], 0)
  76  	d1, borrow = bits.Sub64(r[1], pLimbs[1], borrow)
  77  	d2, borrow = bits.Sub64(r[2], pLimbs[2], borrow)
  78  	d3, borrow = bits.Sub64(r[3], pLimbs[3], borrow)
  79  
  80  	underflow = borrow &^ carry
  81  	mask := uint64(0) - underflow
  82  	r[0] = (r[0] & mask) | (d0 & ^mask)
  83  	r[1] = (r[1] & mask) | (d1 & ^mask)
  84  	r[2] = (r[2] & mask) | (d2 & ^mask)
  85  	r[3] = (r[3] & mask) | (d3 & ^mask)
  86  }
  87  
  88  func feSub(r, a, b *fe) {
  89  	var borrow, carry uint64
  90  	r[0], borrow = bits.Sub64(a[0], b[0], 0)
  91  	r[1], borrow = bits.Sub64(a[1], b[1], borrow)
  92  	r[2], borrow = bits.Sub64(a[2], b[2], borrow)
  93  	r[3], borrow = bits.Sub64(a[3], b[3], borrow)
  94  
  95  	mask := uint64(0) - borrow
  96  	r[0], carry = bits.Add64(r[0], pLimbs[0]&mask, 0)
  97  	r[1], carry = bits.Add64(r[1], pLimbs[1]&mask, carry)
  98  	r[2], carry = bits.Add64(r[2], pLimbs[2]&mask, carry)
  99  	r[3], _ = bits.Add64(r[3], pLimbs[3]&mask, carry)
 100  }
 101  
 102  func feNeg(r, a *fe) {
 103  	z := a[0] | a[1] | a[2] | a[3]
 104  	nonzero := z | (0 - z)
 105  	mask := 0 - (nonzero >> 63)
 106  
 107  	var borrow uint64
 108  	r[0], borrow = bits.Sub64(pLimbs[0], a[0], 0)
 109  	r[1], borrow = bits.Sub64(pLimbs[1], a[1], borrow)
 110  	r[2], borrow = bits.Sub64(pLimbs[2], a[2], borrow)
 111  	r[3], _ = bits.Sub64(pLimbs[3], a[3], borrow)
 112  	r[0] &= mask
 113  	r[1] &= mask
 114  	r[2] &= mask
 115  	r[3] &= mask
 116  }
 117  
 118  // montMul computes r = a * b * R^{-1} mod P using fully-unrolled CIOS.
 119  func montMul(r, a, b *fe) {
 120  	var t0, t1, t2, t3, t4 uint64
 121  	var hi, lo, c0, c1, m, carry uint64
 122  
 123  	// i = 0
 124  	hi, lo = bits.Mul64(a[0], b[0])
 125  	t0, c0 = bits.Add64(lo, t0, 0)
 126  	carry, _ = bits.Add64(hi, 0, c0)
 127  	hi, lo = bits.Mul64(a[0], b[1])
 128  	lo, c0 = bits.Add64(lo, t1, 0)
 129  	hi, c1 = bits.Add64(hi, 0, c0)
 130  	t1, c0 = bits.Add64(lo, carry, 0)
 131  	carry, _ = bits.Add64(hi, c1, c0)
 132  	hi, lo = bits.Mul64(a[0], b[2])
 133  	lo, c0 = bits.Add64(lo, t2, 0)
 134  	hi, c1 = bits.Add64(hi, 0, c0)
 135  	t2, c0 = bits.Add64(lo, carry, 0)
 136  	carry, _ = bits.Add64(hi, c1, c0)
 137  	hi, lo = bits.Mul64(a[0], b[3])
 138  	lo, c0 = bits.Add64(lo, t3, 0)
 139  	hi, c1 = bits.Add64(hi, 0, c0)
 140  	t3, c0 = bits.Add64(lo, carry, 0)
 141  	t4, _ = bits.Add64(hi, c1, c0)
 142  	m = t0 * pPrime
 143  	hi, lo = bits.Mul64(m, pLimbs[0])
 144  	_, c0 = bits.Add64(t0, lo, 0)
 145  	carry, _ = bits.Add64(hi, 0, c0)
 146  	hi, lo = bits.Mul64(m, pLimbs[1])
 147  	lo, c0 = bits.Add64(lo, t1, 0)
 148  	hi, c1 = bits.Add64(hi, 0, c0)
 149  	t0, c0 = bits.Add64(lo, carry, 0)
 150  	carry, _ = bits.Add64(hi, c1, c0)
 151  	hi, lo = bits.Mul64(m, pLimbs[2])
 152  	lo, c0 = bits.Add64(lo, t2, 0)
 153  	hi, c1 = bits.Add64(hi, 0, c0)
 154  	t1, c0 = bits.Add64(lo, carry, 0)
 155  	carry, _ = bits.Add64(hi, c1, c0)
 156  	hi, lo = bits.Mul64(m, pLimbs[3])
 157  	lo, c0 = bits.Add64(lo, t3, 0)
 158  	hi, c1 = bits.Add64(hi, 0, c0)
 159  	t2, c0 = bits.Add64(lo, carry, 0)
 160  	carry, _ = bits.Add64(hi, c1, c0)
 161  	t3, c0 = bits.Add64(t4, carry, 0)
 162  	t4 = c0
 163  
 164  	// i = 1
 165  	hi, lo = bits.Mul64(a[1], b[0])
 166  	t0, c0 = bits.Add64(lo, t0, 0)
 167  	carry, _ = bits.Add64(hi, 0, c0)
 168  	hi, lo = bits.Mul64(a[1], b[1])
 169  	lo, c0 = bits.Add64(lo, t1, 0)
 170  	hi, c1 = bits.Add64(hi, 0, c0)
 171  	t1, c0 = bits.Add64(lo, carry, 0)
 172  	carry, _ = bits.Add64(hi, c1, c0)
 173  	hi, lo = bits.Mul64(a[1], b[2])
 174  	lo, c0 = bits.Add64(lo, t2, 0)
 175  	hi, c1 = bits.Add64(hi, 0, c0)
 176  	t2, c0 = bits.Add64(lo, carry, 0)
 177  	carry, _ = bits.Add64(hi, c1, c0)
 178  	hi, lo = bits.Mul64(a[1], b[3])
 179  	lo, c0 = bits.Add64(lo, t3, 0)
 180  	hi, c1 = bits.Add64(hi, 0, c0)
 181  	t3, c0 = bits.Add64(lo, carry, 0)
 182  	hi, _ = bits.Add64(hi, c1, c0)
 183  	t4 += hi
 184  	m = t0 * pPrime
 185  	hi, lo = bits.Mul64(m, pLimbs[0])
 186  	_, c0 = bits.Add64(t0, lo, 0)
 187  	carry, _ = bits.Add64(hi, 0, c0)
 188  	hi, lo = bits.Mul64(m, pLimbs[1])
 189  	lo, c0 = bits.Add64(lo, t1, 0)
 190  	hi, c1 = bits.Add64(hi, 0, c0)
 191  	t0, c0 = bits.Add64(lo, carry, 0)
 192  	carry, _ = bits.Add64(hi, c1, c0)
 193  	hi, lo = bits.Mul64(m, pLimbs[2])
 194  	lo, c0 = bits.Add64(lo, t2, 0)
 195  	hi, c1 = bits.Add64(hi, 0, c0)
 196  	t1, c0 = bits.Add64(lo, carry, 0)
 197  	carry, _ = bits.Add64(hi, c1, c0)
 198  	hi, lo = bits.Mul64(m, pLimbs[3])
 199  	lo, c0 = bits.Add64(lo, t3, 0)
 200  	hi, c1 = bits.Add64(hi, 0, c0)
 201  	t2, c0 = bits.Add64(lo, carry, 0)
 202  	carry, _ = bits.Add64(hi, c1, c0)
 203  	t3, c0 = bits.Add64(t4, carry, 0)
 204  	t4 = c0
 205  
 206  	// i = 2
 207  	hi, lo = bits.Mul64(a[2], b[0])
 208  	t0, c0 = bits.Add64(lo, t0, 0)
 209  	carry, _ = bits.Add64(hi, 0, c0)
 210  	hi, lo = bits.Mul64(a[2], b[1])
 211  	lo, c0 = bits.Add64(lo, t1, 0)
 212  	hi, c1 = bits.Add64(hi, 0, c0)
 213  	t1, c0 = bits.Add64(lo, carry, 0)
 214  	carry, _ = bits.Add64(hi, c1, c0)
 215  	hi, lo = bits.Mul64(a[2], b[2])
 216  	lo, c0 = bits.Add64(lo, t2, 0)
 217  	hi, c1 = bits.Add64(hi, 0, c0)
 218  	t2, c0 = bits.Add64(lo, carry, 0)
 219  	carry, _ = bits.Add64(hi, c1, c0)
 220  	hi, lo = bits.Mul64(a[2], b[3])
 221  	lo, c0 = bits.Add64(lo, t3, 0)
 222  	hi, c1 = bits.Add64(hi, 0, c0)
 223  	t3, c0 = bits.Add64(lo, carry, 0)
 224  	hi, _ = bits.Add64(hi, c1, c0)
 225  	t4 += hi
 226  	m = t0 * pPrime
 227  	hi, lo = bits.Mul64(m, pLimbs[0])
 228  	_, c0 = bits.Add64(t0, lo, 0)
 229  	carry, _ = bits.Add64(hi, 0, c0)
 230  	hi, lo = bits.Mul64(m, pLimbs[1])
 231  	lo, c0 = bits.Add64(lo, t1, 0)
 232  	hi, c1 = bits.Add64(hi, 0, c0)
 233  	t0, c0 = bits.Add64(lo, carry, 0)
 234  	carry, _ = bits.Add64(hi, c1, c0)
 235  	hi, lo = bits.Mul64(m, pLimbs[2])
 236  	lo, c0 = bits.Add64(lo, t2, 0)
 237  	hi, c1 = bits.Add64(hi, 0, c0)
 238  	t1, c0 = bits.Add64(lo, carry, 0)
 239  	carry, _ = bits.Add64(hi, c1, c0)
 240  	hi, lo = bits.Mul64(m, pLimbs[3])
 241  	lo, c0 = bits.Add64(lo, t3, 0)
 242  	hi, c1 = bits.Add64(hi, 0, c0)
 243  	t2, c0 = bits.Add64(lo, carry, 0)
 244  	carry, _ = bits.Add64(hi, c1, c0)
 245  	t3, c0 = bits.Add64(t4, carry, 0)
 246  	t4 = c0
 247  
 248  	// i = 3
 249  	hi, lo = bits.Mul64(a[3], b[0])
 250  	t0, c0 = bits.Add64(lo, t0, 0)
 251  	carry, _ = bits.Add64(hi, 0, c0)
 252  	hi, lo = bits.Mul64(a[3], b[1])
 253  	lo, c0 = bits.Add64(lo, t1, 0)
 254  	hi, c1 = bits.Add64(hi, 0, c0)
 255  	t1, c0 = bits.Add64(lo, carry, 0)
 256  	carry, _ = bits.Add64(hi, c1, c0)
 257  	hi, lo = bits.Mul64(a[3], b[2])
 258  	lo, c0 = bits.Add64(lo, t2, 0)
 259  	hi, c1 = bits.Add64(hi, 0, c0)
 260  	t2, c0 = bits.Add64(lo, carry, 0)
 261  	carry, _ = bits.Add64(hi, c1, c0)
 262  	hi, lo = bits.Mul64(a[3], b[3])
 263  	lo, c0 = bits.Add64(lo, t3, 0)
 264  	hi, c1 = bits.Add64(hi, 0, c0)
 265  	t3, c0 = bits.Add64(lo, carry, 0)
 266  	hi, _ = bits.Add64(hi, c1, c0)
 267  	t4 += hi
 268  	m = t0 * pPrime
 269  	hi, lo = bits.Mul64(m, pLimbs[0])
 270  	_, c0 = bits.Add64(t0, lo, 0)
 271  	carry, _ = bits.Add64(hi, 0, c0)
 272  	hi, lo = bits.Mul64(m, pLimbs[1])
 273  	lo, c0 = bits.Add64(lo, t1, 0)
 274  	hi, c1 = bits.Add64(hi, 0, c0)
 275  	t0, c0 = bits.Add64(lo, carry, 0)
 276  	carry, _ = bits.Add64(hi, c1, c0)
 277  	hi, lo = bits.Mul64(m, pLimbs[2])
 278  	lo, c0 = bits.Add64(lo, t2, 0)
 279  	hi, c1 = bits.Add64(hi, 0, c0)
 280  	t1, c0 = bits.Add64(lo, carry, 0)
 281  	carry, _ = bits.Add64(hi, c1, c0)
 282  	hi, lo = bits.Mul64(m, pLimbs[3])
 283  	lo, c0 = bits.Add64(lo, t3, 0)
 284  	hi, c1 = bits.Add64(hi, 0, c0)
 285  	t2, c0 = bits.Add64(lo, carry, 0)
 286  	carry, _ = bits.Add64(hi, c1, c0)
 287  	t3, c0 = bits.Add64(t4, carry, 0)
 288  	t4 = c0
 289  
 290  	r[0] = t0
 291  	r[1] = t1
 292  	r[2] = t2
 293  	r[3] = t3
 294  	feReduce(r)
 295  }
 296  
 297  func montSquare(r, a *fe) {
 298  	montMul(r, a, a)
 299  }
 300  
 301  func feToMont(r, a *fe) {
 302  	montMul(r, a, &feR2)
 303  }
 304  
 305  func feFromMont(r, a *fe) {
 306  	one := fe{1, 0, 0, 0}
 307  	montMul(r, a, &one)
 308  }
 309  
 310  func feFromSmall(r *fe, v int64) {
 311  	if v >= 0 {
 312  		*r = fe{uint64(v), 0, 0, 0}
 313  	} else {
 314  		*r = pLimbs
 315  		var borrow uint64
 316  		r[0], borrow = bits.Sub64(r[0], uint64(-v), 0)
 317  		r[1], borrow = bits.Sub64(r[1], 0, borrow)
 318  		r[2], borrow = bits.Sub64(r[2], 0, borrow)
 319  		r[3], _ = bits.Sub64(r[3], 0, borrow)
 320  	}
 321  	feToMont(r, r)
 322  }
 323  
 324  func feFromBytes27(r *fe, b []byte) (ok bool) {
 325  	if int32(len(b)) < 27 {
 326  		*r = feZero
 327  		return false
 328  	}
 329  	r[3] = uint64(b[0])<<16 | uint64(b[1])<<8 | uint64(b[2])
 330  	r[2] = uint64(b[3])<<56 | uint64(b[4])<<48 | uint64(b[5])<<40 | uint64(b[6])<<32 |
 331  		uint64(b[7])<<24 | uint64(b[8])<<16 | uint64(b[9])<<8 | uint64(b[10])
 332  	r[1] = uint64(b[11])<<56 | uint64(b[12])<<48 | uint64(b[13])<<40 | uint64(b[14])<<32 |
 333  		uint64(b[15])<<24 | uint64(b[16])<<16 | uint64(b[17])<<8 | uint64(b[18])
 334  	r[0] = uint64(b[19])<<56 | uint64(b[20])<<48 | uint64(b[21])<<40 | uint64(b[22])<<32 |
 335  		uint64(b[23])<<24 | uint64(b[24])<<16 | uint64(b[25])<<8 | uint64(b[26])
 336  
 337  	var d0, d1, borrow uint64
 338  	d0, borrow = bits.Sub64(r[0], pLimbs[0], 0)
 339  	d1, borrow = bits.Sub64(r[1], pLimbs[1], borrow)
 340  	_, borrow = bits.Sub64(r[2], pLimbs[2], borrow)
 341  	_, borrow = bits.Sub64(r[3], pLimbs[3], borrow)
 342  	_ = d0
 343  	_ = d1
 344  	if borrow == 0 {
 345  		*r = feZero
 346  		return false
 347  	}
 348  
 349  	feToMont(r, r)
 350  	return true
 351  }
 352  
 353  func feToBytes27(b []byte, a *fe) {
 354  	var norm fe
 355  	feFromMont(&norm, a)
 356  	b[0] = byte(norm[3] >> 16)
 357  	b[1] = byte(norm[3] >> 8)
 358  	b[2] = byte(norm[3])
 359  	b[3] = byte(norm[2] >> 56)
 360  	b[4] = byte(norm[2] >> 48)
 361  	b[5] = byte(norm[2] >> 40)
 362  	b[6] = byte(norm[2] >> 32)
 363  	b[7] = byte(norm[2] >> 24)
 364  	b[8] = byte(norm[2] >> 16)
 365  	b[9] = byte(norm[2] >> 8)
 366  	b[10] = byte(norm[2])
 367  	b[11] = byte(norm[1] >> 56)
 368  	b[12] = byte(norm[1] >> 48)
 369  	b[13] = byte(norm[1] >> 40)
 370  	b[14] = byte(norm[1] >> 32)
 371  	b[15] = byte(norm[1] >> 24)
 372  	b[16] = byte(norm[1] >> 16)
 373  	b[17] = byte(norm[1] >> 8)
 374  	b[18] = byte(norm[1])
 375  	b[19] = byte(norm[0] >> 56)
 376  	b[20] = byte(norm[0] >> 48)
 377  	b[21] = byte(norm[0] >> 40)
 378  	b[22] = byte(norm[0] >> 32)
 379  	b[23] = byte(norm[0] >> 24)
 380  	b[24] = byte(norm[0] >> 16)
 381  	b[25] = byte(norm[0] >> 8)
 382  	b[26] = byte(norm[0])
 383  }
 384  
 385  func fePow(r, base, exp *fe) {
 386  	var result fe
 387  	feSet(&result, &feOne)
 388  	var b fe
 389  	feSet(&b, base)
 390  
 391  	for i := int32(0); i < 4; i++ {
 392  		word := exp[i]
 393  		for bit := int32(0); bit < 64; bit++ {
 394  			if word&1 == 1 {
 395  				montMul(&result, &result, &b)
 396  			}
 397  			montSquare(&b, &b)
 398  			word >>= 1
 399  		}
 400  	}
 401  	feSet(r, &result)
 402  }
 403  
 404  func feInv(r, a *fe) {
 405  	fePow(r, a, &pMinus2)
 406  }
 407  
 408  func feSqrt(r, a *fe) (ok bool) {
 409  	if feIsZero(a) == 1 {
 410  		*r = feZero
 411  		return true
 412  	}
 413  
 414  	var feFive fe
 415  	feFive[0] = 5
 416  	feToMont(&feFive, &feFive)
 417  
 418  	var t fe
 419  	fePow(&t, a, &tsOddQ)
 420  
 421  	var rr, tmp fe
 422  	fePow(&tmp, a, &tsQm1h)
 423  	montMul(&rr, &tmp, a)
 424  
 425  	var c fe
 426  	fePow(&c, &feFive, &tsOddQ)
 427  
 428  	m := int32(tsS)
 429  
 430  	for {
 431  		if feEqual(&t, &feOne) == 1 {
 432  			feSet(r, &rr)
 433  			return true
 434  		}
 435  
 436  		var bv fe
 437  		feSet(&bv, &t)
 438  		iv := int32(0)
 439  		for j := int32(1); j < m; j++ {
 440  			montSquare(&bv, &bv)
 441  			if feEqual(&bv, &feOne) == 1 {
 442  				iv = j
 443  				break
 444  			}
 445  		}
 446  		if iv == 0 {
 447  			return false
 448  		}
 449  
 450  		feSet(&bv, &c)
 451  		for j := int32(0); j < m-iv-1; j++ {
 452  			montSquare(&bv, &bv)
 453  		}
 454  
 455  		montMul(&rr, &rr, &bv)
 456  		montSquare(&c, &bv)
 457  		montMul(&t, &t, &c)
 458  		m = iv
 459  	}
 460  }
 461  
 462  func feIsGtThirdP(a *fe) (result int32) {
 463  	for i := int32(3); i >= 0; i-- {
 464  		if a[i] > thirdPNorm[i] {
 465  			return 1
 466  		}
 467  		if a[i] < thirdPNorm[i] {
 468  			return 0
 469  		}
 470  	}
 471  	return 0
 472  }
 473  
 474  func feIsGtHalfP(a *fe) (result int32) {
 475  	for i := int32(3); i >= 0; i-- {
 476  		if a[i] > halfPNorm[i] {
 477  			return 1
 478  		}
 479  		if a[i] < halfPNorm[i] {
 480  			return 0
 481  		}
 482  	}
 483  	return 0
 484  }
 485