From ee5b2ae776328cdf87d3e479ec5ce5df987a91c0 Mon Sep 17 00:00:00 2001
From: Reiner Martin <reinermartin@hotmail.com>
Date: Wed, 21 Aug 2019 23:28:12 +0200
Subject: [PATCH 01/11] Update api.md

Added Primes.totient
---
 docs/src/api.md | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/docs/src/api.md b/docs/src/api.md
index d82168b..16f3430 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -9,7 +9,6 @@ DocTestSetup = :(using Primes)
 ```@docs
 Primes.factor
 Primes.prodfactors
-Primes.radical
 ```
 
 ## Generating prime numbers
@@ -28,3 +27,10 @@ Primes.isprime
 Primes.ismersenneprime
 Primes.primesmask
 ```
+
+## Multiplicative functions
+
+```@docs
+Primes.radical
+Primes.totient
+```

From a345d167fad2e3d59c7a2b163e1e2ccacb88a26d Mon Sep 17 00:00:00 2001
From: Reiner Martin <reinermartin@hotmail.com>
Date: Wed, 21 Aug 2019 23:59:50 +0200
Subject: [PATCH 02/11] Update api.md

---
 docs/src/api.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/api.md b/docs/src/api.md
index 16f3430..5e08842 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -28,7 +28,7 @@ Primes.ismersenneprime
 Primes.primesmask
 ```
 
-## Multiplicative functions
+## Number-theoretic functions
 
 ```@docs
 Primes.radical

From 83c66dd88cdc05e77b3467387c3ba7754ae3d975 Mon Sep 17 00:00:00 2001
From: Reiner Martin <reinermartin@hotmail.com>
Date: Thu, 22 Aug 2019 13:58:23 +0200
Subject: [PATCH 03/11] Added moebiusmu

Added `moebiusmu` with two signatures (integers and factorizations as inputs), and added docstring
---
 src/Primes.jl | 22 ++++++++++++++++++++++
 1 file changed, 22 insertions(+)

diff --git a/src/Primes.jl b/src/Primes.jl
index e7ad2c9..656038c 100644
--- a/src/Primes.jl
+++ b/src/Primes.jl
@@ -646,4 +646,26 @@ julia> prime(3)
 prime(::Type{T}, i::Integer) where {T<:Integer} = i < 0 ? throw(DomainError(i)) : nextprime(T(2), i)
 prime(i::Integer) = prime(Int, i)
 
+"""
+    moebiusmu(n::Integer) -> Int
+    moebiusmu(f::Factorization) -> Int
+
+Compute the Moebius function of the integer `n`, which is ``(-1)`^k`` if `n`has `k` *distinct* prme factors,
+and 0 if it has a multiple prime factor.
+
+If the factorization of `n` is already known, it can passed into the function directly. This is
+useful, as finding the factorization can be expensive.
+
+# Example
+```jldoctest
+julia> map(moebiusmu, -2:10)
+[-1, 1, -1, 1, -1, -1, 0, -1, 1, -1, 0, 0, 1]
+
+julia> moebiusmu(factor(12))
+0
+```
+"""
+moebiusmu(f::Factorization{T}) where T <: Integer = reduce(*, e == 1 ? -1 : 0 for (p, e) in f if p ≥ 0; init=1)
+moebiusmu(n::Integer) = moebiusmu(factor(n))
+                     
 end # module

From c2d060f0afab5497a573f91d08085d571eb0b532 Mon Sep 17 00:00:00 2001
From: Reiner Martin <reinermartin@hotmail.com>
Date: Thu, 22 Aug 2019 13:59:14 +0200
Subject: [PATCH 04/11] Added Primes.moebiusmu

---
 docs/src/api.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/src/api.md b/docs/src/api.md
index 5e08842..63af3e2 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -33,4 +33,5 @@ Primes.primesmask
 ```@docs
 Primes.radical
 Primes.totient
+Primes.moebiusmu
 ```

From aeb72c1cb478fb8961d3729b3091dda550e334f8 Mon Sep 17 00:00:00 2001
From: Reiner Martin <reinermartin@hotmail.com>
Date: Thu, 22 Aug 2019 14:52:24 +0200
Subject: [PATCH 05/11] Extended doc for moebiusmu

---
 src/Primes.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/Primes.jl b/src/Primes.jl
index 656038c..93d0620 100644
--- a/src/Primes.jl
+++ b/src/Primes.jl
@@ -650,8 +650,8 @@ prime(i::Integer) = prime(Int, i)
     moebiusmu(n::Integer) -> Int
     moebiusmu(f::Factorization) -> Int
 
-Compute the Moebius function of the integer `n`, which is ``(-1)`^k`` if `n`has `k` *distinct* prme factors,
-and 0 if it has a multiple prime factor.
+Compute the Moebius function of the integer `n`, which is ``(-1)^k`` if `n` has `k` *distinct* prime factors,
+and 0 if it has a multiple prime factor. Also, if `n` is negative, we apply it to `abs(n)`, and `moebiusmu(0)=0`.
 
 If the factorization of `n` is already known, it can passed into the function directly. This is
 useful, as finding the factorization can be expensive.

From dfbe2f520be5b217094a21a841848e31ed70364c Mon Sep 17 00:00:00 2001
From: Reiner Martin <reinermartin@hotmail.com>
Date: Thu, 22 Aug 2019 15:44:30 +0200
Subject: [PATCH 06/11] Added liouvillelambda

---
 src/Primes.jl | 28 +++++++++++++++++++++++++++-
 1 file changed, 27 insertions(+), 1 deletion(-)

diff --git a/src/Primes.jl b/src/Primes.jl
index 93d0620..4a3abe3 100644
--- a/src/Primes.jl
+++ b/src/Primes.jl
@@ -667,5 +667,31 @@ julia> moebiusmu(factor(12))
 """
 moebiusmu(f::Factorization{T}) where T <: Integer = reduce(*, e == 1 ? -1 : 0 for (p, e) in f if p ≥ 0; init=1)
 moebiusmu(n::Integer) = moebiusmu(factor(n))
