diff --git a/src/Groebner.jl b/src/Groebner.jl
index 3803aa93..2ab1ce61 100644
--- a/src/Groebner.jl
+++ b/src/Groebner.jl
@@ -182,6 +182,7 @@ include("f4/matrix.jl")
 # Linear algebra backends
 include("f4/linalg/linalg.jl")
 include("f4/linalg/backend.jl")
+include("f4/linalg/backend_changematrix.jl")
 include("f4/linalg/backend_threaded.jl")
 include("f4/linalg/backend_randomized.jl")
 include("f4/linalg/backend_randomized_threaded.jl")
@@ -210,6 +211,7 @@ include("groebner/state.jl")
 include("groebner/correctness.jl")
 # `groebner` backend
 include("groebner/groebner.jl")
+include("groebner/groebner_with_change_matrix.jl")
 # `groebner_learn` and `groebner_apply` backend
 include("groebner/learn-apply.jl")
 # `isgroebner` backend
@@ -238,6 +240,7 @@ include("precompile.jl")
 # Exports
 
 export groebner, groebner_learn, groebner_apply!
+export groebner_with_change_matrix
 export isgroebner
 export normalform
 
diff --git a/src/f4/basis.jl b/src/f4/basis.jl
index 91c730d9..697555de 100644
--- a/src/f4/basis.jl
+++ b/src/f4/basis.jl
@@ -139,23 +139,25 @@ mutable struct Basis{C <: Coeff}
 
     # Sugar degrees of basis polynomials
     sugar_cubes::Vector{SugarCube}
-end
-
-# Initialize basis for `ngens` elements with coefficient of type T
-function basis_initialize(ring::PolyRing, ngens::Int, ::Type{T}) where {T <: Coeff}
-    sz = ngens # * 2
-    ndone = 0
-    nfilled = 0
-    nlead = 0
 
-    monoms = Vector{Vector{MonomId}}(undef, sz)
-    coeffs = Vector{Vector{T}}(undef, sz)
-    isred = zeros(Bool, sz)
-    nonred = Vector{Int}(undef, sz)
-    lead = Vector{DivisionMask}(undef, sz)
-    sugar_cubes = Vector{Int}(undef, sz)
+    changematrix::Vector{Dict{Int, Dict{MonomId, C}}}
+end
 
-    Basis(monoms, coeffs, sz, ndone, nfilled, isred, nonred, lead, nlead, sugar_cubes)
+# Initialize basis for `sz` elements with coefficient of type T
+function basis_initialize(ring::PolyRing, sz::Int, ::Type{T}) where {T <: Coeff}
+    Basis(
+        Vector{Vector{MonomId}}(undef, sz),
+        Vector{Vector{T}}(undef, sz),
+        sz,
+        0,
+        0,
+        zeros(Bool, sz),
+        Vector{Int}(undef, sz),
+        Vector{DivisionMask}(undef, sz),
+        0,
+        Vector{Int}(undef, 0),
+        Vector{Dict{Int, Dict{MonomId, T}}}(undef, 0)
+    )
 end
 
 # initialize basis with the given (already hashed) monomials and coefficients.
@@ -165,16 +167,19 @@ function basis_initialize(
     coeffs::Vector{Vector{T}}
 ) where {T <: Coeff}
     sz = length(hashedexps)
-    ndone = 0
-    nfilled = sz
-    nlead = 0
-
-    isred = zeros(Bool, sz)
-    nonred = Vector{Int}(undef, sz)
-    lead = Vector{DivisionMask}(undef, sz)
-    sugar_cubes = Vector{Int}(undef, sz)
-
-    Basis(hashedexps, coeffs, sz, ndone, nfilled, isred, nonred, lead, nlead, sugar_cubes)
+    Basis(
+        hashedexps,
+        coeffs,
+        sz,
+        0,
+        sz,
+        zeros(Bool, sz),
+        Vector{Int}(undef, sz),
+        Vector{DivisionMask}(undef, sz),
+        0,
+        Vector{Int}(undef, 0),
+        Vector{Dict{Int, Dict{MonomId, T}}}(undef, 0)
+    )
 end
 
 # Same as f4_initialize_structs, but uses an existing hashtable
