diff --git a/src/PropertyT.jl b/src/PropertyT.jl
index 72e9247..2e3d9bf 100644
--- a/src/PropertyT.jl
+++ b/src/PropertyT.jl
@@ -5,6 +5,7 @@ using LinearAlgebra
 using SparseArrays
 using Dates
 
+using IntervalArithmetic
 using JuMP
 
 using Groups
@@ -14,7 +15,7 @@ using SymbolicWedderburn
 include("laplacians.jl")
 include("constraint_matrix.jl")
 include("sos_sdps.jl")
-include("checksolution.jl")
+include("certify.jl")
 
 include("1712.07167.jl")
 include("1812.03456.jl")
diff --git a/src/certify.jl b/src/certify.jl
new file mode 100644
index 0000000..bca18c6
--- /dev/null
+++ b/src/certify.jl
@@ -0,0 +1,168 @@
+function augment_columns!(Q::AbstractMatrix)
+    for c in eachcol(Q)
+        c .-= sum(c) ./ length(c)
+    end
+    return Q
+end
+
+function _fma_SOS_thr!(
+    result::AbstractVector{T},
+    mstructure::AbstractMatrix{<:Integer},
+    Q::AbstractMatrix{T},
+    acc_matrix=zeros(T, size(mstructure)...),
+) where {T}
+
+    s1, s2 = size(mstructure)
+
+    @inbounds for k = 1:s2
+        let k = k, s1 = s1, s2 = s2, Q = Q, acc_matrix = acc_matrix
+            Threads.@threads for j = 1:s2
+                for i = 1:s1
+                    @inbounds acc_matrix[i, j] =
+                        muladd(Q[i, k], Q[j, k], acc_matrix[i, j])
+                end
+            end
+        end
+    end
+
+    @inbounds for j = 1:s2
+        for i = 1:s1
+            result[mstructure[i, j]] += acc_matrix[i, j]
+        end
+    end
+
+    return result
+end
+
+function _cnstr_sos!(res::AlgebraElement, Q::AbstractMatrix, cnstrs)
+    StarAlgebras.zero!(res)
+    Q² = Q' * Q
+    for (g, A_g) in cnstrs
+        res[g] = dot(A_g, Q²)
+    end
+    return res
+end
+
+function _augmented_sos!(res::AlgebraElement, Q::AbstractMatrix)
+    A = parent(res)
+    StarAlgebras.zero!(res)
+    Q² = Q' * Q
+
+    N = LinearAlgebra.checksquare(A.mstructure)
+    augmented_basis = [A(1) - A(b) for b in @view basis(A)[1:N]]
+    tmp = zero(res)
+
+    for (j, y) in enumerate(augmented_basis)
+        for (i, x) in enumerate(augmented_basis)
+            # res += Q²[i, j] * x * y
+
+            StarAlgebras.mul!(tmp, x, y)
+            StarAlgebras.mul!(tmp, tmp, Q²[i, j])
+            StarAlgebras.add!(res, res, tmp)
+        end
+    end
+    return res
+end
+
+function compute_sos(A::StarAlgebra, Q::AbstractMatrix; augmented::Bool)
+    if augmented
+        z = zeros(eltype(Q), length(basis(A)))
+        res = AlgebraElement(z, A)
+        return _augmented_sos!(res, Q)
+        cnstrs = constraints(basis(A), A.mstructure; augmented=true)
+        return _cnstr_sos!(res, Q, cnstrs)
+    else
+        @assert size(A.mstructure) == size(Q)
+        z = zeros(eltype(Q), length(basis(A)))
+
+        _fma_SOS_thr!(z, A.mstructure, Q)
+
+        return AlgebraElement(z, A)
+    end
+end
+
+function sufficient_λ(residual::AlgebraElement, λ; halfradius)
+    L1_norm = norm(residual, 1)
+    suff_λ = λ - 2.0^(2ceil(log2(halfradius))) * L1_norm
+
+    eq_sign = let T = eltype(residual)
+        if T <: Interval
+            "∈"
+        elseif T <: Union{Rational,Integer}
+            "="
+        else # if T <: AbstractFloat
+            "≈"
+        end
+    end
+
+    info_strs = [
+        "Numerical metrics of the obtained SOS:",
+        "ɛ(elt - λu - ∑ξᵢ*ξᵢ) $eq_sign $(aug(residual))",
+        "‖elt - λu - ∑ξᵢ*ξᵢ‖₁ $eq_sign $(L1_norm)",
+        " λ $eq_sign $suff_λ",
+    ]
+    @info join(info_strs, "\n")
+
+    return suff_λ
+end
+
+function sufficient_λ(
+    elt::AlgebraElement,
+    order_unit::AlgebraElement,
+    λ,
+    sos::AlgebraElement;
+    halfradius
+)
+
+    @assert parent(elt) === parent(order_unit) == parent(sos)
+    residual = (elt - λ * order_unit) - sos
+
+    return sufficient_λ(residual, λ; halfradius=halfradius)
+end
+
+function certify_solution(
+    elt::AlgebraElement,
+    orderunit::AlgebraElement,
+    λ,
+    Q::AbstractMatrix{<:AbstractFloat};
+    halfradius,
+    augmented=iszero(aug(elt)) && iszero(aug(orderunit))
+)
+
+    should_we_augment = !augmented && aug(elt) == aug(orderunit) == 0
+
+    Q = should_we_augment ? augment_columns!(Q) : Q
+    @time sos = compute_sos(parent(elt), Q, augmented=augmented)
+
+    @info "Checking in $(eltype(sos)) arithmetic with" λ
+
+    λ_flpoint = sufficient_λ(elt, orderunit, λ, sos, halfradius=halfradius)
+
+    if λ_flpoint ≤ 0
+        return false, λ_flpoint
+    end
+
+    λ_int = @interval(λ)
+    Q_int = [@interval(q) for q in Q]
+
+    check, sos_int = @time if should_we_augment
+        @info("Projecting columns of Q to the augmentation ideal...")
+        Q_int = augment_columns!(Q_int)
+        @info "Checking that sum of every column contains 0.0..."
+        check_augmented = all(0 ∈ sum(c) for c in eachcol(Q_int))
+        check_augmented || @error(
+            "Augmentation failed! The following numbers are not certified!"
+        )
+        sos_int = compute_sos(parent(elt), Q_int; augmented=augmented)
+        check_augmented, sos_int
+    else
+        true, compute_sos(parent(elt), Q_int, augmented=augmented)
+    end
+
+    @info "Checking in $(eltype(sos_int)) arithmetic with" λ
+
+    λ_certified =
+        sufficient_λ(elt, orderunit, λ_int, sos_int, halfradius=halfradius)
+
+    return check && inf(λ_certified) > 0.0, inf(λ_certified)
+end
diff --git a/src/checksolution.jl b/src/checksolution.jl
deleted file mode 100644
index 8a4a166..0000000
--- a/src/checksolution.jl
+++ /dev/null
@@ -1,77 +0,0 @@
-using IntervalArithmetic
-
-IntervalArithmetic.setrounding(Interval, :tight)
-IntervalArithmetic.setformat(sigfigs=12)
-
-function fma_SOS_thr!(result::AbstractVector{T}, pm::AbstractMatrix{<:Integer},
-    Q::AbstractMatrix{T}, acc_matrix=zeros(T, size(pm)...)) where T
-
-    s1, s2 = size(pm)
-
-    @inbounds for k in 1:s2
-        let k=k, s1=s1, s2=s2, Q=Q, acc_matrix=acc_matrix
-            Threads.@threads for j in 1:s2
-                for i in 1:s1
-                    @inbounds acc_matrix[i,j] = muladd(Q[i, k], Q[j, k], acc_matrix[i,j])
-                end
-            end
-        end
-    end
-
-    @inbounds for j in 1:s2
-        for i in 1:s1
-            result[pm[i,j]] += acc_matrix[i,j]
-        end
-    end
-
-    return result
-end
-
-function compute_SOS(pm::AbstractMatrix{<:Integer}, Q::AbstractMatrix)
-    result = zeros(eltype(Q), maximum(pm));
-    return fma_SOS_thr!(result, pm, Q)
-end
-
-function compute_SOS(RG::GroupRing, Q::AbstractMatrix{<:Real})
-    result = compute_SOS(RG.pm, Q)
-    return GroupRingElem(result, RG)
-end
-
-function compute_SOS_square(pm::AbstractMatrix{<:Integer}, Q::AbstractMatrix)
-    result = zeros(eltype(Q), maximum(pm));
-
-    for i in 1:size(Q,2)
-        GroupRings.fmac!(result, view(Q,:,i), view(Q,:,i), pm)
-    end
-
-    return result
-end
-
-function compute_SOS_square(RG::GroupRing, Q::AbstractMatrix)
-    return GroupRingElem(compute_SOS_square(RG.pm, Q), RG)
-end
-
-function augIdproj(Q::AbstractMatrix{T}) where T
-    result = zeros(T, size(Q))
-    l = size(Q, 2)
-    Threads.@threads for j in 1:l
-        col = sum(view(Q, :,j))/l
-        for i in 1:size(Q, 1)
-            result[i,j] = Q[i,j] - col
-        end
-    end
-    return result
-end
-
-function augIdproj(::Type{Interval}, Q::AbstractMatrix{T}) where {T<:Real}
-    result = zeros(Interval{T}, size(Q))
-    l = size(Q, 2)
-    Threads.@threads for j in 1:l
-        col = sum(view(Q, :,j))/l
-        for i in 1:size(Q, 1)
-            result[i,j] = @interval(Q[i,j] - col)
-        end
-    end
-    check = all([zero(T) in sum(view(result, :, i)) for i in 1:size(result, 2)])
-    return result, check
-end