130 lines
2.8 KiB
Go
130 lines
2.8 KiB
Go
package main
|
||
|
||
import (
|
||
"fmt"
|
||
big "github.com/ncw/gmp"
|
||
"log"
|
||
)
|
||
|
||
var (
|
||
zero = big.NewInt(0)
|
||
one = big.NewInt(1)
|
||
two = big.NewInt(2)
|
||
three = big.NewInt(3)
|
||
four = big.NewInt(4)
|
||
five = big.NewInt(5)
|
||
six = big.NewInt(6)
|
||
ten = big.NewInt(10)
|
||
k100 = big.NewInt(100000)
|
||
)
|
||
|
||
func pollardRho(n, c *big.Int) *big.Int {
|
||
g := func(x, y *big.Int) *big.Int {
|
||
x2 := new(big.Int)
|
||
x2.Mul(x, x)
|
||
x2.Add(x2, c)
|
||
return x2.Mod(x2, y)
|
||
}
|
||
x, y, z := big.NewInt(2), big.NewInt(2), big.NewInt(1)
|
||
d := new(big.Int)
|
||
count := 0
|
||
for {
|
||
x = g(x, n)
|
||
y = g(g(y, n), n)
|
||
d.Sub(x, y)
|
||
d.Abs(d)
|
||
d.Mod(d, n)
|
||
z.Mul(z, d)
|
||
count++
|
||
if count == 100 {
|
||
d.GCD(nil, nil, z, n)
|
||
if d.Cmp(one) != 0 {
|
||
break
|
||
}
|
||
z.Set(one)
|
||
count = 0
|
||
}
|
||
}
|
||
if d.Cmp(n) == 0 {
|
||
return zero
|
||
}
|
||
return d
|
||
}
|
||
|
||
func smallestPrimeFactorWheel(n, max *big.Int) *big.Int {
|
||
if n.ProbablyPrime(15) {
|
||
return n
|
||
}
|
||
z := new(big.Int)
|
||
if z.Rem(n, two).Cmp(zero) == 0 {
|
||
return two
|
||
}
|
||
if z.Rem(n, three).Cmp(zero) == 0 {
|
||
return three
|
||
}
|
||
if z.Rem(n, five).Cmp(zero) == 0 {
|
||
return five
|
||
}
|
||
k := big.NewInt(7)
|
||
i := 0
|
||
inc := []*big.Int{four, two, four, two, four, six, two, six}
|
||
for z.Mul(k, k).Cmp(n) <= 0 {
|
||
if z.Rem(n, k).Cmp(zero) == 0 {
|
||
return k
|
||
}
|
||
k.Add(k, inc[i])
|
||
if k.Cmp(max) > 0 {
|
||
break
|
||
}
|
||
i = (i + 1) % 8
|
||
}
|
||
return nil
|
||
}
|
||
|
||
func smallestPrimeFactor(n *big.Int) *big.Int {
|
||
s := smallestPrimeFactorWheel(n, k100)
|
||
if s != nil {
|
||
return s
|
||
}
|
||
c := big.NewInt(1)
|
||
s = new(big.Int).Set(n)
|
||
for {
|
||
d := pollardRho(n, c)
|
||
if d.Cmp(zero) == 0 {
|
||
if c.Cmp(ten) == 0 {
|
||
log.Fatal("Pollard Rho doesn't appear to be working.")
|
||
}
|
||
c.Add(c, one)
|
||
} else {
|
||
// get the smallest prime factor of 'd'
|
||
factor := smallestPrimeFactorWheel(d, d)
|
||
// check whether n/d has a smaller prime factor
|
||
s = smallestPrimeFactorWheel(n.Quo(n, d), factor)
|
||
if s != nil {
|
||
if s.Cmp(factor) < 0 {
|
||
return s
|
||
} else {
|
||
return factor
|
||
}
|
||
} else {
|
||
return factor
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
func main() {
|
||
k := 19
|
||
fmt.Println("First", k, "terms of the Euclid–Mullin sequence:")
|
||
fmt.Println(2)
|
||
prod := big.NewInt(2)
|
||
z := new(big.Int)
|
||
count := 1
|
||
for count < k {
|
||
z.Add(prod, one)
|
||
t := smallestPrimeFactor(z)
|
||
fmt.Println(t)
|
||
prod.Mul(prod, t)
|
||
count++
|
||
}
|
||
}
|