-                     
+
+"""
+    liouvillelambda(n::Integer) -> Int
+    liouvillelambda(f::Factorization) -> Int
+
+Compute Liouville's λ function, which gives ``(-1)^k`` where `k` is the number of prime factors of `n`
+counting multiplicity. We assume `n` to be non-zero. Also, if `n` is negative, we use the absolute
+value instead.
+
+If the factorization of `n` is already known, it can passed into the function directly. This is
+useful, as finding the factorization can be expensive.
+
+# Example
+```jldoctest
+julia> map(liouvillelambda, 1:10)
+[1, -1, -1, 1, -1, 1, -1, -1, 1, 1]
+
+julia> liouvillelambda(factor(30))
+-1
+```
+
+# External links
+* [Liouville function](https://en.wikipedia.org/wiki/Liouville_function) on Wikipedia.
+"""
+liouvillelambda(f::Dict{T,Int}) where T <: Integer = (-1)^reduce(+, e for (p, e) in f if p ≥ 0; init=0)
+liouvillelambda(n::Integer) = liouvillelambda(factor(n))
+                                          
 end # module

From da28472d25ce9f69b291f02216dea027fa90a0fe Mon Sep 17 00:00:00 2001
From: Reiner Martin <reinermartin@hotmail.com>
Date: Thu, 22 Aug 2019 18:41:31 +0200
Subject: [PATCH 07/11] Added divisorcount, divisorsum

---
 src/Primes.jl | 42 ++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 40 insertions(+), 2 deletions(-)

diff --git a/src/Primes.jl b/src/Primes.jl
index 4a3abe3..c7d8382 100644
--- a/src/Primes.jl
+++ b/src/Primes.jl
@@ -646,6 +646,9 @@ julia> prime(3)
 prime(::Type{T}, i::Integer) where {T<:Integer} = i < 0 ? throw(DomainError(i)) : nextprime(T(2), i)
 prime(i::Integer) = prime(Int, i)
 
