possible number theory functions to add
jmichel7 opened this issue · 6 comments
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
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).
Then you might want to move totient
and radical
to this new package, since they are no more basic than the two above...
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?
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.
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
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),