有没有一种优化这段代码以找到一个数的因数的方法？

Question

有没有一种优化这段代码以找到一个数的因数的方法？

algorithmjuliainteger-divisionsieve-of-eratosthenesnumber-theory

4

我用Julia编写了一个程序，可以高效地计算一个数n的因数。这个算法是原创的（据我所知），并且松散地基于Eratosthenes筛法。它的工作原理如下：

对于给定的质数p，令p^k || n；列表中满足p^{k+1} | m的每个数字m都被删除，并且对于每个小于n的质数p，重复此过程。

质数使用传统的Eratosthenes筛法在原地计算。

function ν(p, n)     #returns the smallest power of p that does not divide n
    q = 1

    for i = 0:n
        if mod(n, q) != 0
            return (i, q) 
        end

        q *= p
    end
end

function divisors(n)    #returns a list of divisors of n
    dsieve, psieve = BitArray([true for i = 1:n]), BitArray([true for i = 1:n])
    psieve[1] = false

    for i = 1:n
        if psieve[i] && dsieve[i]
            #sieving out the non-primes
            for j = i^2:i:n
                psieve[j] = false
            end

            #sieving out the non-divisors
            v = ν(i, n)[2]
            for j = v:v:n
                dsieve[j] = false
            end
        end
    end
    return dsieve #the code for converting this BitArray to an array of divisors has been omitted for clarity
end

虽然这个方法完全没问题，但同时使用两个筛子效率不高。我认为，通过允许筛选器数组中的每个元素取三个不同的值（对应于“未检查”、“除数”和“非除数”），可以解决这个问题，但是这不能再实现为“BitArray”。我还尝试修改函数“ν”，使其更加有效。

function ν₀(p, n)      #the same as ν, but implemented differently
    q = p
    while mod(n, q) == 0
        q = q^2
    end

    q = floor(Int64, √q)
    q < p ? 1 : q * ν₀(p, n÷q)    #change 1 to p to get the smallest power of p that does not divide n
end

尽管这种方法更加复杂，但在p除n的幂较大时，它比之前的算法要快一些。注意：我知道有更好的算法可以找到一个数的因子。我只是好奇上述算法可以被优化到什么程度。正如我之前提到的，使用两个筛子相当麻烦，如果能找到一种消除传统素数筛子而不影响效率的方法，那就太好了。

- Art

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Miles Lucas · Accepted Answer

我可以指出几点：

dsieve, psieve = BitArray([true for i = 1:n]), BitArray([true for i = 1:n])

每个数组都分配了两次内存（列表推导式，然后是转换）。下面的代码就可以很好地解决这个问题：（编辑：@DNF 指出Vector{Bool} 在这里更优秀）

dsieve = fill(true, n)
psieve = fill(true, n)

接下来，我们可以通过使用任何类型的更快索引来提高利用率。

for i in eachindex(psieve)

可以使用自动化范围（range）代替手动指定范围。然后，在for循环之前添加以下内容：

@inbounds for i in eachindex(psieve)

如果你正在使用 Julia 1.3 或更高版本，而且已经在运行之前设置了 JULIA_NUM_THREADS，那么你甚至可以将其进行多线程处理。

@inbounds Threads.@threads for i in eachindex(psieve)