+multfunction(a::Function, n::Integer) = multfunction(a, factor(n))
+multfunction(a::Function, f::Factorization) = reduce(*, a(p, e) for (p, e) in f if p ≥ 0; init=1)
+
 """
     moebiusmu(n::Integer) -> Int
     moebiusmu(f::Factorization) -> Int
@@ -665,8 +668,14 @@ julia> moebiusmu(factor(12))
 0
 ```
 """
-moebiusmu(f::Factorization{T}) where T <: Integer = reduce(*, e == 1 ? -1 : 0 for (p, e) in f if p ≥ 0; init=1)
-moebiusmu(n::Integer) = moebiusmu(factor(n))
+# moebiusmu(f::Dict{T,Int}) where T <: Integer = reduce(*, e == 1 ? -1 : 0 for (p, e) in f if p ≥ 0; init=1)
+# moebiusmu(n::Integer) = moebiusmu(factor(n))
+
+moebiusmu_aux(p, e) = e == 1 ? -1 : 0
+moebiusmu(n::Integer) = multfunction(moebiusmu_aux, n)
+moebiusmu(f::Factorization) = multfunction(moebiusmu_aux, f)
+
+
 
 """
     liouvillelambda(n::Integer) -> Int
@@ -693,5 +702,34 @@ julia> liouvillelambda(factor(30))
 """
 liouvillelambda(f::Dict{T,Int}) where T <: Integer = (-1)^reduce(+, e for (p, e) in f if p ≥ 0; init=0)
 liouvillelambda(n::Integer) = liouvillelambda(factor(n))
+
+
+
+"""
+    divisorcount(n::Integer)
+    divisorcount(f::Dict{T,Int})
+
+Compute the number of positive divisors of `n`.
+
+If the factorization of `n` is already known, it can passed into the function directly. This is
+useful, as finding the factorization can be expensive.
+"""
+
+divisorcount(n::Integer) = multfunction((p, e) -> e, n)
+divisorcount(f::Factorization) = multfunction((p, e) -> e, f)
+
+
+"""
+    divisorsum(n::Integer)
+    divisorsum(f::Dict{T,Int})
+    
+The sum of all positive divisors of `n`.
+
+If the factorization of `n` is already known, it can passed into the function directly. This is
+useful, as finding the factorization can be expensive.
+"""
+divisorsum(n::Integer) = multfunction((p, e) -> (p^(e+1)-1) ÷ (p-1), n)
+divisorsum(f::Factorization) = multfunction((p, e) -> (p^(e+1)-1) ÷ (p-1), f)
+
                                           
 end # module

From 5ffd4c88a62958983fc715355a8f11e0a9838f24 Mon Sep 17 00:00:00 2001
From: Reiner Martin <reinermartin@hotmail.com>
Date: Fri, 23 Aug 2019 07:25:29 +0200
Subject: [PATCH 08/11] Added liouville, divisorcount, divisorsum

---
 docs/src/api.md |  5 +++-
 src/Primes.jl   | 76 ++++++++++++++++++++++++++++++-------------------
 2 files changed, 50 insertions(+), 31 deletions(-)

diff --git a/docs/src/api.md b/docs/src/api.md
index 63af3e2..d02cde0 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -33,5 +33,8 @@ Primes.primesmask
 ```@docs
 Primes.radical
 Primes.totient
-Primes.moebiusmu
+Primes.moebius
+Primes.liouville
+Primes.divisorcount
+Primes.divisorsum
 ```
diff --git a/src/Primes.jl b/src/Primes.jl
index c7d8382..eb0290d 100644
--- a/src/Primes.jl
+++ b/src/Primes.jl
@@ -8,7 +8,8 @@ using Base: BitSigned
 using Base.Checked: checked_neg
 
 export isprime, primes, primesmask, factor, ismersenneprime, isrieselprime,
-       nextprime, prevprime, prime, prodfactors, radical, totient
+       nextprime, prevprime, prime, prodfactors, radical, totient,
+       moebius, liouville, divisorcount, divisorsum
 
 include("factorization.jl")
 
@@ -646,90 +647,105 @@ julia> prime(3)
 prime(::Type{T}, i::Integer) where {T<:Integer} = i < 0 ? throw(DomainError(i)) : nextprime(T(2), i)
 prime(i::Integer) = prime(Int, i)
 
