我怎样才能加速这个Julia代码？

Question

该代码实现了一个Pollard rho()函数的示例,用于查找正整数n的因子.我已经检查了Julia"Primes"包中的一些代码,这些代码快速运行以试图加速pollard_rho()函数,但都无济于事.代码应该在大约100毫秒到30秒(Erlang,Haskell,Mercury,SWI Prolog)中执行n = 1524157897241274137,但在JuliaBox,IJulia和Julia REPL上需要大约3到4分钟.我怎样才能让它变得快速？

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

这里的类型不稳定存在很多问题.

1. 不要返回字符串"error"或结果; 而是显式调用error().

2. 正如克里斯提到,x并且r应该被注解为类型T,否则,他们将是不稳定的.

溢出似乎也存在潜在问题.一种解决方案是在截断回到类型之前在平方步骤中加宽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

进行这些更改后,该功能将以您预期的速度运行.

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

要查找类型不稳定的这些问题,请使用@code_warntype宏.