Source file src/crypto/internal/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 "encoding/binary" 9 "errors" 10 "math/big" 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 // choice represents a constant-time boolean. The value of choice is always 22 // either 1 or 0. We use an int instead of bool in order to make decisions in 23 // constant time by turning it into a mask. 24 type choice uint 25 26 func not(c choice) choice { return 1 ^ c } 27 28 const yes = choice(1) 29 const no = choice(0) 30 31 // ctMask is all 1s if on is yes, and all 0s otherwise. 32 func ctMask(on choice) uint { return -uint(on) } 33 34 // ctEq returns 1 if x == y, and 0 otherwise. The execution time of this 35 // function does not depend on its inputs. 36 func ctEq(x, y uint) choice { 37 // If x != y, then either x - y or y - x will generate a carry. 38 _, c1 := bits.Sub(x, y, 0) 39 _, c2 := bits.Sub(y, x, 0) 40 return not(choice(c1 | c2)) 41 } 42 43 // ctGeq returns 1 if x >= y, and 0 otherwise. The execution time of this 44 // function does not depend on its inputs. 45 func ctGeq(x, y uint) choice { 46 // If x < y, then x - y generates a carry. 47 _, carry := bits.Sub(x, y, 0) 48 return not(choice(carry)) 49 } 50 51 // Nat represents an arbitrary natural number 52 // 53 // Each Nat has an announced length, which is the number of limbs it has stored. 54 // Operations on this number are allowed to leak this length, but will not leak 55 // any information about the values contained in those limbs. 56 type Nat struct { 57 // limbs is little-endian in base 2^W with W = bits.UintSize. 58 limbs []uint 59 } 60 61 // preallocTarget is the size in bits of the numbers used to implement the most 62 // common and most performant RSA key size. It's also enough to cover some of 63 // the operations of key sizes up to 4096. 64 const preallocTarget = 2048 65 const preallocLimbs = (preallocTarget + _W - 1) / _W 66 67 // NewNat returns a new nat with a size of zero, just like new(Nat), but with 68 // the preallocated capacity to hold a number of up to preallocTarget bits. 69 // NewNat inlines, so the allocation can live on the stack. 70 func NewNat() *Nat { 71 limbs := make([]uint, 0, preallocLimbs) 72 return &Nat{limbs} 73 } 74 75 // expand expands x to n limbs, leaving its value unchanged. 76 func (x *Nat) expand(n int) *Nat { 77 if len(x.limbs) > n { 78 panic("bigmod: internal error: shrinking nat") 79 } 80 if cap(x.limbs) < n { 81 newLimbs := make([]uint, n) 82 copy(newLimbs, x.limbs) 83 x.limbs = newLimbs 84 return x 85 } 86 extraLimbs := x.limbs[len(x.limbs):n] 87 for i := range extraLimbs { 88 extraLimbs[i] = 0 89 } 90 x.limbs = x.limbs[:n] 91 return x 92 } 93 94 // reset returns a zero nat of n limbs, reusing x's storage if n <= cap(x.limbs). 95 func (x *Nat) reset(n int) *Nat { 96 if cap(x.limbs) < n { 97 x.limbs = make([]uint, n) 98 return x 99 } 100 for i := range x.limbs { 101 x.limbs[i] = 0 102 } 103 x.limbs = x.limbs[:n] 104 return x 105 } 106 107 // set assigns x = y, optionally resizing x to the appropriate size. 108 func (x *Nat) set(y *Nat) *Nat { 109 x.reset(len(y.limbs)) 110 copy(x.limbs, y.limbs) 111 return x 112 } 113 114 // setBig assigns x = n, optionally resizing n to the appropriate size. 115 // 116 // The announced length of x is set based on the actual bit size of the input, 117 // ignoring leading zeroes. 118 func (x *Nat) setBig(n *big.Int) *Nat { 119 limbs := n.Bits() 120 x.reset(len(limbs)) 121 for i := range limbs { 122 x.limbs[i] = uint(limbs[i]) 123 } 124 return x 125 } 126 127 // Bytes returns x as a zero-extended big-endian byte slice. The size of the 128 // slice will match the size of m. 129 // 130 // x must have the same size as m and it must be reduced modulo m. 131 func (x *Nat) Bytes(m *Modulus) []byte { 132 i := m.Size() 133 bytes := make([]byte, i) 134 for _, limb := range x.limbs { 135 for j := 0; j < _S; j++ { 136 i-- 137 if i < 0 { 138 if limb == 0 { 139 break 140 } 141 panic("bigmod: modulus is smaller than nat") 142 } 143 bytes[i] = byte(limb) 144 limb >>= 8 145 } 146 } 147 return bytes 148 } 149 150 // SetBytes assigns x = b, where b is a slice of big-endian bytes. 151 // SetBytes returns an error if b >= m. 152 // 153 // The output will be resized to the size of m and overwritten. 154 func (x *Nat) SetBytes(b []byte, m *Modulus) (*Nat, error) { 155 if err := x.setBytes(b, m); err != nil { 156 return nil, err 157 } 158 if x.cmpGeq(m.nat) == yes { 159 return nil, errors.New("input overflows the modulus") 160 } 161 return x, nil 162 } 163 164 // SetOverflowingBytes assigns x = b, where b is a slice of big-endian bytes. 165 // SetOverflowingBytes returns an error if b has a longer bit length than m, but 166 // reduces overflowing values up to 2^⌈log2(m)⌉ - 1. 167 // 168 // The output will be resized to the size of m and overwritten. 169 func (x *Nat) SetOverflowingBytes(b []byte, m *Modulus) (*Nat, error) { 170 if err := x.setBytes(b, m); err != nil { 171 return nil, err 172 } 173 leading := _W - bitLen(x.limbs[len(x.limbs)-1]) 174 if leading < m.leading { 175 return nil, errors.New("input overflows the modulus size") 176 } 177 x.maybeSubtractModulus(no, m) 178 return x, nil 179 } 180 181 // bigEndianUint returns the contents of buf interpreted as a 182 // big-endian encoded uint value. 183 func bigEndianUint(buf []byte) uint { 184 if _W == 64 { 185 return uint(binary.BigEndian.Uint64(buf)) 186 } 187 return uint(binary.BigEndian.Uint32(buf)) 188 } 189 190 func (x *Nat) setBytes(b []byte, m *Modulus) error { 191 x.resetFor(m) 192 i, k := len(b), 0 193 for k < len(x.limbs) && i >= _S { 194 x.limbs[k] = bigEndianUint(b[i-_S : i]) 195 i -= _S 196 k++ 197 } 198 for s := 0; s < _W && k < len(x.limbs) && i > 0; s += 8 { 199 x.limbs[k] |= uint(b[i-1]) << s 200 i-- 201 } 202 if i > 0 { 203 return errors.New("input overflows the modulus size") 204 } 205 return nil 206 } 207 208 // Equal returns 1 if x == y, and 0 otherwise. 209 // 210 // Both operands must have the same announced length. 211 func (x *Nat) Equal(y *Nat) choice { 212 // Eliminate bounds checks in the loop. 213 size := len(x.limbs) 214 xLimbs := x.limbs[:size] 215 yLimbs := y.limbs[:size] 216 217 equal := yes 218 for i := 0; i < size; i++ { 219 equal &= ctEq(xLimbs[i], yLimbs[i]) 220 } 221 return equal 222 } 223 224 // IsZero returns 1 if x == 0, and 0 otherwise. 225 func (x *Nat) IsZero() choice { 226 // Eliminate bounds checks in the loop. 227 size := len(x.limbs) 228 xLimbs := x.limbs[:size] 229 230 zero := yes 231 for i := 0; i < size; i++ { 232 zero &= ctEq(xLimbs[i], 0) 233 } 234 return zero 235 } 236 237 // cmpGeq returns 1 if x >= y, and 0 otherwise. 238 // 239 // Both operands must have the same announced length. 240 func (x *Nat) cmpGeq(y *Nat) choice { 241 // Eliminate bounds checks in the loop. 242 size := len(x.limbs) 243 xLimbs := x.limbs[:size] 244 yLimbs := y.limbs[:size] 245 246 var c uint 247 for i := 0; i < size; i++ { 248 _, c = bits.Sub(xLimbs[i], yLimbs[i], c) 249 } 250 // If there was a carry, then subtracting y underflowed, so 251 // x is not greater than or equal to y. 252 return not(choice(c)) 253 } 254 255 // assign sets x <- y if on == 1, and does nothing otherwise. 256 // 257 // Both operands must have the same announced length. 258 func (x *Nat) assign(on choice, y *Nat) *Nat { 259 // Eliminate bounds checks in the loop. 260 size := len(x.limbs) 261 xLimbs := x.limbs[:size] 262 yLimbs := y.limbs[:size] 263 264 mask := ctMask(on) 265 for i := 0; i < size; i++ { 266 xLimbs[i] ^= mask & (xLimbs[i] ^ yLimbs[i]) 267 } 268 return x 269 } 270 271 // add computes x += y and returns the carry. 272 // 273 // Both operands must have the same announced length. 274 func (x *Nat) add(y *Nat) (c uint) { 275 // Eliminate bounds checks in the loop. 276 size := len(x.limbs) 277 xLimbs := x.limbs[:size] 278 yLimbs := y.limbs[:size] 279 280 for i := 0; i < size; i++ { 281 xLimbs[i], c = bits.Add(xLimbs[i], yLimbs[i], c) 282 } 283 return 284 } 285 286 // sub computes x -= y. It returns the borrow of the subtraction. 287 // 288 // Both operands must have the same announced length. 289 func (x *Nat) sub(y *Nat) (c uint) { 290 // Eliminate bounds checks in the loop. 291 size := len(x.limbs) 292 xLimbs := x.limbs[:size] 293 yLimbs := y.limbs[:size] 294 295 for i := 0; i < size; i++ { 296 xLimbs[i], c = bits.Sub(xLimbs[i], yLimbs[i], c) 297 } 298 return 299 } 300 301 // Modulus is used for modular arithmetic, precomputing relevant constants. 302 // 303 // Moduli are assumed to be odd numbers. Moduli can also leak the exact 304 // number of bits needed to store their value, and are stored without padding. 305 // 306 // Their actual value is still kept secret. 307 type Modulus struct { 308 // The underlying natural number for this modulus. 309 // 310 // This will be stored without any padding, and shouldn't alias with any 311 // other natural number being used. 312 nat *Nat 313 leading int // number of leading zeros in the modulus 314 m0inv uint // -nat.limbs[0]⁻¹ mod _W 315 rr *Nat // R*R for montgomeryRepresentation 316 } 317 318 // rr returns R*R with R = 2^(_W * n) and n = len(m.nat.limbs). 319 func rr(m *Modulus) *Nat { 320 rr := NewNat().ExpandFor(m) 321 n := uint(len(rr.limbs)) 322 mLen := uint(m.BitLen()) 323 logR := _W * n 324 325 // We start by computing R = 2^(_W * n) mod m. We can get pretty close, to 326 // 2^⌊log₂m⌋, by setting the highest bit we can without having to reduce. 327 rr.limbs[n-1] = 1 << ((mLen - 1) % _W) 328 // Then we double until we reach 2^(_W * n). 329 for i := mLen - 1; i < logR; i++ { 330 rr.Add(rr, m) 331 } 332 333 // Next we need to get from R to 2^(_W * n) R mod m (aka from one to R in 334 // the Montgomery domain, meaning we can use Montgomery multiplication now). 335 // We could do that by doubling _W * n times, or with a square-and-double 336 // chain log2(_W * n) long. Turns out the fastest thing is to start out with 337 // doublings, and switch to square-and-double once the exponent is large 338 // enough to justify the cost of the multiplications. 339 340 // The threshold is selected experimentally as a linear function of n. 341 threshold := n / 4 342 343 // We calculate how many of the most-significant bits of the exponent we can 344 // compute before crossing the threshold, and we do it with doublings. 345 i := bits.UintSize 346 for logR>>i <= threshold { 347 i-- 348 } 349 for k := uint(0); k < logR>>i; k++ { 350 rr.Add(rr, m) 351 } 352 353 // Then we process the remaining bits of the exponent with a 354 // square-and-double chain. 355 for i > 0 { 356 rr.montgomeryMul(rr, rr, m) 357 i-- 358 if logR>>i&1 != 0 { 359 rr.Add(rr, m) 360 } 361 } 362 363 return rr 364 } 365 366 // minusInverseModW computes -x⁻¹ mod _W with x odd. 367 // 368 // This operation is used to precompute a constant involved in Montgomery 369 // multiplication. 370 func minusInverseModW(x uint) uint { 371 // Every iteration of this loop doubles the least-significant bits of 372 // correct inverse in y. The first three bits are already correct (1⁻¹ = 1, 373 // 3⁻¹ = 3, 5⁻¹ = 5, and 7⁻¹ = 7 mod 8), so doubling five times is enough 374 // for 64 bits (and wastes only one iteration for 32 bits). 375 // 376 // See https://crypto.stackexchange.com/a/47496. 377 y := x 378 for i := 0; i < 5; i++ { 379 y = y * (2 - x*y) 380 } 381 return -y 382 } 383 384 // NewModulusFromBig creates a new Modulus from a [big.Int]. 385 // 386 // The Int must be odd. The number of significant bits (and nothing else) is 387 // leaked through timing side-channels. 388 func NewModulusFromBig(n *big.Int) (*Modulus, error) { 389 if b := n.Bits(); len(b) == 0 { 390 return nil, errors.New("modulus must be >= 0") 391 } else if b[0]&1 != 1 { 392 return nil, errors.New("modulus must be odd") 393 } 394 m := &Modulus{} 395 m.nat = NewNat().setBig(n) 396 m.leading = _W - bitLen(m.nat.limbs[len(m.nat.limbs)-1]) 397 m.m0inv = minusInverseModW(m.nat.limbs[0]) 398 m.rr = rr(m) 399 return m, nil 400 } 401 402 // bitLen is a version of bits.Len that only leaks the bit length of n, but not 403 // its value. bits.Len and bits.LeadingZeros use a lookup table for the 404 // low-order bits on some architectures. 405 func bitLen(n uint) int { 406 var len int 407 // We assume, here and elsewhere, that comparison to zero is constant time 408 // with respect to different non-zero values. 409 for n != 0 { 410 len++ 411 n >>= 1 412 } 413 return len 414 } 415 416 // Size returns the size of m in bytes. 417 func (m *Modulus) Size() int { 418 return (m.BitLen() + 7) / 8 419 } 420 421 // BitLen returns the size of m in bits. 422 func (m *Modulus) BitLen() int { 423 return len(m.nat.limbs)*_W - int(m.leading) 424 } 425 426 // Nat returns m as a Nat. The return value must not be written to. 427 func (m *Modulus) Nat() *Nat { 428 return m.nat 429 } 430 431 // shiftIn calculates x = x << _W + y mod m. 432 // 433 // This assumes that x is already reduced mod m. 434 func (x *Nat) shiftIn(y uint, m *Modulus) *Nat { 435 d := NewNat().resetFor(m) 436 437 // Eliminate bounds checks in the loop. 438 size := len(m.nat.limbs) 439 xLimbs := x.limbs[:size] 440 dLimbs := d.limbs[:size] 441 mLimbs := m.nat.limbs[:size] 442 443 // Each iteration of this loop computes x = 2x + b mod m, where b is a bit 444 // from y. Effectively, it left-shifts x and adds y one bit at a time, 445 // reducing it every time. 446 // 447 // To do the reduction, each iteration computes both 2x + b and 2x + b - m. 448 // The next iteration (and finally the return line) will use either result 449 // based on whether 2x + b overflows m. 450 needSubtraction := no 451 for i := _W - 1; i >= 0; i-- { 452 carry := (y >> i) & 1 453 var borrow uint 454 mask := ctMask(needSubtraction) 455 for i := 0; i < size; i++ { 456 l := xLimbs[i] ^ (mask & (xLimbs[i] ^ dLimbs[i])) 457 xLimbs[i], carry = bits.Add(l, l, carry) 458 dLimbs[i], borrow = bits.Sub(xLimbs[i], mLimbs[i], borrow) 459 } 460 // Like in maybeSubtractModulus, we need the subtraction if either it 461 // didn't underflow (meaning 2x + b > m) or if computing 2x + b 462 // overflowed (meaning 2x + b > 2^_W*n > m). 463 needSubtraction = not(choice(borrow)) | choice(carry) 464 } 465 return x.assign(needSubtraction, d) 466 } 467 468 // Mod calculates out = x mod m. 469 // 470 // This works regardless how large the value of x is. 471 // 472 // The output will be resized to the size of m and overwritten. 473 func (out *Nat) Mod(x *Nat, m *Modulus) *Nat { 474 out.resetFor(m) 475 // Working our way from the most significant to the least significant limb, 476 // we can insert each limb at the least significant position, shifting all 477 // previous limbs left by _W. This way each limb will get shifted by the 478 // correct number of bits. We can insert at least N - 1 limbs without 479 // overflowing m. After that, we need to reduce every time we shift. 480 i := len(x.limbs) - 1 481 // For the first N - 1 limbs we can skip the actual shifting and position 482 // them at the shifted position, which starts at min(N - 2, i). 483 start := len(m.nat.limbs) - 2 484 if i < start { 485 start = i 486 } 487 for j := start; j >= 0; j-- { 488 out.limbs[j] = x.limbs[i] 489 i-- 490 } 491 // We shift in the remaining limbs, reducing modulo m each time. 492 for i >= 0 { 493 out.shiftIn(x.limbs[i], m) 494 i-- 495 } 496 return out 497 } 498 499 // ExpandFor ensures x has the right size to work with operations modulo m. 500 // 501 // The announced size of x must be smaller than or equal to that of m. 502 func (x *Nat) ExpandFor(m *Modulus) *Nat { 503 return x.expand(len(m.nat.limbs)) 504 } 505 506 // resetFor ensures out has the right size to work with operations modulo m. 507 // 508 // out is zeroed and may start at any size. 509 func (out *Nat) resetFor(m *Modulus) *Nat { 510 return out.reset(len(m.nat.limbs)) 511 } 512 513 // maybeSubtractModulus computes x -= m if and only if x >= m or if "always" is yes. 514 // 515 // It can be used to reduce modulo m a value up to 2m - 1, which is a common 516 // range for results computed by higher level operations. 517 // 518 // always is usually a carry that indicates that the operation that produced x 519 // overflowed its size, meaning abstractly x > 2^_W*n > m even if x < m. 520 // 521 // x and m operands must have the same announced length. 522 func (x *Nat) maybeSubtractModulus(always choice, m *Modulus) { 523 t := NewNat().set(x) 524 underflow := t.sub(m.nat) 525 // We keep the result if x - m didn't underflow (meaning x >= m) 526 // or if always was set. 527 keep := not(choice(underflow)) | choice(always) 528 x.assign(keep, t) 529 } 530 531 // Sub computes x = x - y mod m. 532 // 533 // The length of both operands must be the same as the modulus. Both operands 534 // must already be reduced modulo m. 535 func (x *Nat) Sub(y *Nat, m *Modulus) *Nat { 536 underflow := x.sub(y) 537 // If the subtraction underflowed, add m. 538 t := NewNat().set(x) 539 t.add(m.nat) 540 x.assign(choice(underflow), t) 541 return x 542 } 543 544 // Add computes x = x + y mod m. 545 // 546 // The length of both operands must be the same as the modulus. Both operands 547 // must already be reduced modulo m. 548 func (x *Nat) Add(y *Nat, m *Modulus) *Nat { 549 overflow := x.add(y) 550 x.maybeSubtractModulus(choice(overflow), m) 551 return x 552 } 553 554 // montgomeryRepresentation calculates x = x * R mod m, with R = 2^(_W * n) and 555 // n = len(m.nat.limbs). 556 // 557 // Faster Montgomery multiplication replaces standard modular multiplication for 558 // numbers in this representation. 559 // 560 // This assumes that x is already reduced mod m. 561 func (x *Nat) montgomeryRepresentation(m *Modulus) *Nat { 562 // A Montgomery multiplication (which computes a * b / R) by R * R works out 563 // to a multiplication by R, which takes the value out of the Montgomery domain. 564 return x.montgomeryMul(x, m.rr, m) 565 } 566 567 // montgomeryReduction calculates x = x / R mod m, with R = 2^(_W * n) and 568 // n = len(m.nat.limbs). 569 // 570 // This assumes that x is already reduced mod m. 571 func (x *Nat) montgomeryReduction(m *Modulus) *Nat { 572 // By Montgomery multiplying with 1 not in Montgomery representation, we 573 // convert out back from Montgomery representation, because it works out to 574 // dividing by R. 575 one := NewNat().ExpandFor(m) 576 one.limbs[0] = 1 577 return x.montgomeryMul(x, one, m) 578 } 579 580 // montgomeryMul calculates x = a * b / R mod m, with R = 2^(_W * n) and 581 // n = len(m.nat.limbs), also known as a Montgomery multiplication. 582 // 583 // All inputs should be the same length and already reduced modulo m. 584 // x will be resized to the size of m and overwritten. 585 func (x *Nat) montgomeryMul(a *Nat, b *Nat, m *Modulus) *Nat { 586 n := len(m.nat.limbs) 587 mLimbs := m.nat.limbs[:n] 588 aLimbs := a.limbs[:n] 589 bLimbs := b.limbs[:n] 590 591 switch n { 592 default: 593 // Attempt to use a stack-allocated backing array. 594 T := make([]uint, 0, preallocLimbs*2) 595 if cap(T) < n*2 { 596 T = make([]uint, 0, n*2) 597 } 598 T = T[:n*2] 599 600 // This loop implements Word-by-Word Montgomery Multiplication, as 601 // described in Algorithm 4 (Fig. 3) of "Efficient Software 602 // Implementations of Modular Exponentiation" by Shay Gueron 603 // [https://eprint.iacr.org/2011/239.pdf]. 604 var c uint 605 for i := 0; i < n; i++ { 606 _ = T[n+i] // bounds check elimination hint 607 608 // Step 1 (T = a × b) is computed as a large pen-and-paper column 609 // multiplication of two numbers with n base-2^_W digits. If we just 610 // wanted to produce 2n-wide T, we would do 611 // 612 // for i := 0; i < n; i++ { 613 // d := bLimbs[i] 614 // T[n+i] = addMulVVW(T[i:n+i], aLimbs, d) 615 // } 616 // 617 // where d is a digit of the multiplier, T[i:n+i] is the shifted 618 // position of the product of that digit, and T[n+i] is the final carry. 619 // Note that T[i] isn't modified after processing the i-th digit. 620 // 621 // Instead of running two loops, one for Step 1 and one for Steps 2–6, 622 // the result of Step 1 is computed during the next loop. This is 623 // possible because each iteration only uses T[i] in Step 2 and then 624 // discards it in Step 6. 625 d := bLimbs[i] 626 c1 := addMulVVW(T[i:n+i], aLimbs, d) 627 628 // Step 6 is replaced by shifting the virtual window we operate 629 // over: T of the algorithm is T[i:] for us. That means that T1 in 630 // Step 2 (T mod 2^_W) is simply T[i]. k0 in Step 3 is our m0inv. 631 Y := T[i] * m.m0inv 632 633 // Step 4 and 5 add Y × m to T, which as mentioned above is stored 634 // at T[i:]. The two carries (from a × d and Y × m) are added up in 635 // the next word T[n+i], and the carry bit from that addition is 636 // brought forward to the next iteration. 637 c2 := addMulVVW(T[i:n+i], mLimbs, Y) 638 T[n+i], c = bits.Add(c1, c2, c) 639 } 640 641 // Finally for Step 7 we copy the final T window into x, and subtract m 642 // if necessary (which as explained in maybeSubtractModulus can be the 643 // case both if x >= m, or if x overflowed). 644 // 645 // The paper suggests in Section 4 that we can do an "Almost Montgomery 646 // Multiplication" by subtracting only in the overflow case, but the 647 // cost is very similar since the constant time subtraction tells us if 648 // x >= m as a side effect, and taking care of the broken invariant is 649 // highly undesirable (see https://go.dev/issue/13907). 650 copy(x.reset(n).limbs, T[n:]) 651 x.maybeSubtractModulus(choice(c), m) 652 653 // The following specialized cases follow the exact same algorithm, but 654 // optimized for the sizes most used in RSA. addMulVVW is implemented in 655 // assembly with loop unrolling depending on the architecture and bounds 656 // checks are removed by the compiler thanks to the constant size. 657 case 1024 / _W: 658 const n = 1024 / _W // compiler hint 659 T := make([]uint, n*2) 660 var c uint 661 for i := 0; i < n; i++ { 662 d := bLimbs[i] 663 c1 := addMulVVW1024(&T[i], &aLimbs[0], d) 664 Y := T[i] * m.m0inv 665 c2 := addMulVVW1024(&T[i], &mLimbs[0], Y) 666 T[n+i], c = bits.Add(c1, c2, c) 667 } 668 copy(x.reset(n).limbs, T[n:]) 669 x.maybeSubtractModulus(choice(c), m) 670 671 case 1536 / _W: 672 const n = 1536 / _W // compiler hint 673 T := make([]uint, n*2) 674 var c uint 675 for i := 0; i < n; i++ { 676 d := bLimbs[i] 677 c1 := addMulVVW1536(&T[i], &aLimbs[0], d) 678 Y := T[i] * m.m0inv 679 c2 := addMulVVW1536(&T[i], &mLimbs[0], Y) 680 T[n+i], c = bits.Add(c1, c2, c) 681 } 682 copy(x.reset(n).limbs, T[n:]) 683 x.maybeSubtractModulus(choice(c), m) 684 685 case 2048 / _W: 686 const n = 2048 / _W // compiler hint 687 T := make([]uint, n*2) 688 var c uint 689 for i := 0; i < n; i++ { 690 d := bLimbs[i] 691 c1 := addMulVVW2048(&T[i], &aLimbs[0], d) 692 Y := T[i] * m.m0inv 693 c2 := addMulVVW2048(&T[i], &mLimbs[0], Y) 694 T[n+i], c = bits.Add(c1, c2, c) 695 } 696 copy(x.reset(n).limbs, T[n:]) 697 x.maybeSubtractModulus(choice(c), m) 698 } 699 700 return x 701 } 702 703 // addMulVVW multiplies the multi-word value x by the single-word value y, 704 // adding the result to the multi-word value z and returning the final carry. 705 // It can be thought of as one row of a pen-and-paper column multiplication. 706 func addMulVVW(z, x []uint, y uint) (carry uint) { 707 _ = x[len(z)-1] // bounds check elimination hint 708 for i := range z { 709 hi, lo := bits.Mul(x[i], y) 710 lo, c := bits.Add(lo, z[i], 0) 711 // We use bits.Add with zero to get an add-with-carry instruction that 712 // absorbs the carry from the previous bits.Add. 713 hi, _ = bits.Add(hi, 0, c) 714 lo, c = bits.Add(lo, carry, 0) 715 hi, _ = bits.Add(hi, 0, c) 716 carry = hi 717 z[i] = lo 718 } 719 return carry 720 } 721 722 // Mul calculates x = x * y mod m. 723 // 724 // The length of both operands must be the same as the modulus. Both operands 725 // must already be reduced modulo m. 726 func (x *Nat) Mul(y *Nat, m *Modulus) *Nat { 727 // A Montgomery multiplication by a value out of the Montgomery domain 728 // takes the result out of Montgomery representation. 729 xR := NewNat().set(x).montgomeryRepresentation(m) // xR = x * R mod m 730 return x.montgomeryMul(xR, y, m) // x = xR * y / R mod m 731 } 732 733 // Exp calculates out = x^e mod m. 734 // 735 // The exponent e is represented in big-endian order. The output will be resized 736 // to the size of m and overwritten. x must already be reduced modulo m. 737 func (out *Nat) Exp(x *Nat, e []byte, m *Modulus) *Nat { 738 // We use a 4 bit window. For our RSA workload, 4 bit windows are faster 739 // than 2 bit windows, but use an extra 12 nats worth of scratch space. 740 // Using bit sizes that don't divide 8 are more complex to implement, but 741 // are likely to be more efficient if necessary. 742 743 table := [(1 << 4) - 1]*Nat{ // table[i] = x ^ (i+1) 744 // newNat calls are unrolled so they are allocated on the stack. 745 NewNat(), NewNat(), NewNat(), NewNat(), NewNat(), 746 NewNat(), NewNat(), NewNat(), NewNat(), NewNat(), 747 NewNat(), NewNat(), NewNat(), NewNat(), NewNat(), 748 } 749 table[0].set(x).montgomeryRepresentation(m) 750 for i := 1; i < len(table); i++ { 751 table[i].montgomeryMul(table[i-1], table[0], m) 752 } 753 754 out.resetFor(m) 755 out.limbs[0] = 1 756 out.montgomeryRepresentation(m) 757 tmp := NewNat().ExpandFor(m) 758 for _, b := range e { 759 for _, j := range []int{4, 0} { 760 // Square four times. Optimization note: this can be implemented 761 // more efficiently than with generic Montgomery multiplication. 762 out.montgomeryMul(out, out, m) 763 out.montgomeryMul(out, out, m) 764 out.montgomeryMul(out, out, m) 765 out.montgomeryMul(out, out, m) 766 767 // Select x^k in constant time from the table. 768 k := uint((b >> j) & 0b1111) 769 for i := range table { 770 tmp.assign(ctEq(k, uint(i+1)), table[i]) 771 } 772 773 // Multiply by x^k, discarding the result if k = 0. 774 tmp.montgomeryMul(out, tmp, m) 775 out.assign(not(ctEq(k, 0)), tmp) 776 } 777 } 778 779 return out.montgomeryReduction(m) 780 } 781 782 // ExpShortVarTime calculates out = x^e mod m. 783 // 784 // The output will be resized to the size of m and overwritten. x must already 785 // be reduced modulo m. This leaks the exponent through timing side-channels. 786 func (out *Nat) ExpShortVarTime(x *Nat, e uint, m *Modulus) *Nat { 787 // For short exponents, precomputing a table and using a window like in Exp 788 // doesn't pay off. Instead, we do a simple conditional square-and-multiply 789 // chain, skipping the initial run of zeroes. 790 xR := NewNat().set(x).montgomeryRepresentation(m) 791 out.set(xR) 792 for i := bits.UintSize - bitLen(e) + 1; i < bits.UintSize; i++ { 793 out.montgomeryMul(out, out, m) 794 if k := (e >> (bits.UintSize - i - 1)) & 1; k != 0 { 795 out.montgomeryMul(out, xR, m) 796 } 797 } 798 return out.montgomeryReduction(m) 799 } 800