-multfunction(a::Function, n::Integer) = multfunction(a, factor(n))
-multfunction(a::Function, f::Factorization) = reduce(*, a(p, e) for (p, e) in f if p ≥ 0; init=1)
+
+
+assurenonzero(n::Integer) = begin if n == 0 throw(ArgumentError("Argument must be non-zero")) end; n end
+
+function assurenonzero(f::Factorization{T}) where T <: Integer
+    if !isempty(f) && first(first(f)) == 0
+        throw(ArgumentError("Argument must be non-zero"))
+    end
+    f
+end
+
+multfunct(a::Function, n::Integer) = multfunct(a, factor(n))
+multfunct(a::Function, f::Factorization{T}) where T <: Integer =
+    reduce(*, a(p, e) for (p, e) in f if p ≥ 0; init=1)
 
 """
-    moebiusmu(n::Integer) -> Int
-    moebiusmu(f::Factorization) -> Int
+    moebius(n::Integer) -> Int
+    moebius(f::Factorization) -> Int
 
-Compute the Moebius function of the integer `n`, which is ``(-1)^k`` if `n` has `k` *distinct* prime factors,
-and 0 if it has a multiple prime factor. Also, if `n` is negative, we apply it to `abs(n)`, and `moebiusmu(0)=0`.
+Compute the Moebius function of the integer `n`, which is ``(-1)^k`` if `n` has `k` prime factors which
+are all distinct, and 0 if it has a multiple prime factor. Also, `moebius(0)=0`.
 
 If the factorization of `n` is already known, it can passed into the function directly. This is
 useful, as finding the factorization can be expensive.
 
 # Example
 ```jldoctest
-julia> map(moebiusmu, -2:10)
+julia> map(moebius, -2:10)
 [-1, 1, -1, 1, -1, -1, 0, -1, 1, -1, 0, 0, 1]
 
-julia> moebiusmu(factor(12))
+julia> moebius(factor(12))
 0
 ```
+
+# External links
+* [Möbius function](https://en.wikipedia.org/wiki/Möbius_function) on Wikipedia.
 """
 # moebiusmu(f::Dict{T,Int}) where T <: Integer = reduce(*, e == 1 ? -1 : 0 for (p, e) in f if p ≥ 0; init=1)
 # moebiusmu(n::Integer) = moebiusmu(factor(n))
 
-moebiusmu_aux(p, e) = e == 1 ? -1 : 0
-moebiusmu(n::Integer) = multfunction(moebiusmu_aux, n)
-moebiusmu(f::Factorization) = multfunction(moebiusmu_aux, f)
+m(p, e) = e == 1 ? -1 : 0
+moebius(n::Integer) = multfunct(m, n)
+moebius(f::Factorization{T}) where T <: Integer = multfunct(m, f)
 
 
 
 """
-    liouvillelambda(n::Integer) -> Int
-    liouvillelambda(f::Factorization) -> Int
+    liouville(n::Integer) -> Int
+    liouville(f::Factorization) -> Int
 
-Compute Liouville's λ function, which gives ``(-1)^k`` where `k` is the number of prime factors of `n`
-counting multiplicity. We assume `n` to be non-zero. Also, if `n` is negative, we use the absolute
-value instead.
+Compute Liouville's λ function, which gives ``(-1)^k`` where `k` is the number of prime factors
+of the non-zero integer `n`, counting multiplicity. 
 
 If the factorization of `n` is already known, it can passed into the function directly. This is
 useful, as finding the factorization can be expensive.
 
 # Example
 ```jldoctest
-julia> map(liouvillelambda, 1:10)
+julia> map(liouville, 1:10)
 [1, -1, -1, 1, -1, 1, -1, -1, 1, 1]
 
-julia> liouvillelambda(factor(30))
+julia> liouville(factor(-30))
 -1
 ```
 
 # External links
 * [Liouville function](https://en.wikipedia.org/wiki/Liouville_function) on Wikipedia.
 """
