diff --git a/README.md b/README.md
index 73715f4..f7dc5a8 100644
--- a/README.md
+++ b/README.md
@@ -260,9 +260,9 @@ julia> ones(GaussMod{5},4)
 
 ### Iteration
 
-The `Mod{m}` type can be used as an iterator (in `for` statements and list comprehension):
+`allmods(Mod{m})` is an iterator over all elements of `Mod{m}` (usable in `for` statements and list comprehensions):
 ```
-julia> for a in Mod{5}
+julia> for a in allmods(Mod{5})
        println(a)
        end
 Mod{5}(0)
@@ -271,7 +271,7 @@ Mod{5}(2)
 Mod{5}(3)
 Mod{5}(4)
 
-julia> collect(Mod{6})
+julia> collect(allmods(Mod{6}))
 6-element Vector{Mod{6}}:
  Mod{6}(0)
  Mod{6}(1)
@@ -280,7 +280,7 @@ julia> collect(Mod{6})
  Mod{6}(4)
  Mod{6}(5)
 
-julia> [k*k for k ∈ Mod{7}]
+julia> [k*k for k ∈ allmods(Mod{7})]
 7-element Vector{Mod{7}}:
  Mod{7}(0)
  Mod{7}(1)
@@ -290,13 +290,13 @@ julia> [k*k for k ∈ Mod{7}]
  Mod{7}(4)
  Mod{7}(1)
 
-julia> prod(k for k in Mod{5} if k!=0) == -1  # Wilson's theorem
+julia> prod(k for k in allmods(Mod{5}) if k!=0) == -1  # Wilson's theorem
 true
 ```
 
-One can also use `GaussMod` as an iterator:
+Likewise `allmods(GaussMod{m})` iterates over all Gaussian residues:
 ```
-julia> for z in GaussMod{3}
+julia> for z in allmods(GaussMod{3})
        println(z)
        end
 GaussMod{3}(0 + 0im)
diff --git a/docs/src/index.md b/docs/src/index.md
index 73715f4..f7dc5a8 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -260,9 +260,9 @@ julia> ones(GaussMod{5},4)
 
 ### Iteration
 
-The `Mod{m}` type can be used as an iterator (in `for` statements and list comprehension):
+`allmods(Mod{m})` is an iterator over all elements of `Mod{m}` (usable in `for` statements and list comprehensions):
 ```
-julia> for a in Mod{5}
+julia> for a in allmods(Mod{5})
        println(a)
        end
 Mod{5}(0)
@@ -271,7 +271,7 @@ Mod{5}(2)
 Mod{5}(3)
 Mod{5}(4)
 
-julia> collect(Mod{6})
+julia> collect(allmods(Mod{6}))
 6-element Vector{Mod{6}}:
  Mod{6}(0)
  Mod{6}(1)
@@ -280,7 +280,7 @@ julia> collect(Mod{6})
  Mod{6}(4)
  Mod{6}(5)
 
-julia> [k*k for k ∈ Mod{7}]
+julia> [k*k for k ∈ allmods(Mod{7})]
 7-element Vector{Mod{7}}:
  Mod{7}(0)
  Mod{7}(1)
@@ -290,13 +290,13 @@ julia> [k*k for k ∈ Mod{7}]
  Mod{7}(4)
  Mod{7}(1)
 
-julia> prod(k for k in Mod{5} if k!=0) == -1  # Wilson's theorem
+julia> prod(k for k in allmods(Mod{5}) if k!=0) == -1  # Wilson's theorem
 true
 ```
 
-One can also use `GaussMod` as an iterator:
+Likewise `allmods(GaussMod{m})` iterates over all Gaussian residues:
 ```
-julia> for z in GaussMod{3}
+julia> for z in allmods(GaussMod{3})
        println(z)
        end
 GaussMod{3}(0 + 0im)
