JuliaMath / Primes.jl

Prime numbers in Julia
Other
99 stars 32 forks source link

Divisors function #124

Closed adrsm108 closed 2 years ago

adrsm108 commented 2 years ago

This PR introduces the function divisors, which takes an integer (or integer factorization), and returns a vector of its positive divisors.

Notes

Methods

The function has two methods: divisors(n::T)::Vector{T} where {T<:Integer}, and divisors(f::Factorization{T})::Vector{T} where {T<:Integer}. The former calls the latter directly. I played with the new eachfactor iterator looking for a more efficient technique in the integer case, but was ultimately unsuccessful. Knowing all prime exponents beforehand allows you to preallocate your output array, and it appears the performance benefits of doing so far outweigh the costs of creating a new Factorization.

Divisor sorting

As I've tried to make clear in its documentation, the output of this function is not numerically sorted. Instead, divisors come in lexicographic order: $$n=p_1^{k_1} p_2^{k_2} \dots \mapsto (1, p_1, p_1^2, \dots, p_1^{k_1}, p_2, p_1 p_2, p_1^2 p_2, \dots, p_1^{k_1} p_2, p_2^2, \dots ).$$ The algorithm builds the array in this order, and I think returning it without sorting is a suitable design choice because it makes fewer assumptions about the user's purpose. Numerical order, if desired, can be readily obtained by calling sort!. Alternatively, reshape(divisors(p1^k1 * p2^k2 * ... * pm^km), (k1 + 1, k2 + 1, ..., km + 1)) results in an m-dimensional array where dimension $i$ corresponds to powers of $p_i$—potentially a useful structure. For yet other uses (Dirichlet convolution comes to mind) order is irrelevant, and speed is the number one priority.

Zero

What should divisors(0) do? I throw an ArgumentError, same as totient, but perhaps it would be more convenient to quietly return an empty array?

adrsm108 commented 2 years ago

Related to #122

oscardssmith commented 2 years ago

I think an empty array would be a reasonable behavior for divisors (similar to factor). I think lexicographic order makes sense, but I would want to hear from others on that topic if possible. One other possibility worth considering: Can the ith divisior be determined lazily given the list of factors, and if so, is the right return a lazy object that computes divisiors on the fly? This would have the advantage of letting users do things like pick 5 random divisors without paying for the cost of computing all of them.

adrsm108 commented 2 years ago

I can change the behavior for divisors(0), no problem.

As for your second question, yeah, that would absolutely be possible. I don't know if a lazy iterator is really the best choice for divisors, though. If you will be using all of them (as is often the case for me), it's considerably more efficient to generate them at once, taking advantage of the fact that every divisor greater than 1 is a prime multiple of some previous divisor. I'd also point out that you'd have to be working with some truly enormous inputs for keeping them all in memory to become an issue, since the number-of-divisors function grows really slowly. (If you're curious like I was, you can look at the highly composite numbers OEIS A002182 to find that the upper bound on $\sigma_0(n)$ is 1600 (at n=2095133040) for an Int32, and 161280 (at n=9200527969062830400) for an Int64).

For cases where you only need a few, perhaps an nth_divisor function would be a better alternative? Something like

nth_divisor(f::Factorization, n::Int) = prod(
        ((p, s),) -> p^(s - 1),
        zip(keys(f), Base._ind2sub(Tuple(values(f) .+ 1), n));
        init = one(keytype(f))
    )

You'd still have to take the special cases like 0 and negative factorizations into account, but that's the basic idea.

oscardssmith commented 2 years ago

good point. In that case, I think this looks good as is.