rref_finite_field.jl

function rref!(A::Matrix{Int},
    n::Integer;
    colswap::Bool=false,
    verbose::Bool=false)::Matrix{Int}
    
    nrows, ncols = size(A)
    i = j = 1
    
    while i ≤ nrows && j ≤ ncols
        # Ignore zero rows.
        # if isnothing(findfirst(!iszero, A[i, :]))
        #     i += 1
        #     continue
        # end
        for k in i+1:nrows
            if isnothing(findfirst(!iszero, A[k, :]))
                continue
            end
            A[i, :], A[k, :] = A[k, :], A[i, :]
            if verbose
                println("r$(i) ⟷  r$(k)")
                __displaymatrix(A);println()
            end
            break
        end
        
        # Rule 1: Swap with the row above if out of order.
        # if i > 1
        #     if isnothing(findfirst(!iszero, A[i - 1, :]))
        #         A[i,:], A[i - 1,:] = A[i - 1,:], A[i,:]
        #         if verbose
        #             println("r$(i) ⟷  r$(i - 1)")
        #             __displaymatrix(A);println()
        #         end
        #         continue
        #     end
        # end
        
        # Rule 2: Normalize each row
        s = findfirst(!iszero, A[i,:])
        α = invmod(A[i, s], n)
        A[i,:] .= mod.(A[i,:] * α, n)
        if verbose
            if ! iszero(α)
                isone(α) || (println("r$(i) ⟵  r$(i) × $(α)"); __displaymatrix(A);println())
            end
        end
        
        # Rule 3: Subtract it from the others
        s = findfirst(!iszero, A[i,:])
        isnothing(s) && break
        for k in 1:nrows
            if i ≠ k
                β = A[k, s]
                A[k,:] .= mod.(A[k,:] - β * A[i,:], n)
                if verbose
                    if ! iszero(β)
                        if isone(β)
                            println("r$(k) ⟵  r$(k) - r$(i)")
                            __displaymatrix(A);println()
                        else
                            println("r$(k) ⟵  r$(k) - $(β)r$(i)")
                            __displaymatrix(A);println()
                        end
                    end
                end
            end
        end
        # i += 1
        
        # Rule 4: Swap columns if needed (and allowed)
        if colswap
            if iszero(A[i, j])
                s = findfirst(!iszero, A[i,:])
                A[:,s], A[:,j] = A[:,j], A[:,s]
                if verbose
                    println("c$(j) ⟷ c$(s)")
                end
            end
        end
        
        # increment counters
        i += 1
        j += 1
    end
    
    if verbose
        println()
    end

    return A
end