diff --git a/src/Mods.jl b/src/Mods.jl
index 78a1b0a..1c8f884 100644
--- a/src/Mods.jl
+++ b/src/Mods.jl
@@ -6,6 +6,7 @@ import Base: hash, iszero, isone, mod, abs, conj, real, imag
 
 export GaussMod, Mod, AbstractMod
 export modulus, value, is_invertible
+export allmods
 
 QZ = Union{<:Rational,<:Integer}
 CZ = Union{Complex,<:Integer}
diff --git a/src/iterate.jl b/src/iterate.jl
index a14b5f4..9e0a0d2 100644
--- a/src/iterate.jl
+++ b/src/iterate.jl
@@ -1,36 +1,35 @@
-function Base.iterate(::Type{Mod{m}}) where {m}
-    return Mod{m}(0), 0
+"""
+    allmods(Mod{m})
+    allmods(GaussMod{m})
+
+An iterator over all elements of `Mod{m}` (the residues `0, 1, …, m-1`) or of
+`GaussMod{m}` (`a + b*im` for `a, b` in `0:m-1`, `b` varying fastest).
+
+```
+julia> collect(allmods(Mod{4}))
+4-element Vector{Mod{4}}:
+ Mod{4}(0)
+ Mod{4}(1)
+ Mod{4}(2)
+ Mod{4}(3)
+```
+"""
+allmods(::Type{T}) where {T<:AbstractMod} = AllMods{T}()
+
+struct AllMods{T<:AbstractMod} end
+
+Base.eltype(::Type{AllMods{T}}) where {T} = T
+Base.IteratorSize(::Type{<:AllMods}) = Base.HasLength()
+Base.length(::AllMods{Mod{m}}) where {m} = m
+Base.length(::AllMods{GaussMod{m}}) where {m} = m * m
+
+function Base.iterate(::AllMods{Mod{m}}, s = 0) where {m}
+    s < m || return nothing
+    return Mod{m}(s), s + 1
 end
 
-function Base.iterate(::Type{Mod{m}}, s) where {m}
-    if s == m - 1
-        return nothing
-    end
-    s += 1
-    return Mod{m}(s), s
-end
-
-Base.IteratorSize(::Type{Mod{m}}) where {m} = Base.HasLength()
-
-Base.length(::Type{Mod{m}}) where {m} = m
-
-
-
-
-function Base.iterate(::Type{GaussMod{m}}) where {m}
-    return GaussMod{m}(0), 0
-end
-
-function Base.iterate(::Type{GaussMod{m}}, s) where {m}
-    if s == m * m - 1
-        return nothing
-    end
-    s += 1
+function Base.iterate(::AllMods{GaussMod{m}}, s = 0) where {m}
+    s < m * m || return nothing
     a, b = divrem(s, m)
-    return GaussMod{m}(a + b * im), s
+    return GaussMod{m}(a + b * im), s + 1
 end
-
-
-Base.IteratorSize(::Type{GaussMod{m}}) where {m} = Base.HasLength()
-
-Base.length(::Type{GaussMod{m}}) where {m} = m * m
diff --git a/test/runtests.jl b/test/runtests.jl
index 737d94b..a9b9977 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -119,12 +119,14 @@ end
 
 
 @testset "Iteration" begin
-    @test sum(k for k in Mod{7}) == zero(Mod{7})
-    @test collect(Mod{7}) == Mod{7}.(0:6)
-    @test prod(k for k in Mod{7} if k != 0) == -1  # Wilson's theorem
+    @test sum(k for k in allmods(Mod{7})) == zero(Mod{7})
+    @test collect(allmods(Mod{7})) == Mod{7}.(0:6)
+    @test prod(k for k in allmods(Mod{7}) if k != 0) == -1  # Wilson's theorem
 
     # complex interation
-    @test sum(GaussMod{7}) == 0
+    @test sum(allmods(GaussMod{7})) == 0
+    @test length(allmods(GaussMod{7})) == 49
+    @test eltype(allmods(Mod{7})) == Mod{7}
 end