Adding some multiplicative functions

[reinermartin]

I am considering adding some basic prime-factorization-related functions to 'Primes.jl', in particular some multiplicative functions (in the sense of number theory, i.e, such that f(a*b) = f(a)f(b) if a and b are relative prime).

Primes.jl already has totient (but it does not appear to be documented; I might fix this). I would add at least the Möbius function, the divisor function, and the Liouville function. Some of these functions I have used to help solve Project Euler problems. In Python, such functions are provided by SmyPy .

Similar to the totient function I would add two versions each, to have efficiency in case the factorization is already known, e.g.:

moebiusmu(f::Factorization{T}) where T <: Integer
moebiusmu(n::Integer)

Also, I would add divisors(n::Integer) and divisors(f::Factorization{T}) which return an array (or an iterator) giving all the positive divisors of n.

What do you think about this?

[reinermartin]

Also, I would add a generic version, where you pass in a function f(p, e) which is applied to each prime power p^e in the factorization. So for the Moebius function, for example, one would use f(p, e) = -1 if e is 1, and zero otherwise.

[reinermartin]

To add a bit more color, the implementation can be as simple as

moebiusmu(f::Factorization{T}) where T <: Integer = reduce(*, e == 1 ? -1 : 0 for (p, e) in f; init=1)
moebiusmu(n::Integer) = moebiusmu(factor(n))

in the example of the Möbius function. For a function like divisorsigma, which can be large, there should be a choice whther to return a BigInt, or whether to return the result modulo some number.

[rfourquet]

Sounds like a good idea to me. But I would just call the functions moebius and divisors (i.e. not moebiusmu and divisorsigma).

[scheinerman]

I was looking for a divisors function, so I wrote one. It'd be great if this were already included in Primes. My preference would be for the function to return a Set of the positive divisors. In the mean time, here's the code I'm using:

"""
`divisors(n)` returns a `Set` of all the positive divisors of the integer `n`.
"""
function divisors(n::Integer)
    @assert n != 0 "Can't list all the divisors of zero"
    f = collect(factor(abs(n)))
    powers = tuple([x[2] for x in f]...)
    plist = [x[1] for x in f]
    np = length(plist)
    master_idx = tuple([0:t for t in powers]...)

    return Set(
        prod(plist[k]^idx[k] for k = 1:np) for idx in Iterators.product(master_idx...)
    )
end

I suspect a Julia expert could make this cleaner.

[scheinerman]

Also, if one only wants the number of divisors, the following is more efficient:

"""
`num_divisors(n)` returns the number of positive divisors of `n`.
"""
function num_divisors(n::Integer)
    @assert n != 0 "Can't count the number of divisors of zero"
    f = factor(abs(n))
    v = values(f)
    return prod(t + 1 for t in v)
end

[scheinerman]

I was looking for a divisors function, so I wrote one. It'd be great if this were already included in Primes. My preference would be for the function to return a Set of the positive divisors. In the mean time, here's the code I'm using:

"""
`divisors(n)` returns a `Set` of all the positive divisors of the integer `n`.
"""
function divisors(n::Integer)
    @assert n != 0 "Can't list all the divisors of zero"
    f = collect(factor(abs(n)))
    powers = tuple([x[2] for x in f]...)
    plist = [x[1] for x in f]
    np = length(plist)
    master_idx = tuple([0:t for t in powers]...)

    return Set(
        prod(plist[k]^idx[k] for k = 1:np) for idx in Iterators.product(master_idx...)
    )
end

I suspect a Julia expert could make this cleaner.

I have a slight bug in this code if n==±1; easy to fix with an if statement.