I've been trying to implement Pollard's Rho algorithm (https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm) in Nim, using the BigNum module (which in turn uses GMP).
Up to the 6th Fermat number (18446744073709551617) it works fine. And it works fine for F(8) as well.
However, it keeps crashing for F(7). Basically, it consumes close to 100GB of disk space and then the process gets killed.
I know there is an issue with F(7) factorization, given that it doesn't have small factors.
The question is: do you see anything wrong with my code or anything that could be written more efficiently?
proc pollardRho*(n: Int): Int =
var x = newInt(2)
var y = newInt(2)
var d = newInt(1)
var z = newInt(1)
var count = 0
while true:
x = (x*x + 1) mod n
y = (((y*y+1) mod n)*((y*y+1) mod n) + 1) mod n
var t = abs(x-y)
t = t mod n
z = z * t
inc(count)
if count==100:
d = gcd(z,n)
if d!=1:
break
z = newInt(1)
count = 0
if cmp(d,n)==0:
return newInt(0)
else:
return d
proc primeFactors*(num: Int): seq[Int] =
result = @[]
var n = num
if n.probablyPrime(10)!=0:
result.add(n)
let factor1 = pollardRho(num)
if factor1==0:
return @[]
if factor1.probablyPrime(10)==0:
return @[]
let factor2 = n div factor1
if factor2.probablyPrime(10)==0:
return @[]
result.add(factor1)
result.add(factor2)
The issue is probably not in you algorithm but in memory management.
You are using out-of-place routines like x = (x*x + 1) mod n, that means an extra allocation at each loop. F(7) is 2^128 + 1 so it does take some space and for F(7) it seems like the GC cannot reclaim that memory. And you have a dozens of out-of-place operand per-loop.
Instead you should refactor your program to allocate the minimum possible by using the in-place operators that would re-use one of your bigint buffer instead of allocating a new one.
That will make it uglier obviously.
With the new Nim semantics of destructors, especially the sink parameters, we could recreate a wrapper that would detect when a parameter is actually not used later and its buffer could be reused without allocating. I leave such an exercise to the reader ;).
You're a star! So simple, but I was missing it. I did it with the in-place routines and it worked beautifully!
No memory issues, and SUPER fast (relatively).
Here's my updated code in case somebody needs it:
proc pollardG*(n: var Int, m: Int) {.inline.} =
discard mul(n,n,n)
discard add(n,n,1)
discard `mod`(n,n,m)
proc pollardRho*(n: Int): Int {.noSideEffect.} =
var x = newInt(2)
var y = newInt(2)
var d = newInt(1)
var z = newInt(1)
var count = 0
var t = newInt(0)
while true:
pollardG(x,n) # x = (x*x + 1) mod n
pollardG(y,n)
pollardG(y,n) # y = (((y*y+1) mod n)*((y*y+1) mod n) + 1) mod n
discard abs(t,sub(t,x,y)) # t = abs(x-y)
#var t = abs(x-y)
discard `mod`(t,t,n) # t = t mod n
discard mul(z,z,t) # z = z * t
inc(count)
if count==100:
discard gcd(d,z,n) # d = gcd(z,n)
if cmp(d,1)!=0: # d!=1:
break
discard set(z,1) # z = newInt(1)
count = 0
if cmp(d,n)==0:
return newInt(0)
else:
return d
proc primeFactors*(num: Int): seq[Int] {.noSideEffect.} =
result = @[]
var n = num
if n.probablyPrime(10)!=0:
result.add(n)
let factor1 = pollardRho(num)
if factor1==0:
return @[]
if factor1.probablyPrime(10)==0:
return @[]
let factor2 = n div factor1
if factor2.probablyPrime(10)==0:
return @[]
result.add(factor1)
result.add(factor2)
Thanks again!
And... without the one million ugly comments:
proc pollardG*(n: var Int, m: Int) {.inline.} =
discard mul(n,n,n)
discard add(n,n,1)
discard `mod`(n,n,m)
proc pollardRho*(n: Int): Int {.noSideEffect.} =
var x = newInt(2)
var y = newInt(2)
var d = newInt(1)
var z = newInt(1)
var count = 0
var t = newInt(0)
while true:
pollardG(x,n)
pollardG(y,n)
pollardG(y,n)
discard abs(t,sub(t,x,y))
discard `mod`(t,t,n)
discard mul(z,z,t)
inc(count)
if count==100:
discard gcd(d,z,n)
if cmp(d,1)!=0:
break
discard set(z,1)
count = 0
if cmp(d,n)==0:
return newInt(0)
else:
return d
proc primeFactors*(num: Int): seq[Int] {.noSideEffect.} =
result = @[]
var n = num
if n.probablyPrime(10)!=0:
result.add(n)
let factor1 = pollardRho(num)
if factor1==0:
return @[]
if factor1.probablyPrime(10)==0:
return @[]
let factor2 = n div factor1
if factor2.probablyPrime(10)==0:
return @[]
result.add(factor1)
result.add(factor2)