Is there a way to optimize this code for finding the divisors of a number?

Question

I've written a program in Julia to compute the divisors of a number n efficiently. The algorithm is original (as far as I know), and is loosely based on the Sieve of Eratosthenes. It essentially works like this:

For a given prime p, let p^k || n; every number m in the list satisfying p^{k+1} | m is removed, and this process is repeated for every prime p < n.

The primes are calculated in-situ using a traditional Sieve of 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

While this works perfectly fine, I find it inefficient to use two sieves simultaneously. I think this problem can be fixed by allowing each element in the sieve array to take three different values (corresponding to unchecked, divisor, and not divisor), but then this can no longer be implemented as a BitArray.

I've also tried modifying the function ν to make it more efficient:

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

Although this is more complicated, it's a little faster than the previous algorithm - especially when the power of p dividing n is large.

Note: I'm aware that there are much better algorithms for finding the divisors of a number. I'm just curious to see the extent to which the above algorithm can be optimised. As I mentioned earlier, using two sieves is rather cumbersome, and it would be nice to find a way to eliminate the traditional sieve for prime numbers without affecting the efficiency.

Answer 1

There's a couple things I can point out-

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

allocates twice for each array (the list comp and then the conversion). This will do it just fine: (Edit: @DNF points out superiority of Vector{Bool} here)

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

Next we can make sure to leverage any kind of faster indexing by using

for i in eachindex(psieve)

instead of the manual range. Then you can prepend the for-loop with

@inbounds for i in eachindex(psieve)

Or go even further, if you're on Julia 1.3 or later, and multi-thread it (assuming you've set JULIA_NUM_THREADS before running it)

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