-liouvillelambda(f::Dict{T,Int}) where T <: Integer = (-1)^reduce(+, e for (p, e) in f if p ≥ 0; init=0)
-liouvillelambda(n::Integer) = liouvillelambda(factor(n))
-
+liouville(n::Integer) = liouville(factor(assurenonzero(n)))
+liouville(f::Factorization{T}) where T <: Integer =
+    (-1)^reduce(+, e for (p, e) in assurenonzero(f) if p ≥ 0; init=0)
 
 
 """
     divisorcount(n::Integer)
-    divisorcount(f::Dict{T,Int})
+    divisorcount(f::Factorization)
 
-Compute the number of positive divisors of `n`.
+Compute the number of positive divisors of the nonzero integer `n`.
 
 If the factorization of `n` is already known, it can passed into the function directly. This is
 useful, as finding the factorization can be expensive.
 """
 
-divisorcount(n::Integer) = multfunction((p, e) -> e, n)
-divisorcount(f::Factorization) = multfunction((p, e) -> e, f)
+divisorcount(n::Integer) = multfunct((p, e) -> e, assurenonzero(n))
+divisorcount(f::Factorization{T}) where T <: Integer = multfunct((p, e) -> e, assurenonzero(f))
 
 
 """
     divisorsum(n::Integer)
-    divisorsum(f::Dict{T,Int})
+    divisorsum(f::Factorization)
     
-The sum of all positive divisors of `n`.
+The sum of all positive divisors of the nonzero integer `n`.
 
 If the factorization of `n` is already known, it can passed into the function directly. This is
 useful, as finding the factorization can be expensive.
 """
-divisorsum(n::Integer) = multfunction((p, e) -> (p^(e+1)-1) ÷ (p-1), n)
-divisorsum(f::Factorization) = multfunction((p, e) -> (p^(e+1)-1) ÷ (p-1), f)
+δ(p, e) = (p^(e+1)-1) ÷ (p-1)
+divisorsum(n::Integer) = multfunct(δ, assurenonzero(n))
+divisorsum(f::Factorization{T}) where T <: Integer = multfunct(δ, assurenonzero(f))
 
                                           
 end # module

From 20b2e853196935fc521e1ca5d94db1be93bdd3eb Mon Sep 17 00:00:00 2001
From: Reiner Martin <reinermartin@hotmail.com>
Date: Fri, 23 Aug 2019 18:08:34 +0200
Subject: [PATCH 09/11] Added tests

Added tests for moebius(), liouville(), divisorcount(), divisorsum()
---
 test/runtests.jl | 80 +++++++++++++++++++++++++++++++++++++++---------
 1 file changed, 65 insertions(+), 15 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 6fd28e6..384c3d5 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -290,21 +290,6 @@ end
 @testset "totient() correctness" begin
     @test totient(2^4 * 3^4 * 5^4) == 216000
     @test totient(big"2"^1000) == big"2"^999
-
-    some_coprime_numbers = BigInt[
-        450000000, 1099427429702334733, 200252151854851, 1416976291499, 7504637909,
-        1368701327204614490999, 662333585807659, 340557446329, 1009091
-    ]
-
-    for i in some_coprime_numbers
-        for j in some_coprime_numbers
-            if i ≠ j
-                @test totient(i*j) == totient(i) * totient(j)
-            end
-        end
-        # can use directly with Factorization
-        @test totient(i) == totient(factor(i))
-    end
 end
 
 # check copy property for big primes relied upon in nextprime/prevprime
@@ -386,3 +371,68 @@ for T in (Int, UInt, BigInt)
     @test prodfactors(factor(Set, T(123456))) == 3858
     @test prod(factor(T(123456))) == 123456
 end
