diff --git a/docs/src/api.md b/docs/src/api.md
index 5e08842..d02cde0 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -33,4 +33,8 @@ Primes.primesmask
 ```@docs
 Primes.radical
 Primes.totient
+Primes.moebius
+Primes.liouville
+Primes.divisorcount
+Primes.divisorsum
 ```
diff --git a/src/Primes.jl b/src/Primes.jl
index 2aa77b0..219a6e8 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")
 
@@ -331,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])
@@ -346,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
@@ -375,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
@@ -440,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
@@ -466,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
@@ -552,6 +558,8 @@ 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
+=======
 If `interval` is provided, primes are sought in increments of `interval`.
 This can be useful to ensure the presence of certain divisors in `p-1`.
 The selected interval should be even.
@@ -614,6 +622,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
@@ -652,16 +661,135 @@ 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)
 prime(i::Integer) = prime(Int, i)
 
+
+
+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)
+
+"""
+    moebius(n::Integer) -> Int
+    moebius(f::Factorization) -> Int
+
+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.
+
+# Examples
+```jldoctest
+julia> map(moebius, 0:12)
+[-1, 1, -1, -1, 0, -1, 1, -1, 0, 0, 1, -1, 0]
+
+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))
+
+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)
+
+
+
+"""
+    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 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.
+
+# Examples
+```jldoctest
+julia> map(liouville, 1:12)
+[1, -1, -1, 1, -1, 1, -1, -1, 1, 1, -1, -1]
+
+julia> liouville(factor(12))
+-1
+```
+
+# External links
+* [Liouville function](https://en.wikipedia.org/wiki/Liouville_function) on Wikipedia.
+"""
+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::Factorization)
+
+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+1, assurenonzero(n))
+divisorcount(f::Factorization{T}) where T <: Integer = multfunct((p, e) -> e, assurenonzero(f))
+
+
+"""
+    divisorsum(n::Integer)
+    divisorsum(f::Factorization)
+    
+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::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
diff --git a/test/runtests.jl b/test/runtests.jl
index df6312e..9eb68ac 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
@@ -390,3 +375,52 @@ 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