diff --git a/src/numbers.jl b/src/numbers.jl
index 4039f9f..f1c2bb3 100644
--- a/src/numbers.jl
+++ b/src/numbers.jl
@@ -151,6 +151,11 @@ function stirlings1(n::Int, k::Int, signed::Bool=false)
     return (n - 1) * stirlings1(n - 1, k) + stirlings1(n - 1, k - 1)
 end
 
+"""
+    stirlings2(n::Int, k::Int)
+
+Compute the Stirling number of the second kind, `S(n,k)`.
+"""
 function stirlings2(n::Int, k::Int)
     if n < 0
         throw(DomainError(n, "n must be nonnegative"))
diff --git a/src/partitions.jl b/src/partitions.jl
index 305e23b..36ffb79 100644
--- a/src/partitions.jl
+++ b/src/partitions.jl
@@ -8,7 +8,7 @@ export
     #nextprod,
 
 
-#integer partitions
+# integer partitions
 
 struct IntegerPartitions
     n::Int
@@ -24,17 +24,56 @@ function Base.iterate(p::IntegerPartitions, xs = Int[])
 end
 
 """
-    partitions(n)
-
-Generate all integer arrays that sum to `n`. Because the number of partitions can be very
-large, this function returns an iterator object. Use `collect(partitions(n))` to get an
-array of all partitions. The number of partitions to generate can be efficiently computed
-using `length(partitions(n))`.
+    partitions(n::Integer)
+
+Generate all integer arrays that sum to `n`.
+
+Because the number of partitions can be very large,
+this function returns an iterator object.
+Use `collect(partitions(n))` to get an array of all partitions.
+
+The number of partitions to generate can be efficiently computed using
+`length(partitions(n))`.
+
+See also:
+- [`integer_partitions(n::Integer)`](@ref)
+    for a non-iterator version that returns all partitions as a array
+- [`partitions(n::Integer, m::Integer)`](@ref)
+    for partitions with exactly `m` parts.
+
+## Examples
+```jldoctest
+julia> collect(partitions(2))
+2-element Vector{Vector{Int64}}:
+ [2]
+ [1, 1]
+
+julia> collect(partitions(3))
+3-element Vector{Vector{Int64}}:
+ [3]
+ [2, 1]
+ [1, 1, 1]
+
+julia> integer_partitions(3)
+3-element Vector{Vector{Int64}}:
+ [1, 1, 1]
+ [2, 1]
+ [3]
+
+julia> first(partitions(10))
+1-element Vector{Int64}:
+ 10
+
+julia> length(partitions(10))
+42
+```
+
+# References
+- [Integer partition - Wikipedia](https://en.wikipedia.org/wiki/Integer_partition)
 """
 partitions(n::Integer) = IntegerPartitions(n)
 
 
-
 function nextpartition(n, as)
     isempty(as) && return Int[n]
 
@@ -97,12 +136,47 @@ Base.length(f::FixedPartitions) = npartitions(f.n,f.m)
 Base.eltype(f::FixedPartitions) = Vector{Int}
 
 """
-    partitions(n, m)
-
-Generate all arrays of `m` integers that sum to `n`. Because the number of partitions can
-be very large, this function returns an iterator object. Use `collect(partitions(n, m))` to
-get an array of all partitions. The number of partitions to generate can be efficiently
-computed using `length(partitions(n, m))`.
+    partitions(n::Integer, m::Integer)
+
+Generate all integer partitions of `n` into exactly `m` parts, that sum to `n`.
+
+Because the number of partitions can be very large,
+this function returns an iterator object.
+Use `collect(partitions(n, m))` to get an array of all partitions.
+
+The number of partitions to generate can be efficiently computed using
+`length(partitions(n, m))`.
+
+See also: [`partitions(n::Integer)`](@ref)
+
+## Examples
+```jldoctest
+julia> collect(partitions(4))
+5-element Vector{Vector{Int64}}:
+ [4]
+ [3, 1]
+ [2, 2]
+ [2, 1, 1]
+ [1, 1, 1, 1]
+
+julia> collect(partitions(4, 2))
+2-element Vector{Vector{Int64}}:
+ [3, 1]
+ [2, 2]
+
+julia> collect(partitions(4, 4))
+1-element Vector{Vector{Int64}}:
+ [1, 1, 1, 1]
+
+julia> collect(partitions(4, 5))
+Vector{Int64}[]
+
+julia> partitions(4, 0)
+ERROR: DomainError with (4, 0):
+n and m must be positive
+Stacktrace:
+[...]
+```
 """
 partitions(n::Integer, m::Integer) =
     n >= 1 && m >= 1 ?
