use IntegerMathUtils

Author: oscardssmith

src/Primes.jl

@@ -320,11 +320,10 @@ Base.isempty(f::FactorIterator) = f.n == 1
 # Uses a variety of algorithms depending on the size of n to find a factor.
 #     https://en.algorithmica.org/hpc/algorithms/factorization
 # Cache of small factors for small numbers (n < 2^16)
-# Trial division of small (< 2^16) precomputed primes
-# Pollard rho's algorithm with Richard P. Brent optimizations
+# Trial division of small (< 2^16) precomputed primes and
+# Lenstra elliptic curve algorithm
 #     https://en.wikipedia.org/wiki/Trial_division
-#     https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm
-#     http://maths-people.anu.edu.au/~brent/pub/pub051.html
+#     https://en.wikipedia.org/wiki/Lenstra_elliptic-curve_factorization
 #
 
 """
@@ -668,7 +667,7 @@ given by `f`. This method may be preferable to [`totient(::Integer)`](@ref)
 when the factorization can be reused for other purposes.
 """
 function totient(f::Factorization{T}) where T <: Integer
-    if iszero(sign(f))
+    if !isempty(f) && first(first(f)) == 0
         throw(ArgumentError("ϕ(0) is not defined"))
     end
     result = one(T)
@@ -687,16 +686,8 @@ positive integers less than or equal to ``n`` that are relatively prime to
 The totient function of `n` when `n` is negative is defined to be
 `totient(abs(n))`.
 """
-function totient(n::T) where T<:Integer
-    n = abs(n)
-    if n == 0
-        throw(ArgumentError("ϕ(0) is not defined"))
-    end
-    result = one(T)
-    for (p, k) in eachfactor(n)
-        result *= p^(k-1) * (p - 1)
-    end
-    result
+function totient(n::Integer)
+    totient(factor(abs(n)))
 end
 
 # add: checked add (when makes sense), result of same type as first argument
@@ -984,90 +975,4 @@ julia> prevprimes(10, 10)
 prevprimes(start::T, n::Integer) where {T<:Integer} =
     collect(T, Iterators.take(prevprimes(start), n))
 
-"""
-    divisors(n::Integer) -> Vector
-
-Return a vector with the positive divisors of `n`.
-
-For a nonzero 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
-lexicographic (rather than numerical) order.
-
-`divisors(-n)` is equivalent to `divisors(n)`.
-
-For convenience, `divisors(0)` returns `[]`.
-
-# Example
-
-```jldoctest; filter = r"(\\s+#.*)?"
-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
-
-julia> divisors(-10)
-4-element Vector{Int64}:
-  1
-  2
-  5
- 10
-
-julia> divisors(0)
-Int64[]
-```
-"""
-function divisors(n::T) where {T<:Integer}
-    n = abs(n)
-    if iszero(n)
-        return T[]
-    elseif isone(n)
-        return [n]
-    else
-        return divisors(factor(n))
-    end
-end
-
-"""
-    divisors(f::Factorization) -> Vector
-
-Return a vector with the positive divisors of the number whose factorization is `f`.
-Divisors are sorted lexicographically, rather than numerically.
-"""
-function divisors(f::Factorization{T}) where T
-    sgn = sign(f)
-    if iszero(sgn) # n == 0
-        return T[]
-    elseif isempty(f) || length(f) == 1 && sgn < 0 # n == 1 or n == -1
-        return [one(T)]
-    end
-
-    i = m = 1
-    fs = rest(f, 1 + (sgn < 0))
-    divs = Vector{T}(undef, prod(x -> x.second + 1, fs))
-    divs[i] = one(T)
-
-    for (p, k) in fs
-        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
-
 end # module