diff --git a/.github/workflows/Format.yml b/.github/workflows/Format.yml
index be66417b..1d7384a6 100644
--- a/.github/workflows/Format.yml
+++ b/.github/workflows/Format.yml
@@ -8,5 +8,3 @@ jobs:
     steps:
     - uses: actions/checkout@master
     - uses: julia-actions/julia-format@master
-      with:
-        args: -v .
diff --git a/CHANGELOG.md b/CHANGELOG.md
index a9788a21..2e7d8a00 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,7 @@
+# v0.17
+
+The homology algorithm can now find infinite intervals and is a bit more efficient.
+
 # v0.16.9
 
 Replace LightGraps with Graphs.
diff --git a/Project.toml b/Project.toml
index 86605ebe..5ec92f5d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Ripserer"
 uuid = "aa79e827-bd0b-42a8-9f10-2b302677a641"
 authors = ["mtsch <matijacufar@gmail.com>"]
-version = "0.16.12"
+version = "0.17.0"
 
 [deps]
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
diff --git a/src/Ripserer.jl b/src/Ripserer.jl
index 93d8d752..6b17f919 100644
--- a/src/Ripserer.jl
+++ b/src/Ripserer.jl
@@ -70,7 +70,10 @@ include("base/simplexrecipes.jl")
 include("computation/utils.jl")
 include("computation/zerodimensional.jl")
 include("computation/reducedmatrix.jl")
+include("computation/persistence.jl")
 include("computation/cohomology.jl")
+include("computation/homology.jl")
+include("computation/involuted.jl")
 include("computation/ripserer.jl")
 
 include("filtrations/utils.jl")
diff --git a/src/base/abstractfiltration.jl b/src/base/abstractfiltration.jl
index 07257617..9e78ca2d 100644
--- a/src/base/abstractfiltration.jl
+++ b/src/base/abstractfiltration.jl
@@ -162,10 +162,7 @@ julia> vertices(Rips([0 1 1; 1 0 1; 1 1 0]))
 Base.OneTo(3)
 
 julia> vertices(Cubical([0 1 1; 1 0 1; 1 1 0]))