@@ -176,11 +250,45 @@ Base.eltype(p::SetPartitions) = Vector{Vector{eltype(p.s)}}
 """
     partitions(s::AbstractVector)
 
-Generate all set partitions of the elements of an array `s`, represented as arrays of
-arrays. Because the number of partitions can be very large, this function returns an
-iterator object. Use `collect(partitions(s))` to get an array of all partitions. The
-number of partitions to generate can be efficiently computed using
-`length(partitions(s))`.
+Generate all set partitions of the elements of an array `s`,
+represented as arrays of arrays.
+
+Because the number of partitions can be very large,
+this function returns an iterator object.
+Use `collect(partitions(s))` to get an array of all partitions.
+
+The number of partitions of an `n`-element set
+is given by the n-th Bell number `Bn`:
+`length(partitions(s)) == bellnum(length(s))`.
+
+See also: [`bellnum`](@ref)
+
+# Examples
+```jldoctest
+julia> collect(partitions([1, 1]))
+2-element Vector{Vector{Vector{Int64}}}:
+ [[1, 1]]
+ [[1], [1]]
+
+julia> collect(partitions(-1:-1:-2))
+2-element Vector{Vector{Vector{Int64}}}:
+ [[-1, -2]]
+ [[-1], [-2]]
+
+julia> collect(partitions('a':'c'))
+5-element Vector{Vector{Vector{Char}}}:
+ [['a', 'b', 'c']]
+ [['a', 'b'], ['c']]
+ [['a', 'c'], ['b']]
+ [['a'], ['b', 'c']]
+ [['a'], ['b'], ['c']]
+
+julia> length(partitions(1:10)) == bellnum(10)
+true
+```
+
+# References
+- [Partition of a set - Wikipedia](https://en.wikipedia.org/wiki/Partition_of_a_set)
 """
 partitions(s::AbstractVector) = SetPartitions(s)
 
@@ -257,11 +365,45 @@ Base.eltype(p::FixedSetPartitions) = Vector{Vector{eltype(p.s)}}
     partitions(s::AbstractVector, m::Int)
 
 Generate all set partitions of the elements of an array `s` into exactly `m` subsets,
