Source file src/crypto/internal/fips140/bigmod/nat.go
1 // Copyright 2021 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 bigmod 6 7 import ( 8 _ "crypto/internal/fips140/check" 9 "crypto/internal/fips140deps/byteorder" 10 "errors" 11 "math/bits" 12 ) 13 14 const ( 15 // _W is the size in bits of our limbs. 16 _W = bits.UintSize 17 // _S is the size in bytes of our limbs. 18 _S = _W / 8 19 ) 20 21 // Note: These functions make many loops over all the words in a Nat. 22 // These loops used to be in assembly, invisible to -race, -asan, and -msan, 23 // but now they are in Go and incur significant overhead in those modes. 24 // To bring the old performance back, we mark all functions that loop 25 // over Nat words with //go:norace. Because //go:norace does not 26 // propagate across inlining, we must also mark functions that inline 27 // //go:norace functions - specifically, those that inline add, addMulVVW, 28 // assign, cmpGeq, rshift1, and sub. 29 30 // choice represents a constant-time boolean. The value of choice is always 31 // either 1 or 0. We use an int instead of bool in order to make decisions in 32 // constant time by turning it into a mask. 33 type choice uint 34 35 func not(c choice) choice { return 1 ^ c } 36 37 const yes = choice(1) 38 const no = choice(0) 39 40 // ctMask is all 1s if on is yes, and all 0s otherwise. 41 func ctMask(on choice) uint { return -uint(on) } 42 43 // ctEq returns 1 if x == y, and 0 otherwise. The execution time of this 44 // function does not depend on its inputs. 45 func ctEq(x, y uint) choice { 46 // If x != y, then either x - y or y - x will generate a carry. 47 _, c1 := bits.Sub(x, y, 0) 48 _, c2 := bits.Sub(y, x, 0) 49 return not(choice(c1 | c2)) 50 } 51 52 // Nat represents an arbitrary natural number 53 // 54 // Each Nat has an announced length, which is the number of limbs it has stored. 55 // Operations on this number are allowed to leak this length, but will not leak 56 // any information about the values contained in those limbs. 57 type Nat struct { 58 // limbs is little-endian in base 2^W with W = bits.UintSize. 59 limbs []uint 60 } 61 62 // preallocTarget is the size in bits of the numbers used to implement the most 63 // common and most performant RSA key size. It's also enough to cover some of 64 // the operations of key sizes up to 4096. 65 const preallocTarget = 2048 66 const preallocLimbs = (preallocTarget + _W - 1) / _W 67 68 // NewNat returns a new nat with a size of zero, just like new(Nat), but with 69 // the preallocated capacity to hold a number of up to preallocTarget bits. 70 // NewNat inlines, so the allocation can live on the stack. 71 func NewNat() *Nat { 72 limbs := make([]uint, 0, preallocLimbs) 73 return &Nat{limbs} 74 } 75 76 // expand expands x to n limbs, leaving its value unchanged. 77 func (x *Nat) expand(n int) *Nat { 78 if len(x.limbs) > n { 79 panic("bigmod: internal error: shrinking nat") 80 } 81 if cap(x.limbs) < n { 82 newLimbs := make([]uint, n) 83 copy(newLimbs, x.limbs) 84 x.limbs = newLimbs 85 return x 86 } 87 extraLimbs := x.limbs[len(x.limbs):n] 88 clear(extraLimbs) 89 x.limbs = x.limbs[:n] 90 return x 91 } 92 93 // reset returns a zero nat of n limbs, reusing x's storage if n <= cap(x.limbs). 94 func (x *Nat) reset(n int) *Nat { 95 if cap(x.limbs) < n { 96 x.limbs = make([]uint, n) 97 return x 98 } 99 // Clear both the returned limbs and the previously used ones. 100 clear(x.limbs[:max(n, len(x.limbs))]) 101 x.limbs = x.limbs[:n] 102 return x 103 } 104 105 // resetToBytes assigns x = b, where b is a slice of big-endian bytes, resizing 106 // n to the appropriate size. 107 // 108 // The announced length of x is set based on the actual bit size of the input, 109 // ignoring leading zeroes. 110 func (x *Nat) resetToBytes(b []byte) *Nat { 111 x.reset((len(b) + _S - 1) / _S) 112 if err := x.setBytes(b); err != nil { 113 panic("bigmod: internal error: bad arithmetic") 114 } 115 return x.trim() 116 } 117 118 // trim reduces the size of x to match its value. 119 func (x *Nat) trim() *Nat { 120 // Trim most significant (trailing in little-endian) zero limbs. 121 // We assume comparison with zero (but not the branch) is constant time. 122 for i := len(x.limbs) - 1; i >= 0; i-- { 123 if x.limbs[i] != 0 { 124 break 125 } 126 x.limbs = x.limbs[:i] 127 } 128 return x 129 } 130 131 // set assigns x = y, optionally resizing x to the appropriate size. 132 func (x *Nat) set(y *Nat) *Nat { 133 x.reset(len(y.limbs)) 134 copy(x.limbs, y.limbs) 135 return x 136 } 137 138 // Bits returns x as a little-endian slice of uint. The length of the slice 139 // matches the announced length of x. The result and x share the same underlying 140 // array. 141 func (x *Nat) Bits() []uint { 142 return x.limbs 143 } 144 145 // SetBits assigns x = y, where y is a slice of little-endian uint. x is resized 146 // to the length of y. 147 func (x *Nat) SetBits(y []uint) *Nat { 148 x.reset(len(y)) 149 copy(x.limbs, y) 150 return x 151 } 152 153 // Bytes returns x as a zero-extended big-endian byte slice. The size of the 154 // slice will match the size of m. 155 // 156 // x must have the same size as m and it must be less than or equal to m. 157 func (x *Nat) Bytes(m *Modulus) []byte { 158 i := m.Size() 159 bytes := make([]byte, i) 160 for _, limb := range x.limbs { 161 for j := 0; j < _S; j++ { 162 i-- 163 if i < 0 { 164 if limb == 0 { 165 break 166 } 167 panic("bigmod: modulus is smaller than nat") 168 } 169 bytes[i] = byte(limb) 170 limb >>= 8 171 } 172 } 173 return bytes 174 } 175 176 // SetBytes assigns x = b, where b is a slice of big-endian bytes. 177 // SetBytes returns an error if b >= m. 178 // 179 // The output will be resized to the size of m and overwritten. 180 // 181 //go:norace 182 func (x *Nat) SetBytes(b []byte, m *Modulus) (*Nat, error) { 183 x.resetFor(m) 184 if err := x.setBytes(b); err != nil { 185 return nil, err 186 } 187 if x.cmpGeq(m.nat) == yes { 188 return nil, errors.New("input overflows the modulus") 189 } 190 return x, nil 191 } 192 193 // SetOverflowingBytes assigns x = b, where b is a slice of big-endian bytes. 194 // SetOverflowingBytes returns an error if b has a longer bit length than m, but 195 // reduces overflowing values up to 2^⌈log2(m)⌉ - 1. 196 // 197 // The output will be resized to the size of m and overwritten. 198 func (x *Nat) SetOverflowingBytes(b []byte, m *Modulus) (*Nat, error) { 199 x.resetFor(m) 200 if err := x.setBytes(b); err != nil { 201 return nil, err 202 } 203 // setBytes would have returned an error if the input overflowed the limb 204 // size of the modulus, so now we only need to check if the most significant 205 // limb of x has more bits than the most significant limb of the modulus. 206 if bitLen(x.limbs[len(x.limbs)-1]) > bitLen(m.nat.limbs[len(m.nat.limbs)-1]) { 207 return nil, errors.New("input overflows the modulus size") 208 } 209 x.maybeSubtractModulus(no, m) 210 return x, nil 211 } 212 213 // bigEndianUint returns the contents of buf interpreted as a 214 // big-endian encoded uint value. 215 func bigEndianUint(buf []byte) uint { 216 if _W == 64 { 217 return uint(byteorder.BEUint64(buf)) 218 } 219 return uint(byteorder.BEUint32(buf)) 220 } 221 222 func (x *Nat) setBytes(b []byte) error { 223 i, k := len(b), 0 224 for k < len(x.limbs) && i >= _S { 225 x.limbs[k] = bigEndianUint(b[i-_S : i]) 226 i -= _S 227 k++ 228 } 229 for s := 0; s < _W && k < len(x.limbs) && i > 0; s += 8 { 230 x.limbs[k] |= uint(b[i-1]) << s 231 i-- 232 } 233 if i > 0 { 234 return errors.New("input overflows the modulus size") 235 } 236 return nil 237 } 238 239 // SetUint assigns x = y. 240 // 241 // The output will be resized to a single limb and overwritten. 242 func (x *Nat) SetUint(y uint) *Nat { 243 x.reset(1) 244 x.limbs[0] = y 245 return x 246 } 247 248 // Equal returns 1 if x == y, and 0 otherwise. 249 // 250 // Both operands must have the same announced length. 251 // 252 //go:norace 253 func (x *Nat) Equal(y *Nat) choice { 254 // Eliminate bounds checks in the loop. 255 size := len(x.limbs) 256 xLimbs := x.limbs[:size] 257 yLimbs := y.limbs[:size] 258 259 equal := yes 260 for i := 0; i < size; i++ { 261 equal &= ctEq(xLimbs[i], yLimbs[i]) 262 } 263 return equal 264 } 265 266 // IsZero returns 1 if x == 0, and 0 otherwise. 267 // 268 //go:norace 269 func (x *Nat) IsZero() choice { 270 // Eliminate bounds checks in the loop. 271 size := len(x.limbs) 272 xLimbs := x.limbs[:size] 273 274 zero := yes 275 for i := 0; i < size; i++ { 276 zero &= ctEq(xLimbs[i], 0) 277 } 278 return zero 279 } 280 281 // IsOne returns 1 if x == 1, and 0 otherwise. 282 // 283 //go:norace 284 func (x *Nat) IsOne() choice { 285 // Eliminate bounds checks in the loop. 286 size := len(x.limbs) 287 xLimbs := x.limbs[:size] 288 289 if len(xLimbs) == 0 { 290 return no 291 } 292 293 one := ctEq(xLimbs[0], 1) 294 for i := 1; i < size; i++ { 295 one &= ctEq(xLimbs[i], 0) 296 } 297 return one 298 } 299 300 // IsMinusOne returns 1 if x == -1 mod m, and 0 otherwise. 301 // 302 // The length of x must be the same as the modulus. x must already be reduced 303 // modulo m. 304 // 305 //go:norace 306 func (x *Nat) IsMinusOne(m *Modulus) choice { 307 minusOne := m.Nat() 308 minusOne.SubOne(m) 309 return x.Equal(minusOne) 310 } 311 312 // IsOdd returns 1 if x is odd, and 0 otherwise. 313 func (x *Nat) IsOdd() choice { 314 if len(x.limbs) == 0 { 315 return no 316 } 317 return choice(x.limbs[0] & 1) 318 } 319 320 // TrailingZeroBitsVarTime returns the number of trailing zero bits in x. 321 func (x *Nat) TrailingZeroBitsVarTime() uint { 322 var t uint 323 limbs := x.limbs 324 for _, l := range limbs { 325 if l == 0 { 326 t += _W 327 continue 328 } 329 t += uint(bits.TrailingZeros(l)) 330 break 331 } 332 return t 333 } 334 335 // cmpGeq returns 1 if x >= y, and 0 otherwise. 336 // 337 // Both operands must have the same announced length. 338 // 339 //go:norace 340 func (x *Nat) cmpGeq(y *Nat) choice { 341 // Eliminate bounds checks in the loop. 342 size := len(x.limbs) 343 xLimbs := x.limbs[:size] 344 yLimbs := y.limbs[:size] 345 346 var c uint 347 for i := 0; i < size; i++ { 348 _, c = bits.Sub(xLimbs[i], yLimbs[i], c) 349 } 350 // If there was a carry, then subtracting y underflowed, so 351 // x is not greater than or equal to y. 352 return not(choice(c)) 353 } 354 355 // assign sets x <- y if on == 1, and does nothing otherwise. 356 // 357 // Both operands must have the same announced length. 358 // 359 //go:norace 360 func (x *Nat) assign(on choice, y *Nat) *Nat { 361 // Eliminate bounds checks in the loop. 362 size := len(x.limbs) 363 xLimbs := x.limbs[:size] 364 yLimbs := y.limbs[:size] 365 366 mask := ctMask(on) 367 for i := 0; i < size; i++ { 368 xLimbs[i] ^= mask & (xLimbs[i] ^ yLimbs[i]) 369 } 370 return x 371 } 372 373 // add computes x += y and returns the carry. 374 // 375 // Both operands must have the same announced length. 376 // 377 //go:norace 378 func (x *Nat) add(y *Nat) (c uint) { 379 // Eliminate bounds checks in the loop. 380 size := len(x.limbs) 381 xLimbs := x.limbs[:size] 382 yLimbs := y.limbs[:size] 383 384 for i := 0; i < size; i++ { 385 xLimbs[i], c = bits.Add(xLimbs[i], yLimbs[i], c) 386 } 387 return 388 } 389 390 // sub computes x -= y. It returns the borrow of the subtraction. 391 // 392 // Both operands must have the same announced length. 393 // 394 //go:norace 395 func (x *Nat) sub(y *Nat) (c uint) { 396 // Eliminate bounds checks in the loop. 397 size := len(x.limbs) 398 xLimbs := x.limbs[:size] 399 yLimbs := y.limbs[:size] 400 401 for i := 0; i < size; i++ { 402 xLimbs[i], c = bits.Sub(xLimbs[i], yLimbs[i], c) 403 } 404 return 405 } 406 407 // ShiftRightVarTime sets x = x >> n. 408 // 409 // The announced length of x is unchanged. 410 // 411 //go:norace 412 func (x *Nat) ShiftRightVarTime(n uint) *Nat { 413 // Eliminate bounds checks in the loop. 414 size := len(x.limbs) 415 xLimbs := x.limbs[:size] 416 417 shift := int(n % _W) 418 shiftLimbs := int(n / _W) 419 420 var shiftedLimbs []uint 421 if shiftLimbs < size { 422 shiftedLimbs = xLimbs[shiftLimbs:] 423 } 424 425 for i := range xLimbs { 426 if i >= len(shiftedLimbs) { 427 xLimbs[i] = 0 428 continue 429 } 430 431 xLimbs[i] = shiftedLimbs[i] >> shift 432 if i+1 < len(shiftedLimbs) { 433 xLimbs[i] |= shiftedLimbs[i+1] << (_W - shift) 434 } 435 } 436 437 return x 438 } 439 440 // BitLenVarTime returns the actual size of x in bits. 441 // 442 // The actual size of x (but nothing more) leaks through timing side-channels. 443 // Note that this is ordinarily secret, as opposed to the announced size of x. 444 func (x *Nat) BitLenVarTime() int { 445 // Eliminate bounds checks in the loop. 446 size := len(x.limbs) 447 xLimbs := x.limbs[:size] 448 449 for i := size - 1; i >= 0; i-- { 450 if xLimbs[i] != 0 { 451 return i*_W + bitLen(xLimbs[i]) 452 } 453 } 454 return 0 455 } 456 457 // bitLen is a version of bits.Len that only leaks the bit length of n, but not 458 // its value. bits.Len and bits.LeadingZeros use a lookup table for the 459 // low-order bits on some architectures. 460 func bitLen(n uint) int { 461 len := 0 462 // We assume, here and elsewhere, that comparison to zero is constant time 463 // with respect to different non-zero values. 464 for n != 0 { 465 len++ 466 n >>= 1 467 } 468 return len 469 } 470 471 // Modulus is used for modular arithmetic, precomputing relevant constants. 472 // 473 // A Modulus can leak the exact number of bits needed to store its value 474 // and is stored without padding. Its actual value is still kept secret. 475 type Modulus struct { 476 // The underlying natural number for this modulus. 477 // 478 // This will be stored without any padding, and shouldn't alias with any 479 // other natural number being used. 480 nat *Nat 481 482 // If m is even, the following fields are not set. 483 odd bool 484 m0inv uint // -nat.limbs[0]⁻¹ mod _W 485 rr *Nat // R*R for montgomeryRepresentation 486 } 487 488 // rr returns R*R with R = 2^(_W * n) and n = len(m.nat.limbs). 489 func rr(m *Modulus) *Nat { 490 rr := NewNat().ExpandFor(m) 491 n := uint(len(rr.limbs)) 492 mLen := uint(m.BitLen()) 493 logR := _W * n 494 495 // We start by computing R = 2^(_W * n) mod m. We can get pretty close, to 496 // 2^⌊log₂m⌋, by setting the highest bit we can without having to reduce. 497 rr.limbs[n-1] = 1 << ((mLen - 1) % _W) 498 // Then we double until we reach 2^(_W * n). 499 for i := mLen - 1; i < logR; i++ { 500 rr.Add(rr, m) 501 } 502 503 // Next we need to get from R to 2^(_W * n) R mod m (aka from one to R in 504 // the Montgomery domain, meaning we can use Montgomery multiplication now). 505 // We could do that by doubling _W * n times, or with a square-and-double 506 // chain log2(_W * n) long. Turns out the fastest thing is to start out with 507 // doublings, and switch to square-and-double once the exponent is large 508 // enough to justify the cost of the multiplications. 509 510 // The threshold is selected experimentally as a linear function of n. 511 threshold := n / 4 512 513 // We calculate how many of the most-significant bits of the exponent we can 514 // compute before crossing the threshold, and we do it with doublings. 515 i := bits.UintSize 516 for logR>>i <= threshold { 517 i-- 518 } 519 for k := uint(0); k < logR>>i; k++ { 520 rr.Add(rr, m) 521 } 522 523 // Then we process the remaining bits of the exponent with a 524 // square-and-double chain. 525 for i > 0 { 526 rr.montgomeryMul(rr, rr, m) 527 i-- 528 if logR>>i&1 != 0 { 529 rr.Add(rr, m) 530 } 531 } 532 533 return rr 534 } 535 536 // minusInverseModW computes -x⁻¹ mod _W with x odd. 537 // 538 // This operation is used to precompute a constant involved in Montgomery 539 // multiplication. 540 func minusInverseModW(x uint) uint { 541 // Every iteration of this loop doubles the least-significant bits of 542 // correct inverse in y. The first three bits are already correct (1⁻¹ = 1, 543 // 3⁻¹ = 3, 5⁻¹ = 5, and 7⁻¹ = 7 mod 8), so doubling five times is enough 544 // for 64 bits (and wastes only one iteration for 32 bits). 545 // 546 // See https://crypto.stackexchange.com/a/47496. 547 y := x 548 for i := 0; i < 5; i++ { 549 y = y * (2 - x*y) 550 } 551 return -y 552 } 553 554 // NewModulus creates a new Modulus from a slice of big-endian bytes. The 555 // modulus must be greater than one. 556 // 557 // The number of significant bits and whether the modulus is even is leaked 558 // through timing side-channels. 559 func NewModulus(b []byte) (*Modulus, error) { 560 n := NewNat().resetToBytes(b) 561 return newModulus(n) 562 } 563 564 // NewModulusProduct creates a new Modulus from the product of two numbers 565 // represented as big-endian byte slices. The result must be greater than one. 566 // 567 //go:norace 568 func NewModulusProduct(a, b []byte) (*Modulus, error) { 569 x := NewNat().resetToBytes(a) 570 y := NewNat().resetToBytes(b) 571 n := NewNat().reset(len(x.limbs) + len(y.limbs)) 572 for i := range y.limbs { 573 n.limbs[i+len(x.limbs)] = addMulVVW(n.limbs[i:i+len(x.limbs)], x.limbs, y.limbs[i]) 574 } 575 return newModulus(n.trim()) 576 } 577 578 func newModulus(n *Nat) (*Modulus, error) { 579 m := &Modulus{nat: n} 580 if m.nat.IsZero() == yes || m.nat.IsOne() == yes { 581 return nil, errors.New("modulus must be > 1") 582 } 583 if m.nat.IsOdd() == 1 { 584 m.odd = true 585 m.m0inv = minusInverseModW(m.nat.limbs[0]) 586 m.rr = rr(m) 587 } 588 return m, nil 589 } 590 591 // Size returns the size of m in bytes. 592 func (m *Modulus) Size() int { 593 return (m.BitLen() + 7) / 8 594 } 595 596 // BitLen returns the size of m in bits. 597 func (m *Modulus) BitLen() int { 598 return m.nat.BitLenVarTime() 599 } 600 601 // Nat returns m as a Nat. 602 func (m *Modulus) Nat() *Nat { 603 // Make a copy so that the caller can't modify m.nat or alias it with 604 // another Nat in a modulus operation. 605 n := NewNat() 606 n.set(m.nat) 607 return n 608 } 609 610 // shiftIn calculates x = x << _W + y mod m. 611 // 612 // This assumes that x is already reduced mod m. 613 // 614 //go:norace 615 func (x *Nat) shiftIn(y uint, m *Modulus) *Nat { 616 d := NewNat().resetFor(m) 617 618 // Eliminate bounds checks in the loop. 619 size := len(m.nat.limbs) 620 xLimbs := x.limbs[:size] 621 dLimbs := d.limbs[:size] 622 mLimbs := m.nat.limbs[:size] 623 624 // Each iteration of this loop computes x = 2x + b mod m, where b is a bit 625 // from y. Effectively, it left-shifts x and adds y one bit at a time, 626 // reducing it every time. 627 // 628 // To do the reduction, each iteration computes both 2x + b and 2x + b - m. 629 // The next iteration (and finally the return line) will use either result 630 // based on whether 2x + b overflows m. 631 needSubtraction := no 632 for i := _W - 1; i >= 0; i-- { 633 carry := (y >> i) & 1 634 var borrow uint 635 mask := ctMask(needSubtraction) 636 for i := 0; i < size; i++ { 637 l := xLimbs[i] ^ (mask & (xLimbs[i] ^ dLimbs[i])) 638 xLimbs[i], carry = bits.Add(l, l, carry) 639 dLimbs[i], borrow = bits.Sub(xLimbs[i], mLimbs[i], borrow) 640 } 641 // Like in maybeSubtractModulus, we need the subtraction if either it 642 // didn't underflow (meaning 2x + b > m) or if computing 2x + b 643 // overflowed (meaning 2x + b > 2^(_W * n) > m). 644 needSubtraction = not(choice(borrow)) | choice(carry) 645 } 646 return x.assign(needSubtraction, d) 647 } 648 649 // Mod calculates out = x mod m. 650 // 651 // This works regardless how large the value of x is. 652 // 653 // The output will be resized to the size of m and overwritten. 654 // 655 //go:norace 656 func (out *Nat) Mod(x *Nat, m *Modulus) *Nat { 657 out.resetFor(m) 658 // Working our way from the most significant to the least significant limb, 659 // we can insert each limb at the least significant position, shifting all 660 // previous limbs left by _W. This way each limb will get shifted by the 661 // correct number of bits. We can insert at least N - 1 limbs without 662 // overflowing m. After that, we need to reduce every time we shift. 663 i := len(x.limbs) - 1 664 // For the first N - 1 limbs we can skip the actual shifting and position 665 // them at the shifted position, which starts at min(N - 2, i). 666 start := len(m.nat.limbs) - 2 667 if i < start { 668 start = i 669 } 670 for j := start; j >= 0; j-- { 671 out.limbs[j] = x.limbs[i] 672 i-- 673 } 674 // We shift in the remaining limbs, reducing modulo m each time. 675 for i >= 0 { 676 out.shiftIn(x.limbs[i], m) 677 i-- 678 } 679 return out 680 } 681 682 // ExpandFor ensures x has the right size to work with operations modulo m. 683 // 684 // The announced size of x must be smaller than or equal to that of m. 685 func (x *Nat) ExpandFor(m *Modulus) *Nat { 686 return x.expand(len(m.nat.limbs)) 687 } 688 689 // resetFor ensures out has the right size to work with operations modulo m. 690 // 691 // out is zeroed and may start at any size. 692 func (out *Nat) resetFor(m *Modulus) *Nat { 693 return out.reset(len(m.nat.limbs)) 694 } 695 696 // maybeSubtractModulus computes x -= m if and only if x >= m or if "always" is yes. 697 // 698 // It can be used to reduce modulo m a value up to 2m - 1, which is a common 699 // range for results computed by higher level operations. 700 // 701 // always is usually a carry that indicates that the operation that produced x 702 // overflowed its size, meaning abstractly x > 2^(_W * n) > m even if x < m. 703 // 704 // x and m operands must have the same announced length. 705 // 706 //go:norace 707 func (x *Nat) maybeSubtractModulus(always choice, m *Modulus) { 708 t := NewNat().set(x) 709 underflow := t.sub(m.nat) 710 // We keep the result if x - m didn't underflow (meaning x >= m) 711 // or if always was set. 712 keep := not(choice(underflow)) | choice(always) 713 x.assign(keep, t) 714 } 715 716 // Sub computes x = x - y mod m. 717 // 718 // The length of both operands must be the same as the modulus. Both operands 719 // must already be reduced modulo m. 720 // 721 //go:norace 722 func (x *Nat) Sub(y *Nat, m *Modulus) *Nat { 723 underflow := x.sub(y) 724 // If the subtraction underflowed, add m. 725 t := NewNat().set(x) 726 t.add(m.nat) 727 x.assign(choice(underflow), t) 728 return x 729 } 730 731 // SubOne computes x = x - 1 mod m. 732 // 733 // The length of x must be the same as the modulus. 734 func (x *Nat) SubOne(m *Modulus) *Nat { 735 one := NewNat().ExpandFor(m) 736 one.limbs[0] = 1 737 // Sub asks for x to be reduced modulo m, while SubOne doesn't, but when 738 // y = 1, it works, and this is an internal use. 739 return x.Sub(one, m) 740 } 741 742 // Add computes x = x + y mod m. 743 // 744 // The length of both operands must be the same as the modulus. Both operands 745 // must already be reduced modulo m. 746 // 747 //go:norace 748 func (x *Nat) Add(y *Nat, m *Modulus) *Nat { 749 overflow := x.add(y) 750 x.maybeSubtractModulus(choice(overflow), m) 751 return x 752 } 753 754 // montgomeryRepresentation calculates x = x * R mod m, with R = 2^(_W * n) and 755 // n = len(m.nat.limbs). 756 // 757 // Faster Montgomery multiplication replaces standard modular multiplication for 758 // numbers in this representation. 759 // 760 // This assumes that x is already reduced mod m. 761 func (x *Nat) montgomeryRepresentation(m *Modulus) *Nat { 762 // A Montgomery multiplication (which computes a * b / R) by R * R works out 763 // to a multiplication by R, which takes the value out of the Montgomery domain. 764 return x.montgomeryMul(x, m.rr, m) 765 } 766 767 // montgomeryReduction calculates x = x / R mod m, with R = 2^(_W * n) and 768 // n = len(m.nat.limbs). 769 // 770 // This assumes that x is already reduced mod m. 771 func (x *Nat) montgomeryReduction(m *Modulus) *Nat { 772 // By Montgomery multiplying with 1 not in Montgomery representation, we 773 // convert out back from Montgomery representation, because it works out to 774 // dividing by R. 775 one := NewNat().ExpandFor(m) 776 one.limbs[0] = 1 777 return x.montgomeryMul(x, one, m) 778 } 779 780 // montgomeryMul calculates x = a * b / R mod m, with R = 2^(_W * n) and 781 // n = len(m.nat.limbs), also known as a Montgomery multiplication. 782 // 783 // All inputs should be the same length and already reduced modulo m. 784 // x will be resized to the size of m and overwritten. 785 // 786 //go:norace 787 func (x *Nat) montgomeryMul(a *Nat, b *Nat, m *Modulus) *Nat { 788 n := len(m.nat.limbs) 789 mLimbs := m.nat.limbs[:n] 790 aLimbs := a.limbs[:n] 791 bLimbs := b.limbs[:n] 792 793 switch n { 794 default: 795 // Attempt to use a stack-allocated backing array. 796 T := make([]uint, 0, preallocLimbs*2) 797 if cap(T) < n*2 { 798 T = make([]uint, 0, n*2) 799 } 800 T = T[:n*2] 801 802 // This loop implements Word-by-Word Montgomery Multiplication, as 803 // described in Algorithm 4 (Fig. 3) of "Efficient Software 804 // Implementations of Modular Exponentiation" by Shay Gueron 805 // [https://eprint.iacr.org/2011/239.pdf]. 806 var c uint 807 for i := 0; i < n; i++ { 808 _ = T[n+i] // bounds check elimination hint 809 810 // Step 1 (T = a × b) is computed as a large pen-and-paper column 811 // multiplication of two numbers with n base-2^_W digits. If we just 812 // wanted to produce 2n-wide T, we would do 813 // 814 // for i := 0; i < n; i++ { 815 // d := bLimbs[i] 816 // T[n+i] = addMulVVW(T[i:n+i], aLimbs, d) 817 // } 818 // 819 // where d is a digit of the multiplier, T[i:n+i] is the shifted 820 // position of the product of that digit, and T[n+i] is the final carry. 821 // Note that T[i] isn't modified after processing the i-th digit. 822 // 823 // Instead of running two loops, one for Step 1 and one for Steps 2–6, 824 // the result of Step 1 is computed during the next loop. This is 825 // possible because each iteration only uses T[i] in Step 2 and then 826 // discards it in Step 6. 827 d := bLimbs[i] 828 c1 := addMulVVW(T[i:n+i], aLimbs, d) 829 830 // Step 6 is replaced by shifting the virtual window we operate 831 // over: T of the algorithm is T[i:] for us. That means that T1 in 832 // Step 2 (T mod 2^_W) is simply T[i]. k0 in Step 3 is our m0inv. 833 Y := T[i] * m.m0inv 834 835 // Step 4 and 5 add Y × m to T, which as mentioned above is stored 836 // at T[i:]. The two carries (from a × d and Y × m) are added up in 837 // the next word T[n+i], and the carry bit from that addition is 838 // brought forward to the next iteration. 839 c2 := addMulVVW(T[i:n+i], mLimbs, Y) 840 T[n+i], c = bits.Add(c1, c2, c) 841 } 842 843 // Finally for Step 7 we copy the final T window into x, and subtract m 844 // if necessary (which as explained in maybeSubtractModulus can be the 845 // case both if x >= m, or if x overflowed). 846 // 847 // The paper suggests in Section 4 that we can do an "Almost Montgomery 848 // Multiplication" by subtracting only in the overflow case, but the 849 // cost is very similar since the constant time subtraction tells us if 850 // x >= m as a side effect, and taking care of the broken invariant is 851 // highly undesirable (see https://go.dev/issue/13907). 852 copy(x.reset(n).limbs, T[n:]) 853 x.maybeSubtractModulus(choice(c), m) 854 855 // The following specialized cases follow the exact same algorithm, but 856 // optimized for the sizes most used in RSA. addMulVVW is implemented in 857 // assembly with loop unrolling depending on the architecture and bounds 858 // checks are removed by the compiler thanks to the constant size. 859 case 1024 / _W: 860 const n = 1024 / _W // compiler hint 861 T := make([]uint, n*2) 862 var c uint 863 for i := 0; i < n; i++ { 864 d := bLimbs[i] 865 c1 := addMulVVW1024(&T[i], &aLimbs[0], d) 866 Y := T[i] * m.m0inv 867 c2 := addMulVVW1024(&T[i], &mLimbs[0], Y) 868 T[n+i], c = bits.Add(c1, c2, c) 869 } 870 copy(x.reset(n).limbs, T[n:]) 871 x.maybeSubtractModulus(choice(c), m) 872 873 case 1536 / _W: 874 const n = 1536 / _W // compiler hint 875 T := make([]uint, n*2) 876 var c uint 877 for i := 0; i < n; i++ { 878 d := bLimbs[i] 879 c1 := addMulVVW1536(&T[i], &aLimbs[0], d) 880 Y := T[i] * m.m0inv 881 c2 := addMulVVW1536(&T[i], &mLimbs[0], Y) 882 T[n+i], c = bits.Add(c1, c2, c) 883 } 884 copy(x.reset(n).limbs, T[n:]) 885 x.maybeSubtractModulus(choice(c), m) 886 887 case 2048 / _W: 888 const n = 2048 / _W // compiler hint 889 T := make([]uint, n*2) 890 var c uint 891 for i := 0; i < n; i++ { 892 d := bLimbs[i] 893 c1 := addMulVVW2048(&T[i], &aLimbs[0], d) 894 Y := T[i] * m.m0inv 895 c2 := addMulVVW2048(&T[i], &mLimbs[0], Y) 896 T[n+i], c = bits.Add(c1, c2, c) 897 } 898 copy(x.reset(n).limbs, T[n:]) 899 x.maybeSubtractModulus(choice(c), m) 900 } 901 902 return x 903 } 904 905 // addMulVVW multiplies the multi-word value x by the single-word value y, 906 // adding the result to the multi-word value z and returning the final carry. 907 // It can be thought of as one row of a pen-and-paper column multiplication. 908 // 909 //go:norace 910 func addMulVVW(z, x []uint, y uint) (carry uint) { 911 _ = x[len(z)-1] // bounds check elimination hint 912 for i := range z { 913 hi, lo := bits.Mul(x[i], y) 914 lo, c := bits.Add(lo, z[i], 0) 915 // We use bits.Add with zero to get an add-with-carry instruction that 916 // absorbs the carry from the previous bits.Add. 917 hi, _ = bits.Add(hi, 0, c) 918 lo, c = bits.Add(lo, carry, 0) 919 hi, _ = bits.Add(hi, 0, c) 920 carry = hi 921 z[i] = lo 922 } 923 return carry 924 } 925 926 // Mul calculates x = x * y mod m. 927 // 928 // The length of both operands must be the same as the modulus. Both operands 929 // must already be reduced modulo m. 930 // 931 //go:norace 932 func (x *Nat) Mul(y *Nat, m *Modulus) *Nat { 933 if m.odd { 934 // A Montgomery multiplication by a value out of the Montgomery domain 935 // takes the result out of Montgomery representation. 936 xR := NewNat().set(x).montgomeryRepresentation(m) // xR = x * R mod m 937 return x.montgomeryMul(xR, y, m) // x = xR * y / R mod m 938 } 939 940 n := len(m.nat.limbs) 941 xLimbs := x.limbs[:n] 942 yLimbs := y.limbs[:n] 943 944 switch n { 945 default: 946 // Attempt to use a stack-allocated backing array. 947 T := make([]uint, 0, preallocLimbs*2) 948 if cap(T) < n*2 { 949 T = make([]uint, 0, n*2) 950 } 951 T = T[:n*2] 952 953 // T = x * y 954 for i := 0; i < n; i++ { 955 T[n+i] = addMulVVW(T[i:n+i], xLimbs, yLimbs[i]) 956 } 957 958 // x = T mod m 959 return x.Mod(&Nat{limbs: T}, m) 960 961 // The following specialized cases follow the exact same algorithm, but 962 // optimized for the sizes most used in RSA. See montgomeryMul for details. 963 case 1024 / _W: 964 const n = 1024 / _W // compiler hint 965 T := make([]uint, n*2) 966 for i := 0; i < n; i++ { 967 T[n+i] = addMulVVW1024(&T[i], &xLimbs[0], yLimbs[i]) 968 } 969 return x.Mod(&Nat{limbs: T}, m) 970 case 1536 / _W: 971 const n = 1536 / _W // compiler hint 972 T := make([]uint, n*2) 973 for i := 0; i < n; i++ { 974 T[n+i] = addMulVVW1536(&T[i], &xLimbs[0], yLimbs[i]) 975 } 976 return x.Mod(&Nat{limbs: T}, m) 977 case 2048 / _W: 978 const n = 2048 / _W // compiler hint 979 T := make([]uint, n*2) 980 for i := 0; i < n; i++ { 981 T[n+i] = addMulVVW2048(&T[i], &xLimbs[0], yLimbs[i]) 982 } 983 return x.Mod(&Nat{limbs: T}, m) 984 } 985 } 986 987 // Exp calculates out = x^e mod m. 988 // 989 // The exponent e is represented in big-endian order. The output will be resized 990 // to the size of m and overwritten. x must already be reduced modulo m. 991 // 992 // m must be odd, or Exp will panic. 993 // 994 //go:norace 995 func (out *Nat) Exp(x *Nat, e []byte, m *Modulus) *Nat { 996 if !m.odd { 997 panic("bigmod: modulus for Exp must be odd") 998 } 999 1000 // We use a 4 bit window. For our RSA workload, 4 bit windows are faster 1001 // than 2 bit windows, but use an extra 12 nats worth of scratch space. 1002 // Using bit sizes that don't divide 8 are more complex to implement, but 1003 // are likely to be more efficient if necessary. 1004 1005 table := [(1 << 4) - 1]*Nat{ // table[i] = x ^ (i+1) 1006 // newNat calls are unrolled so they are allocated on the stack. 1007 NewNat(), NewNat(), NewNat(), NewNat(), NewNat(), 1008 NewNat(), NewNat(), NewNat(), NewNat(), NewNat(), 1009 NewNat(), NewNat(), NewNat(), NewNat(), NewNat(), 1010 } 1011 table[0].set(x).montgomeryRepresentation(m) 1012 for i := 1; i < len(table); i++ { 1013 table[i].montgomeryMul(table[i-1], table[0], m) 1014 } 1015 1016 out.resetFor(m) 1017 out.limbs[0] = 1 1018 out.montgomeryRepresentation(m) 1019 tmp := NewNat().ExpandFor(m) 1020 for _, b := range e { 1021 for _, j := range []int{4, 0} { 1022 // Square four times. Optimization note: this can be implemented 1023 // more efficiently than with generic Montgomery multiplication. 1024 out.montgomeryMul(out, out, m) 1025 out.montgomeryMul(out, out, m) 1026 out.montgomeryMul(out, out, m) 1027 out.montgomeryMul(out, out, m) 1028 1029 // Select x^k in constant time from the table. 1030 k := uint((b >> j) & 0b1111) 1031 for i := range table { 1032 tmp.assign(ctEq(k, uint(i+1)), table[i]) 1033 } 1034 1035 // Multiply by x^k, discarding the result if k = 0. 1036 tmp.montgomeryMul(out, tmp, m) 1037 out.assign(not(ctEq(k, 0)), tmp) 1038 } 1039 } 1040 1041 return out.montgomeryReduction(m) 1042 } 1043 1044 // ExpShortVarTime calculates out = x^e mod m. 1045 // 1046 // The output will be resized to the size of m and overwritten. x must already 1047 // be reduced modulo m. This leaks the exponent through timing side-channels. 1048 // 1049 // m must be odd, or ExpShortVarTime will panic. 1050 func (out *Nat) ExpShortVarTime(x *Nat, e uint, m *Modulus) *Nat { 1051 if !m.odd { 1052 panic("bigmod: modulus for ExpShortVarTime must be odd") 1053 } 1054 // For short exponents, precomputing a table and using a window like in Exp 1055 // doesn't pay off. Instead, we do a simple conditional square-and-multiply 1056 // chain, skipping the initial run of zeroes. 1057 xR := NewNat().set(x).montgomeryRepresentation(m) 1058 out.set(xR) 1059 for i := bits.UintSize - bits.Len(e) + 1; i < bits.UintSize; i++ { 1060 out.montgomeryMul(out, out, m) 1061 if k := (e >> (bits.UintSize - i - 1)) & 1; k != 0 { 1062 out.montgomeryMul(out, xR, m) 1063 } 1064 } 1065 return out.montgomeryReduction(m) 1066 } 1067 1068 // InverseVarTime calculates x = a⁻¹ mod m and returns (x, true) if a is 1069 // invertible. Otherwise, InverseVarTime returns (x, false) and x is not 1070 // modified. 1071 // 1072 // a must be reduced modulo m, but doesn't need to have the same size. The 1073 // output will be resized to the size of m and overwritten. 1074 // 1075 //go:norace 1076 func (x *Nat) InverseVarTime(a *Nat, m *Modulus) (*Nat, bool) { 1077 u, A, err := extendedGCD(a, m.nat) 1078 if err != nil { 1079 return x, false 1080 } 1081 if u.IsOne() == no { 1082 return x, false 1083 } 1084 return x.set(A), true 1085 } 1086 1087 // GCDVarTime calculates x = GCD(a, b) where at least one of a or b is odd, and 1088 // both are non-zero. If GCDVarTime returns an error, x is not modified. 1089 // 1090 // The output will be resized to the size of the larger of a and b. 1091 func (x *Nat) GCDVarTime(a, b *Nat) (*Nat, error) { 1092 u, _, err := extendedGCD(a, b) 1093 if err != nil { 1094 return nil, err 1095 } 1096 return x.set(u), nil 1097 } 1098 1099 // extendedGCD computes u = GCD(a, m). Additionally, if a < m, it computes A 1100 // such that u = A*a - B*m. If a >= m, A is undefined. 1101 // 1102 // u will have the size of the larger of a and m, and A will have the size of m. 1103 // 1104 // It is an error if either a or m is zero, or if they are both even. 1105 func extendedGCD(a, m *Nat) (u, A *Nat, err error) { 1106 // This is the extended binary GCD algorithm described in the Handbook of 1107 // Applied Cryptography, Algorithm 14.61, adapted by BoringSSL to bound 1108 // coefficients and avoid negative numbers. For more details and proof of 1109 // correctness, see https://github.com/mit-plv/fiat-crypto/pull/333/files. 1110 // 1111 // Following the proof linked in the PR above, the changes are: 1112 // 1113 // 1. Negate [B] and [C] so they are positive. The invariant now involves a 1114 // subtraction. 1115 // 2. If step 2 (both [x] and [y] are even) runs, abort immediately. This 1116 // case needs to be handled by the caller. 1117 // 3. Subtract copies of [x] and [y] as needed in step 6 (both [u] and [v] 1118 // are odd) so coefficients stay in bounds. 1119 // 4. Replace the [u >= v] check with [u > v]. This changes the end 1120 // condition to [v = 0] rather than [u = 0]. This saves an extra 1121 // subtraction due to which coefficients were negated. 1122 // 5. Rename x and y to a and n, to capture that one is a modulus. 1123 // 6. Rearrange steps 4 through 6 slightly. Merge the loops in steps 4 and 1124 // 5 into the main loop (step 7's goto), and move step 6 to the start of 1125 // the loop iteration, ensuring each loop iteration halves at least one 1126 // value. 1127 // 1128 // Note this algorithm does not handle either input being zero. 1129 // 1130 // See https://go.dev/issue/78218 for a Gobra proof of this implementation. 1131 1132 if a.IsZero() == yes || m.IsZero() == yes { 1133 return nil, nil, errors.New("extendedGCD: a or m is zero") 1134 } 1135 if a.IsOdd() == no && m.IsOdd() == no { 1136 return nil, nil, errors.New("extendedGCD: both a and m are even") 1137 } 1138 1139 size := max(len(a.limbs), len(m.limbs)) 1140 u = NewNat().set(a).expand(size) 1141 v := NewNat().set(m).expand(size) 1142 1143 A = NewNat().reset(len(m.limbs)) 1144 A.limbs[0] = 1 1145 B := NewNat().reset(len(a.limbs)) 1146 C := NewNat().reset(len(m.limbs)) 1147 D := NewNat().reset(len(a.limbs)) 1148 D.limbs[0] = 1 1149 1150 // Before and after each loop iteration, the following hold: 1151 // 1152 // 0 < m 1153 // 0 < u <= a 1154 // 0 <= v <= m 1155 // a or m is odd 1156 // u or v is odd 1157 // gcd(u, v) = gcd(a, m) 1158 // 1159 // If a < m, then the following also hold: 1160 // 1161 // 0 <= A < m 1162 // 0 <= B < a 1163 // 0 <= C < m 1164 // 0 <= D <= a 1165 // u = A*a - B*m 1166 // v = D*m - C*a 1167 // 1168 // After each loop iteration, u + v only gets smaller, and at least one of 1169 // u and v shrinks by at least a factor of two. 1170 for { 1171 // If both u and v are odd, subtract the smaller from the larger. 1172 // If u = v, we need to subtract from v to hit the modified exit condition. 1173 if u.IsOdd() == yes && v.IsOdd() == yes { 1174 if v.cmpGeq(u) == no { 1175 u.sub(v) 1176 syncAdd(A, C, B, D, m, a) 1177 } else { 1178 v.sub(u) 1179 syncAdd(C, A, D, B, m, a) 1180 } 1181 } 1182 1183 // Exactly one of u and v is now even. 1184 if u.IsOdd() == v.IsOdd() { 1185 panic("bigmod: internal error: u and v are not in the expected state") 1186 } 1187 1188 // Halve the even one and adjust the corresponding coefficient. 1189 if u.IsOdd() == no { 1190 rshift1(u, 0) 1191 if A.IsOdd() == yes || B.IsOdd() == yes { 1192 rshift1(A, A.add(m)) 1193 rshift1(B, B.add(a)) 1194 } else { 1195 rshift1(A, 0) 1196 rshift1(B, 0) 1197 } 1198 } else { // v.IsOdd() == no 1199 rshift1(v, 0) 1200 if C.IsOdd() == yes || D.IsOdd() == yes { 1201 rshift1(C, C.add(m)) 1202 rshift1(D, D.add(a)) 1203 } else { 1204 rshift1(C, 0) 1205 rshift1(D, 0) 1206 } 1207 } 1208 1209 if v.IsZero() == yes { 1210 // Base case: v = 0 -> gcd(a, m) = gcd(u, 0) = u. 1211 return u, A, nil 1212 } 1213 } 1214 } 1215 1216 // syncAdd adds Y to X and W to Z, then subtracts m from X and a from Z if 1217 // X + Y >= m. This is synchronized single-subtraction modular reduction: 1218 // X = (X + Y) mod m, with Z tracking the same wrap/no-wrap. 1219 // 1220 //go:norace 1221 func syncAdd(X, Y, Z, W, m, a *Nat) { 1222 c := X.add(Y) 1223 Z.add(W) 1224 1225 // Like in maybeSubtractModulus, we need the subtraction if either 1226 // X + Y >= m, or if X + Y overflowed (meaning X + Y >= 2^(_W * n) > m). 1227 if choice(c) == yes || X.cmpGeq(m) == yes { 1228 X.sub(m) 1229 Z.sub(a) 1230 } 1231 } 1232 1233 //go:norace 1234 func rshift1(a *Nat, carry uint) { 1235 size := len(a.limbs) 1236 aLimbs := a.limbs[:size] 1237 1238 for i := range size { 1239 aLimbs[i] >>= 1 1240 if i+1 < size { 1241 aLimbs[i] |= aLimbs[i+1] << (_W - 1) 1242 } else { 1243 aLimbs[i] |= carry << (_W - 1) 1244 } 1245 } 1246 } 1247 1248 // ShiftRightByOne sets x = x >> 1. 1249 // 1250 // The announced length of x is unchanged. 1251 // 1252 //go:norace 1253 func (x *Nat) ShiftRightByOne() *Nat { 1254 rshift1(x, 0) 1255 return x 1256 } 1257 1258 // DivShortVarTime calculates x = x / y and returns the remainder. 1259 // 1260 // It panics if y is zero. 1261 // 1262 //go:norace 1263 func (x *Nat) DivShortVarTime(y uint) uint { 1264 if y == 0 { 1265 panic("bigmod: division by zero") 1266 } 1267 1268 var r uint 1269 for i := len(x.limbs) - 1; i >= 0; i-- { 1270 x.limbs[i], r = bits.Div(r, x.limbs[i], y) 1271 } 1272 return r 1273 } 1274