+
+# check multiplicativity for multiplicative functions, and whether f(factor(n)) = f(n)
+@testset "multiplicativity and consistency" begin
+    some_coprime_numbers = BigInt[
+        450000000, 1099427429702334733, 200252151854851, 1416976291499, 7504637909,
+        1368701327204614490999, 662333585807659, 340557446329, 1009091
+    ]
+
+    for i in some_coprime_numbers
+        for j in some_coprime_numbers
+            if i ≠ j
+                for fct in [totient, moebius, liouville, divisorcount, divisorsum]
+                    @test fct(i*j) == fct(i) * fct(j)
+                end
+            end
+        end
+        # can use directly with Factorization
+        for fct in [totient, moebius, liouville, divisorcount, divisorsum]
+            @test totient(i) == totient(factor(i))
+        end
+    end
+end
+
+@testset "Int moebius(), liouville(), divisorcount(), divisorsum()" begin
+    ns = [3^n+7^n+1 for n in 1:22]
+    moebius_res = [-1, -1, 1, 1, 0, 1, 1, -1, -1, 1, 1, 1, -1, 1, 1, 1, -1, 1, 1, -1, 0, 1]
+    liouville_res = [-1, -1, 1, 1, -1, 1, 1, -1, -1, 1, 1, 1, -1, 1, 1, 1, -1, 1, 1, -1, 1, 1]
+    divisorcount_res = [2, 2, 4, 4, 6, 4, 4, 8, 8, 4, 4, 4, 2, 4, 4, 4, 2, 4, 4, 8, 48, 4]
+    divisorsum_res = [12, 60, 432, 2688, 18420, 121176, 828072, 6416256, 46199040, 283101000, 2157276984, 13850631048, 96890604732, 694000596696, 5425800981552, 35789356202208, 232630643127372, 1628414668930224, 11475399007686000, 86013321721581408, 750825696336150240, 3950128513778110104]
+    @test map(moebius, ns) == map(moebius, -ns) == moebius_res
+    @test map(liouville, ns) == map(liouville, -ns) == liouville_res
+    @test map(divisorcount, ns) == map(divisorcount, -ns) == divisorcount_res
+    @test map(divisorsum, ns) == map(divisorsum, -ns) == divisorsum_res
+end
+
+@testset "BigInt moebius(), liouville(), divisorcount(), divisorsum()" begin
+    ns = BigInt[
+        450000000, 1099427429702334733, 200252151854851, 1416976291499, 7504637909,
+        1368701327204614490999, 662333585807659, 340557446329, 1009091
+    ]
+    moebius_res = [0, 0, 0, 0, 0, 0, 0, 0, -1]
+    liouville_res = [-1, -1, -1, 1, 1, 1, 1, 1, -1]
+    divisorcount_res = [216, 308, 90, 40, 24, 120, 40, 27, 8]
+    divisorsum_res = [1618651515, 1528455927220372220, 234749079860820, 1557002958720, 8042136960, 1442725739861582095920, 691292021934800, 353095503663, 1039584]
+    @test map(moebius, ns) == map(moebius, -ns) == moebius_res
+    @test map(liouville, ns) == map(liouville, -ns) == liouville_res
+    @test map(divisorcount, ns) == map(divisorcount, -ns) == divisorcount_res
+    @test map(divisorsum, ns) == map(divisorsum, -ns) == divisorsum_res
+end
+
+@testset "moebius(), liouville(), divisorcount(), divisorsum()" begin
+    some_coprime_numbers = BigInt[
+        450000000, 1099427429702334733, 200252151854851, 1416976291499, 7504637909,
+        1368701327204614490999, 662333585807659, 340557446329, 1009091
+    ]
+    moebius_res = [0, 0, 0, 0, 0, 0, 0, 0, -1]
+    liouville_res = [-1, -1, -1, 1, 1, 1, 1, 1, -1]
+    divisorcount_res = [216, 308, 90, 40, 24, 120, 40, 27, 8]
+    divisorsum_res = [1618651515, 1528455927220372220, 234749079860820, 1557002958720, \
+8042136960, 1442725739861582095920, 691292021934800, 353095503663, \
+1039584]
+    for n in some_coprime_numbers
+        @test moebius(n) == 
+    end
+end

From 6824af7efe07a1efe857cfd96ef7a9359b57e330 Mon Sep 17 00:00:00 2001
From: Reiner Martin <reinermartin@hotmail.com>
Date: Fri, 23 Aug 2019 18:11:50 +0200
Subject: [PATCH 10/11] Bugfix divisorcount, improved docs