-represented as arrays of arrays. Because the number of partitions can be very large,
-this function returns an iterator object. Use `collect(partitions(s, m))` to get
-an array of all partitions. The number of partitions into `m` subsets is equal to the
-Stirling number of the second kind, and can be efficiently computed using
-`length(partitions(s, m))`.
+represented as arrays of arrays.
+
+Because the number of partitions can be very large,
+this function returns an iterator object.
+Use `collect(partitions(s, m))` to get an array of all partitions.
+
+The number of partitions into `m` subsets is equal to
+the Stirling number of the second kind,
+and can be efficiently computed using `length(partitions(s, m))`.
+
+See also: [`stirlings2(n::Int, k::Int)`](@ref)
+
+# Examples
+```jldoctest
+julia> collect(partitions('a':'c', 3))
+1-element Vector{Vector{Vector{Char}}}:
+ [['a'], ['b'], ['c']]
+
+julia> collect(partitions([1, 1, 1], 2))
+3-element Vector{Vector{Vector{Int64}}}:
+ [[1, 1], [1]]
+ [[1, 1], [1]]
+ [[1], [1, 1]]
+
+julia> collect(partitions(1:3, 2))
+3-element Vector{Vector{Vector{Int64}}}:
+ [[1, 2], [3]]
+ [[1, 3], [2]]
+ [[1], [2, 3]]
+
+julia> stirlings2(3, 2)
+3
+
+julia> length(partitions(1:10, 3)) == stirlings2(10, 3)
+true
+```
+
+# References
+- [Partition of a set - Wikipedia](https://en.wikipedia.org/wiki/Partition_of_a_set)
 """
 partitions(s::AbstractVector, m::Int) =
     length(s) >= 1 && m >= 1 ?
@@ -394,15 +536,31 @@ end
 """
     prevprod(a::Vector{Int}, x)
 
-Previous integer not greater than `x` that can be written as ``\\prod k_i^{p_i}`` for
-integers ``p_1``, ``p_2``, etc.
+Find the largest integer not greater than `x`
+that can be expressed as a product of powers of the elements in `a`.
+
+This function computes the largest value `y ≤ x` that can be written as:
+```math
+y = \\prod a_i^{n_i}
+  = a_1^{n_1} a_2^{n_2} \\cdots a_k^{n_k}
+  \\leq x
+```
+where ``n_i`` is a non-negative integer, `k` is the length of Vector `a`.
+
+# Examples
+```jldoctest
+julia> prevprod([10], 1000)   # 1000 = 10^3
+1000
 
-For integers ``i_1``, ``i_2``, ``i_3``, this is equivalent to finding the largest ``x``
-such that
+julia> prevprod([2, 5], 30)   # 25 = 2^0 * 5^2
+25
 
-``i_1^{n_1} i_2^{n_2} i_3^{n_3} \\leq x``
+julia> prevprod([2, 3], 100)  # 96 = 2^5 * 3^1
+96
 
-for integers ``n_1``, ``n_2``, ``n_3``.
+julia> prevprod([2, 3, 5], 1) # 1 = 2^0 * 3^0 * 5^0
+1
+```
 """
 function prevprod(a::Vector{Int}, x)
     if x > typemax(Int)
@@ -446,11 +604,43 @@ end
 """
     integer_partitions(n)
 
-List the partitions of the integer `n`.
+Generates all partitions of the integer `n` as a list of integer arrays,
+where each partition represents a way to write `n` as a sum of positive integers.
+
+See also: [`partitions(n::Integer)`](@ref)
 
 !!! note
     The order of the resulting array is consistent with that produced by the computational
     discrete algebra software GAP.
+
+# Examples
+```jldoctest
+julia> integer_partitions(2)
+2-element Vector{Vector{Int64}}:
+ [1, 1]
+ [2]
+
+julia> integer_partitions(3)
+3-element Vector{Vector{Int64}}:
+ [1, 1, 1]
+ [2, 1]
+ [3]
+
+julia> collect(partitions(3))
+3-element Vector{Vector{Int64}}:
+ [3]
+ [2, 1]
+ [1, 1, 1]
+
+julia> integer_partitions(-1)
+ERROR: DomainError with -1:
+n must be nonnegative
+Stacktrace:
+[...]
+```
+
+# References
+- [Integer partition - Wikipedia](https://en.wikipedia.org/wiki/Integer_partition)
 """
 function integer_partitions(n::Integer)
     if n < 0
@@ -474,12 +664,46 @@ function integer_partitions(n::Integer)
 end
 
 
-
-#Noncrossing partitions
+# Noncrossing partitions
 
 const _cmp = cmp
 
-#Produces noncrossing partitions of length n
+"""
+    ncpartitions(n::Int)
+
+Generates all noncrossing partitions of a set of `n` elements,
+returning them as a `Vector` of partition representations.
+
+The number of noncrossing partitions of an `n`-element set
+is given by the n-th Catalan number `Cn`:
+`length(ncpartitions(n)) == catalannum(n)`.
+
+See also: [`catalannum`](@ref)
+
+# Examples
+```jldoctest
+julia> ncpartitions(1)
+1-element Vector{Vector{Vector{Int64}}}:
+ [[1]]
+
+julia> ncpartitions(3)
+5-element Vector{Vector{Vector{Int64}}}:
+ [[1], [2], [3]]
+ [[1], [2, 3]]
+ [[1, 2], [3]]
+ [[1, 3], [2]]
+ [[1, 2, 3]]
+
+julia> catalannum(3)
+5
+
+julia> length(ncpartitions(6)) == catalannum(6)
+true
+```
+
+# References
+- [Noncrossing partition - Wikipedia](https://en.wikipedia.org/wiki/Noncrossing_partition)
+"""
 function ncpartitions(n::Int)
     partitions = Vector{Vector{Int}}[]
     _ncpart!(1,n,n,Vector{Int}[], partitions)
diff --git a/test/partitions.jl b/test/partitions.jl
index 4e42ffb..fba1eb6 100644
--- a/test/partitions.jl
+++ b/test/partitions.jl
@@ -1,59 +1,164 @@
 @testset "partitions" begin
-    @test collect(partitions(4)) == Any[[4], [3, 1], [2, 2], [2, 1, 1], [1, 1, 1, 1]]
-    @test collect(partitions(8, 3)) == Any[[6, 1, 1], [5, 2, 1], [4, 3, 1], [4, 2, 2], [3, 3, 2]]
-    @test collect(partitions(8, 1)) == Any[[8]]
-    @test collect(partitions(8, 9)) == []
-    @test collect(partitions([1, 2, 3])) == Any[Any[[1, 2, 3]], Any[[1, 2], [3]], Any[[1, 3], [2]], Any[[1], [2, 3]], Any[[1], [2], [3]]]
-    @test collect(partitions([1, 2, 3, 4], 3)) == Any[Any[[1, 2], [3], [4]], Any[[1, 3], [2], [4]], Any[[1], [2, 3], [4]],
-        Any[[1, 4], [2], [3]], Any[[1], [2, 4], [3]], Any[[1], [2], [3, 4]]]
-    @test collect(partitions([1, 2, 3, 4], 1)) == Any[Any[[1, 2, 3, 4]]]
-    @test collect(partitions([1, 2, 3, 4], 5)) == []
-
-    @inferred first(partitions(4))
-    @inferred first(partitions(8, 3))
-    @inferred first(partitions([1, 2, 3]))
-    @inferred first(partitions([1, 2, 3, 4], 3))
-
-    @test isa(collect(partitions(4)), Vector{Vector{Int}})
-    @test isa(collect(partitions(8, 3)), Vector{Vector{Int}})
-    @test isa(collect(partitions([1, 2, 3])), Vector{Vector{Vector{Int}}})
-    @test isa(collect(partitions([1, 2, 3, 4], 3)), Vector{Vector{Vector{Int}}})
-
-    @test length(partitions(0)) == 1
-    @test length(partitions(-1)) == 0
-    @test length(collect(partitions(30))) == length(partitions(30))
-    @test length(collect(partitions(90, 4))) == length(partitions(90, 4))
-    @test length(collect(partitions('a':'h'))) == length(partitions('a':'h'))
-    @test length(collect(partitions('a':'h', 5))) == length(partitions('a':'h', 5))
-
-    # integer_partitions
-    @test integer_partitions(0) == []
-    @test integer_partitions(5) == Any[[1, 1, 1, 1, 1], [2, 1, 1, 1], [2, 2, 1], [3, 1, 1], [3, 2], [4, 1], [5]]
-    @test_throws DomainError integer_partitions(-1)
-
-    @test_throws ArgumentError prevprod([2, 3, 5], Int128(typemax(Int)) + 1)
-    @test prevprod([2, 3, 5], 30) == 30
-    @test prevprod([2, 3, 5], 33) == 32
-
-    # noncrossing partitions
-    let nc4 = ncpartitions(4)
-        @test nc4 == Any[
-            Any[[1], [2], [3], [4]],
-            Any[[1], [2], [3, 4]],
-            Any[[1], [2, 3], [4]],
-            Any[[1], [2, 4], [3]],
-            Any[[1], [2, 3, 4]],
-            Any[[1, 2], [3], [4]],
-            Any[[1, 2], [3, 4]],
-            Any[[1, 3], [2], [4]],
-            Any[[1, 4], [2], [3]],
-            Any[[1, 4], [2, 3]],
-            Any[[1, 2, 3], [4]],
-            Any[[1, 3, 4], [2]],
-            Any[[1, 2, 4], [3]],
-            Any[[1, 2, 3, 4]],
+
+    @testset "partitions(n::Integer)" begin
+        @test_broken collect(partitions(0)) == [[]]
+        @test collect(partitions(1)) == [[1]]
+        @test collect(partitions(2)) == [[2], [1, 1]]
+        @test collect(partitions(3)) == [[3], [2, 1], [1, 1, 1]]
+        @test collect(partitions(4)) == [[4], [3, 1], [2, 2], [2, 1, 1], [1, 1, 1, 1]]
+
+        # Partition function p(n):  https://oeis.org/A000041
+        @test length(partitions(-1)) == 0
+        @test length(partitions(0)) == 1
+        @test length(partitions(1)) == 1
+        @test length(partitions(40)) == 37338
+        @test length(partitions(49)) == 173525
+        @test length(collect(partitions(30))) == length(partitions(30))
+
+        # Type stable
+        @inferred first(partitions(4))
+        @test isa(collect(partitions(4)), Vector{Vector{Int}})
+    end
+
+    @testset "partitions(n::Integer, m::Integer)" begin
+        @test collect(partitions(8, 1)) == [[8]]
+        @test collect(partitions(8, 3)) == [[6, 1, 1], [5, 2, 1], [4, 3, 1], [4, 2, 2], [3, 3, 2]]
+        @test collect(partitions(8, 8)) == [ones(Int, 8)]
+        @test collect(partitions(8, 9)) == []
+        @test collect(partitions(90, 90)) == [ones(Int, 90)]
+        @test collect(partitions(90, 91)) == []
+
+        # legnth
+        @test length(partitions(8, 1)) == 1
+        @test length(partitions(8, 3)) == 5
+        @test length(partitions(8, 8)) == 1
+        @test length(partitions(8, 9)) == 0
+        @test length(partitions(90, 1)) == 1
+        # gap> NrPartitions(90, 4);
+        @test length(partitions(90, 4)) == 5231
+        @test length(collect(partitions(90, 4))) == length(partitions(90, 4))
+        @test length(partitions(90, 90)) == 1
+
+        # Type stable
+        @inferred first(partitions(8, 3))
+        @test isa(collect(partitions(8, 3)), Vector{Vector{Int}})
+
+        # errors
+        @test_throws DomainError partitions(-1, 1)
+        @test_throws DomainError partitions(8, 0)
+        @test_throws DomainError partitions(8, -1)
+    end
+
+    @testset "partitions(s::AbstractVector)" begin
+        @test collect(partitions([1])) == [[[1]]]
+        @test collect(partitions([1, 2])) == [[[1, 2]], [[1], [2]]]
+        @test collect(partitions([1, 2, 3])) == [[[1, 2, 3]], [[1, 2], [3]], [[1, 3], [2]], [[1], [2, 3]], [[1], [2], [3]]]
+        @test collect(partitions(1:3)) == collect(partitions([1, 2, 3]))
+        @test collect(partitions('a':'b')) == [[['a', 'b']], [['a'], ['b']]]
+
+        # length:  https://oeis.org/A000110
+        @test length(partitions(1:10)) == 115975
+        @test length(partitions(1:20)) == 51724158235372
+        @test length(collect(partitions('a':'h'))) == length(partitions('a':'h'))
+
+        # Type stable
+        @inferred first(partitions([1, 2, 3]))
+        @test isa(collect(partitions([1, 2, 3])), Vector{Vector{Vector{Int}}})
+    end
+
+    @testset "partitions(s::AbstractVector, m::Int)" begin
+        @test collect(partitions(1:3, 2)) == [
+            [[1, 2], [3]],
+            [[1, 3], [2]],
+            [[1], [2, 3]],
+        ]
+        @test collect(partitions([1, 2, 3, 4], 1)) == [[[1, 2, 3, 4]]]
+        @test collect(partitions([1, 2, 3, 4], 3)) == [
+            [[1, 2], [3], [4]], [[1, 3], [2], [4]], [[1], [2, 3], [4]],
+            [[1, 4], [2], [3]], [[1], [2, 4], [3]], [[1], [2], [3, 4]]
         ]
-        @test length(nc4) == catalannum(4)
+        @test collect(partitions([1, 2, 3, 4], 4)) == [[[1], [2], [3], [4]]]
+        @test collect(partitions([1, 2, 3, 4], 5)) == []
+
+        # length: Stirling numbers of the second kind
+        #   https://dlmf.nist.gov/26.8#T2
+        @test length(partitions(1:3, 2)) == 3
+        @test length(partitions([1, 2, 3, 4], 3)) == 6
+        @test length(partitions(1:10, 4)) == 34105
+        @test length(partitions(1:10, 5)) == 42525
+        @test length(partitions(1:10, 5)) == stirlings2(10, 5)
+        @test length(partitions(1:10, 6)) == 22827
+        @test length(collect(partitions('a':'h', 5))) == length(partitions('a':'h', 5))
+
+        # Type stable
+        @inferred first(partitions([1, 2, 3, 4], 3))
+        @test isa(collect(partitions([1, 2, 3, 4], 3)), Vector{Vector{Vector{Int}}})
+    end
+
+    @testset "integer partitions" begin
+        @test_broken integer_partitions(0) == [[]]
+        @test integer_partitions(1) == [[1]]
+        @test integer_partitions(2) == [[1, 1], [2]]
+        @test integer_partitions(3) == [[1, 1, 1], [2, 1], [3]]
+        # gap> Partitions( 5 );
+        @test integer_partitions(5) == [
+            [1, 1, 1, 1, 1],
+            [2, 1, 1, 1],
+            [2, 2, 1],
+            [3, 1, 1],
+            [3, 2],
+            [4, 1],
+            [5]
+        ]
+        # integer_partitions <--> partitions(::Integer)
+        @test Set(integer_partitions(5)) == Set(partitions(5))
+
+        @test_throws DomainError integer_partitions(-1)
+    end
+
+    @testset "prevprod" begin
+        @test prevprod([2, 3, 5], 30) == 30         # 30 = 2 * 3 * 5
+        @test prevprod([2, 3, 5], 33) == 32         # 32 = 2^5
+        @test prevprod([2, 3, 5, 7], 420) == 420    # 420 = 2^2 * 3 * 5 * 7
+
+        # prime factor:  https://en.wikipedia.org/wiki/Table_of_prime_factors
+        prime_factors = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]
+        @test prevprod(prime_factors, 37) == 37
+        @test prevprod(prime_factors, 97) == 96     # 96 = 2^5 * 3
+        @test prevprod(prime_factors, 149) == 148   # 148 = 2^2 * 37
+        @test prevprod(prime_factors, 911) == 910   # 910 = 2 * 5 * 7 * 13
+        @test prevprod(prime_factors, 4999) == 4998 # 4998 = 2 * 3 * 7^2 * 17
+
+        # errors
+        @test_throws ArgumentError prevprod([2, 3, 5], Int128(typemax(Int)) + 1)
+    end
+
+    @testset "noncrossing partitions" begin
+        @test ncpartitions(0) == []
+        @test ncpartitions(1) == [[[1]]]
+        @test ncpartitions(2) == [[[1], [2]], [[1, 2]]]
+        # The 14 noncrossing partitions of a 4-element set ordered in a Hasse diagram
+        #   https://commons.wikimedia.org/wiki/File:Noncrossing_partitions_4;_Hasse.svg
+        @test ncpartitions(4) == [
+            [[1], [2], [3], [4]],
+            [[1], [2], [3, 4]],
+            [[1], [2, 3], [4]],
+            [[1], [2, 4], [3]],
+            [[1], [2, 3, 4]],
+            [[1, 2], [3], [4]],
+            [[1, 2], [3, 4]],
+            [[1, 3], [2], [4]],
+            [[1, 4], [2], [3]],
+            [[1, 4], [2, 3]],
+            [[1, 2, 3], [4]],
+            [[1, 3, 4], [2]],
+            [[1, 2, 4], [3]],
+            [[1, 2, 3, 4]],
+        ]
+
+        for n in 1:8
+            @test length(ncpartitions(n)) == catalannum(n)
+        end
     end
 
 end