@@ -189,6 +194,136 @@ function basis_initialize_using_existing_hashtable(
     basis
 end
 
+###
+# Change matrix
+
+function basis_changematrix_initialize!(
+    basis::Basis{C},
+    hashtable::MonomialHashtable{M, Ord}
+) where {C <: Coeff, M <: Monom, Ord}
+    resize!(basis.changematrix, basis.nfilled)
+    id_of_1 = hashtable_insert!(hashtable, monom_construct_const_monom(M, hashtable.nvars))
+    for i in 1:(basis.nfilled)
+        basis.changematrix[i] = Dict{Int, Dict{MonomId, C}}()
+        basis.changematrix[i][i] = Dict{MonomId, C}(id_of_1 => one(C))
+    end
+    nothing
+end
+
+function basis_changematrix_mul!(
+    basis::Basis,
+    idx::Int,
+    cf::C,
+    arithmetic::AbstractArithmetic
+) where {C <: Coeff}
+    row = basis.changematrix[idx]
+    new_row = empty(row)
+    for (pos, poly) in row
+        new_poly = empty(poly)
+        for (monom_id, monom_cf) in poly
+            new_poly[monom_id] = mod_p(monom_cf * cf, arithmetic)
+        end
+        new_row[pos] = new_poly
+    end
+    basis.changematrix[idx] = new_row
+    nothing
+end
+
+function basis_changematrix_addmul!(
+    basis::Basis,
+    ht::MonomialHashtable,
+    symbol_ht,
+    idx::Int,
+    poly_idx,
+    poly_mult,
+    cf::C,
+    arithmetic::AbstractArithmeticZp{AccumType, CoeffType}
+) where {C <: Coeff, AccumType, CoeffType}
+    row = basis.changematrix[idx]
+    if true
+        @inbounds for (ref_poly_idx, ref_quo) in basis.changematrix[poly_idx]
+            if !haskey(row, ref_poly_idx)
+                row[ref_poly_idx] = Dict{MonomId, C}()
+            end
+            poly = row[ref_poly_idx]
+            hashtable_resize_if_needed!(ht, length(ref_quo))
+            for (ref_mult, ref_cf) in ref_quo
+                ref_monom = ht.monoms[ref_mult]
+                quo_monom = ht.monoms[poly_mult]
+                new_monom = monom_copy(ref_monom)
+                new_monom = monom_product!(new_monom, ref_monom, quo_monom)
+                new_monom_id = hashtable_insert!(ht, new_monom)
+                if !haskey(poly, new_monom_id)
+                    poly[new_monom_id] = zero(CoeffType)
+                end
+                poly[new_monom_id] =
+                    mod_p(poly[new_monom_id] + ref_cf * AccumType(cf), arithmetic)
+            end
+        end
+    else
+        if !haskey(row, poly_idx)
+            row[poly_idx] = Dict{MonomId, C}()
+        end
+        poly = row[poly_idx]
+        if !haskey(poly, poly_mult)
+            poly[poly_mult] = zero(CoeffType)
+        end
+        poly[poly_mult] = mod_p(poly[poly_mult] + AccumType(cf), arithmetic)
+    end
+    nothing
+end
+
+function basis_changematrix_deep_copy_with_new_type(
+    changematrix,
+    new_sparse_row_coeffs::Vector{Vector{C}}
+) where {C}
+    new_changematrix = Vector{Dict{Int, Dict{MonomId, C}}}(undef, length(changematrix))
+    for i in 1:length(changematrix)
+        !isassigned(changematrix, i) && continue
+        obj = changematrix[i]
+        new_obj = Dict{Int, Dict{MonomId, C}}()
+        for (k, v) in obj
+            new_obj[k] = Dict{MonomId, C}()
+            for (k2, v2) in v
+                new_obj[k][k2] = C(v2)
+            end
+        end
+        new_changematrix[i] = new_obj
+    end
+    new_changematrix
+end
+
+function basis_changematrix_shallow_copy_with_new_type(changematrix, new_sparse_row_coeffs)
+    basis_changematrix_deep_copy_with_new_type(changematrix, new_sparse_row_coeffs)
+end
+
+function basis_changematrix_export(
+    basis::Basis{C},
+    ht::MonomialHashtable{M},
+    npolys
+) where {C <: Coeff, M <: Monom}
+    @log :debug basis.changematrix
+    matrix_monoms = Vector{Vector{Vector{M}}}(undef, basis.nnonredundant)
+    matrix_coeffs = Vector{Vector{Vector{C}}}(undef, basis.nnonredundant)
+    @inbounds for i in 1:(basis.nnonredundant)
+        matrix_monoms[i] = Vector{Vector{M}}(undef, npolys)
+        matrix_coeffs[i] = Vector{Vector{M}}(undef, npolys)
+        idx = basis.nonredundant[i]
+        row = basis.changematrix[idx]
+        for poly_idx in 1:npolys
+            matrix_monoms[i][poly_idx] = Vector{M}()
+            matrix_coeffs[i][poly_idx] = Vector{C}()
+            !haskey(row, poly_idx) && continue
+            for (monom_id, monom_cf) in row[poly_idx]
+                push!(matrix_monoms[i][poly_idx], ht.monoms[monom_id])
+                push!(matrix_coeffs[i][poly_idx], monom_cf)
+            end
+        end
+        sort_input_terms_to_change_ordering!(matrix_monoms[i], matrix_coeffs[i], ht.ord)
+    end
+    matrix_monoms, matrix_coeffs
+end
+
 ###
 # Basis utils
 
@@ -206,7 +341,11 @@ function basis_shallow_copy_with_new_coeffs(
         basis.nonredundant,
         basis.divmasks,
         basis.nnonredundant,
-        basis.sugar_cubes
+        basis.sugar_cubes,
+        basis_changematrix_shallow_copy_with_new_type(
+            basis.changematrix,
+            new_sparse_row_coeffs
+        )
     )
 end
 
@@ -236,7 +375,11 @@ function basis_deep_copy_with_new_coeffs(
         copy(basis.nonredundant),
         copy(basis.divmasks),
         basis.nnonredundant,
-        copy(basis.sugar_cubes)
+        copy(basis.sugar_cubes),
+        basis_changematrix_deep_copy_with_new_type(
+            basis.changematrix,
+            new_sparse_row_coeffs
+        )
     )
 end
 
@@ -274,16 +417,17 @@ function basis_resize_if_needed!(basis::Basis{T}, to_add::Int) where {T}
         @inbounds basis.isredundant[(basis.nprocessed + 1):end] .= false
         resize!(basis.nonredundant, basis.size)
         resize!(basis.divmasks, basis.size)
-        resize!(basis.sugar_cubes, basis.size)
+        # resize!(basis.sugar_cubes, basis.size)
     end
     @invariant basis.size >= basis.nprocessed + to_add
     nothing
 end
 
 # Normalize each element of the basis to have leading coefficient equal to 1
-@timeit function basis_normalize!(
+function basis_normalize!(
     basis::Basis{C},
-    arithmetic::AbstractArithmeticZp{A, C}
+    arithmetic::AbstractArithmeticZp{A, C},
+    changematrix::Bool
 ) where {A <: Union{CoeffZp, CompositeCoeffZp}, C <: Union{CoeffZp, CompositeCoeffZp}}
     @log :debug "Normalizing polynomials in the basis"
     cfs = basis.coeffs
@@ -295,6 +439,9 @@ end
             cfs[i][j] = mod_p(A(cfs[i][j]) * A(mul), arithmetic) % C
         end
         @invariant isone(cfs[i][1])
+        if changematrix
+            basis_changematrix_mul!(basis, i, A(mul), arithmetic)
+        end
     end
     basis
 end
@@ -302,7 +449,8 @@ end
 # Normalize each element of the basis by dividing it by its leading coefficient
 function basis_normalize!(
     basis::Basis{C},
-    arithmetic::AbstractArithmeticQQ
+    arithmetic::AbstractArithmeticQQ,
+    changematrix::Bool
 ) where {C <: CoeffQQ}
     @log :debug "Normalizing polynomials in the basis"
     cfs = basis.coeffs
@@ -560,7 +708,8 @@ end
     basis::Basis,
     ht::MonomialHashtable,
     ord,
-    arithmetic
+    arithmetic,
+    changematrix::Bool
 )
     @inbounds for i in 1:(basis.nnonredundant)
         idx = basis.nonredundant[i]
@@ -568,6 +717,9 @@ end
         basis.isredundant[i] = false
         basis.coeffs[i] = basis.coeffs[idx]
         basis.monoms[i] = basis.monoms[idx]
+        if changematrix
+            basis.changematrix[i] = basis.changematrix[idx]
+        end
     end
     basis.size = basis.nprocessed = basis.nfilled = basis.nnonredundant
     resize!(basis.coeffs, basis.nprocessed)
@@ -575,9 +727,10 @@ end
     resize!(basis.divmasks, basis.nprocessed)
     resize!(basis.nonredundant, basis.nprocessed)
     resize!(basis.isredundant, basis.nprocessed)
-    resize!(basis.sugar_cubes, basis.nprocessed)
-    sort_polys_by_lead_increasing!(basis, ht, ord=ord)
-    basis_normalize!(basis, arithmetic)
+    resize!(basis.changematrix, basis.nprocessed)
+    # resize!(basis.sugar_cubes, basis.nprocessed)
+    sort_polys_by_lead_increasing!(basis, ht, changematrix, ord=ord)
+    basis_normalize!(basis, arithmetic, changematrix)
 end
 
 # Returns the monomials of the polynomials in the basis
diff --git a/src/f4/f4.jl b/src/f4/f4.jl
index 457a582f..cd8e16a3 100644
--- a/src/f4/f4.jl
+++ b/src/f4/f4.jl
@@ -48,8 +48,12 @@
     basis_fill_data!(basis, hashtable, monoms, coeffs)
     hashtable_fill_divmasks!(hashtable)
 
+    if params.changematrix
+        basis_changematrix_initialize!(basis, hashtable)
+    end
+
     if sort_input
-        permutation = sort_polys_by_lead_increasing!(basis, hashtable)
+        permutation = sort_polys_by_lead_increasing!(basis, hashtable, params.changematrix)
     else
         permutation = collect(1:(basis.nfilled))
     end
@@ -58,7 +62,7 @@
     # We do not need normalization in some cases, e.g., when computing the
     # normal forms
     if normalize_input
-        basis_normalize!(basis, params.arithmetic)
+        basis_normalize!(basis, params.arithmetic, params.changematrix)
     end
 
     basis, pairset, hashtable, permutation
@@ -78,7 +82,7 @@ end
     # Call the linear algebra backend
     linalg_main!(matrix, basis, params)
     # Extract nonzero rows from the matrix into the basis
-    matrix_convert_rows_to_basis_elements!(matrix, basis, ht, symbol_ht)
+    matrix_convert_rows_to_basis_elements!(matrix, basis, ht, symbol_ht, params)
 end
 
 # F4 update
@@ -208,7 +212,7 @@ function f4_autoreduce!(
 
     linalg_autoreduce!(matrix, basis, params)
 
-    matrix_convert_rows_to_basis_elements!(matrix, basis, ht, symbol_ht)
+    matrix_convert_rows_to_basis_elements!(matrix, basis, ht, symbol_ht, params)
 
     basis.nfilled = matrix.npivots + basis.nprocessed
     basis.nprocessed = matrix.npivots
@@ -575,6 +579,8 @@ function f4_add_critical_pairs_to_matrix!(
         # over all polys with same lcm,
         # add them to the lower part of matrix
         for k in 1:npolys
+            # TODO: some lower rows can be discarded as useless on this stage.
+
             # duplicate generator,
             # we can do so as long as generators are sorted
             if polys[k] == prev
@@ -696,7 +702,7 @@ end
     # @invariant pairset_well_formed(:input_f4!, pairset, basis, ht)
 
     @log :debug "Entering F4."
-    basis_normalize!(basis, params.arithmetic)
+    basis_normalize!(basis, params.arithmetic, params.changematrix)
 
     matrix = matrix_initialize(ring, C)
 
@@ -730,16 +736,13 @@ end
                 symbol_ht,
                 maxpairs=params.maxpairs
             )
-            # matrix    = matrix_initialize(ring, C)
-            # symbol_ht = hashtable_initialize_secondary(hashtable)
             matrix_reinitialize!(matrix, 0)
             hashtable_reinitialize!(symbol_ht)
             continue
         end
 
-        # selects pairs for reduction from pairset following normal strategy
-        # (minimal lcm degrees are selected),
-        # and puts these into the matrix rows
+        # selects pairs for reduction from the pairset and puts these pairs into
+        # the matrix rows
         f4_select_critical_pairs!(
             pairset,
             basis,
@@ -794,7 +797,14 @@ end
         f4_autoreduce!(ring, basis, matrix, hashtable, symbol_ht, params)
     end
 
-    basis_standardize!(ring, basis, hashtable, hashtable.ord, params.arithmetic)
+    basis_standardize!(
+        ring,
+        basis,
+        hashtable,
+        hashtable.ord,
+        params.arithmetic,
+        params.changematrix
+    )
 
     # @invariant hashtable_well_formed(:output_f4!, ring, hashtable)
     @invariant basis_well_formed(:output_f4!, ring, basis, hashtable)
diff --git a/src/f4/learn-apply.jl b/src/f4/learn-apply.jl
index ce2af9af..eacb6e3f 100644
--- a/src/f4/learn-apply.jl
+++ b/src/f4/learn-apply.jl
@@ -55,9 +55,9 @@ function standardize_basis_in_learn!(
     resize!(basis.divmasks, basis.nprocessed)
     resize!(basis.nonredundant, basis.nprocessed)
     resize!(basis.isredundant, basis.nprocessed)
-    perm = sort_polys_by_lead_increasing!(basis, ht, ord=ord)
+    perm = sort_polys_by_lead_increasing!(basis, ht, false, ord=ord)
     trace.output_sort_indices = perm
-    basis_normalize!(basis, arithmetic)
+    basis_normalize!(basis, arithmetic, false)
 end
 
 function matrix_compute_pivot_signature(pivots::Vector{Vector{MonomId}}, from::Int, sz::Int)
@@ -84,7 +84,7 @@ function reduction_learn!(
 )
     matrix_fill_column_to_monom_map!(matrix, symbol_ht)
     linalg_main!(matrix, basis, params, trace, linalg=LinearAlgebra(:learn, :sparse))
-    matrix_convert_rows_to_basis_elements!(matrix, basis, hashtable, symbol_ht)
+    matrix_convert_rows_to_basis_elements!(matrix, basis, hashtable, symbol_ht, params)
     pivot_indices =
         map(i -> Int32(basis.monoms[basis.nprocessed + i][1]), 1:(matrix.npivots))
     push!(trace.matrix_pivot_indices, pivot_indices)
@@ -144,7 +144,7 @@ function f4_reducegb_learn!(
 
     linalg_autoreduce!(matrix, basis, params, trace, linalg=LinearAlgebra(:learn, :sparse))
 
-    matrix_convert_rows_to_basis_elements!(matrix, basis, ht, symbol_ht)
+    matrix_convert_rows_to_basis_elements!(matrix, basis, ht, symbol_ht, params)
 
     basis.nfilled = matrix.npivots + basis.nprocessed
     basis.nprocessed = matrix.npivots
@@ -191,7 +191,7 @@ end
     @invariant params.reduced
 
     @log :debug "Entering F4 Learn phase."
-    basis_normalize!(basis, params.arithmetic)
+    basis_normalize!(basis, params.arithmetic, params.changematrix)
 
     matrix = matrix_initialize(ring, C)
 
@@ -355,7 +355,7 @@ function reduction_apply!(
         return (false, false)
     end
 
-    matrix_convert_rows_to_basis_elements!(matrix, basis, ht, symbol_ht)
+    matrix_convert_rows_to_basis_elements!(matrix, basis, ht, symbol_ht, params)
 
     # Check that the leading terms were not reduced to zero accidentally
     pivot_indices = trace.matrix_pivot_indices[f4_iteration]
@@ -557,7 +557,7 @@ function autoreduce_f4_apply!(
         @log :warn "In apply, the final autoreduction of the basis failed"
         return false
     end
-    matrix_convert_rows_to_basis_elements!(matrix, basis, hashtable, symbol_ht)
+    matrix_convert_rows_to_basis_elements!(matrix, basis, hashtable, symbol_ht, params)
 
     basis.nfilled = matrix.npivots + basis.nprocessed
     basis.nprocessed = matrix.npivots
@@ -585,7 +585,7 @@ function standardize_basis_in_apply!(ring::PolyRing, trace::TraceF4, arithmetic)
     end
     buf.nprocessed = buf.nnonredundant = 0
     buf.nfilled = trace.input_basis.nfilled
-    basis_normalize!(basis, arithmetic)
+    basis_normalize!(basis, arithmetic, false)
 end
 
 @timeit function f4_apply!(
@@ -597,7 +597,7 @@ end
     @invariant basis_well_formed(:input_f4_apply!, ring, basis, trace.hashtable)
     @invariant params.reduced
 
-    basis_normalize!(basis, params.arithmetic)
+    basis_normalize!(basis, params.arithmetic, params.changematrix)
 
     iters_total = length(trace.matrix_infos) - 1
     iters = 0
diff --git a/src/f4/linalg/backend.jl b/src/f4/linalg/backend.jl
index ed9ce4d2..91040da2 100644
--- a/src/f4/linalg/backend.jl
+++ b/src/f4/linalg/backend.jl
@@ -32,6 +32,7 @@ function linalg_deterministic_sparse_interreduction!(
     basis::Basis,
     arithmetic::AbstractArithmetic
 )
+    sort_matrix_upper_rows!(matrix)
     @log :matrix "linalg_deterministic_sparse_interreduction!"
     @log :matrix matrix_string_repr(matrix)
 
@@ -825,7 +826,7 @@ function linalg_normalize_row!(
     @invariant !iszero(row[first_nnz_index])
 
     lead = row[first_nnz_index]
-    isone(lead) && return row
+    isone(lead) && return lead
 
     @inbounds pinv = inv_mod_p(A(lead), arithmetic) % T
     @inbounds row[first_nnz_index] = one(T)
@@ -833,7 +834,7 @@ function linalg_normalize_row!(
         row[i] = mod_p(A(row[i]) * A(pinv), arithmetic) % T
     end
 
-    row
+    pinv
 end
 
 # Normalize the row to have the leading coefficient equal to 1
@@ -845,7 +846,7 @@ function linalg_normalize_row!(
     @invariant !iszero(row[first_nnz_index])
 
     lead = row[first_nnz_index]
-    isone(lead) && return row
+    isone(lead) && return lead
 
     @inbounds pinv = inv(lead)
     @inbounds row[1] = one(T)
@@ -853,7 +854,7 @@ function linalg_normalize_row!(
         row[i] = row[i] * pinv
     end
 
-    row
+    pinv
 end
 
 # Linear combination of dense vector and sparse vector modulo a prime.
@@ -876,7 +877,7 @@ function linalg_vector_addmul_sparsedense_mod_p!(
         row[idx] = mod_p(a, arithmetic)
     end
 
-    nothing
+    mul
 end
 
 # Linear combination of dense vector and sparse vector.
diff --git a/src/f4/linalg/backend_changematrix.jl b/src/f4/linalg/backend_changematrix.jl
new file mode 100644
index 00000000..89db09d0
--- /dev/null
+++ b/src/f4/linalg/backend_changematrix.jl
@@ -0,0 +1,346 @@
+# This file is a part of Groebner.jl. License is GNU GPL v2.
+
+###
+# High level
+
+function linalg_deterministic_sparse_changematrix!(
+    matrix::MacaulayMatrix,
+    basis::Basis,
+    linalg::LinearAlgebra,
+    arithmetic::AbstractArithmetic
+)
+    sort_matrix_upper_rows!(matrix) # for the AB part
+    sort_matrix_lower_rows!(matrix) # for the CD part
+
+    @log :matrix "linalg_deterministic_sparse_changematrix!"
+    @log :matrix matrix_string_repr(matrix)
+
+    matrix_changematrix_initialize!(matrix, matrix.nrows_filled_lower)
+
+    # Reduce CD with AB
+    linalg_reduce_matrix_lower_part_with_changematrix!(matrix, basis, arithmetic)
+    # Interreduce CD
+    linalg_interreduce_matrix_pivots_with_changematrix!(matrix, basis, arithmetic)
+    true
+end
+
+function linalg_deterministic_sparse_interreduction_changematrix!(
+    matrix::MacaulayMatrix,
+    basis::Basis{C},
+    arithmetic::AbstractArithmetic
+) where {C}
+    sort_matrix_upper_rows!(matrix)
+    @log :matrix "linalg_deterministic_sparse_interreduction_changematrix!"
+    @log :matrix matrix_string_repr(matrix)
+
+    matrix_changematrix_initialize!(matrix, matrix.nrows_filled_upper)
+    for i in 1:(matrix.nrows_filled_upper)
+        poly_id = matrix.upper_to_coeffs[i]
+        mult_id = matrix.upper_to_mult[i]
+        matrix.changematrix[i] = Dict{Tuple{Int, MonomId}, C}((poly_id, mult_id) => one(C))
+    end
+
+    # Prepare the matrix
+    linalg_prepare_matrix_pivots_in_interreduction!(matrix, basis)
+    # Interreduce AB
+    linalg_interreduce_matrix_pivots_with_changematrix!(
+        matrix,
+        basis,
+        arithmetic,
+        reversed_rows=true
+    )
+    true
+end
+
+###
+# Low level
+
+function linalg_reduce_matrix_lower_part_with_changematrix!(
+    matrix::MacaulayMatrix{CoeffType},
+    basis::Basis{CoeffType},
+    arithmetic::AbstractArithmetic{AccumType, CoeffType}
+) where {CoeffType <: Coeff, AccumType <: Coeff}
+    _, ncols = size(matrix)
+    nleft, _ = matrix_ncols_filled(matrix)
+    nup, nlow = matrix_nrows_filled(matrix)
+
+    # Prepare the matrix
+    pivots, row_idx_to_coeffs = linalg_prepare_matrix_pivots!(matrix)
+    resize!(matrix.some_coeffs, nlow)
+    changematrix = matrix.changematrix
+
+    # Allocate the buffers
+    reducer_rows = Vector{Tuple{Int, CoeffType}}()
+    row = zeros(AccumType, ncols)
+    new_sparse_row_support, new_sparse_row_coeffs = linalg_new_empty_sparse_row(CoeffType)
+
+    @inbounds for i in 1:nlow
+        sparse_row_support = matrix.lower_rows[i]
+        sparse_row_coeffs = basis.coeffs[row_idx_to_coeffs[i]]
+
+        linalg_load_sparse_row!(row, sparse_row_support, sparse_row_coeffs)
+
+        # Additionally record the indices of rows that participated in reduction
+        # of the given row
+        empty!(reducer_rows)
+        first_nnz_column = sparse_row_support[1]
+        zeroed = linalg_reduce_dense_row_by_pivots_sparse_changematrix!(
+            new_sparse_row_support,
+            new_sparse_row_coeffs,
+            row,
+            matrix,
+            basis,
+            pivots,
+            first_nnz_column,
+            ncols,
+            arithmetic,
+            reducer_rows,
+            tmp_pos=-1
+        )
+
+        # if fully reduced
+        zeroed && continue
+
+        # NOTE: we are not recording reducers from the lower part of the matrix
+        cm = changematrix[i]
+        cm[(row_idx_to_coeffs[i], matrix.lower_to_mult[i])] = one(CoeffType)
+        for (idx, quo) in reducer_rows
+            if idx <= nleft
+                poly_idx, poly_mult = matrix.upper_to_coeffs[idx], matrix.upper_to_mult[idx]
+                if !haskey(cm, (poly_idx, poly_mult))
+                    cm[(poly_idx, poly_mult)] = zero(CoeffType)
+                end
+                cm[(poly_idx, poly_mult)] =
+                    mod_p(cm[(poly_idx, poly_mult)] + AccumType(quo), arithmetic)
+            else
+                row_idx = matrix.lower_to_coeffs[idx]
+                for ((poly_idx, poly_mult), cf) in changematrix[row_idx]
+                    if !haskey(cm, (poly_idx, poly_mult))
+                        cm[(poly_idx, poly_mult)] = zero(CoeffType)
+                    end
+                    cm[(poly_idx, poly_mult)] =
+                        mod_p(cm[(poly_idx, poly_mult)] + quo * AccumType(cf), arithmetic)
+                end
+            end
+        end
+
+        @invariant length(new_sparse_row_support) == length(new_sparse_row_coeffs)
+        pinv = linalg_normalize_row!(new_sparse_row_coeffs, arithmetic)
+
+        cm2 = empty(cm)
+        for ((poly_idx, quo_idx), quo_cf) in cm
+            cm2[(poly_idx, quo_idx)] = mod_p(AccumType(pinv) * quo_cf, arithmetic)
+        end
+        changematrix[i] = cm2
+
+        matrix.some_coeffs[i] = new_sparse_row_coeffs
+        pivots[new_sparse_row_support[1]] = new_sparse_row_support
+        matrix.lower_to_coeffs[new_sparse_row_support[1]] = i
+
+        new_sparse_row_support, new_sparse_row_coeffs =
+            linalg_new_empty_sparse_row(CoeffType)
+    end
+
+    true
+end
+
+function linalg_interreduce_matrix_pivots_with_changematrix!(
+    matrix::MacaulayMatrix{CoeffType},
+    basis::Basis{CoeffType},
+    arithmetic::AbstractArithmetic{AccumType, CoeffType};
+    reversed_rows::Bool=false
+) where {CoeffType <: Coeff, AccumType <: Coeff}
+    _, ncols = size(matrix)
+    nleft, nright = matrix_ncols_filled(matrix)
+    nupper, _ = matrix_nrows_filled(matrix)
+
+    # Prepare the matrix
+    resize!(matrix.lower_rows, nright)
+    pivots = matrix.pivots
+    new_pivots = 0
+    any_zeroed = false
+    changematrix = matrix.changematrix
+    compact_changematrix = similar(changematrix)
+
+    # Allocate the buffers
+    reducer_rows = Vector{Tuple{Int, CoeffType}}()
+    row = zeros(AccumType, ncols)
+    # Indices of rows that did no reduce to zero
+    not_reduced_to_zero = Vector{Int}(undef, nright)
+
+    # for each column in the block D..
+    @inbounds for i in 1:nright
+        abs_column_idx = ncols - i + 1
+        # Check if there is a row that starts at `abs_column_idx`
+        !isassigned(pivots, abs_column_idx) && continue
+
+        # Locate the row support and coefficients
+        sparse_row_support = pivots[abs_column_idx]
+        if abs_column_idx <= nleft
+            # upper part of matrix
+            sparse_row_coeffs = basis.coeffs[matrix.upper_to_coeffs[abs_column_idx]]
+        else
+            # lower part of matrix
+            sparse_row_coeffs = matrix.some_coeffs[matrix.lower_to_coeffs[abs_column_idx]]
+        end
+        @invariant length(sparse_row_support) == length(sparse_row_coeffs)
+
+        # Load the row into a dense array
+        linalg_load_sparse_row!(row, sparse_row_support, sparse_row_coeffs)
+
+        empty!(reducer_rows)
+        new_sparse_row_support, new_sparse_row_coeffs =
+            linalg_new_empty_sparse_row(CoeffType)
+        first_nnz_column = sparse_row_support[1]
+        zeroed = linalg_reduce_dense_row_by_pivots_sparse_changematrix!(
+            new_sparse_row_support,
+            new_sparse_row_coeffs,
+            row,
+            matrix,
+            basis,
+            pivots,
+            first_nnz_column,
+            ncols,
+            arithmetic,
+            tmp_pos=first_nnz_column,
+            reducer_rows
+        )
+
+        # If the row reduced to zero
+        if zeroed
+            any_zeroed = true
+            continue
+        end
+
+        current_row_idx = matrix.lower_to_coeffs[abs_column_idx]
+        cm = changematrix[current_row_idx]
+        for (idx, quo) in reducer_rows
+            if idx <= nleft
+                @assert false
+                poly_idx, poly_mult = matrix.upper_to_coeffs[idx], matrix.upper_to_mult[idx]
+                cm[(poly_idx, poly_mult)] = quo
+            else
+                row_idx = matrix.lower_to_coeffs[idx]
+                for ((poly_idx, poly_mult), cf) in changematrix[row_idx]
+                    if !haskey(cm, (poly_idx, poly_mult))
+                        cm[(poly_idx, poly_mult)] = zero(CoeffType)
+                    end
+                    cm[(poly_idx, poly_mult)] =
+                        mod_p(cm[(poly_idx, poly_mult)] + quo * AccumType(cf), arithmetic)
+                end
+            end
+        end
+
+        new_pivots += 1
+        not_reduced_to_zero[new_pivots] = i
+
+        # Update the row support and coefficients.
+        # TODO: maybe get rid of the reversed_rows parameter?
+        if !reversed_rows
+            compact_changematrix[new_pivots] = cm
+            matrix.lower_rows[new_pivots] = new_sparse_row_support
+            matrix.some_coeffs[matrix.lower_to_coeffs[abs_column_idx]] =
+                new_sparse_row_coeffs
+            pivots[abs_column_idx] = matrix.lower_rows[new_pivots]
+        else
+            compact_changematrix[nupper - new_pivots + 1] = cm
+            matrix.lower_rows[nupper - new_pivots + 1] = new_sparse_row_support
+            matrix.some_coeffs[matrix.lower_to_coeffs[abs_column_idx]] =
+                new_sparse_row_coeffs
+            pivots[abs_column_idx] = matrix.lower_rows[nupper - new_pivots + 1]
+        end
+    end
+
+    matrix.npivots = new_pivots
+    matrix.changematrix = compact_changematrix
+    resize!(matrix.changematrix, new_pivots)
+    resize!(matrix.lower_rows, new_pivots)
+    resize!(not_reduced_to_zero, new_pivots)
+
+    true, any_zeroed, not_reduced_to_zero
+end
+
+function linalg_reduce_dense_row_by_pivots_sparse_changematrix!(
+    new_sparse_row_support::Vector{I},
+    new_sparse_row_coeffs::Vector{C},
+    row::Vector{A},
+    matrix::MacaulayMatrix{C},
+    basis::Basis{C},
+    pivots::Vector{Vector{I}},
+    start_column::Integer,
+    end_column::Integer,
+    arithmetic::AbstractArithmeticZp{A, C},
+    active_reducers;
+    tmp_pos::Integer=-1,
+    exact_column_mapping::Bool=false,
+    computing_rref::Bool=false
+) where {I, C <: Union{CoeffZp, CompositeCoeffZp}, A <: Union{CoeffZp, CompositeCoeffZp}}
+    _, ncols = size(matrix)
+    nleft, _ = matrix_ncols_filled(matrix)
+
+    n_nonzeros = 0
+    new_pivot_column = -1
+
+    @inbounds for i in start_column:end_column
+        # if the element is zero -- no reduction is needed
+        if iszero(row[i])
+            continue
+        end
+
+        # if there is no pivot with the leading column equal to i
+        if !isassigned(pivots, i) || (tmp_pos != -1 && tmp_pos == i)
+            if new_pivot_column == -1
+                new_pivot_column = i
+            end
+            n_nonzeros += 1
+            continue
+        end
+        # At this point, we have determined that pivots[i] is a valid pivot.
+
+        # Locate the support and the coefficients of the pivot row
+        pivot_support = pivots[i]
+        if exact_column_mapping
+            # if pivot comes from the new pivots
+            pivot_coeffs = matrix.some_coeffs[tmp_pos]
+        elseif i <= nleft
+            # if pivot comes from the original pivots
+            if matrix.upper_part_is_rref || computing_rref
+                pivot_coeffs = matrix.upper_coeffs[i]
+            else
+                pivot_coeffs = basis.coeffs[matrix.upper_to_coeffs[i]]
+            end
+        else
+            pivot_coeffs = matrix.some_coeffs[matrix.lower_to_coeffs[i]]
+        end
+        @invariant length(pivot_support) == length(pivot_coeffs)
+
+        mul = linalg_vector_addmul_sparsedense_mod_p!(
+            row,
+            pivot_support,
+            pivot_coeffs,
+            arithmetic
+        )
+
+        push!(active_reducers, (i, mul))
+
+        @invariant iszero(row[i])
+    end
+
+    # all row elements reduced to zero!
+    if n_nonzeros == 0
+        return true
+    end
+
+    # form the resulting row in sparse format
+    resize!(new_sparse_row_support, n_nonzeros)
+    resize!(new_sparse_row_coeffs, n_nonzeros)
+    linalg_extract_sparse_row!(
+        new_sparse_row_support,
+        new_sparse_row_coeffs,
+        row,
+        convert(Int, start_column),
+        ncols
+    )
+
+    false
+end
diff --git a/src/f4/linalg/linalg.jl b/src/f4/linalg/linalg.jl
index f5cb9a71..45dc2a0a 100644
--- a/src/f4/linalg/linalg.jl
+++ b/src/f4/linalg/linalg.jl
@@ -44,6 +44,7 @@ Returns `true` if successful and `false` otherwise.
     if isnothing(linalg)
         linalg = params.linalg
     end
+    changematrix = params.changematrix
 
     # Decide if multi-threading should be used. Further in dispatch, the backend
     # may opt to NOT use multi-threading regardless of this setting.
@@ -65,9 +66,18 @@ Returns `true` if successful and `false` otherwise.
     end
 
     flag = if !isnothing(trace)
-        _linalg_main_with_trace!(trace, matrix, basis, linalg, threaded, arithmetic, rng)
+        _linalg_main_with_trace!(
+            trace,
+            matrix,
+            basis,
+            linalg,
+            threaded,
+            changematrix,
+            arithmetic,
+            rng
+        )
     else
-        _linalg_main!(matrix, basis, linalg, threaded, arithmetic, rng)
+        _linalg_main!(matrix, basis, linalg, threaded, changematrix, arithmetic, rng)
     end
 
     flag
@@ -93,11 +103,12 @@ function linalg_autoreduce!(
     if isnothing(linalg)
         linalg = params.linalg
     end
+    changematrix = params.changematrix
 
     flag = if !isnothing(trace)
-        _linalg_autoreduce_with_trace!(trace, matrix, basis, linalg, arithmetic)
+        _linalg_autoreduce_with_trace!(trace, matrix, basis, linalg, changematrix, arithmetic)
     else
-        _linalg_autoreduce!(matrix, basis, linalg, arithmetic)
+        _linalg_autoreduce!(matrix, basis, linalg, changematrix, arithmetic)
     end
 
     flag
@@ -175,10 +186,13 @@ function _linalg_main!(
     basis::Basis,
     linalg::LinearAlgebra,
     threaded::Symbol,
+    changematrix::Bool,
     arithmetic::AbstractArithmetic,
     rng::AbstractRNG
 )
-    flag = if linalg.algorithm === :deterministic
+    flag = if changematrix
+        linalg_deterministic_sparse_changematrix!(matrix, basis, linalg, arithmetic)
+    elseif linalg.algorithm === :deterministic
         if linalg_should_use_threading(matrix, linalg, threaded)
             linalg_deterministic_sparse_threaded!(matrix, basis, linalg, arithmetic)
         else
@@ -200,9 +214,10 @@ function _linalg_main!(
         linalg_randomized_hashcolumns_sparse!(matrix, basis, linalg, arithmetic, rng)
     else
         __throw_linalg_error("Cannot pick a suitable linear algebra backend for parameters
-                             linalg     = $linalg
-                             threaded   = $threaded
-                             arithmetic = $arithmetic")
+                             linalg         = $linalg
+                             threaded       = $threaded
+                             arithmetic     = $arithmetic
+                             changematrix   = $changematrix")
     end
 
     flag
@@ -214,6 +229,7 @@ function _linalg_main_with_trace!(
     basis::Basis,
     linalg::LinearAlgebra,
     threaded::Symbol,
+    changematrix::Bool,
     arithmetic::AbstractArithmetic,
     rng::AbstractRNG
 )
@@ -236,11 +252,14 @@ function _linalg_autoreduce!(
     matrix::MacaulayMatrix,
     basis::Basis,
     linalg::LinearAlgebra,
+    changematrix::Bool,
     arithmetic::AbstractArithmetic
 )
-    sort_matrix_upper_rows!(matrix)
-    linalg_deterministic_sparse_interreduction!(matrix, basis, arithmetic)
-    true
+    if changematrix
+        linalg_deterministic_sparse_interreduction_changematrix!(matrix, basis, arithmetic)
+    else
+        linalg_deterministic_sparse_interreduction!(matrix, basis, arithmetic)
+    end
 end
 
 function _linalg_autoreduce_with_trace!(
@@ -248,8 +267,10 @@ function _linalg_autoreduce_with_trace!(
     matrix::MacaulayMatrix,
     basis::Basis,
     linalg::LinearAlgebra,
+    changematrix::Bool,
     arithmetic::AbstractArithmetic
 )
+    @assert !changematrix
     sort_matrix_upper_rows!(matrix)
 
     if linalg.algorithm === :learn
diff --git a/src/f4/matrix.jl b/src/f4/matrix.jl
index a37256a5..9108b407 100644
--- a/src/f4/matrix.jl
+++ b/src/f4/matrix.jl
@@ -109,6 +109,8 @@ mutable struct MacaulayMatrix{T <: Coeff}
     use_preallocated_buffers::Bool
     preallocated_buffer1_load::Int
     preallocated_buffer1::Vector{Vector{ColumnLabel}}
+
+    changematrix::Vector{Dict{Tuple{Int, MonomId}, T}}
 end
 
 # The number of allocated (not necessarily filled) rows and columns in the
@@ -164,7 +166,8 @@ function matrix_initialize(ring::PolyRing, ::Type{T}) where {T <: Coeff}
         Vector{Int8}(),
         false,
         0,
-        Vector{Vector{ColumnLabel}}()
+        Vector{Vector{ColumnLabel}}(),
+        Vector{Dict{Tuple{Int, MonomId}, T}}()
     )
 end
 
@@ -306,14 +309,28 @@ end
 # end
 
 ###
+# Change matrix
+
+function matrix_changematrix_initialize!(matrix::MacaulayMatrix{T}, n) where {T <: Coeff}
+    resize!(matrix.changematrix, n)
+    for i in 1:length(matrix.changematrix)
+        matrix.changematrix[i] = Dict{Tuple{Int, MonomId}, T}()
+    end
+    nothing
+end
+
+###
+# Stuff
 
 function matrix_convert_rows_to_basis_elements!(
     matrix::MacaulayMatrix,
-    basis::Basis,
+    basis::Basis{C},
     ht::MonomialHashtable,
-    symbol_ht::MonomialHashtable
-)
+    symbol_ht::MonomialHashtable,
+    params::AlgorithmParameters
+) where {C}
     # We mutate the basis directly by adding new elements
+    @log :debug "# new pivots: $(matrix.npivots)"
 
     basis_resize_if_needed!(basis, matrix.npivots)
     rows = matrix.lower_rows
@@ -332,6 +349,25 @@ function matrix_convert_rows_to_basis_elements!(
         @invariant length(basis.coeffs[crs + i]) == length(basis.monoms[crs + i])
     end
 
+    if params.changematrix
+        resize!(basis.changematrix, basis.nfilled + matrix.npivots)
+        for i in 1:(matrix.npivots)
+            basis.changematrix[crs + i] = Dict{Int, Dict{MonomId, C}}()
+            for ((poly_idx, poly_mult), cf) in matrix.changematrix[i]
+                basis_changematrix_addmul!(
+                    basis,
+                    ht,
+                    symbol_ht,
+                    crs + i,
+                    poly_idx,
+                    poly_mult,
+                    cf,
+                    params.arithmetic
+                )
+            end
+        end
+    end
+
     basis.nfilled += matrix.npivots
 end
 
diff --git a/src/f4/sort.jl b/src/f4/sort.jl
index 30b6ce3e..eb42e1af 100644
--- a/src/f4/sort.jl
+++ b/src/f4/sort.jl
@@ -51,6 +51,7 @@ end
 function sort_polys_by_lead_increasing!(
     basis::Basis,
     hashtable::MonomialHashtable,
+    changematrix::Bool,
     abc...;
     ord::Ord=hashtable.ord
 ) where {Ord <: AbstractInternalOrdering}
@@ -78,6 +79,11 @@ function sort_polys_by_lead_increasing!(
         a[1:(basis.nfilled)] = a[permutation]
     end
 
+    if changematrix
+        @invariant length(basis.changematrix) >= basis.nfilled
+        @inbounds basis.changematrix[1:(basis.nfilled)] = basis.changematrix[permutation]
+    end
+
     permutation
 end
 
diff --git a/src/fglm/fglm_internal.jl b/src/fglm/fglm_internal.jl
index addbeb5b..93cf89e2 100644
--- a/src/fglm/fglm_internal.jl
+++ b/src/fglm/fglm_internal.jl
@@ -134,7 +134,7 @@ end
         enumerator_produce_next_monomials!(monom_enumerator, ht, monom)
     end
 
-    basis_standardize!(ring, new_basis, ht, ord, params.arithmetic)
+    basis_standardize!(ring, new_basis, ht, ord, params.arithmetic, params.changematrix)
 
     linbasis = extract_linear_basis(ring, matrix)
     new_basis, linbasis, ht
@@ -223,7 +223,7 @@ function _fglm_residuals_in_batch!(
 
     @label End
 
-    basis_standardize!(ring, new_basis, ht, ord, params.arithmetic)
+    basis_standardize!(ring, new_basis, ht, ord, params.arithmetic, params.changematrix)
 
     linbasis = extract_linear_basis(ring, matrix)
     new_basis, linbasis, ht
diff --git a/src/groebner/autoreduce.jl b/src/groebner/autoreduce.jl
index 9c207540..49c9d59a 100644
--- a/src/groebner/autoreduce.jl
+++ b/src/groebner/autoreduce.jl
@@ -14,6 +14,13 @@ function _autoreduce1(
     matrix = matrix_initialize(ring, C)
     symbol_ht = hashtable_initialize_secondary(hashtable)
     f4_autoreduce!(ring, basis, matrix, hashtable, symbol_ht, params)
-    basis_standardize!(ring, basis, hashtable, hashtable.ord, params.arithmetic)
+    basis_standardize!(
+        ring,
+        basis,
+        hashtable,
+        hashtable.ord,
+        params.arithmetic,
+        params.changematrix
+    )
     basis_export_data(basis, hashtable)
 end
diff --git a/src/groebner/correctness.jl b/src/groebner/correctness.jl
index 3a325455..98252dd0 100644
--- a/src/groebner/correctness.jl
+++ b/src/groebner/correctness.jl
@@ -129,7 +129,7 @@ function randomized_correctness_check!(
     # random prime
     arithmetic = select_arithmetic(CoeffModular, prime, :auto, false)
     # TODO: Why is this here? F4 normalizes the basis on entry
-    basis_normalize!(gb_ff, arithmetic)
+    basis_normalize!(gb_ff, arithmetic, params.changematrix)
     f4_normalform!(ring_ff, gb_ff, input_ff, hashtable, arithmetic)
     for i in 1:(input_ff.nprocessed)
         # meaning that something is not reduced
diff --git a/src/groebner/groebner.jl b/src/groebner/groebner.jl
index 9e3d8125..b5095f7e 100644
--- a/src/groebner/groebner.jl
+++ b/src/groebner/groebner.jl
@@ -4,7 +4,7 @@
 # Backend for `groebner`
 
 # Proxy function for handling exceptions.
-function _groebner0(polynomials, kws::KeywordsHandler)
+function _groebner0(polynomials, kws::KeywordArguments)
     # We try to select an efficient internal polynomial representation, i.e., a
     # suitable representation of monomials and coefficients.
     polynomial_repr = io_select_polynomial_representation(polynomials, kws)
@@ -28,7 +28,7 @@ function _groebner0(polynomials, kws::KeywordsHandler)
     end
 end
 
-function _groebner1(polynomials, kws::KeywordsHandler, representation)
+function _groebner1(polynomials, kws::KeywordArguments, representation)
     # Extract ring information, exponents, and coefficients from input
     # polynomials. Convert these to an internal polynomial representation. 
     # NOTE: This must copy the input, so that input `polynomials` is never
diff --git a/src/groebner/groebner_with_change_matrix.jl b/src/groebner/groebner_with_change_matrix.jl
new file mode 100644
index 00000000..3c97f647
--- /dev/null
+++ b/src/groebner/groebner_with_change_matrix.jl
@@ -0,0 +1,371 @@
+# This file is a part of Groebner.jl. License is GNU GPL v2.
+
+### 
+# Backend for `groebner_with_change_matrix`
+
+# Proxy function for handling exceptions.
+function _groebner_with_change_matrix0(polynomials, kws::KeywordArguments)
+    # We try to select an efficient internal polynomial representation, i.e., a
+    # suitable representation of monomials and coefficients.
+    polynomial_repr = io_select_polynomial_representation(polynomials, kws)
+    try
+        # The backend is wrapped in a try/catch to catch exceptions that one can
+        # hope to recover from (and, perhaps, restart the computation with safer
+        # parameters).
+        return _groebner_with_change_matrix1(polynomials, kws, polynomial_repr)
+    catch err
+        if isa(err, MonomialDegreeOverflow)
+            @log :info """
+            Possible overflow of exponent vector detected. 
+            Restarting with at least $(32) bits per exponent."""
+            polynomial_repr =
+                io_select_polynomial_representation(polynomials, kws, hint=:large_exponents)
+            return _groebner_with_change_matrix1(polynomials, kws, polynomial_repr)
+        else
+            # Something bad happened.
+            rethrow(err)
+        end
+    end
+end
+
+function _groebner_with_change_matrix1(polynomials, kws::KeywordArguments, representation)
+    # Extract ring information, exponents, and coefficients from input
+    # polynomials. Convert these to an internal polynomial representation. 
+    # NOTE: This must copy the input, so that input `polynomials` is never
+    # modified.
+    # NOTE: The body of this function is type-unstable (by design)
+    ring, var_to_index, monoms, coeffs =
+        io_convert_to_internal(representation, polynomials, kws)
+
+    # Check and set parameters and monomial ordering
+    params = AlgorithmParameters(ring, representation, kws)
+    ring, _ = io_set_monomial_ordering!(ring, var_to_index, monoms, coeffs, params)
+
+    # Fast path for the input of zeros
+    if isempty(monoms)
+        @log :misc "Input consisting of zero polynomials. Returning zero."
+        changematrix = io_convert_changematrix_to_output(
+            ring,
+            polynomials,
+            0,
+            [empty(monoms)],
+            [empty(coeffs)],
+            params
+        )
+        return io_convert_to_output(ring, polynomials, monoms, coeffs, params), changematrix
+    end
+
+    if params.homogenize
+        # this also performs saturation w.r.t. the homogenizing variable
+        _, ring, monoms, coeffs = homogenize_generators!(ring, monoms, coeffs, params)
+    end
+
+    # Compute a groebner basis!
+    gbmonoms, gbcoeffs, matrix_monoms, matrix_coeffs =
+        _groebner_with_change_matrix2(ring, monoms, coeffs, params)
+
+    if params.homogenize
+        ring, gbmonoms, gbcoeffs =
+            dehomogenize_generators!(ring, gbmonoms, gbcoeffs, params)
+    end
+
+    # Convert result back to the representation of input
+    basis = io_convert_to_output(ring, polynomials, gbmonoms, gbcoeffs, params)
+    changematrix = io_convert_changematrix_to_output(
+        ring,
+        polynomials,
+        length(monoms),
+        matrix_monoms,
+        matrix_coeffs,
+        params
+    )
+
+    basis, changematrix
+end
+
+###
+# Groebner basis over Z_p.
+# Just calls f4 directly.
+
+@timeit function _groebner_with_change_matrix2(
+    ring::PolyRing,
+    monoms::Vector{Vector{M}},
+    coeffs::Vector{Vector{C}},
+    params::AlgorithmParameters
+) where {M <: Monom, C <: CoeffZp}
+    # NOTE: we can mutate ring, monoms, and coeffs here.
+    @log :misc "Backend: F4 over Z_$(ring.ch) (with change matrix)"
+    # NOTE: the sorting of input polynomials is not deterministic across
+    # different Julia versions when sorting only w.r.t. the leading term
+    basis, pairset, hashtable = f4_initialize_structs(ring, monoms, coeffs, params)
+    tracer = TinyTraceF4()
+    f4!(ring, basis, pairset, hashtable, tracer, params)
+    # Extract monomials and coefficients from basis and hashtable
+    gbmonoms, gbcoeffs = basis_export_data(basis, hashtable)
+    matrix_monoms, matrix_coeffs =
+        basis_changematrix_export(basis, hashtable, length(monoms))
+    gbmonoms, gbcoeffs, matrix_monoms, matrix_coeffs
+end
+
+###
+# Groebner basis over Q.
+
+# GB over the rationals uses modular computation.
+@timeit function _groebner_with_change_matrix2(
+    ring::PolyRing,
+    monoms::Vector{Vector{M}},
+    coeffs::Vector{Vector{C}},
+    params::AlgorithmParameters
+) where {M <: Monom, C <: CoeffQQ}
+    _groebner_with_change_classic_modular(ring, monoms, coeffs, params)
+end
+
+###
+# Classic multi-modular strategy
+
+function _groebner_with_change_classic_modular(
+    ring::PolyRing,
+    monoms::Vector{Vector{M}},
+    coeffs::Vector{Vector{C}},
+    params::AlgorithmParameters
+) where {M <: Monom, C <: CoeffQQ}
+    # NOTE: we can mutate ring, monoms, and coeffs here.
+    @log :misc "Backend: classic multi-modular F4 (with change matrix)"
+
+    # Initialize supporting structs
+    state = GroebnerState{BigInt, C, CoeffModular}(params)
+    # Initialize F4 structs
+    basis, pairset, hashtable =
+        f4_initialize_structs(ring, monoms, coeffs, params, normalize_input=false)
+    tracer = TinyTraceF4(params)
+    crt_algorithm = params.crt_algorithm
+
+    # Scale the input coefficients to integers to speed up the subsequent search
+    # for lucky primes
+    @log :debug "Input polynomials" basis
+    @log :misc "Clearing the denominators of the input polynomials"
+    basis_zz = clear_denominators!(state.buffer, basis, deepcopy=false)
+    @log :debug "Integer coefficients are" basis_zz.coeffs
+
+    # Handler for lucky primes
+    luckyprimes = LuckyPrimes(basis_zz.coeffs)
+    prime = next_lucky_prime!(luckyprimes)
+    @log :misc "The first lucky prime is $prime"
+    @log :misc "Reducing input generators modulo $prime"
+
+    # Perform reduction modulo prime and store result in basis_ff
+    ring_ff, basis_ff = reduce_modulo_p!(state.buffer, ring, basis_zz, prime, deepcopy=true)
+    @log :debug "Reduced coefficients are" basis_ff.coeffs
+
+    @log :debug "Before F4" basis_ff
+    params_zp = params_mod_p(params, prime)
+    f4!(ring_ff, basis_ff, pairset, hashtable, tracer, params_zp)
+    @log :debug "After F4:" basis_ff
+    # NOTE: basis_ff may not own its coefficients, one should not mutate it
+    # directly further in the code
+
+    push!(state.gb_coeffs_ff_all, basis_ff.coeffs)
+
+    changematrix_monoms, changematrix_coeffs =
+        basis_changematrix_export(basis_ff, hashtable, length(monoms))
+    push!(state.changematrix_coeffs_ff_all, changematrix_coeffs)
+
+    # Reconstruct coefficients and write results to the accumulator.
+    # CRT reconstrction is trivial here.
+    @log :misc "Reconstructing coefficients from Z_$prime to QQ"
+    full_simultaneous_crt_reconstruct!(state, luckyprimes)
+    full_simultaneous_crt_reconstruct_changematrix!(state, luckyprimes)
+
+    success_reconstruct = full_rational_reconstruct!(state, luckyprimes, params.use_flint)
+    @log :debug "Reconstructed coefficients" state.gb_coeffs_qq
+    @log :misc "Successfull reconstruction: $success_reconstruct"
+
+    changematrix_success_reconstruct =
+        full_rational_reconstruct_changematrix!(state, luckyprimes, params.use_flint)
+    @log :misc "Successfull change matrix reconstruction: $success_reconstruct"
+
+    correct_basis = false
+    if success_reconstruct && changematrix_success_reconstruct
+        @log :misc "Verifying the correctness of reconstruction"
+        correct_basis = correctness_check!(
+            state,
+            luckyprimes,
+            ring_ff,
+            basis,
+            basis_zz,
+            basis_ff,
+            hashtable,
+            params
+        )
+        @log :misc "Passed correctness check: $correct_basis"
+        # At this point, if the constructed basis is correct, we return it.
+        if correct_basis
+            # take monomials from the basis modulo a prime
+            gb_monoms, _ = basis_export_data(basis_ff, hashtable)
+            # take coefficients from the reconstrcted basis
+            gb_coeffs_qq = state.gb_coeffs_qq
+            changematrix_coeffs_qq = state.changematrix_coeffs_qq
+            return gb_monoms, gb_coeffs_qq, changematrix_monoms, changematrix_coeffs_qq
+        end
+    end
+
+    # At this point, either the reconstruction or the correctness check failed.
+    # Continue to compute Groebner bases modulo different primes in batches. 
+    primes_used = 1
+    batchsize = 1
+    batchsize_scaling = 0.10
+    @log :misc """
+    Preparing to compute bases in batches.. 
+    The initial size of the batch is $batchsize. 
+    The batch scale factor is $batchsize_scaling."""
+
+    if !tracer.ready_to_use
+        @log :misc """
+          The tracer is disabled until the shape of the basis is not determined via majority vote.
+          The threshold for the majority vote is $(params.majority_threshold)
+          """
+    end
+
+    # CRT and rational reconstrction settings
+    indices_selection = Vector{Tuple{Int, Int}}(undef, length(state.gb_coeffs_zz))
+    k = 1
+    for i in 1:length(state.gb_coeffs_zz)
+        l = length(state.gb_coeffs_zz[i])
+        nl = max(isqrt(l) - 1, 1)
+        if isone(l)
+            continue
+        end
+        for j in 1:nl
+            if k > length(indices_selection)
+                resize!(indices_selection, 2 * length(indices_selection))
+            end
+            indices_selection[k] = (i, rand(2:l))
+            k += 1
+        end
+    end
+    resize!(indices_selection, k - 1)
+    unique!(indices_selection)
+
+    partial_simultaneous_crt_reconstruct!(state, luckyprimes, indices_selection)
+
+    iters = 0
+    while !correct_basis
+        @log :misc "Iteration # $iters of modular Groebner"
+
+        for j in 1:batchsize
+            prime = next_lucky_prime!(luckyprimes)
+            @log :debug "The lucky prime is $prime"
+            @log :debug "Reducing input generators modulo $prime"
+
+            # Perform reduction modulo prime and store result in basis_ff
+            ring_ff, basis_ff =
+                reduce_modulo_p!(state.buffer, ring, basis_zz, prime, deepcopy=true)
+            params_zp = params_mod_p(params, prime)
+
+            f4!(ring_ff, basis_ff, pairset, hashtable, tracer, params_zp)
+
+            push!(state.gb_coeffs_ff_all, basis_ff.coeffs)
+            _, changematrix_coeffs =
+                basis_changematrix_export(basis_ff, hashtable, length(monoms))
+            push!(state.changematrix_coeffs_ff_all, changematrix_coeffs)
+
+            if !majority_vote!(state, basis_ff, tracer, params)
+                @log :debug "Majority vote is not conclusive, aborting reconstruction!"
+                continue
+            end
+
+            if crt_algorithm === :incremental
+                partial_incremental_crt_reconstruct!(state, luckyprimes, indices_selection)
+            end
+            primes_used += 1
+        end
+        if crt_algorithm === :simultaneous
+            partial_simultaneous_crt_reconstruct!(state, luckyprimes, indices_selection)
+        end
+
+        @log :misc "Partially reconstructing coefficients to QQ"
+        @log :debug "Partially reconstructing coefficients from Z_$(luckyprimes.modulo * prime) to QQ"
+        success_reconstruct = partial_rational_reconstruct!(
+            state,
+            luckyprimes,
+            indices_selection,
+            params.use_flint
+        )
+        @log :misc "Partial reconstruction successfull: $success_reconstruct"
+        @log :misc """
+          Used $(length(luckyprimes.primes)) primes in total over $(iters + 1) iterations.
+          The current batch size is $batchsize.
+          """
+
+        if !success_reconstruct
+            @log :misc "Partial rational reconstruction failed"
+            iters += 1
+            batchsize = get_next_batchsize(primes_used, batchsize, batchsize_scaling)
+            continue
+        end
+
+        if params.heuristic_check
+            success_check =
+                heuristic_correctness_check(state.selected_coeffs_qq, luckyprimes.modulo)
+            if !success_check
+                @log :misc "Heuristic check failed for partial reconstruction"
+                iters += 1
+                batchsize = get_next_batchsize(primes_used, batchsize, batchsize_scaling)
+                continue
+            end
+        end
+
+        # Perform full reconstruction
+        @log :misc "Performing full CRT and rational reconstruction.."
+        if crt_algorithm === :incremental
+            full_incremental_crt_reconstruct!(state, luckyprimes)
+        else
+            full_simultaneous_crt_reconstruct!(state, luckyprimes)
+        end
+        success_reconstruct =
+            full_rational_reconstruct!(state, luckyprimes, params.use_flint)
+
+        if !success_reconstruct
+            @log :misc "Full reconstruction failed"
+            iters += 1
+            batchsize = get_next_batchsize(primes_used, batchsize, batchsize_scaling)
+            continue
+        end
+
+        full_simultaneous_crt_reconstruct_changematrix!(state, luckyprimes)
+        success_reconstruct_changematrix =
+            full_rational_reconstruct_changematrix!(state, luckyprimes, params.use_flint)
+        if !success_reconstruct_changematrix
+            @log :misc "Failed to reconstruct the change matrix"
+            iters += 1
+            batchsize = get_next_batchsize(primes_used, batchsize, batchsize_scaling)
+            continue
+        end
+
+        correct_basis = correctness_check!(
+            state,
+            luckyprimes,
+            ring_ff,
+            basis,
+            basis_zz,
+            basis_ff,
+            hashtable,
+            params
+        )
+
+        iters += 1
+        batchsize = get_next_batchsize(primes_used, batchsize, batchsize_scaling)
+    end
+
+    @log :misc "Correctness check passed!"
+    @log :misc "Used $(length(luckyprimes.primes)) primes in total over $(iters) iterations"
+
+    # Construct the output basis.
+    # Take monomials from the basis modulo a prime
+    gb_monoms, _ = basis_export_data(basis_ff, hashtable)
+    # Take coefficients from the reconstructed basis
+    gb_coeffs_qq = state.gb_coeffs_qq
+    changematrix_coeffs_qq = state.changematrix_coeffs_qq
+
+    return gb_monoms, gb_coeffs_qq, changematrix_monoms, changematrix_coeffs_qq
+end
diff --git a/src/groebner/isgroebner.jl b/src/groebner/isgroebner.jl
index 830c49d6..185bde2d 100644
--- a/src/groebner/isgroebner.jl
+++ b/src/groebner/isgroebner.jl
@@ -3,7 +3,7 @@
 ### 
 # Backend for `isgroebner`
 
-function _isgroebner0(polynomials, kws::KeywordsHandler)
+function _isgroebner0(polynomials, kws::KeywordArguments)
     polynomial_repr =
         io_select_polynomial_representation(polynomials, kws, hint=:large_exponents)
     ring, var_to_index, monoms, coeffs =
diff --git a/src/groebner/learn-apply.jl b/src/groebner/learn-apply.jl
index 33616d04..f0487aa9 100644
--- a/src/groebner/learn-apply.jl
+++ b/src/groebner/learn-apply.jl
@@ -9,7 +9,7 @@
 # Proxy function for handling exceptions.
 # NOTE: probably at some point we'd want to merge this with error handling in
 # _groebner. But for now, we keep it simple.
-function _groebner_learn0(polynomials, kws::KeywordsHandler)
+function _groebner_learn0(polynomials, kws::KeywordArguments)
     # We try to select an efficient internal polynomial representation, i.e., a
     # suitable representation of monomials and coefficients.
     polynomial_repr = io_select_polynomial_representation(polynomials, kws)
@@ -92,7 +92,7 @@ end
 function _groebner_apply0!(
     wrapped_trace::WrappedTraceF4,
     polynomials::AbstractVector,
-    kws::KeywordsHandler
+    kws::KeywordArguments
 )
     trace = get_trace!(wrapped_trace, polynomials, kws)
     @log :debug "Selected trace" trace.representation.coefftype
@@ -124,7 +124,7 @@ end
 function _groebner_apply0!(
     wrapped_trace::WrappedTraceF4,
     batch::NTuple{N, T},
-    kws::KeywordsHandler
+    kws::KeywordArguments
 ) where {N, T <: AbstractVector}
     trace = get_trace!(wrapped_trace, batch, kws)
     @log :debug "Selected trace" trace.representation.coefftype
@@ -185,7 +185,7 @@ function groebner_applyX!(
     modulo::UInt32;
     options...
 )
-    kws = KeywordsHandler(:groebner_apply!, options)
+    kws = KeywordArguments(:groebner_apply!, options)
 
     logging_setup(kws)
     statistics_setup(kws)
diff --git a/src/groebner/normalform.jl b/src/groebner/normalform.jl
index 95998860..e783befc 100644
--- a/src/groebner/normalform.jl
+++ b/src/groebner/normalform.jl
@@ -3,7 +3,7 @@
 ###
 # Backend for `normalform`
 
-function _normalform0(polynomials, to_be_reduced, kws::KeywordsHandler)
+function _normalform0(polynomials, to_be_reduced, kws::KeywordArguments)
     polynomial_repr =
         io_select_polynomial_representation(polynomials, kws, hint=:large_exponents)
     ring, var_to_index1, monoms, coeffs =
diff --git a/src/groebner/parameters.jl b/src/groebner/parameters.jl
index 7f656bcf..1f41b22d 100644
--- a/src/groebner/parameters.jl
+++ b/src/groebner/parameters.jl
@@ -3,7 +3,7 @@
 ### 
 # Select parameters in Groebner basis computation
 
-# It seems there is no Xoshiro rng in Julia v < 1.8.
+# There is no Xoshiro rng in Julia v < 1.8.
 # Use Random.Xoshiro, if available, as it is a bit faster.
 const _default_rng_type = @static if VERSION >= v"1.8.0"
     Random.Xoshiro
@@ -11,7 +11,7 @@ else
     Random.MersenneTwister
 end
 
-# Specifies linear backend algorithm
+# Specifies linear algebra backend algorithm
 struct LinearAlgebra
     # One of :deterministic, :randomized, :experimental_1, :experimental_2,
     # :experimental_3,
@@ -25,14 +25,14 @@ struct LinearAlgebra
 end
 
 # Stores parameters for a single GB computation.
-# NOTE: in principle, MonomOrd1, ..., MonomOrd3 can be subtypes of any type
-# besides the usual Groebner.AbstractInternalOrdering
 mutable struct AlgorithmParameters{
     MonomOrd1,
     MonomOrd2,
     MonomOrd3,
     Arithmetic <: AbstractArithmetic
 }
+    # NOTE: in principle, MonomOrd1, ..., MonomOrd3 can be subtypes of any type
+
     # Desired monomial ordering of output polynomials
     target_ord::MonomOrd1
     # Monomial ordering for the actual computation
@@ -109,12 +109,14 @@ mutable struct AlgorithmParameters{
     statistics::Symbol
 
     use_flint::Bool
+
+    changematrix::Bool
 end
 
 function AlgorithmParameters(
     ring,
     representation,
-    kwargs::KeywordsHandler;
+    kwargs::KeywordArguments;
     orderings=nothing
 )
     # TODO: we should probably document this better
@@ -250,6 +252,16 @@ function AlgorithmParameters(
 
     use_flint = kwargs.use_flint
 
+    changematrix = kwargs.changematrix
+    if changematrix
+        if !(target_ord isa DegRevLex)
+            __throw_input_not_supported(
+                "Only DegRevLex is supported with changematrix = true.",
+                target_ord
+            )
+        end
+    end
+
     @log :misc """
     Selected parameters:
     target_ord = $target_ord
@@ -277,7 +289,8 @@ function AlgorithmParameters(
     rng = $rng
     sweep = $sweep
     statistics = $statistics
-    use_flint = $use_flint"""
+    use_flint = $use_flint
+    changematrix = $changematrix"""
 
     AlgorithmParameters(
         target_ord,
@@ -305,7 +318,8 @@ function AlgorithmParameters(
         rng,
         sweep,
         statistics,
-        use_flint
+        use_flint,
+        changematrix
     )
 end
 
@@ -345,6 +359,7 @@ function params_mod_p(
         params.rng,
         params.sweep,
         params.statistics,
-        params.use_flint
+        params.use_flint,
+        params.changematrix
     )
 end
diff --git a/src/groebner/state.jl b/src/groebner/state.jl
index 7c4ae251..ef9b5330 100644
--- a/src/groebner/state.jl
+++ b/src/groebner/state.jl
@@ -48,11 +48,9 @@ mutable struct GroebnerState{T1 <: CoeffZZ, T2 <: CoeffQQ, T3}
     prev_gb_coeffs_zz::Vector{Vector{T1}}
     gb_coeffs_qq::Vector{Vector{T2}}
 
-    #
     gb_coeffs_ff_all::Vector{Vector{Vector{T3}}}
     prev_index::Int
 
-    #
     selected_coeffs_zz::Vector{T1}
     selected_prev_coeffs_zz::Vector{T1}
     selected_coeffs_qq::Vector{T2}
@@ -60,6 +58,11 @@ mutable struct GroebnerState{T1 <: CoeffZZ, T2 <: CoeffQQ, T3}
     is_rational_reconstructed_mask::Vector{BitVector}
 
     buffer::CoefficientBuffer
+
+    changematrix_coeffs_ff_all::Vector{Vector{Vector{Vector{T3}}}}
+    changematrix_coeffs_zz::Vector{Vector{Vector{T1}}}
+    changematrix_coeffs_qq::Vector{Vector{Vector{T2}}}
+
     function GroebnerState{T1, T2, T3}(
         params
     ) where {T1 <: CoeffZZ, T2 <: CoeffQQ, T3 <: CoeffZp}
@@ -74,7 +77,10 @@ mutable struct GroebnerState{T1 <: CoeffZZ, T2 <: CoeffQQ, T3}
             Vector{T2}(),
             Vector{BitVector}(),
             Vector{BitVector}(),
-            CoefficientBuffer()
+            CoefficientBuffer(),
+            Vector{Vector{Vector{Vector{T3}}}}(),
+            Vector{Vector{Vector{T1}}}(),
+            Vector{Vector{Vector{T2}}}()
         )
     end
 end
@@ -441,6 +447,89 @@ function full_simultaneous_crt_reconstruct!(state::GroebnerState, lucky::LuckyPr
     nothing
 end
 
+function full_simultaneous_crt_reconstruct_changematrix!(
+    state::GroebnerState,
+    lucky::LuckyPrimes
+)
+    if isempty(state.changematrix_coeffs_zz)
+        @log :misc "Using full trivial CRT reconstruction"
+        coeffs_ff = state.changematrix_coeffs_ff_all[1]
+
+        resize!(state.changematrix_coeffs_zz, length(coeffs_ff))
+        resize!(state.changematrix_coeffs_qq, length(coeffs_ff))
+
+        @inbounds for i in 1:length(coeffs_ff)
+            state.changematrix_coeffs_zz[i] =
+                Vector{Vector{BigInt}}(undef, length(coeffs_ff[i]))
+            state.changematrix_coeffs_qq[i] =
+                Vector{Vector{Rational{BigInt}}}(undef, length(coeffs_ff[i]))
+            for j in 1:length(coeffs_ff[i])
+                state.changematrix_coeffs_zz[i][j] =
+                    [BigInt(0) for _ in 1:length(coeffs_ff[i][j])]
+                state.changematrix_coeffs_qq[i][j] =
+                    [Rational{BigInt}(1) for _ in 1:length(coeffs_ff[i][j])]
+            end
+        end
+
+        changematrix_coeffs_zz = state.changematrix_coeffs_zz
+        @inbounds for i in 1:length(coeffs_ff)
+            for j in 1:length(coeffs_ff[i])
+                for k in 1:length(coeffs_ff[i][j])
+                    Base.GMP.MPZ.set_ui!(
+                        changematrix_coeffs_zz[i][j][k],
+                        coeffs_ff[i][j][k]
+                    )
+                end
+            end
+        end
+        # Base.GMP.MPZ.mul_ui!(lucky.modulo, lucky.primes[1])
+        return nothing
+    end
+    # @invariant length(coeffs_ff) == length(state.changematrix_coeffs_zz)
+
+    @log :misc "Using full CRT reconstruction"
+    changematrix_coeffs_zz = state.changematrix_coeffs_zz
+
+    # Takes the lock..
+    @invariant length(state.changematrix_coeffs_ff_all) == length(lucky.primes)
+
+    # Set the buffers for CRT and precompute some values
+    buffer = state.buffer
+    buf = buffer.reconstructbuf1
+    n1, n2 = buffer.reconstructbuf2, buffer.reconstructbuf3
+    M = buffer.reconstructbuf4
+    invm1, invm2 = buffer.reconstructbuf6, buffer.reconstructbuf7
+    M0 = buffer.reconstructbuf8
+    MM0 = buffer.reconstructbuf9
+
+    n = length(lucky.primes)
+    rems = Vector{CoeffModular}(undef, n)
+    mults = Vector{BigInt}(undef, n)
+    for i in 1:length(mults)
+        mults[i] = BigInt(0)
+    end
+    moduli = lucky.primes
+    crt_precompute!(M, n1, n2, mults, moduli)
+
+    @log :debug "Using simultaneous CRT with moduli $moduli"
+
+    @inbounds for i in 1:length(changematrix_coeffs_zz)
+        for j in 1:length(changematrix_coeffs_zz[i])
+            for k in 1:length(changematrix_coeffs_zz[i][j])
+                for ell in 1:length(lucky.primes)
+                    rems[ell] = state.changematrix_coeffs_ff_all[ell][i][j][k]
+                end
+                crt!(M, buf, n1, n2, rems, mults)
+                Base.GMP.MPZ.set!(changematrix_coeffs_zz[i][j][k], buf)
+            end
+        end
+    end
+
+    # Base.GMP.MPZ.set!(lucky.modulo, M)
+
+    nothing
+end
+
 @timeit function partial_incremental_crt_reconstruct!(
     state::GroebnerState,
     lucky::LuckyPrimes,
@@ -679,6 +768,48 @@ end
     true
 end
 
+function full_rational_reconstruct_changematrix!(
+    state::GroebnerState,
+    lucky::LuckyPrimes,
+    use_flint::Bool
+)
+    modulo = lucky.modulo
+    @invariant modulo == prod(BigInt, lucky.primes)
+
+    buffer = state.buffer
+    bnd = ratrec_reconstruction_bound(modulo)
+
+    buf, buf1 = buffer.reconstructbuf1, buffer.reconstructbuf2
+    buf2, buf3 = buffer.reconstructbuf3, buffer.reconstructbuf4
+    u1, u2 = buffer.reconstructbuf5, buffer.reconstructbuf6
+    u3, v1 = buffer.reconstructbuf7, buffer.reconstructbuf8
+    v2, v3 = buffer.reconstructbuf9, buffer.reconstructbuf10
+    changematrix_coeffs_zz = state.changematrix_coeffs_zz
+    changematrix_coeffs_qq = state.changematrix_coeffs_qq
+    # is_rational_reconstructed_mask = state.is_rational_reconstructed_mask
+
+    @invariant length(changematrix_coeffs_zz) == length(changematrix_coeffs_qq)
+    @assert use_flint
+    nemo_modulo = Nemo.ZZRingElem(modulo)
+
+    @inbounds for i in 1:length(changematrix_coeffs_zz)
+        @invariant length(changematrix_coeffs_zz[i]) == length(changematrix_coeffs_qq[i])
+        # !!!
+        for j in 1:length(changematrix_coeffs_zz[i])
+            for k in 1:length(changematrix_coeffs_zz[i][j])
+                cz = changematrix_coeffs_zz[i][j][k]
+                nemo_rem = Nemo.ZZRingElem(cz)
+                success, (num, den) = ratrec_nemo(nemo_rem, nemo_modulo)
+                changematrix_coeffs_qq[i][j][k] = Base.unsafe_rational(num, den)
+
+                !success && return false
+            end
+        end
+    end
+
+    true
+end
+
 @timeit function partial_rational_reconstruct!(
     state::GroebnerState,
     lucky::LuckyPrimes,
diff --git a/src/input-output/AbstractAlgebra.jl b/src/input-output/AbstractAlgebra.jl
index d76448dc..fdcf885b 100644
--- a/src/input-output/AbstractAlgebra.jl
+++ b/src/input-output/AbstractAlgebra.jl
@@ -239,7 +239,7 @@ function extract_coeffs_raw!(
     trace,
     representation::PolynomialRepresentation,
     polys::Vector{T},
-    kws::KeywordsHandler
+    kws::KeywordArguments
 ) where {T}
     ring = extract_ring(polys)
     !is_input_compatible_in_apply(trace, ring, polys, kws) && __throw_input_not_supported(
@@ -290,7 +290,7 @@ function extract_coeffs_raw_X!(
     representation::PolynomialRepresentation,
     coeffs_zp,
     modulo,
-    kws::KeywordsHandler
+    kws::KeywordArguments
 )
     ring = PolyRing(trace.ring.nvars, trace.ring.ord, UInt64(modulo))
 
@@ -380,7 +380,7 @@ function io_extract_coeffs_raw_batched!(
     trace,
     representation::PolynomialRepresentation,
     batch::NTuple{N, T},
-    kws::KeywordsHandler
+    kws::KeywordArguments
 ) where {N, T <: AbstractVector}
     rings = map(extract_ring, batch)
     chars = (representation.coefftype)(map(ring -> ring.ch, rings))
diff --git a/src/input-output/input-output.jl b/src/input-output/input-output.jl
index 09a6e282..ae3fb805 100644
--- a/src/input-output/input-output.jl
+++ b/src/input-output/input-output.jl
@@ -78,7 +78,7 @@ Additionally, `hint` can be specified to one of the following:
 """
 function io_select_polynomial_representation(
     polynomials::AbstractVector,
-    kws::KeywordsHandler;
+    kws::KeywordArguments;
     hint::Symbol=:none
 )
     if !(hint in (:none, :large_exponents))
@@ -288,7 +288,7 @@ representation specified by the given `representation`.
 @timeit function io_convert_to_internal(
     representation::PolynomialRepresentation,
     polynomials,
-    kws::KeywordsHandler;
+    kws::KeywordArguments;
     dropzeros=true
 )
     io_check_input(polynomials, kws)
@@ -406,6 +406,34 @@ reference).
     _io_convert_to_output(ring, polynomials, monoms, coeffs, params)
 end
 
+function io_convert_changematrix_to_output(
+    ring::PolyRing,
+    polynomials,
+    npolys::Int,
+    monoms::Vector{Vector{Vector{M}}},
+    coeffs::Vector{Vector{Vector{C}}},
+    params::AlgorithmParameters
+) where {M <: Monom, C <: Coeff}
+    @assert !isempty(polynomials)
+    @log :misc "Converting polynomials from internal representation to output format"
+    changematrix = Matrix{eltype(polynomials)}(undef, length(monoms), length(polynomials))
+    for i in 1:length(monoms)
+        matrix_row = io_convert_to_output(ring, polynomials, monoms[i], coeffs[i], params)
+        matrix_row_full = Vector{eltype(matrix_row)}(undef, length(polynomials))
+        k = 1
+        for j in 1:length(polynomials)
+            if iszero(polynomials[j])
+                matrix_row_full[j] = zero(polynomials[j])
+            else
+                matrix_row_full[j] = matrix_row[k]
+                k += 1
+            end
+        end
+        changematrix[i, :] .= matrix_row_full
+    end
+    changematrix
+end
+
 @timeit function io_convert_to_output_batched(
     ring::PolyRing,
     batch::NTuple{N, V},
diff --git a/src/interface.jl b/src/interface.jl
index 36844103..b3a6d8bd 100644
--- a/src/interface.jl
+++ b/src/interface.jl
@@ -117,7 +117,7 @@ The function `groebner` is **thread-safe**.
 function groebner(polynomials::AbstractVector; options...)
     Base.require_one_based_indexing(polynomials)
 
-    keywords = KeywordsHandler(:groebner, options)
+    keywords = KeywordArguments(:groebner, options)
 
     logging_setup(keywords)
     statistics_setup(keywords)
@@ -132,6 +132,22 @@ function groebner(polynomials::AbstractVector; options...)
     result
 end
 
+function groebner_with_change_matrix(polynomials::AbstractVector; options...)
+    Base.require_one_based_indexing(polynomials)
+
+    keywords = KeywordArguments(:groebner_with_change_matrix, options)
+
+    logging_setup(keywords)
+    statistics_setup(keywords)
+
+    result = _groebner_with_change_matrix0(polynomials, keywords)
+
+    performance_counters_print(keywords)
+    statistics_print(keywords)
+
+    result
+end
+
 """
     groebner_learn(polynomials; options...)
 
@@ -250,7 +266,7 @@ The function `groebner_learn` is **thread-safe**.
 function groebner_learn(polynomials::AbstractVector; options...)
     Base.require_one_based_indexing(polynomials)
 
-    keywords = KeywordsHandler(:groebner_learn, options)
+    keywords = KeywordArguments(:groebner_learn, options)
 
     logging_setup(keywords)
     statistics_setup(keywords)
@@ -298,7 +314,7 @@ function groebner_apply! end
 function groebner_apply!(trace, polynomials::AbstractVector; options...)
     Base.require_one_based_indexing(polynomials)
 
-    keywords = KeywordsHandler(:groebner_apply!, options)
+    keywords = KeywordArguments(:groebner_apply!, options)
 
     logging_setup(keywords)
     statistics_setup(keywords)
@@ -320,7 +336,7 @@ function groebner_apply!(
     @assert N in (1, 2, 4, 8, 16) "The batch size must be one of the following: 1, 2, 4, 8, 16"
     all(Base.require_one_based_indexing, batch)
 
-    keywords = KeywordsHandler(:groebner_apply!, options)
+    keywords = KeywordArguments(:groebner_apply!, options)
 
     logging_setup(keywords)
     statistics_setup(keywords)
@@ -396,7 +412,7 @@ The function `isgroebner` is **thread-safe**.
 function isgroebner(polynomials::AbstractVector; options...)
     Base.require_one_based_indexing(polynomials)
 
-    keywords = KeywordsHandler(:isgroebner, options)
+    keywords = KeywordArguments(:isgroebner, options)
 
     logging_setup(keywords)
     statistics_setup(keywords)
@@ -469,7 +485,7 @@ function normalform(basis::AbstractVector, to_be_reduced::AbstractVector; option
     Base.require_one_based_indexing(basis)
     Base.require_one_based_indexing(to_be_reduced)
 
-    keywords = KeywordsHandler(:normalform, options)
+    keywords = KeywordArguments(:normalform, options)
 
     logging_setup(keywords)
     statistics_setup(keywords)
@@ -523,7 +539,7 @@ The function `kbase` is **thread-safe**.
 function kbase(basis::AbstractVector; options...)
     Base.require_one_based_indexing(basis)
 
-    keywords = KeywordsHandler(:kbase, options)
+    keywords = KeywordArguments(:kbase, options)
 
     logging_setup(keywords)
     statistics_setup(keywords)
@@ -551,7 +567,7 @@ function fglm(
 )
     Base.require_one_based_indexing(basis)
 
-    keywords = KeywordsHandler(:fglm, options)
+    keywords = KeywordArguments(:fglm, options)
 
     logging_setup(keywords)
     statistics_setup(keywords)
diff --git a/src/utils/keywords.jl b/src/utils/keywords.jl
index bed67f78..8d92cc1d 100644
--- a/src/utils/keywords.jl
+++ b/src/utils/keywords.jl
@@ -9,23 +9,24 @@
 # arguments with their corresponding default values.
 const _supported_kw_args = (
     groebner = (
-        reduced     = true,
-        ordering    = InputOrdering(),
-        certify     = false,
-        linalg      = :auto,
-        monoms      = :auto,
-        arithmetic  = :auto,
-        seed        = 42,
-        loglevel    = _loglevel_default,
-        maxpairs    = typemax(Int),   # NOTE: maybe use Inf?
-        selection   = :auto,
-        modular     = :auto,
-        threaded    = :auto,
-        sweep       = false,
-        homogenize  = :auto,
-        statistics  = :no,
-        batched     = true,
-        use_flint   = true
+        reduced      = true,
+        ordering     = InputOrdering(),
+        certify      = false,
+        linalg       = :auto,
+        monoms       = :auto,
+        arithmetic   = :auto,
+        seed         = 42,
+        loglevel     = _loglevel_default,
+        maxpairs     = typemax(Int),   # NOTE: maybe use Inf?
+        selection    = :auto,
+        modular      = :auto,
+        threaded     = :auto,
+        sweep        = false,
+        homogenize   = :auto,
+        statistics   = :no,
+        batched      = true,
+        use_flint    = true,
+        changematrix = false
     ),
     normalform = (
         check       = false,
@@ -75,15 +76,42 @@ const _supported_kw_args = (
         sweep       = true,
         statistics  = :no,
         threaded    = :auto,
-    )
+    ),
+    groebner_with_change_matrix = (
+        reduced      = true,
+        ordering     = InputOrdering(),
+        certify      = false,
+        linalg       = :auto,
+        monoms       = :auto,
+        arithmetic   = :auto,
+        seed         = 42,
+        loglevel     = _loglevel_default,
+        maxpairs     = typemax(Int),   # NOTE: maybe use Inf?
+        selection    = :auto,
+        modular      = :auto,
+        threaded     = :auto,
+        sweep        = false,
+        homogenize   = :auto,
+        statistics   = :no,
+        batched      = true,
+        use_flint    = true,
+        changematrix = true
+    ),
 )
 #! format: on
 
 """
-    KeywordsHandler
+    KeywordArguments
+
+Stores keyword arguments passed to a function in the interface in Groebner.jl.
+    
+Can be manually created with 
 
-Stores keyword arguments passed to one of the functions in the interface."""
-struct KeywordsHandler{Ord}
+```julia
+kwargs = Groebner.KeywordArguments(:groebner, seed = 99, reduced = false)
+```
+"""
+struct KeywordArguments{Ord}
     reduced::Bool
     ordering::Ord
     certify::Bool
@@ -102,8 +130,12 @@ struct KeywordsHandler{Ord}
     homogenize::Symbol
     statistics::Symbol
     use_flint::Bool
+    changematrix::Bool
+
+    KeywordArguments(function_key::Symbol; passthrough_options...) =
+        KeywordArguments(function_key, passthrough_options)
 
-    function KeywordsHandler(function_key, kws)
+    function KeywordArguments(function_key::Symbol, kws)
         @assert haskey(_supported_kw_args, function_key)
         default_kw_args = _supported_kw_args[function_key]
         for (key, _) in kws
@@ -189,6 +221,8 @@ struct KeywordsHandler{Ord}
         Possible choices for keyword "statistics" are:
         `:no`, `:timings`, `:stats`, `:all`"""
 
+        changematrix = get(kws, :changematrix, get(default_kw_args, :changematrix, false))
+
         new{typeof(ordering)}(
             reduced,
             ordering,
@@ -207,17 +241,18 @@ struct KeywordsHandler{Ord}
             sweep,
             homogenize,
             statistics,
-            use_flint
+            use_flint,
+            changematrix
         )
     end
 end
 
-function logging_setup(keywords::KeywordsHandler)
+function logging_setup(keywords::KeywordArguments)
     logger_update(loglevel=keywords.loglevel)
     nothing
 end
 
-function statistics_setup(keywords::KeywordsHandler)
+function statistics_setup(keywords::KeywordArguments)
     log_simdinfo()
     if keywords.loglevel <= 0
         performance_counters_refresh()
diff --git a/test/groebner/groebner_with_change_matrix.jl b/test/groebner/groebner_with_change_matrix.jl
new file mode 100644
index 00000000..c4b7090d
--- /dev/null
+++ b/test/groebner/groebner_with_change_matrix.jl
@@ -0,0 +1,44 @@
+
+@testset "groebner, change matrix" begin
+    R, (x, y, z) =
+        polynomial_ring(GF(2^30 + 3), ["x", "y", "z"], internal_ordering=:degrevlex)
+    f = [x * y * z - 1, x * y + x * z + y * z, x + y + z]
+    g, m = Groebner.groebner_with_change_matrix(f)
+    @test m * f == g
+    @test size(m) == (length(g), length(f))
+    @test Groebner.isgroebner(g)
+
+    for ground in [GF(2^30 + 3), QQ]
+        R, (x, y) = polynomial_ring(ground, ["x", "y"], internal_ordering=:degrevlex)
+
+        cases = [
+            [R(0), R(0), R(0)],
+            [R(0), R(1), R(0)],
+            [x],
+            [y],
+            [x, x, x, x, R(0), R(0), x, x, x, x],
+            [x + y, R(1), x^5 - y^5],
+            [x^200 + y^250, x^100 * y^200 + 1],
+            Groebner.katsuran(4, k=ground),
+            Groebner.noonn(4, k=ground),
+            Groebner.cyclicn(4, k=ground),
+            [x^i * y + x for i in 1:500]
+        ]
+
+        for f in cases
+            g, m = Groebner.groebner_with_change_matrix(f)
+            @test m * f == g
+            @test size(m) == (length(g), length(f))
+            @test Groebner.isgroebner(g)
+        end
+    end
+
+    f = Groebner.katsuran(4, k=QQ, internal_ordering=:lex)
+    g, m = Groebner.groebner_with_change_matrix(f, ordering=Groebner.DegRevLex())
+    @test m * f == g
+    @test_throws DomainError Groebner.groebner_with_change_matrix(
+        f,
+        ordering=Groebner.DegLex()
+    )
+    @test_throws DomainError Groebner.groebner_with_change_matrix(f)
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index aa9ee1b3..972c81a0 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -55,6 +55,8 @@ end
         "groebner/multi_threading"
     ]
 
+    @time @includetests ["groebner/groebner_with_change_matrix"]
+
     @time @includetests [
         "learn_and_apply/learn_and_apply",
         "learn_and_apply/apply_in_batches"