如何加速这个 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. 正如 Chris 所说，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 宏。