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

possible number theory functions to add #122

Open
jmichel7 opened this issue Jun 27, 2022 · 6 comments
Open

possible number theory functions to add #122

jmichel7 opened this issue Jun 27, 2022 · 6 comments

Comments

@jmichel7
Copy link
Contributor

jmichel7 commented Jun 27, 2022

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

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

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

@oscardssmith
Copy link
Member

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

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

adrsm108 commented Aug 4, 2022

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

oscardssmith commented Aug 4, 2022

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),

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

No branches or pull requests

4 participants