Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Divisors function #124

Merged
merged 6 commits into from
Aug 21, 2022
Merged

Divisors function #124

merged 6 commits into from
Aug 21, 2022

Conversation

adrsm108
Copy link
Contributor

@adrsm108 adrsm108 commented Aug 5, 2022

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
Copy link
Contributor Author

adrsm108 commented Aug 5, 2022

Related to #122

@oscardssmith
Copy link
Member

oscardssmith commented Aug 5, 2022

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.

Implement `sign(::Factorization)` to easily check if the number
represented by a Factorization is lt, gt, or eq to 0, and use it when
identifying special cases in `divisors` and `totient`.

In the `divisors` function:
- No more wasteful allocations from pointlessly `collect`ing factors.
- Make `divisors(0)` return an empty vector instead of throwing.
- Remove unnecessary type annotations and let return types be inferred.
@adrsm108
Copy link
Contributor Author

adrsm108 commented Aug 6, 2022

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
Copy link
Member

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

@oscardssmith oscardssmith merged commit fde566c into JuliaMath:master Aug 21, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants