possible number theory functions to add

[jmichel7]

Here is a example of two number theory functions made fast by eachfactor

"""
`prime_residues(n)` the numbers less than `n` and prime to `n`

julia> [prime_residues(24)]
1-element Vector{Vector{Int64}}:
 [1, 5, 7, 11, 13, 17, 19, 23]
"""
function prime_residues(n)
  if n==1 return [0] end
  pp=trues(n) # use a sieve to go fast
  for (p,np) in eachfactor(n)
    pp[p:p:n].=false
  end
  (1:n)[pp]
end

"""
`divisors(n)` the increasing list of divisors of `n`.

julia> [divisors(24)]
1-element Vector{Vector{Int64}}:
 [1, 2, 3, 4, 6, 8, 12, 24]
"""
function divisors(n::Integer)::Vector{Int}
  sort(vec(map(prod,Iterators.product((p.^(0:m) for (p,m) in eachfactor(n))...))))
end

[oscardssmith]

I think it might be good to make a new package for number theory functions (which depends on Primes.jl). There are a bunch of them, and (IMO) Primes.jl should focus on doing the basics well (factoring, primality checks, and sieving).

[jmichel7]

Then you might want to move totientand radical to this new package, since they are no more basic than the two above...

[oscardssmith]

yeah, I'm not really sure what I want here. On the one hand, these are simple relatively optimal implementations that shouldn't take much maintenance, on the other, they do seem somewhat disconnected to the core of Primes. Does anyone else have strong opinions either way on this one?

[wherrera10]

I find I have used this one frequently since Julia 0.7: using Primes

function factors(n)
    f = [one(n)]
    for (p, e) in factor(n)
        f = reduce(vcat, [f * p^j for j in 1:e], init = f)
    end
    return length(f) == 1 ? [one(n), n] : sort!(f)
end

which is more or less divisors above in a different form. So I think divisors belongs in Primes for me more than totient does.

[adrsm108]

I'd also agree that Primes needs a divisors function. I'll throw in my own implementation, which I've tuned a bit for performance since I find myself in need of it surprisingly often. @oscardssmith would you consider a pull request?

"""
    divisors(n::T) -> Vector{T}

Returns all positive divisors of the integer `n`.

For an integer `n` with prime factorization `n = p₁^k₁ ⋯ pₘ^kₘ`, `divisors(n)`
returns a vector of length (k₁ + 1)⋯(kₘ + 1) containing the divisors of `n` in
lexicographical (rather than numerical) order.

```julia
julia> divisors(60)
12-element Vector{Int64}:
  1      # 1
  2      # 2
  4      # 2 * 2
  3      # 3
  6      # 3 * 2
 12      # 3 * 2 * 2
  5      # 5
 10      # 5 * 2
 20      # 5 * 2 * 2
 15      # 5 * 3
 30      # 5 * 3 * 2
 60      # 5 * 3 * 2 * 2

divisors(-n) is equivalent to divisors(n).

divisors(0) returns []. """ function divisors(n::T)::Vector{T} where {T <: Integer} if iszero(n) return T[] elseif isone(n) return T[n] elseif n < 0 return divisors(-n) else return divisors(factor(n)) end end

""" divisors(factors::Factorization{T}) -> Vector{T}

Returns divisors of the number whose factorization is given by factors. """ function divisors(factors::Primes.Factorization{T})::Vector{T} where {T <: Integer} pe = factors.pe

if isempty(pe)
    return T[one(T)] # n == 1
elseif pe[1][1] == 0 # n == 0
    return T[]
elseif pe[1][1] == -1 # n < 0
    if length(pe) == 1 # n == -1
        return T[one(T)]
    else
        pe = pe[2:end]
    end
end

i::Int = 1
m::Int = 1
divs = Vector{T}(undef, prod(x -> x.second + 1, pe))
divs[i] = one(T)

for (p, k) in pe
    i = 1
    for _ in 1:k
        for j in i:(i + m - 1)
            divs[j + m] = divs[j] * p
        end
        i += m
    end
    m += i - 1
end

return divs

end

[oscardssmith]

Given the widespread interest in this function, I think we should probably add it. If you make a PR, I would be happy to merge it (possibly with a few minor changes),