---
 src/Primes.jl | 51 ++++++++++++++++++++++++++++++++++++++-------------
 1 file changed, 38 insertions(+), 13 deletions(-)

diff --git a/src/Primes.jl b/src/Primes.jl
index eb0290d..21475dd 100644
--- a/src/Primes.jl
+++ b/src/Primes.jl
@@ -332,6 +332,7 @@ true
 When `ContainerType == Set`, this returns the distinct prime
 factors as a set.
 
+# Examples
 ```julia
 julia> factor(Set, 100)
 Set([2,5])
@@ -347,6 +348,7 @@ factor(::Type{T}, n) where {T<:Any} = throw(MethodError(factor, (T, n)))
 """
     prodfactors(factors)
 
+# Examples
 Compute `n` (or the radical of `n` when `factors` is of type `Set` or
 `BitSet`) where `factors` is interpreted as the result of
 `factor(typeof(factors), n)`. Note that if `factors` is of type
@@ -376,6 +378,7 @@ Base.prod(factors::Factorization) = prodfactors(factors)
 Compute the radical of `n`, i.e. the largest square-free divisor of `n`.
 This is equal to the product of the distinct prime numbers dividing `n`.
 
+# Example
 ```jldoctest
 julia> radical(2*2*3)
 6
@@ -441,6 +444,7 @@ whether `M` is a valid Mersenne number; to be used with caution.
 Return `true` if the given Mersenne number is prime, and `false`
 otherwise.
 
+# Examples
 ```jldoctest
 julia> ismersenneprime(2^11 - 1)
 false
@@ -467,6 +471,7 @@ with `0 < k < Q`, `Q = 2^n - 1` and `n > 0`, also known as Riesel primes.
 Returns `true` if `R` is prime, and `false` otherwise or
 if the combination of k and n is not supported.
 
+# Examples
 ```jldoctest
 julia> isrieselprime(1, 2^11 - 1)  # == ismersenneprime(2^11 - 1)
 false
@@ -553,6 +558,7 @@ prevprime(n, -i). Note that for `n::BigInt`, the returned number is
 only a pseudo-prime (the function [`isprime`](@ref) is used
 internally). See also [`prevprime`](@ref).
 
+# Examples
 ```jldoctest
 julia> nextprime(4)
 5
@@ -604,6 +610,7 @@ The `i`-th largest prime not greater than `n` (in particular
 only a pseudo-prime (the function [`isprime`](@ref) is used internally). See
 also [`nextprime`](@ref).
 
+# Examples
 ```jldoctest
 julia> prevprime(4)
 3
@@ -635,13 +642,13 @@ end
 
 The `i`-th prime number.
 
+# Examples
 ```jldoctest
 julia> prime(1)
 2
 
 julia> prime(3)
 5
-
 ```
 """
 prime(::Type{T}, i::Integer) where {T<:Integer} = i < 0 ? throw(DomainError(i)) : nextprime(T(2), i)
@@ -672,10 +679,10 @@ are all distinct, and 0 if it has a multiple prime factor. Also, `moebius(0)=0`.
 If the factorization of `n` is already known, it can passed into the function directly. This is
 useful, as finding the factorization can be expensive.
 
-# Example
+# Examples
 ```jldoctest
-julia> map(moebius, -2:10)
-[-1, 1, -1, 1, -1, -1, 0, -1, 1, -1, 0, 0, 1]
+julia> map(moebius, 0:12)
+[-1, 1, -1, -1, 0, -1, 1, -1, 0, 0, 1, -1, 0]
 
 julia> moebius(factor(12))
 0
@@ -687,7 +694,7 @@ julia> moebius(factor(12))
 # moebiusmu(f::Dict{T,Int}) where T <: Integer = reduce(*, e == 1 ? -1 : 0 for (p, e) in f if p ≥ 0; init=1)
 # moebiusmu(n::Integer) = moebiusmu(factor(n))
 
-m(p, e) = e == 1 ? -1 : 0
+m(p, e) = p == 0 ? 0 : e == 1 ? -1 : 0
 moebius(n::Integer) = multfunct(m, n)
 moebius(f::Factorization{T}) where T <: Integer = multfunct(m, f)
 
@@ -703,12 +710,12 @@ of the non-zero integer `n`, counting multiplicity.
 If the factorization of `n` is already known, it can passed into the function directly. This is
 useful, as finding the factorization can be expensive.
 
-# Example
+# Examples
 ```jldoctest
-julia> map(liouville, 1:10)
-[1, -1, -1, 1, -1, 1, -1, -1, 1, 1]
+julia> map(liouville, 1:12)
+[1, -1, -1, 1, -1, 1, -1, -1, 1, 1, -1, -1]
 
-julia> liouville(factor(-30))
+julia> liouville(factor(12))
 -1
 ```
 
@@ -728,9 +735,18 @@ Compute the number of positive divisors of the nonzero integer `n`.
 
 If the factorization of `n` is already known, it can passed into the function directly. This is
 useful, as finding the factorization can be expensive.
+
+# Examples
+```jldoctest
+julia> map(divisorcount, 1:12)
+[1, 1, 1, 2, 1, 1, 1, 3, 2, 1, 1, 2]
+
+julia> divisorcount(factor(12))
+2
+```
 """
 
