summaryrefslogtreecommitdiff
path: root/libgo/go/math/big/nat.go
diff options
context:
space:
mode:
Diffstat (limited to 'libgo/go/math/big/nat.go')
-rw-r--r--libgo/go/math/big/nat.go176
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
+ }
+}