diff options
Diffstat (limited to 'libgo/go/math/big/nat.go')
-rw-r--r-- | libgo/go/math/big/nat.go | 176 |
1 files changed, 65 insertions, 111 deletions
diff --git a/libgo/go/math/big/nat.go b/libgo/go/math/big/nat.go index 2e65d2a7ef7..9b1a626c4cf 100644 --- a/libgo/go/math/big/nat.go +++ b/libgo/go/math/big/nat.go @@ -542,16 +542,21 @@ func (z nat) div(z2, u, v nat) (q, r nat) { return } -// getNat returns a nat of len n. The contents may not be zero. -func getNat(n int) nat { - var z nat +// getNat returns a *nat of len n. The contents may not be zero. +// The pool holds *nat to avoid allocation when converting to interface{}. +func getNat(n int) *nat { + var z *nat if v := natPool.Get(); v != nil { - z = v.(nat) + z = v.(*nat) } - return z.make(n) + if z == nil { + z = new(nat) + } + *z = z.make(n) + return z } -func putNat(x nat) { +func putNat(x *nat) { natPool.Put(x) } @@ -575,7 +580,8 @@ func (z nat) divLarge(u, uIn, v nat) (q, r nat) { } q = z.make(m + 1) - qhatv := getNat(n + 1) + qhatvp := getNat(n + 1) + qhatv := *qhatvp if alias(u, uIn) || alias(u, v) { u = nil // u is an alias for uIn or v - cannot reuse } @@ -583,36 +589,40 @@ func (z nat) divLarge(u, uIn, v nat) (q, r nat) { u.clear() // TODO(gri) no need to clear if we allocated a new u // D1. - var v1 nat + var v1p *nat shift := nlz(v[n-1]) if shift > 0 { // do not modify v, it may be used by another goroutine simultaneously - v1 = getNat(n) + v1p = getNat(n) + v1 := *v1p shlVU(v1, v, shift) v = v1 } u[len(uIn)] = shlVU(u[0:len(uIn)], uIn, shift) // D2. + vn1 := v[n-1] for j := m; j >= 0; j-- { // D3. qhat := Word(_M) - if u[j+n] != v[n-1] { + if ujn := u[j+n]; ujn != vn1 { var rhat Word - qhat, rhat = divWW(u[j+n], u[j+n-1], v[n-1]) + qhat, rhat = divWW(ujn, u[j+n-1], vn1) // x1 | x2 = q̂v_{n-2} - x1, x2 := mulWW(qhat, v[n-2]) + vn2 := v[n-2] + x1, x2 := mulWW(qhat, vn2) // test if q̂v_{n-2} > br̂ + u_{j+n-2} - for greaterThan(x1, x2, rhat, u[j+n-2]) { + ujn2 := u[j+n-2] + for greaterThan(x1, x2, rhat, ujn2) { qhat-- prevRhat := rhat - rhat += v[n-1] + rhat += vn1 // v[n-1] >= 0, so this tests for overflow. if rhat < prevRhat { break } - x1, x2 = mulWW(qhat, v[n-2]) + x1, x2 = mulWW(qhat, vn2) } } @@ -628,10 +638,10 @@ func (z nat) divLarge(u, uIn, v nat) (q, r nat) { q[j] = qhat } - if v1 != nil { - putNat(v1) + if v1p != nil { + putNat(v1p) } - putNat(qhatv) + putNat(qhatvp) q = q.norm() shrVU(u, u, shift) @@ -650,14 +660,14 @@ func (x nat) bitLen() int { const deBruijn32 = 0x077CB531 -var deBruijn32Lookup = []byte{ +var deBruijn32Lookup = [...]byte{ 0, 1, 28, 2, 29, 14, 24, 3, 30, 22, 20, 15, 25, 17, 4, 8, 31, 27, 13, 23, 21, 19, 16, 7, 26, 12, 18, 6, 11, 5, 10, 9, } const deBruijn64 = 0x03f79d71b4ca8b09 -var deBruijn64Lookup = []byte{ +var deBruijn64Lookup = [...]byte{ 0, 1, 56, 2, 57, 49, 28, 3, 61, 58, 42, 50, 38, 29, 17, 4, 62, 47, 59, 36, 45, 43, 51, 22, 53, 39, 33, 30, 24, 18, 12, 5, 63, 55, 48, 27, 60, 41, 37, 16, 46, 35, 44, 21, 52, 32, 23, 11, @@ -950,7 +960,7 @@ func (z nat) expNN(x, y, m nat) nat { // (x^2...x^15) but then reduces the number of multiply-reduces by a // third. Even for a 32-bit exponent, this reduces the number of // operations. Uses Montgomery method for odd moduli. - if len(x) > 1 && len(y) > 1 && len(m) > 0 { + if x.cmp(natOne) > 0 && len(y) > 1 && len(m) > 0 { if m[0]&1 == 1 { return z.expNNMontgomery(x, y, m) } @@ -1169,96 +1179,6 @@ func (z nat) expNNMontgomery(x, y, m nat) nat { return zz.norm() } -// probablyPrime performs n Miller-Rabin tests to check whether x is prime. -// If x is prime, it returns true. -// If x is not prime, it returns false with probability at least 1 - ¼ⁿ. -// -// It is not suitable for judging primes that an adversary may have crafted -// to fool this test. -func (n nat) probablyPrime(reps int) bool { - if len(n) == 0 { - return false - } - - if len(n) == 1 { - if n[0] < 2 { - return false - } - - if n[0]%2 == 0 { - return n[0] == 2 - } - - // We have to exclude these cases because we reject all - // multiples of these numbers below. - switch n[0] { - case 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53: - return true - } - } - - if n[0]&1 == 0 { - return false // n is even - } - - const primesProduct32 = 0xC0CFD797 // Π {p ∈ primes, 2 < p <= 29} - const primesProduct64 = 0xE221F97C30E94E1D // Π {p ∈ primes, 2 < p <= 53} - - var r Word - switch _W { - case 32: - r = n.modW(primesProduct32) - case 64: - r = n.modW(primesProduct64 & _M) - default: - panic("Unknown word size") - } - - if r%3 == 0 || r%5 == 0 || r%7 == 0 || r%11 == 0 || - r%13 == 0 || r%17 == 0 || r%19 == 0 || r%23 == 0 || r%29 == 0 { - return false - } - - if _W == 64 && (r%31 == 0 || r%37 == 0 || r%41 == 0 || - r%43 == 0 || r%47 == 0 || r%53 == 0) { - return false - } - - nm1 := nat(nil).sub(n, natOne) - // determine q, k such that nm1 = q << k - k := nm1.trailingZeroBits() - q := nat(nil).shr(nm1, k) - - nm3 := nat(nil).sub(nm1, natTwo) - rand := rand.New(rand.NewSource(int64(n[0]))) - - var x, y, quotient nat - nm3Len := nm3.bitLen() - -NextRandom: - for i := 0; i < reps; i++ { - x = x.random(rand, nm3, nm3Len) - x = x.add(x, natTwo) - y = y.expNN(x, q, n) - if y.cmp(natOne) == 0 || y.cmp(nm1) == 0 { - continue - } - for j := uint(1); j < k; j++ { - y = y.mul(y, y) - quotient, y = quotient.div(y, y, n) - if y.cmp(nm1) == 0 { - continue NextRandom - } - if y.cmp(natOne) == 0 { - return false - } - } - return false - } - - return true -} - // bytes writes the value of z into buf using big-endian encoding. // len(buf) must be >= len(z)*_S. The value of z is encoded in the // slice buf[i:]. The number i of unused bytes at the beginning of @@ -1303,3 +1223,37 @@ func (z nat) setBytes(buf []byte) nat { return z.norm() } + +// sqrt sets z = ⌊√x⌋ +func (z nat) sqrt(x nat) nat { + if x.cmp(natOne) <= 0 { + return z.set(x) + } + if alias(z, x) { + z = nil + } + + // 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 z1, z2 nat + z1 = z + z1 = z1.setUint64(1) + z1 = z1.shl(z1, uint(x.bitLen()/2+1)) // must be ≥ √x + for n := 0; ; n++ { + z2, _ = z2.div(nil, x, z1) + z2 = z2.add(z2, z1) + z2 = z2.shr(z2, 1) + if z2.cmp(z1) >= 0 { + // z1 is answer. + // Figure out whether z1 or z2 is currently aliased to z by looking at loop count. + if n&1 == 0 { + return z1 + } + return z.set(z1) + } + z1, z2 = z2, z1 + } +} |