// Copyright 2009 The Go Authors. All rights reserved.// Use of this source code is governed by a BSD-style// license that can be found in the LICENSE file.// This file implements unsigned multi-precision integers (natural// numbers). They are the building blocks for the implementation// of signed integers, rationals, and floating-point numbers.//// Caution: This implementation relies on the function "alias"// which assumes that (nat) slice capacities are never// changed (no 3-operand slice expressions). If that// changes, alias needs to be updated for correctness.package bigimport ()// An unsigned integer x of the form//// x = x[n-1]*_B^(n-1) + x[n-2]*_B^(n-2) + ... + x[1]*_B + x[0]//// with 0 <= x[i] < _B and 0 <= i < n is stored in a slice of length n,// with the digits x[i] as the slice elements.//// A number is normalized if the slice contains no leading 0 digits.// During arithmetic operations, denormalized values may occur but are// always normalized before returning the final result. The normalized// representation of 0 is the empty or nil slice (length = 0).type nat []Wordvar ( natOne = nat{1} natTwo = nat{2} natFive = nat{5} natTen = nat{10})func ( nat) () string {return"0x" + string(.itoa(false, 16))}func ( nat) () nat { := len()for > 0 && [-1] == 0 { -- }return [0:]}func ( nat) ( int) nat {if <= cap() {return [:] // reuse z }if == 1 {// Most nats start small and stay that way; don't over-allocate.returnmake(nat, 1) }// Choosing a good value for e has significant performance impact // because it increases the chance that a value can be reused.const = 4// extra capacityreturnmake(nat, , +)}func ( nat) ( Word) nat {if == 0 {return [:0] } = .make(1) [0] = return}func ( nat) ( uint64) nat {// single-word valueif := Word(); uint64() == {return .setWord() }// 2-word value = .make(2) [1] = Word( >> 32) [0] = Word()return}func ( nat) ( nat) nat { = .make(len())copy(, )return}func ( nat) (, nat) nat { := len() := len()switch {case < :return .(, )case == 0:// n == 0 because m >= n; result is 0return [:0]case == 0:// result is xreturn .set() }// m > 0 = .make( + 1) := addVV([:], [:], [:])if > { = addVW([:], [:], ) } [] = return .norm()}func ( nat) (, nat) nat { := len() := len()switch {case < :panic("underflow")case == 0:// n == 0 because m >= n; result is 0return [:0]case == 0:// result is xreturn .set() }// m > 0 = .make() := subVV([:], [:], [:])if > { = subVW([:], [:], ) }if != 0 {panic("underflow") }return .norm()}func ( nat) ( nat) ( int) { := len() := len()if != || == 0 {switch {case < : = -1case > : = 1 }return } := - 1for > 0 && [] == [] { -- }switch {case [] < []: = -1case [] > []: = 1 }return}// montgomery computes z mod m = x*y*2**(-n*_W) mod m,// assuming k = -1/m mod 2**_W.// z is used for storing the result which is returned;// z must not alias x, y or m.// See Gueron, "Efficient Software Implementations of Modular Exponentiation".// https://eprint.iacr.org/2011/239.pdf// In the terminology of that paper, this is an "Almost Montgomery Multiplication":// x and y are required to satisfy 0 <= z < 2**(n*_W) and then the result// z is guaranteed to satisfy 0 <= z < 2**(n*_W), but it may not be < m.func ( nat) (, , nat, Word, int) nat {// This code assumes x, y, m are all the same length, n. // (required by addMulVVW and the for loop). // It also assumes that x, y are already reduced mod m, // or else the result will not be properly reduced.iflen() != || len() != || len() != {panic("math/big: mismatched montgomery number lengths") } = .make( * 2)clear()varWordfor := 0; < ; ++ { := [] := addMulVVWW([:+], [:+], , , 0) := [] * := addMulVVWW([:+], [:+], , , 0) := + := + [+] = if < || < { = 1 } else { = 0 } }if != 0 {subVV([:], [:], ) } else {copy([:], [:]) }return [:]}// alias reports whether x and y share the same base array.//// Note: alias assumes that the capacity of underlying arrays// is never changed for nat values; i.e. that there are// no 3-operand slice expressions in this code (or worse,// reflect-based operations to the same effect).func alias(, nat) bool {returncap() > 0 && cap() > 0 && &[0:cap()][cap()-1] == &[0:cap()][cap()-1]}// addTo implements z += x; z must be long enough.// (we don't use nat.add because we need z to stay the same// slice, and we don't need to normalize z after each addition)func addTo(, nat) {if := len(); > 0 {if := addVV([:], [:], [:]); != 0 {if < len() {addVW([:], [:], ) } } }}// mulRange computes the product of all the unsigned integers in the// range [a, b] inclusively. If a > b (empty range), the result is 1.// The caller may pass stk == nil to request that mulRange obtain and release one itself.func ( nat) ( *stack, , uint64) nat {switch {case == 0:// cut long ranges short (optimization)return .setUint64(0)case > :return .setUint64(1)case == :return .setUint64()case +1 == :return .mul(, nat(nil).setUint64(), nat(nil).setUint64()) }if == nil { = getStack()defer .free() } := + (-)/2// avoid overflowreturn .mul(, nat(nil).(, , ), nat(nil).(, +1, ))}// A stack provides temporary storage for complex calculations// such as multiplication and division.// The stack is a simple slice of words, extended as needed// to hold all the temporary storage for a calculation.// In general, if a function takes a *stack, it expects a non-nil *stack.// However, certain functions may allow passing a nil *stack instead,// so that they can handle trivial stack-free cases without forcing the// caller to obtain and free a stack that will be unused. These functions// document that they accept a nil *stack in their doc comments.type stack struct { w []Word}var stackPool sync.Pool// getStack returns a temporary stack.// The caller must call [stack.free] to give up use of the stack when finished.func getStack() *stack { , := stackPool.Get().(*stack)if == nil { = new(stack) }return}// free returns the stack for use by another calculation.func ( *stack) () { .w = .w[:0]stackPool.Put()}// save returns the current stack pointer.// A future call to restore with the same value// frees any temporaries allocated on the stack after the call to save.func ( *stack) () int {returnlen(.w)}// restore restores the stack pointer to n.// It is almost always invoked as//// defer stk.restore(stk.save())//// which makes sure to pop any temporaries allocated in the current function// from the stack before returning.func ( *stack) ( int) { .w = .w[:]}// nat returns a nat of n words, allocated on the stack.func ( *stack) ( int) nat { := ( + 3) &^ 3// round up to multiple of 4 := len(.w) .w = slices.Grow(.w, ) .w = .w[:+] := .w[ : + : +]if > 0 { [0] = 0xfedcb }return}// bitLen returns the length of x in bits.// Unlike most methods, it works even if x is not normalized.func ( nat) () int {// This function is used in cryptographic operations. It must not leak // anything but the Int's sign and bit size through side-channels. Any // changes must be reviewed by a security expert.if := len() - 1; >= 0 {// bits.Len uses a lookup table for the low-order bits on some // architectures. Neutralize any input-dependent behavior by setting all // bits after the first one bit. := uint([]) |= >> 1 |= >> 2 |= >> 4 |= >> 8 |= >> 16 |= >> 16 >> 16// ">> 32" doesn't compile on 32-bit architecturesreturn *_W + bits.Len() }return0}// trailingZeroBits returns the number of consecutive least significant zero// bits of x.func ( nat) () uint {iflen() == 0 {return0 }varuintfor [] == 0 { ++ }// x[i] != 0return *_W + uint(bits.TrailingZeros(uint([])))}// isPow2 returns i, true when x == 2**i and 0, false otherwise.func ( nat) () (uint, bool) {varuintfor [] == 0 { ++ }if == uint(len())-1 && []&([]-1) == 0 {return *_W + uint(bits.TrailingZeros(uint([]))), true }return0, false}func same(, nat) bool {returnlen() == len() && len() > 0 && &[0] == &[0]}// z = x << sfunc ( nat) ( nat, uint) nat {if == 0 {ifsame(, ) {return }if !alias(, ) {return .set() } } := len()if == 0 {return [:0] }// m > 0 := + int(/_W) = .make( + 1)if %= _W; == 0 {copy([-:], ) [] = 0 } else { [] = lshVU([-:], , ) }clear([0 : -])return .norm()}// z = x >> sfunc ( nat) ( nat, uint) nat {if == 0 {ifsame(, ) {return }if !alias(, ) {return .set() } } := len() := - int(/_W)if <= 0 {return [:0] }// n > 0 = .make()if %= _W; == 0 {copy(, [-:]) } else {rshVU(, [-:], ) }return .norm()}func ( nat) ( nat, uint, uint) nat { := int( / _W) := Word(1) << ( % _W) := len()switch {case0: = .make()copy(, )if >= {// no need to growreturn } [] &^= return .norm()case1:if >= { = .make( + 1)clear([:]) } else { = .make() }copy(, ) [] |= // no need to normalizereturn }panic("set bit is not 0 or 1")}// bit returns the value of the i'th bit, with lsb == bit 0.func ( nat) ( uint) uint { := / _Wif >= uint(len()) {return0 }// 0 <= j < len(x)returnuint([] >> ( % _W) & 1)}// sticky returns 1 if there's a 1 bit within the// i least significant bits, otherwise it returns 0.func ( nat) ( uint) uint { := / _Wif >= uint(len()) {iflen() == 0 {return0 }return1 }// 0 <= j < len(x)for , := range [:] {if != 0 {return1 } }if []<<(_W-%_W) != 0 {return1 }return0}func ( nat) (, nat) nat { := len() := len()if > { = }// m <= n = .make()for := 0; < ; ++ { [] = [] & [] }return .norm()}// trunc returns z = x mod 2ⁿ.func ( nat) ( nat, uint) nat { := ( + _W - 1) / _Wifuint(len()) < {return .set() } = .make(int())copy(, )if %_W != 0 { [len()-1] &= 1<<(%_W) - 1 }return .norm()}func ( nat) (, nat) nat { := len() := len()if > { = }// m >= n = .make()for := 0; < ; ++ { [] = [] &^ [] }copy([:], [:])return .norm()}func ( nat) (, nat) nat { := len() := len() := if < { , = , = }// m >= n = .make()for := 0; < ; ++ { [] = [] | [] }copy([:], [:])return .norm()}func ( nat) (, nat) nat { := len() := len() := if < { , = , = }// m >= n = .make()for := 0; < ; ++ { [] = [] ^ [] }copy([:], [:])return .norm()}// random creates a random integer in [0..limit), using the space in z if// possible. n is the bit length of limit.func ( nat) ( *rand.Rand, nat, int) nat {ifalias(, ) { = nil// z is an alias for limit - cannot reuse } = .make(len()) := uint( % _W)if == 0 { = _W } := Word((1 << ) - 1)for {switch_W {case32:for := range { [] = Word(.Uint32()) }case64:for := range { [] = Word(.Uint32()) | Word(.Uint32())<<32 }default:panic("unknown word size") } [len()-1] &= if .cmp() < 0 {break } }return .norm()}// If m != 0 (i.e., len(m) != 0), expNN sets z to x**y mod m;// otherwise it sets z to x**y. The result is the value of z.// The caller may pass stk == nil to request that expNN obtain and release one itself.func ( nat) ( *stack, , , nat, bool) nat {ifalias(, ) || alias(, ) {// We cannot allow in-place modification of x or y. = nil }// x**y mod 1 == 0iflen() == 1 && [0] == 1 {return .setWord(0) }// m == 0 || m > 1// x**0 == 1iflen() == 0 {return .setWord(1) }// y > 0// 0**y = 0iflen() == 0 {return .setWord(0) }// x > 0// 1**y = 1iflen() == 1 && [0] == 1 {return .setWord(1) }// x > 1// x**1 == xiflen() == 1 && [0] == 1 && len() == 0 {return .set() }if == nil { = getStack()defer .free() }iflen() == 1 && [0] == 1 { // len(m) > 0return .rem(, , ) }// y > 1iflen() != 0 {// We likely end up being as long as the modulus. = .make(len())// If the exponent is large, we use the Montgomery method for odd values, // and a 4-bit, windowed exponentiation for powers of two, // and a CRT-decomposed Montgomery method for the remaining values // (even values times non-trivial odd values, which decompose into one // instance of each of the first two cases).iflen() > 1 && ! {if [0]&1 == 1 {return .expNNMontgomery(, , , ) }if , := .isPow2(); {return .expNNWindowed(, , , ) }return .expNNMontgomeryEven(, , , ) } } = .set() := [len()-1] // v > 0 because y is normalized and y > 0 := nlz() + 1 <<= varnatconst = 1 << (_W - 1)// We walk through the bits of the exponent one by one. Each time we // see a bit, we square, thus doubling the power. If the bit is a one, // we also multiply by x, thus adding one to the power. := _W - int()// zz and r are used to avoid allocating in mul and div as // otherwise the arguments would alias.var , natfor := 0; < ; ++ { = .sqr(, ) , = , if & != 0 { = .mul(, , ) , = , }iflen() != 0 { , = .div(, , , ) , , , = , , , } <<= 1 }for := len() - 2; >= 0; -- { = []for := 0; < _W; ++ { = .sqr(, ) , = , if & != 0 { = .mul(, , ) , = , }iflen() != 0 { , = .div(, , , ) , , , = , , , } <<= 1 } }return .norm()}// expNNMontgomeryEven calculates x**y mod m where m = m1 × m2 for m1 = 2ⁿ and m2 odd.// It uses two recursive calls to expNN for x**y mod m1 and x**y mod m2// and then uses the Chinese Remainder Theorem to combine the results.// The recursive call using m1 will use expNNWindowed,// while the recursive call using m2 will use expNNMontgomery.// For more details, see Ç. K. Koç, “Montgomery Reduction with Even Modulus”,// IEE Proceedings: Computers and Digital Techniques, 141(5) 314-316, September 1994.// http://www.people.vcu.edu/~jwang3/CMSC691/j34monex.pdffunc ( nat) ( *stack, , , nat) nat {// Split m = m₁ × m₂ where m₁ = 2ⁿ := .trailingZeroBits() := nat(nil).lsh(natOne, ) := nat(nil).rsh(, )// We want z = x**y mod m. // z₁ = x**y mod m1 = (x**y mod m) mod m1 = z mod m1 // z₂ = x**y mod m2 = (x**y mod m) mod m2 = z mod m2 // (We are using the math/big convention for names here, // where the computation is z = x**y mod m, so its parts are z1 and z2. // The paper is computing x = a**e mod n; it refers to these as x2 and z1.) := nat(nil).expNN(, , , , false) := nat(nil).expNN(, , , , false)// Reconstruct z from z₁, z₂ using CRT, using algorithm from paper, // which uses only a single modInverse (and an easy one at that). // p = (z₁ - z₂) × m₂⁻¹ (mod m₁) // z = z₂ + p × m₂ // The final addition is in range because: // z = z₂ + p × m₂ // ≤ z₂ + (m₁-1) × m₂ // < m₂ + (m₁-1) × m₂ // = m₁ × m₂ // = m. = .set()// Compute (z₁ - z₂) mod m1 [m1 == 2**n] into z1. = .subMod2N(, , )// Reuse z2 for p = (z₁ - z₂) [in z1] * m2⁻¹ (mod m₁ [= 2ⁿ]). := nat(nil).modInverse(, ) = .mul(, , ) = .trunc(, )// Reuse z1 for p * m2. = .add(, .mul(, , ))return}// expNNWindowed calculates x**y mod m using a fixed, 4-bit window,// where m = 2**logM.func ( nat) ( *stack, , nat, uint) nat {iflen() <= 1 {panic("big: misuse of expNNWindowed") }if [0]&1 == 0 {// len(y) > 1, so y > logM. // x is even, so x**y is a multiple of 2**y which is a multiple of 2**logM.return .setWord(0) }if == 1 {return .setWord(1) }// zz is used to avoid allocating in mul as otherwise // the arguments would alias.defer .restore(.save()) := int(( + _W - 1) / _W) := .nat()const = 4// powers[i] contains x^i.var [1 << ]natfor := range { [] = .nat() } [0] = [0].set(natOne) [1] = [1].trunc(, )for := 2; < 1<<; += 2 { , , := &[/2], &[], &[+1] * = .sqr(, *) * = .trunc(*, ) * = .mul(, *, ) * = .trunc(*, ) }// Because phi(2**logM) = 2**(logM-1), x**(2**(logM-1)) = 1, // so we can compute x**(y mod 2**(logM-1)) instead of x**y. // That is, we can throw away all but the bottom logM-1 bits of y. // Instead of allocating a new y, we start reading y at the right word // and truncate it appropriately at the start of the loop. := len() - 1 := int(( - 2) / _W) // -2 because the top word of N bits is the (N-1)/W'th word. := ^Word(0)if := ( - 1) & (_W - 1); != 0 { = (1 << ) - 1 }if > { = } := false = .setWord(1)for ; >= 0; -- { := []if == { &= }for := 0; < _W; += {if {// Account for use of 4 bits in previous iteration. // Unrolled loop for significant performance // gain. Use go test -bench=".*" in crypto/rsa // to check performance before making changes. = .sqr(, ) , = , = .trunc(, ) = .sqr(, ) , = , = .trunc(, ) = .sqr(, ) , = , = .trunc(, ) = .sqr(, ) , = , = .trunc(, ) } = .mul(, , [>>(_W-)]) , = , = .trunc(, ) <<= = true } }return .norm()}// expNNMontgomery calculates x**y mod m using a fixed, 4-bit window.// Uses Montgomery representation.func ( nat) ( *stack, , , nat) nat { := len()// We want the lengths of x and m to be equal. // It is OK if x >= m as long as len(x) == len(m).iflen() > { _, = nat(nil).div(, nil, , )// Note: now len(x) <= numWords, not guaranteed ==. }iflen() < { := make(nat, )copy(, ) = }// Ideally the precomputations would be performed outside, and reused // k0 = -m**-1 mod 2**_W. Algorithm from: Dumas, J.G. "On Newton–Raphson // Iteration for Multiplicative Inverses Modulo Prime Powers". := 2 - [0] := [0] - 1for := 1; < _W; <<= 1 { *= *= ( + 1) } = -// RR = 2**(2*_W*len(m)) mod m := nat(nil).setWord(1) := nat(nil).lsh(, uint(2**_W)) _, = nat(nil).div(, , , )iflen() < { = .make()copy(, ) = }// one = 1, with equal length to that of m := make(nat, ) [0] = 1const = 4// powers[i] contains x^ivar [1 << ]nat [0] = [0].montgomery(, , , , ) [1] = [1].montgomery(, , , , )for := 2; < 1<<; ++ { [] = [].montgomery([-1], [1], , , ) }// initialize z = 1 (Montgomery 1) = .make()copy(, [0]) = .make()// same windowed exponent, but with Montgomery multiplicationsfor := len() - 1; >= 0; -- { := []for := 0; < _W; += {if != len()-1 || != 0 { = .montgomery(, , , , ) = .montgomery(, , , , ) = .montgomery(, , , , ) = .montgomery(, , , , ) } = .montgomery(, [>>(_W-)], , , ) , = , <<= } }// convert to regular number = .montgomery(, , , , )// One last reduction, just in case. // See golang.org/issue/13907.if .cmp() >= 0 {// Common case is m has high bit set; in that case, // since zz is the same length as m, there can be just // one multiple of m to remove. Just subtract. // We think that the subtract should be sufficient in general, // so do that unconditionally, but double-check, // in case our beliefs are wrong. // The div is not expected to be reached. = .sub(, )if .cmp() >= 0 { _, = nat(nil).div(, nil, , ) } }return .norm()}// bytes writes the value of z into buf using big-endian encoding.// The value of z is encoded in the slice buf[i:]. If the value of z// cannot be represented in buf, bytes panics. The number i of unused// bytes at the beginning of buf is returned as result.func ( nat) ( []byte) ( int) {// This function is used in cryptographic operations. It must not leak // anything but the Int's sign and bit size through side-channels. Any // changes must be reviewed by a security expert. = len()for , := range {for := 0; < _S; ++ { --if >= 0 { [] = byte() } elseifbyte() != 0 {panic("math/big: buffer too small to fit value") } >>= 8 } }if < 0 { = 0 }for < len() && [] == 0 { ++ }return}// bigEndianWord returns the contents of buf interpreted as a big-endian encoded Word value.func bigEndianWord( []byte) Word {if_W == 64 {returnWord(byteorder.BEUint64()) }returnWord(byteorder.BEUint32())}// setBytes interprets buf as the bytes of a big-endian unsigned// integer, sets z to that value, and returns z.func ( nat) ( []byte) nat { = .make((len() + _S - 1) / _S) := len()for := 0; >= _S; ++ { [] = bigEndianWord([-_S : ]) -= _S }if > 0 {varWordfor := uint(0); > 0; += 8 { |= Word([-1]) << -- } [len()-1] = }return .norm()}// sqrt sets z = ⌊√x⌋// The caller may pass stk == nil to request that sqrt obtain and release one itself.func ( nat) ( *stack, nat) nat {if .cmp(natOne) <= 0 {return .set() }ifalias(, ) { = nil }if == nil { = getStack()defer .free() }// Start with value known to be too large and repeat "z = ⌊(z + ⌊x/z⌋)/2⌋" until it stops getting smaller. // See Brent and Zimmermann, Modern Computer Arithmetic, Algorithm 1.13 (SqrtInt). // https://members.loria.fr/PZimmermann/mca/pub226.html // If x is one less than a perfect square, the sequence oscillates between the correct z and z+1; // otherwise it converges to the correct z and stays there.var , nat = = .setUint64(1) = .lsh(, uint(.bitLen()+1)/2) // must be ≥ √xfor := 0; ; ++ { , _ = .div(, nil, , ) = .add(, ) = .rsh(, 1)if .cmp() >= 0 {// z1 is answer. // Figure out whether z1 or z2 is currently aliased to z by looking at loop count.if &1 == 0 {return }return .set() } , = , }}// subMod2N returns z = (x - y) mod 2ⁿ.func ( nat) (, nat, uint) nat {ifuint(.bitLen()) > {ifalias(, ) {// ok to overwrite x in place = .trunc(, ) } else { = nat(nil).trunc(, ) } }ifuint(.bitLen()) > {ifalias(, ) {// ok to overwrite y in place = .trunc(, ) } else { = nat(nil).trunc(, ) } }if .cmp() >= 0 {return .sub(, ) }// x - y < 0; x - y mod 2ⁿ = x - y + 2ⁿ = 2ⁿ - (y - x) = 1 + 2ⁿ-1 - (y - x) = 1 + ^(y - x). = .sub(, )foruint(len())*_W < { = append(, 0) }for := range { [] = ^[] } = .trunc(, )return .add(, natOne)}
The pages are generated with Goldsv0.7.7-preview. (GOOS=linux GOARCH=amd64)
Golds is a Go 101 project developed by Tapir Liu.
PR and bug reports are welcome and can be submitted to the issue list.
Please follow @zigo_101 (reachable from the left QR code) to get the latest news of Golds.