How can I speedup this Julia code?

Question

The code implements an example of a Pollard rho() function for finding a factor of a positive integer, n. I've examined some of the code in the Julia "Primes" package that runs rapidly in an attempt to speedup the pollard_rho() function, all to no avail. The code should execute n = 1524157897241274137 in approximately 100 mSec to 30 Sec (Erlang, Haskell, Mercury, SWI Prolog) but takes about 3 to 4 minutes on JuliaBox, IJulia, and the Julia REPL. How can I make this go fast?

pollard_rho(1524157897241274137) = 1234567891

__precompile__()
module Pollard

export pollard_rho

function pollard_rho{T<:Integer}(n::T)
    f(x::T, r::T, n) = rem(((x ^ T(2)) + r), n)
    r::T = 7; x::T = 2; y::T = 11; y1::T = 11; z::T = 1
    while z == 1
        x  = f(x, r, n)
        y1 = f(y, r, n)
        y  = f(y1, r, n)
        z  = gcd(n, abs(x - y))
    end
    z >= n ? "error" : z
end

end # module

Answer 1

There are quite a few problems with type instability here.

1. Don't return either the string "error" or a result; instead explicitly call error().

2. As Chris mentioned, x and r ought to be annotated to be of type T, else they will be unstable.

There also seems to be a potential problem with overflow. A solution is to widen in the squaring step before truncating back to type T.

function pollard_rho{T<:Integer}(n::T)
    f(x::T, r::T, n) = rem(Base.widemul(x, x) + r, n) % T
    r::T = 7; x::T = 2; y::T = 11; y1::T = 11; z::T = 1
    while z == 1
        x  = f(x, r, n)
        y1 = f(y, r, n)
        y  = f(y1, r, n)
        z  = gcd(n, abs(x - y))
    end
    z >= n ? error() : z
end

After making these changes the function will run as fast as you could expect.

julia> @btime pollard_rho(1524157897241274137)
  4.128 ms (0 allocations: 0 bytes)
1234567891

To find these problems with type instability, use the @code_warntype macro.