-divisorcount(n::Integer) = multfunct((p, e) -> e, assurenonzero(n))
+divisorcount(n::Integer) = multfunct((p, e) -> e+1, assurenonzero(n))
 divisorcount(f::Factorization{T}) where T <: Integer = multfunct((p, e) -> e, assurenonzero(f))
 
 
@@ -742,10 +758,19 @@ The sum of all positive divisors of the nonzero integer `n`.
 
 If the factorization of `n` is already known, it can passed into the function directly. This is
 useful, as finding the factorization can be expensive.
+
+# Examples
+```jldoctest
+julia> map(divisorsum, 1:12)
+[1, 3, 4, 7, 6, 12, 8, 15, 13, 18, 12, 28]
+
+julia> divisorsum(factor(12))
+28
+```
 """
-δ(p, e) = (p^(e+1)-1) ÷ (p-1)
+δ(p::Integer, e) = foldl((x, _) -> p*x+1, 1:e; init=1)
+δ(p::BigInt, e) = e > 4 ? foldl((x, _) -> p*x+1, 1:e; init=1) : (p^(e+1)-1) ÷ (p-1)
 divisorsum(n::Integer) = multfunct(δ, assurenonzero(n))
 divisorsum(f::Factorization{T}) where T <: Integer = multfunct(δ, assurenonzero(f))
-
-                                          
+                              
 end # module

From a2085ab6f7ca065fac0bfe4221276f254c8a3e8f Mon Sep 17 00:00:00 2001
From: Reiner Martin <reinermartin@hotmail.com>
Date: Fri, 23 Aug 2019 19:16:54 +0200
Subject: [PATCH 11/11] Deleted some nonsense

---
 test/runtests.jl | 16 ----------------
 1 file changed, 16 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 384c3d5..4da9a29 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -420,19 +420,3 @@ end
     @test map(divisorcount, ns) == map(divisorcount, -ns) == divisorcount_res
     @test map(divisorsum, ns) == map(divisorsum, -ns) == divisorsum_res
 end
-
-@testset "moebius(), liouville(), divisorcount(), divisorsum()" begin
-    some_coprime_numbers = BigInt[
-        450000000, 1099427429702334733, 200252151854851, 1416976291499, 7504637909,
-        1368701327204614490999, 662333585807659, 340557446329, 1009091
-    ]
-    moebius_res = [0, 0, 0, 0, 0, 0, 0, 0, -1]
-    liouville_res = [-1, -1, -1, 1, 1, 1, 1, 1, -1]
-    divisorcount_res = [216, 308, 90, 40, 24, 120, 40, 27, 8]
-    divisorsum_res = [1618651515, 1528455927220372220, 234749079860820, 1557002958720, \
-8042136960, 1442725739861582095920, 691292021934800, 353095503663, \
-1039584]
-    for n in some_coprime_numbers
-        @test moebius(n) == 
-    end
-end