-3×3 CartesianIndices{2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}:
- CartesianIndex(1, 1)  CartesianIndex(1, 2)  CartesianIndex(1, 3)
- CartesianIndex(2, 1)  CartesianIndex(2, 2)  CartesianIndex(2, 3)
- CartesianIndex(3, 1)  CartesianIndex(3, 2)  CartesianIndex(3, 3)
+CartesianIndices((3, 3))
 
 ```
 """
diff --git a/src/computation/cohomology.jl b/src/computation/cohomology.jl
index 08ced0bf..821eb931 100644
--- a/src/computation/cohomology.jl
+++ b/src/computation/cohomology.jl
@@ -9,8 +9,8 @@ struct CoboundaryMatrix{
 }
     filtration::F
     reduced::R
-    buffer::B
-    chain::C
+    buffer::B     # stores the columns that were added to current chain
+    chain::C      # current column as it's being reduced
     columns_to_reduce::Vector{S}
     columns_to_skip::Vector{S}
 end
@@ -44,6 +44,8 @@ is_implicit(::CoboundaryMatrix{I}) where {I} = I
 is_cohomology(::CoboundaryMatrix) = true
 ordering(::CoboundaryMatrix) = Base.Order.Forward
 
+attempt_early_stop!(::CoboundaryMatrix, _, _) = false
+
 function coboundary(matrix::CoboundaryMatrix, simplex::AbstractCell)
     return coboundary(matrix.filtration, simplex)
 end
@@ -80,362 +82,22 @@ function next_matrix(matrix::CoboundaryMatrix{I}, verbose) where {I}
     )
 end
 
-"""
-    BoundaryMatrix{I}
-
-This `struct` is used to compute homology. The `I` parameter sets whether the implicit
-algoritm is used or not.
-"""
-struct BoundaryMatrix{
-    I,T<:Number,F,S<:AbstractCell,R<:ReducedMatrix,B<:Chain{T},C<:Chain{T}
-}
-    filtration::F
-    reduced::R
-    buffer::B
-    chain::C
-    columns_to_reduce::Vector{S}
-end
-
-function BoundaryMatrix{I}(::Type{T}, filtration, columns_to_reduce) where {I,T}
-    if eltype(columns_to_reduce) === Any
-        S = typeof(first(columns_to_reduce))
-    else
-        S = eltype(columns_to_reduce)
-    end
-    C = simplex_type(filtration, dim(S) - 1)
-
-    columns = S[]
-    foreach(columns_to_reduce) do c
-        push!(columns, abs(c))
-    end
-
-    if !I
-        reduced = ReducedMatrix{C,T,C}()
-        buffer = Chain{T,C}()
-    else
-        reduced = ReducedMatrix{C,T,S}()
-        buffer = Chain{T,S}()
-    end
-    sizehint!(reduced, length(columns))
-    chain = Chain{T,C}()
-
-    return BoundaryMatrix{
-        I,T,typeof(filtration),S,typeof(reduced),typeof(buffer),typeof(chain)
-    }(
-        filtration, reduced, buffer, chain, columns
-    )
-end
-
-field_type(::BoundaryMatrix{<:Any,T}) where {T} = T
-dim(::BoundaryMatrix{<:Any,<:Any,<:Any,S}) where {S} = dim(S) - 1
-ordering(::BoundaryMatrix) = Base.Order.Reverse
-
-is_implicit(::BoundaryMatrix{I}) where {I} = I
-is_cohomology(::BoundaryMatrix) = false
-
-# The naming here is not ideal...
-function coboundary(matrix::BoundaryMatrix, simplex::AbstractCell)
-    return boundary(matrix.filtration, simplex)
-end
-
-"""
-    initialize_coboundary!(matrix, column)
-
-Initialize the column indexed by `column` by adding its (co)boundary to `matrix.chain`. This
-is where the emergent pairs optimization gets triggered for implicit versions of the
-algorithm.
-
-"""
-function initialize_coboundary!(matrix, column)
-    empty!(matrix.chain)
-    # Emergent pairs: we are looking for pairs of simplices (σ, τ) where σ is the youngest
-    # facet of τ and τ is the oldest cofacet of σ. These pairs give birth to persistence
-    # intervals with zero length and can be skipped.
-
-    # This implementation of this optimization only works if (co)boundary simplices are
-    # returned in the correct order and if the birth times of σ and τ are the same.
-    emergent_check = emergent_pairs(matrix.filtration) && is_implicit(matrix)
-    for cofacet in coboundary(matrix, column)
-        if emergent_check && birth(cofacet) == birth(column)
-            emergent_check = false
-            if !haskey(matrix.reduced, cofacet)
-                return ChainElement{typeof(cofacet),field_type(matrix)}(cofacet)
-            end
-        end
-        push!(matrix.chain, cofacet)
-    end
-    if isempty(matrix.chain)
-        return nothing
-    else
-        heapify!(matrix.chain, ordering(matrix))
-        return heappop!(matrix.chain, ordering(matrix))
-    end
-end
-
-"""
-    add!(matrix, column, pivot)
-
-Add already reduced column indexed by `column` multiplied by `-coefficient(pivot)` to
-`matrix.chain`.
-
-"""
-function add!(matrix, column, pivot)
-    return add!(Val(is_implicit(matrix)), matrix, column, pivot)
-end
-# Implicit version
-function add!(::Val{true}, matrix, column, pivot)
-    factor = -coefficient(pivot)
-    for element in column
-        for cofacet in coboundary(matrix, simplex(element))
-            simplex(pivot) == cofacet && continue
-            heappush!(
-                matrix.chain, (cofacet, coefficient(element) * factor), ordering(matrix)
-            )
-        end
-    end
-    return append!(matrix.buffer, (s, c * factor) for (s, c) in column)
-end
-# Explicit version
-function add!(::Val{false}, matrix, column, pivot)
-    factor = -coefficient(pivot)
-    for element in column
-        # The pivot is not stored in the column, so there is no need to check for it.
-        heappush!(matrix.chain, element * factor, ordering(matrix))
-    end
-end
-
-"""
-    finalize!(matrix, column, pivot)
-
-After reduction is done, finalize the column by adding it to the reduced matrix.
-"""
-function finalize!(matrix, column, pivot)
-    return finalize!(Val(is_implicit(matrix)), matrix, column, pivot)
-end
-# Implicit version
-function finalize!(::Val{true}, matrix, column, pivot)
-    push!(matrix.buffer, column)
-    return matrix.reduced[simplex(pivot)] = clean!(
-        matrix.buffer, ordering(matrix), inv(coefficient(pivot))
-    )
-end
-# Explicit version
-function finalize!(::Val{false}, matrix, column, pivot)
-    return matrix.reduced[simplex(pivot)] = clean!(
-        matrix.chain, ordering(matrix), inv(coefficient(pivot))
-    )
-end
-
-"""
-    reduce_column!(matrix, column_to_reduce)
-
-Reduce `column_to_reduce` by repeatedly adding other columns to it. Once nothing more can be
-added, `finalize!` the column.
-
-"""
-function reduce_column!(matrix, column_to_reduce)
-    empty!(matrix.buffer)
-    pivot = initialize_coboundary!(matrix, column_to_reduce)
-    while !isnothing(pivot)
-        column = matrix.reduced[pivot]
-        isempty(column) && break
-
-        add!(matrix, column, pivot)
-        pivot = heappop!(matrix.chain, ordering(matrix))
-    end
-    if !isnothing(pivot)
-        finalize!(matrix, column_to_reduce, pivot)
-    end
-
-    return pivot
-end
-
-"""
-    collect_cocycle!(matrix, pivot, reps)
-
-Collect the representative (co)cycle.
-"""
-function collect_cocycle!(matrix, column, pivot)
-    if is_cohomology(matrix)
-        if !is_implicit(matrix)
-            error("representative cocycles for explicit cohomology not supported")
-        elseif isnothing(pivot)
-            push!(matrix.buffer, column)
-            return copy(clean!(matrix.buffer, ordering(matrix)))
-        else
-            return Chain(matrix.reduced[pivot])
-        end
-    else
-        if is_implicit(matrix)
-            heappush!(matrix.chain, pivot, ordering(matrix))
-            return heapmove!(matrix.chain, ordering(matrix))
-        else
-            rep = Chain(matrix.reduced[pivot])
-            push!(rep, pivot)
-            return clean!(rep, ordering(matrix))
-        end
-    end
-end
-
-"""
-    interval(matrix, column, pivot, cutoff, reps)
-
-Construct a persistence interval.
-"""
-function interval(matrix, column, pivot, cutoff, reps)
-    if is_cohomology(matrix)
-        birth_simplex = column
-        death_simplex = isnothing(pivot) ? nothing : simplex(pivot)
-    elseif isnothing(pivot)
-        # In homology, birth simplex is nothing when column is fully reduced.
-        return nothing
-    else
-        birth_simplex, death_simplex = simplex(pivot), column
-    end
-    birth_time = Float64(birth(birth_simplex))
-    death_time = isnothing(death_simplex) ? Inf : Float64(birth(death_simplex))
-    if death_time - birth_time > cutoff
-        if reps
-            rep = (; representative=collect_cocycle!(matrix, column, pivot))
-        else
-            rep = NamedTuple()
-        end
-        meta = (; birth_simplex=birth_simplex, death_simplex=death_simplex, rep...)
-        return PersistenceInterval(birth_time, death_time, meta)
-    else
-        return nothing
-    end
-end
-
-"""
-    handle_apparent_pairs!(matrix, intervals, cutoff, verbose, reps)
-
-Handle apparent pair precomputation, if defined for `matrix.filtration`. Only does anything
-for implicit cohomology. Resulting intervals (if any) are `push!`ed to `intervals`.
-"""
-function handle_apparent_pairs!(matrix, intervals, cutoff, verbose, reps)
-    coho = Val(is_cohomology(matrix))
-    impl = Val(is_implicit(matrix))
-    return handle_apparent_pairs!(coho, impl, matrix, intervals, cutoff, verbose, reps)
-end
-# Implicit cohomology version
-function handle_apparent_pairs!(
-    ::Val{true}, ::Val{true}, matrix, intervals, cutoff, verbose, reps
+function _ripserer(
+    ::Val{:cohomology}, filtration, cutoff, verbose, field, dim_max, reps, implicit
 )
-    columns, apparent = find_apparent_pairs(
-        matrix.filtration, matrix.columns_to_reduce, verbose
+    result = PersistenceDiagram[]
+    zeroth, to_reduce, to_skip = zeroth_intervals(
+        filtration, cutoff, verbose, field, _reps(reps, 0)
     )
-    bulk_insert!(matrix.reduced, apparent)
-    for (σ, τ) in apparent
-        τ_elem = ChainElement{typeof(τ),field_type(matrix)}(τ)
-        int = interval(matrix, σ, τ_elem, cutoff, reps)
-        !isnothing(int) && push!(intervals, int)
-    end
-    return columns
-end
-# Other versions don't support this
-function handle_apparent_pairs!(::Val, ::Val, matrix, _, _, _, _)
-    return matrix.columns_to_reduce
-end
-
-"""
-    compute_intervals!(matrix, cutoff, verbose, reps)
-
-Compute all intervals by fully reducing `matrix`.
-"""
-function compute_intervals!(matrix, cutoff, verbose, reps)
-    ###
-    ### Set up output.
-    ###
-    intervals = PersistenceInterval[]
-
-    ###
-    ### Apparent pair stuff.
-    ###
-    columns = handle_apparent_pairs!(matrix, intervals, cutoff, verbose, reps)
-
-    ###
-    ### Interval computation.
-    ###
-    @prog_print(
-        verbose,
-        fmt_number(length(columns)),
-        " ",
-        simplex_name(eltype(columns)),
-        " to reduce.",
-    )
-    # One-dimensional columns in cohomology are already sorted.
-    if !is_cohomology(matrix) || dim(matrix) > 1
-        sort_t = time_ns()
-        sort!(columns; rev=is_cohomology(matrix))
-        elapsed = round((time_ns() - sort_t) / 1e9; digits=3)
-        @prog_println verbose " Sorted in " elapsed "s."
-    else
-        @prog_println verbose
-    end
-
-    if verbose
-        progbar = Progress(length(columns); desc="Computing $(dim(matrix))d intervals... ")
-    end
-    for column in columns
-        pivot = reduce_column!(matrix, column)
-        int = interval(matrix, column, pivot, cutoff, reps)
-        !isnothing(int) && push!(intervals, int)
-
-        verbose && next!(progbar; showvalues=((:intervals, length(intervals)),))
-    end
-
-    diagram = PersistenceDiagram(
-        sort!(intervals; by=persistence);
-        threshold=Float64(threshold(matrix.filtration)),
-        dim=dim(matrix),
-        field=field_type(matrix),
-        filtration=matrix.filtration,
-    )
-    return postprocess_diagram(matrix.filtration, diagram)
-end
-
-"""
-    compute_death_simplices!(matrix::CoboundaryMatrix{true}, verbose, cutoff)
-
-Fully reduce `matrix`, but only compute (homological) death simplices. Return all death
-simplices up to the last that produces an interval with persistence greater than `cutoff`.
-
-Used for assisted homology.
-"""
-function compute_death_simplices!(matrix::CoboundaryMatrix{true}, verbose, cutoff)
-    columns, apparent = find_apparent_pairs(
-        matrix.filtration, matrix.columns_to_reduce, verbose
-    )
-    bulk_insert!(matrix.reduced, apparent)
-    deaths = simplex_type(matrix.filtration, dim(matrix) + 1)[]
-    inf_births = simplex_type(matrix.filtration, dim(matrix))[]
-    if isempty(columns)
-        return deaths, inf_births
-    else
-        dim(matrix) > 1 && sort!(columns; rev=true)
-        thresh = typemin(birth(first(columns)))
-        for pair in apparent
-            if birth(pair[2]) - birth(pair[1]) > cutoff
-                thresh = max(thresh, birth(pair[2]))
-            end
-            push!(deaths, pair[2])
-        end
-        if verbose
-            progbar = Progress(length(columns); desc="Precomputing columns...   ")
-        end
-        for column in columns
-            pivot = reduce_column!(matrix, column)
-            if !isnothing(pivot)
-                if birth(pivot) - birth(column) > cutoff
-                    thresh = max(thresh, birth(pivot))
-                end
-                push!(deaths, simplex(pivot))
-            else
-                push!(inf_births, column)
+    push!(result, zeroth)
+    if dim_max > 0
+        matrix = CoboundaryMatrix{implicit}(field, filtration, to_reduce, to_skip)
+        for dim in 1:dim_max
+            push!(result, compute_intervals!(matrix, cutoff, verbose, _reps(reps, dim)))
+            if dim < dim_max
+                matrix = next_matrix(matrix, verbose)
             end
-            verbose && next!(progbar; showvalues=((:simplices, length(deaths)),))
         end
-        return filter!(x -> birth(x) ≤ thresh, deaths), inf_births
     end
+    return result
 end
diff --git a/src/computation/homology.jl b/src/computation/homology.jl
new file mode 100644
index 00000000..ada838b3
--- /dev/null
+++ b/src/computation/homology.jl
@@ -0,0 +1,147 @@
+"""
+Write a blurb here.
+"""
+struct Homology
+    implicit::Bool
+
+    Homology(; implicit=true) = new(implicit)
+end
+
+"""
+    BoundaryMatrix{I}
+
+This `struct` is used to compute homology. The `I` parameter sets whether the implicit
+algoritm is used or not.
+"""
+struct BoundaryMatrix{
+    I,T<:Number,F,S1<:AbstractCell,S2<:AbstractCell,R<:ReducedMatrix,B<:Chain{T},C<:Chain{T}
+}
+    filtration::F
+    reduced::R
+    buffer::B
+    chain::C
+    birth_candidates::Vector{S1} # used to find infinite intervals
+    columns_to_reduce::Vector{S2}
+    zeroed::Set{S2}
+    infinite_intervals::Bool
+end
+
+function BoundaryMatrix{I}(
+    ::Type{T}, filtration, birth_candidates, columns_to_reduce; infinite_intervals=true
+) where {I,T}
+    S2 = eltype(columns_to_reduce)
+    S1 = simplex_type(filtration, dim(S2) - 1)
+
+    if !I
+        reduced = ReducedMatrix{S1,T,S1}()
+        buffer = Chain{T,S1}()
+    else
+        reduced = ReducedMatrix{S1,T,S2}()
+        buffer = Chain{T,S2}()
+    end
+    sizehint!(reduced, length(columns_to_reduce))
+    chain = Chain{T,S1}()
+
+    return BoundaryMatrix{
+        I,T,typeof(filtration),S1,S2,typeof(reduced),typeof(buffer),typeof(chain)
+    }(
+        filtration,
+        reduced,
+        buffer,
+        chain,
+        birth_candidates,
+        columns_to_reduce,
+        Set{S1}(),
+        infinite_intervals,
+    )
+end
+
+field_type(::BoundaryMatrix{<:Any,T}) where {T} = T
+dim(::BoundaryMatrix{<:Any,<:Any,<:Any,S}) where {S} = dim(S)
+ordering(::BoundaryMatrix) = Base.Order.Reverse
+
+is_implicit(::BoundaryMatrix{I}) where {I} = I
+is_cohomology(::BoundaryMatrix) = false
+
+function attempt_early_stop!(matrix::BoundaryMatrix, i, columns)
+    if length(matrix.reduced.column_index) ≥ length(matrix.birth_candidates)
+        # At this point, all potential births have been found. The rest of the columns
+        # should be marked as zeroed.
+        for j in i:length(columns)
+            push!(matrix.zeroed, columns[j])
+        end
+        return true
+    else
+        return false
+    end
+end
+
+# The naming here is not ideal...
+function coboundary(matrix::BoundaryMatrix, simplex::AbstractCell)
+    return boundary(matrix.filtration, simplex)
+end
+
+function mark_zero_column!(matrix::BoundaryMatrix, column_index)
+    if matrix.infinite_intervals
+        push!(matrix.zeroed, column_index)
+    end
+    return nothing
+end
+
+function append_infinite_intervals!(intervals, matrix::BoundaryMatrix)
+    if matrix.infinite_intervals && length(matrix.birth_candidates) ≠ length(intervals)
+        for simplex in matrix.birth_candidates
+            if !haskey(matrix.reduced, simplex)
+                push!(
+                    intervals,
+                    PersistenceInterval(
+                        birth(simplex),
+                        Inf,
+                        (; birth_simplex=simplex, death_simplex=nothing),
+                    ),
+                )
+            end
+        end
+    end
+    return intervals
+end
+
+function next_matrix(matrix::BoundaryMatrix{I}) where {I}
+    birth_candidates = filter(matrix.columns_to_reduce) do sx
+        sx in matrix.zeroed
+    end
+    columns = simplex_type(matrix.filtration, dim(matrix) + 2)[]
+    for col in columns_to_reduce(matrix.filtration, matrix.columns_to_reduce)
+        push!(columns, abs(col))
+    end
+
+    return BoundaryMatrix{I}(
+        field_type(matrix), matrix.filtration, birth_candidates, columns
+    )
+end
+
+function _ripserer(
+    ::Val{:homology}, filtration, cutoff, verbose, field, dim_max, reps, implicit
+)
+    result = PersistenceDiagram[]
+    zeroth, to_reduce, to_skip = zeroth_intervals(
+        filtration, cutoff, verbose, field, _reps(reps, 0)
+    )
+    push!(result, zeroth)
+    if dim_max > 0
+        birth_candidates = to_reduce
+        columns = simplex_type(filtration, 2)[]
+        for col in columns_to_reduce(filtration, Iterators.flatten((to_reduce, to_skip)))
+            push!(columns, abs(col))
+        end
+        matrix = BoundaryMatrix{implicit}(field, filtration, birth_candidates, columns)
+        for dim in 1:dim_max
+            push!(result, compute_intervals!(matrix, cutoff, verbose, _reps(reps, dim)))
+
+            if dim < dim_max
+                matrix = next_matrix(matrix)
+            end
+        end
+    end
+    return result
+end
diff --git a/src/computation/involuted.jl b/src/computation/involuted.jl
new file mode 100644
index 00000000..b8364c36
--- /dev/null
+++ b/src/computation/involuted.jl
@@ -0,0 +1,76 @@
+"""
+    compute_death_simplices!(matrix::CoboundaryMatrix{true}, verbose, cutoff)
+
+Fully reduce `matrix`, but only compute (homological) death simplices. Return all death
+simplices up to the last that produces an interval with persistence greater than `cutoff`.
+
+Used for involuted homology.
+"""
+function compute_death_simplices!(matrix::CoboundaryMatrix{true}, verbose, cutoff)
+    columns, apparent = find_apparent_pairs(
+        matrix.filtration, matrix.columns_to_reduce, verbose
+    )
+    bulk_insert!(matrix.reduced, apparent)
+    deaths = simplex_type(matrix.filtration, dim(matrix) + 1)[]
+    inf_births = simplex_type(matrix.filtration, dim(matrix))[]
+    if isempty(columns)
+        return deaths, inf_births
+    else
+        dim(matrix) > 1 && sort!(columns; rev=true)
+        thresh = typemin(birth(first(columns)))
+        for pair in apparent
+            if birth(pair[2]) - birth(pair[1]) > cutoff
+                thresh = max(thresh, birth(pair[2]))
+            end
+            push!(deaths, pair[2])
+        end
+        if verbose
+            progbar = Progress(length(columns); desc="Precomputing columns...   ")
+        end
+        for column in columns
+            pivot = reduce_column!(matrix, column)
+            if !isnothing(pivot)
+                if birth(pivot) - birth(column) > cutoff
+                    thresh = max(thresh, birth(pivot))
+                end
+                push!(deaths, simplex(pivot))
+            else
+                push!(inf_births, column)
+            end
+            verbose && next!(progbar; showvalues=((:simplices, length(deaths)),))
+        end
+        return filter!(x -> birth(x) ≤ thresh, deaths), inf_births
+    end
+end
+
+function _ripserer(
+    ::Val{:involuted}, filtration, cutoff, verbose, field, dim_max, reps, implicit
+)
+    result = PersistenceDiagram[]
+    zeroth, to_reduce, to_skip = zeroth_intervals(
+        filtration, cutoff, verbose, field, _reps(reps, 0)
+    )
+    push!(result, zeroth)
+    if dim_max > 0
+        comatrix = CoboundaryMatrix{true}(field, filtration, to_reduce, to_skip)
+        for dim in 1:dim_max
+            columns, inf_births = compute_death_simplices!(comatrix, verbose, cutoff)
+            birth_candidates = comatrix.columns_to_reduce
+            matrix = BoundaryMatrix{implicit}(
+                field, filtration, birth_candidates, columns; infinite_intervals=false
+            )
+            diagram = compute_intervals!(matrix, cutoff, verbose, _reps(reps, dim))
+            for birth_simplex in inf_births
+                push!(
+                    diagram.intervals,
+                    interval(comatrix, birth_simplex, nothing, 0, _reps(reps, dim)),
+                )
+            end
+            push!(result, diagram)
+            if dim < dim_max
+                comatrix = next_matrix(comatrix, verbose)
+            end
+        end
+    end
+    return result
+end
diff --git a/src/computation/persistence.jl b/src/computation/persistence.jl
new file mode 100644
index 00000000..b58f486d
--- /dev/null
+++ b/src/computation/persistence.jl
@@ -0,0 +1,270 @@
+"""
+    initialize_column!(matrix, column_index)
+
+Initialize the column indexed by `column_index` by adding its (co)boundary to
+`matrix.chain`. This is where the emergent pairs optimization gets triggered for implicit
+versions of the algorithm.
+
+Return the pivot as a [`ChainElement`](@ref).
+"""
+function initialize_column!(matrix, column_index)
+    empty!(matrix.chain)
+    # Emergent pairs: we are looking for pairs of simplices (σ, τ) where σ is the youngest
+    # facet of τ and τ is the oldest cofacet of σ. These pairs give birth to persistence
+    # intervals with zero length and can be skipped.
+
+    # This implementation of this optimization only works if (co)boundary simplices are
+    # returned in the correct order and if the birth times of σ and τ are the same.
+    emergent_check = emergent_pairs(matrix.filtration) && is_implicit(matrix)
+    for cofacet in coboundary(matrix, column_index)
+        if emergent_check && birth(cofacet) == birth(column_index)
+            emergent_check = false
+            if !haskey(matrix.reduced, cofacet)
+                return ChainElement{typeof(cofacet),field_type(matrix)}(cofacet)
+            end
+        end
+        push!(matrix.chain, cofacet)
+    end
+    if isempty(matrix.chain)
+        return nothing
+    else
+        heapify!(matrix.chain, ordering(matrix))
+        return heappop!(matrix.chain, ordering(matrix))
+    end
+end
+
+"""
+    add!(matrix, column, pivot)
+
+Add already `column` multiplied by `-coefficient(pivot)` to `matrix.chain`.
+
+"""
+function add!(matrix, column, pivot)
+    add!(Val(is_implicit(matrix)), matrix, column, pivot)
+    return nothing
+end
+# Implicit version
+function add!(::Val{true}, matrix, column, pivot)
+    factor = -coefficient(pivot)
+    for element in column
+        for cofacet in coboundary(matrix, simplex(element))
+            simplex(pivot) == cofacet && continue
+            heappush!(
+                matrix.chain, (cofacet, coefficient(element) * factor), ordering(matrix)
+            )
+        end
+    end
+    append!(matrix.buffer, (s, c * factor) for (s, c) in column)
+    return nothing
+end
+# Explicit version
+function add!(::Val{false}, matrix, column, pivot)
+    factor = -coefficient(pivot)
+    for element in column
+        # The pivot is not stored in the column, so there is no need to check for it.
+        heappush!(matrix.chain, element * factor, ordering(matrix))
+    end
+    return nothing
+end
+
+"""
+    finalize!(matrix, column_index, pivot)
+
+After reduction is done, finalize the current column being reduced by adding it to the
+reduced matrix.
+"""
+function finalize!(matrix, column_index, pivot)
+    finalize!(Val(is_implicit(matrix)), matrix, column_index, pivot)
+    return nothing
+end
+# Implicit version
+function finalize!(::Val{true}, matrix, column_index, pivot)
+    push!(matrix.buffer, column_index)
+    matrix.reduced[simplex(pivot)] = clean!(
+        matrix.buffer, ordering(matrix), inv(coefficient(pivot))
+    )
+    return nothing
+end
+# Explicit version
+function finalize!(::Val{false}, matrix, _, pivot)
+    matrix.reduced[simplex(pivot)] = clean!(
+        matrix.chain, ordering(matrix), inv(coefficient(pivot))
+    )
+    return nothing
+end
+
+"""
+    mark_zero_column!(matrix, column_index)
+
+This function is called on a (co)boundary matrix when a column is completely reduced.
+"""
+function mark_zero_column!(_, _)
+    return nothing
+end
+
+"""
+    reduce_column!(matrix, column_index)
+
+Reduce `column_index` by repeatedly adding other columns to it. Once nothing more can be
+added, `finalize!` the column.
+"""
+function reduce_column!(matrix, column_index)
+    empty!(matrix.buffer)
+    pivot = initialize_column!(matrix, column_index)
+    while !isnothing(pivot)
+        column = matrix.reduced[pivot]
+        isempty(column) && break
+
+        add!(matrix, column, pivot)
+        pivot = heappop!(matrix.chain, ordering(matrix))
+    end
+    if !isnothing(pivot)
+        finalize!(matrix, column_index, pivot)
+    else
+        mark_zero_column!(matrix, column_index)
+    end
+
+    return pivot
+end
+
+"""
+    collect_cocycle!(matrix, column, pivot)
+
+Collect the representative (co)cycle.
+"""
+function collect_cocycle!(matrix, column, pivot)
+    if is_cohomology(matrix)
+        if !is_implicit(matrix)
+            error("representative cocycles for explicit cohomology not supported")
+        elseif isnothing(pivot)
+            push!(matrix.buffer, column)
+            return copy(clean!(matrix.buffer, ordering(matrix)))
+        else
+            return Chain(matrix.reduced[pivot])
+        end
+    else
+        if is_implicit(matrix)
+            heappush!(matrix.chain, pivot, ordering(matrix))
+            return heapmove!(matrix.chain, ordering(matrix))
+        else
+            rep = Chain(matrix.reduced[pivot])
+            push!(rep, pivot)
+            return clean!(rep, ordering(matrix))
+        end
+    end
+end
+
+"""
+    interval(matrix, column, pivot, cutoff, reps)
+
+Construct a persistence interval.
+"""
+function interval(matrix, column, pivot, cutoff, reps)
+    if is_cohomology(matrix)
+        birth_simplex = column
+        death_simplex = isnothing(pivot) ? nothing : simplex(pivot)
+    elseif isnothing(pivot)
+        # In homology, birth simplex is nothing when column is fully reduced.
+        return nothing
+    else
+        birth_simplex, death_simplex = simplex(pivot), column
+    end
+    birth_time = Float64(birth(birth_simplex))
+    death_time = isnothing(death_simplex) ? Inf : Float64(birth(death_simplex))
+    if death_time - birth_time > cutoff
+        if reps
+            rep = (; representative=collect_cocycle!(matrix, column, pivot))
+        else
+            rep = NamedTuple()
+        end
+        meta = (; birth_simplex=birth_simplex, death_simplex=death_simplex, rep...)
+        return PersistenceInterval(birth_time, death_time, meta)
+    else
+        return nothing
+    end
+end
+
+"""
+    handle_apparent_pairs!(matrix, intervals, cutoff, verbose, reps)
+
+Handle apparent pair precomputation, if defined for `matrix.filtration`. Only does anything
+for implicit cohomology. Resulting intervals (if any) are `push!`ed to `intervals`.
+"""
+function handle_apparent_pairs!(matrix, intervals, cutoff, verbose, reps)
+    coho = Val(is_cohomology(matrix))
+    impl = Val(is_implicit(matrix))
+    return handle_apparent_pairs!(coho, impl, matrix, intervals, cutoff, verbose, reps)
+end
+# Implicit cohomology version
+function handle_apparent_pairs!(
+    ::Val{true}, ::Val{true}, matrix, intervals, cutoff, verbose, reps
+)
+    columns, apparent = find_apparent_pairs(
+        matrix.filtration, matrix.columns_to_reduce, verbose
+    )
+    bulk_insert!(matrix.reduced, apparent)
+    for (σ, τ) in apparent
+        τ_elem = ChainElement{typeof(τ),field_type(matrix)}(τ)
+        int = interval(matrix, σ, τ_elem, cutoff, reps)
+        !isnothing(int) && push!(intervals, int)
+    end
+    return columns
+end
+# Other versions don't support this
+function handle_apparent_pairs!(::Val, ::Val, matrix, _, _, _, _)
+    return matrix.columns_to_reduce
+end
+
+"""
+    compute_intervals!(matrix, cutoff, verbose, reps)
+
+Compute all intervals by fully reducing `matrix`.
+"""
+function compute_intervals!(matrix, cutoff, verbose, reps)
+    intervals = PersistenceInterval[]
+
+    columns = handle_apparent_pairs!(matrix, intervals, cutoff, verbose, reps)
+
+    @prog_print(
+        verbose,
+        fmt_number(length(columns)),
+        " ",
+        simplex_name(eltype(columns)),
+        " to reduce.",
+    )
+    # One-dimensional columns in cohomology are already sorted.
+    if !is_cohomology(matrix) || dim(matrix) > 1
+        sort_t = time_ns()
+        sort!(columns; rev=is_cohomology(matrix))
+        elapsed = round((time_ns() - sort_t) / 1e9; digits=3)
+        @prog_println verbose " Sorted in " elapsed "s."
+    else
+        @prog_println verbose
+    end
+
+    if verbose
+        progbar = Progress(length(columns); desc="Computing $(dim(matrix))d intervals... ")
+    end
+    for (i, column) in enumerate(columns)
+        if attempt_early_stop!(matrix, i, columns)
+            break
+        end
+        pivot = reduce_column!(matrix, column)
+        int = interval(matrix, column, pivot, cutoff, reps)
+        !isnothing(int) && push!(intervals, int)
+
+        verbose && next!(progbar; showvalues=((:intervals, length(intervals)),))
+    end
+    if !is_cohomology(matrix)
+        append_infinite_intervals!(intervals, matrix)
+    end
+
+    diagram = PersistenceDiagram(
+        sort!(intervals; by=persistence);
+        threshold=Float64(threshold(matrix.filtration)),
+        dim=dim(matrix),
+        field=field_type(matrix),
+        filtration=matrix.filtration,
+    )
+    return postprocess_diagram(matrix.filtration, diagram)
+end
diff --git a/src/computation/ripserer.jl b/src/computation/ripserer.jl
index 546c8dcd..0e34a32e 100644
--- a/src/computation/ripserer.jl
+++ b/src/computation/ripserer.jl
@@ -78,7 +78,7 @@ diagram, and the last is the (`dim_max`)-dimensional diagram.
     with representative _co_cycles.
 
   - `:homology`: Significantly slower than `:cohomology`, but finds representative
-    cycles. Does not find infinite intervals beyond dimension 0.
+    cycles.
 
   - `:involuted`: Use cohomology result to compute representative cycles. Can be extremely
     efficient compared to `:homology`, especially with `Rips` filtrations. See [this
@@ -158,76 +158,6 @@ function _ripserer(
     return result
 end
 
-function _ripserer(
-    ::Val{:cohomology}, filtration, cutoff, verbose, field, dim_max, reps, implicit
-)
-    result = PersistenceDiagram[]
-    zeroth, to_reduce, to_skip = zeroth_intervals(
-        filtration, cutoff, verbose, field, _reps(reps, 0)
-    )
-    push!(result, zeroth)
-    if dim_max > 0
-        matrix = CoboundaryMatrix{implicit}(field, filtration, to_reduce, to_skip)
-        for dim in 1:dim_max
-            push!(result, compute_intervals!(matrix, cutoff, verbose, _reps(reps, dim)))
-            if dim < dim_max
-                matrix = next_matrix(matrix, verbose)
-            end
-        end
-    end
-    return result
-end
-
-function _ripserer(
-    ::Val{:homology}, filtration, cutoff, verbose, field, dim_max, reps, implicit
-)
-    result = PersistenceDiagram[]
-    zeroth, to_reduce, to_skip = zeroth_intervals(
-        filtration, cutoff, verbose, field, _reps(reps, 0)
-    )
-    push!(result, zeroth)
-    if dim_max > 0
-        simplices = columns_to_reduce(filtration, Iterators.flatten((to_reduce, to_skip)))
-        for dim in 1:dim_max
-            matrix = BoundaryMatrix{implicit}(field, filtration, simplices)
-            push!(result, compute_intervals!(matrix, cutoff, verbose, _reps(reps, dim)))
-            if dim < dim_max
-                simplices = columns_to_reduce(filtration, simplices)
-            end
-        end
-    end
-    return result
-end
-
-function _ripserer(
-    ::Val{:involuted}, filtration, cutoff, verbose, field, dim_max, reps, implicit
-)
-    result = PersistenceDiagram[]
-    zeroth, to_reduce, to_skip = zeroth_intervals(
-        filtration, cutoff, verbose, field, _reps(reps, 0)
-    )
-    push!(result, zeroth)
-    if dim_max > 0
-        comatrix = CoboundaryMatrix{true}(field, filtration, to_reduce, to_skip)
-        for dim in 1:dim_max
-            columns, inf_births = compute_death_simplices!(comatrix, verbose, cutoff)
-            matrix = BoundaryMatrix{implicit}(field, filtration, columns)
-            diagram = compute_intervals!(matrix, cutoff, verbose, _reps(reps, dim))
-            for birth_simplex in inf_births
-                push!(
-                    diagram.intervals,
-                    interval(comatrix, birth_simplex, nothing, 0, _reps(reps, dim)),
-                )
-            end
-            push!(result, diagram)
-            if dim < dim_max
-                comatrix = next_matrix(comatrix, verbose)
-            end
-        end
-    end
-    return result
-end
-
 function _ripserer(::Val{A}, args...) where {A}
     return throw(ArgumentError("unsupported alg=$A"))
 end
diff --git a/src/filtrations/rips.jl b/src/filtrations/rips.jl
index 45db713e..9944cae5 100644
--- a/src/filtrations/rips.jl
+++ b/src/filtrations/rips.jl
@@ -219,7 +219,13 @@ function Rips{I}(
     if issparse(adj)
         adj isa SparseMatrixCSC ||
             throw(ArgumentError("only SparseMatrixCSC sparse matrices supported"))
-        SparseArrays.fkeep!(adj, (_, _, v) -> v ≤ thresh)
+
+        # This is janky, but using the deprecated syntax prints depwarns in the docstrings.
+        if VERSION ≥ v"1.9.0-DEV"
+            SparseArrays.fkeep!((_, _, v) -> v ≤ thresh, adj)
+        else
+            SparseArrays.fkeep!(adj, (_, _, v) -> v ≤ thresh)
+        end
     end
     return Rips{I,T,typeof(adj)}(adj, thresh)
 end
diff --git a/test/doctests.jl b/test/doctests.jl
index 9c215520..c1146dd4 100644
--- a/test/doctests.jl
+++ b/test/doctests.jl
@@ -2,8 +2,8 @@ using Documenter
 using Distances
 using Ripserer
 
-if VERSION ≥ v"1.8-DEV" || VERSION < v"1.7-DEV" || !Sys.islinux()
-    @warn "Doctests only run on Linux and Julia 1.7"
+if VERSION ≥ v"1.11-DEV" || VERSION < v"1.10-DEV" || !Sys.islinux()
+    @warn "Doctests only run on Linux and Julia 1.10"
 else
     DocMeta.setdocmeta!(
         Ripserer, :DocTestSetup, :(using Ripserer; using Distances); recursive=true
diff --git a/test/filtrations/interfacetest.jl b/test/filtrations/interfacetest.jl
index 698e943d..2fbfadc7 100644
--- a/test/filtrations/interfacetest.jl
+++ b/test/filtrations/interfacetest.jl
@@ -56,8 +56,13 @@ function test_filtration(F, args...; test_verbose=true, flt_kwargs=(), kwargs...
                 @test coh_imp[i] == coh_exp[i]
                 @test hom_imp[i] == hom_exp[i] == hom_inv[i]
                 @test hom_inv[i] == coh_imp[i]
-                @test representative.(hom_imp[i]) == representative.(hom_exp[i])
-                @test representative.(hom_inv[i]) == representative.(hom_exp[i])
+
+                hom_imp_reps = [representative(int) for int in hom_imp[i] if isfinite(int)]
+                hom_exp_reps = [representative(int) for int in hom_exp[i] if isfinite(int)]
+                hom_inv_reps = [representative(int) for int in hom_inv[i] if isfinite(int)]
+
+                @test hom_imp_reps == hom_exp_reps
+                @test hom_imp_reps == hom_inv_reps
             end
         end
     end
diff --git a/test/filtrations/rips.jl b/test/filtrations/rips.jl
index f37e71a1..703faceb 100644
--- a/test/filtrations/rips.jl
+++ b/test/filtrations/rips.jl
@@ -324,12 +324,15 @@ end
                 sort!(vcat([(i + 1, i) for i in 1:17], [(18, 1)]))
         end
         @testset "Infinite intervals" begin
-            @test_broken ripserer(
+            @test ripserer(
                 Rips(cycle; threshold=2); alg=:homology, implicit=true
             )[2][1] == (1.0, Inf)
-            @test_broken ripserer(
+            @test ripserer(
                 Rips, cycle; alg=:homology, threshold=2, implicit=false
             )[2][1] == (1.0, Inf)
+            @test ripserer(
+                Rips, projective_plane; alg=:homology, threshold=1, dim_max=3
+            )[3][1] == (1.0, Inf)
             @test ripserer(cycle; alg=:involuted, threshold=2)[2][1] == (1.0, Inf)
         end
     end