Optimize a divisors algorithm in Julia?

Question

I'm trying to write an optimized algorithm in Julia to find all divisors of an integer. The heavy part, i.e. factoring, is done using Primes.jl and the algorithm seems quite fast. But I'm thinking of making it more faster, specifically, as I'll be iterating over the divisors at the end, I tried to use a generator to save memory allocations.

How can I optimize this function using generators? I tried some alternatives but hit an instability issue in Iterators.Flatten(). For example, Iterators.Flatten(((0 for i in 1:5), (1 for i in 1:5))) will yield eltype Any. I appreciate your suggestions.

using Primes

function divisors(n)
    d = Int64[1]
    for (p, e) in factor(n)
        r = 1
        l = length(d)
        for i in 1:e
            r *= p
            for j in 1:l
                push!(d, d[j]*r)
            end
        end
    end
    return sort(d) # not strictly required
end

Edit: Preallocating the output vector saves another 25% of time.

function divisors2(n)
    f = factor(n)
    m = prod(e + 1 for (p, e) in f)
    d = Vector{Int64}(undef, m)
    k = 1
    d[k] = 1
    for (p, e) in f
        r = 1
        l = k
        for i in 1:e
            r *= p
            for j in 1:l
                d[k+=1] = d[j]*r
            end
        end
    end
    return sort(d) # not strictly required
end

Answer 1

For the given example you could achieve type stability where flattening iterators created using the Iterators package.

Hence this is type stable (no more Any in the returned list):

julia> collect(Iterators.flatten((Iterators.repeated(0,5), Iterators.repeated(1,5))))
10-element Vector{Int64}:
 0
 0
 0
 0
 0
 1
 1
 1
 1
 1

Answer 2

The following code is an example of using a ton of Iterators to get the desired divisors. This is mainly based on tensorprod construction. Does this help?

_tensorprod(A,B) = Iterators.map(x->(x[2],x[1]),Iterators.product(A,B))

tensorprod(A,B) = Iterators.map(x->tuple(Iterators.flatten(x)...),_tensorprod(B,A))

using Primes
f = factor(3*25*7)
_f = map(x -> [x[1]^i for i=0:x[2]], sort(collect(f); rev=true))
factored_divisors_iter = foldl(tensorprod, _f)

With factored_divisors_iter you can do the following:

julia> vec(map(prod,factored_divisors_iter))
12-element Vector{Int64}:
   1
   3
   5
  15
  25
  75
   7
  21
  35
 105
 175
 525

Note that Iterators.map allows this last prod to be done in the iterator and the divisor list need not be materialized. Along the lines of:

julia> for d in Iterators.map(prod, factored_divisors_iter)
       println("$(prod(f)) ÷ $d = $(prod(f) ÷ d)")
       end
525 ÷ 1 = 525
525 ÷ 3 = 175
525 ÷ 5 = 105
525 ÷ 15 = 35
525 ÷ 25 = 21
525 ÷ 75 = 7
525 ÷ 7 = 75
525 ÷ 21 = 25
525 ÷ 35 = 15
525 ÷ 105 = 5
525 ÷ 175 = 3
525 ÷ 525 = 1

Answer 3

Hopefully, one cluttered line is not too cumbersome:

using Primes
f = factor(3*25*7)

divisors = Iterators.map(x->prod(map(y->y[1]^(y[2]-1),zip(tuple(keys(f)...),Tuple(x)))),CartesianIndices(tuple((values(f).+1)...)))

and now:

julia> for d in divisors
       println("$(d) * $(prod(f)÷d) = $(prod(f))")
       end
1 * 525 = 525
3 * 175 = 525
5 * 105 = 525
15 * 35 = 525
25 * 21 = 525
75 * 7 = 525
7 * 75 = 525
21 * 25 = 525
35 * 15 = 525
105 * 5 = 525
175 * 3 = 525
525 * 1 = 525

Explanation: To go over divisors, a multi-radix counter needs to be used (corresponding to prime powers), and this is provided by CartesianIndexing into high dimensional array.

Answer 4

For about the past 6 months I have been using:

""" Return the factors of n, including 1, n """
function factors(n::T)::Vector{T} where T <: Integer
    sort(vec(map(prod, Iterators.product((p.^(0:m) for (p, m) in eachfactor(n))...))))
end

Which works with any Integer type to return divisors in that type. The type instability you must work around comes in because when Prime.jl's eachfactor returns a Pair of {prime, exponent}, the exponent is always returned as an Int64. The advantage to using eachfactor over factor is exactly your goal of less allocation.