diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml
new file mode 100644
index 0000000..f196f10
--- /dev/null
+++ b/.JuliaFormatter.toml
@@ -0,0 +1,9 @@
+# Configuration file for JuliaFormatter.jl
+# For more information, see: https://domluna.github.io/JuliaFormatter.jl/stable/config/
+
+always_for_in = true
+always_use_return = true
+margin = 80
+remove_extra_newlines = true
+separate_kwargs_with_semicolon = true
+short_to_long_function_def = true
diff --git a/Project.toml b/Project.toml
index 8255f05..c13977e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,6 +6,7 @@ version = "0.4.0"
 [deps]
 Groups = "5d8bd718-bd84-11e8-3b40-ad14f4a32557"
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
+IntervalMatrices = "5c1f47dc-42dd-5697-8aaa-4d102d140ba9"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
@@ -17,11 +18,12 @@ SymbolicWedderburn = "858aa9a9-4c7c-4c62-b466-2421203962a2"
 COSMO = "0.8"
 Groups = "0.7"
 IntervalArithmetic = "0.20"
+IntervalMatrices = "0.8"
 JuMP = "1.3"
 ProgressMeter = "1.7"
-SCS = "1.1.0"
+SCS = "1.1"
 StaticArrays = "1"
-SymbolicWedderburn = "0.3.2"
+SymbolicWedderburn = "0.3.4"
 julia = "1.6"
 
 [extras]
diff --git a/scripts/G₂_Adj.jl b/scripts/G₂_Adj.jl
new file mode 100644
index 0000000..8a499d0
--- /dev/null
+++ b/scripts/G₂_Adj.jl
@@ -0,0 +1,98 @@
+using LinearAlgebra
+BLAS.set_num_threads(8)
+using MKL_jll
+ENV["OMP_NUM_THREADS"] = 4
+
+using Groups
+import Groups.MatrixGroups
+
+include(joinpath(@__DIR__, "../test/optimizers.jl"))
+using PropertyT
+
+using PropertyT.SymbolicWedderburn
+using PropertyT.PermutationGroups
+using PropertyT.StarAlgebras
+
+include(joinpath(@__DIR__, "argparse.jl"))
+include(joinpath(@__DIR__, "utils.jl"))
+
+# const N = parsed_args["N"]
+const HALFRADIUS = parsed_args["halfradius"]
+const UPPER_BOUND = parsed_args["upper_bound"]
+
+include(joinpath(@__DIR__, "./G₂_gens.jl"))
+
+G, roots, Weyl = G₂_roots_weyl()
+@info "Running Adj² - λ·Δ sum of squares decomposition for G₂"
+
+@info "computing group algebra structure"
+RG, S, sizes = @time PropertyT.group_algebra(G, halfradius = HALFRADIUS)
+
+@info "computing WedderburnDecomposition"
+wd = let Σ = Weyl, RG = RG
+    act = PropertyT.AlphabetPermutation{eltype(Σ),Int64}(
+        Dict(g => PermutationGroups.perm(g) for g in Σ),
+    )
+
+    @time SymbolicWedderburn.WedderburnDecomposition(
+        Float64,
+        Σ,
+        act,
+        basis(RG),
+        StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[HALFRADIUS]]),
+        semisimple = false,
+    )
+end
+@info wd
+
+function desubscriptify(symbol::Symbol)
+    digits = [
+        Int(l) - 0x2080 for
+        l in reverse(string(symbol)) if 0 ≤ Int(l) - 0x2080 ≤ 9
+    ]
+    res = 0
+    for (i, d) in enumerate(digits)
+        res += 10^(i - 1) * d
+    end
+    return res
+end
+
+function PropertyT.grading(g::MatrixGroups.MatrixElt, roots = roots)
+    id = desubscriptify(g.id)
+    return roots[id]
+end
+
+Δ = RG(length(S)) - sum(RG(s) for s in S)
+Δs = PropertyT.laplacians(
+    RG,
+    S,
+    x -> (gx = PropertyT.grading(x); Set([gx, -gx])),
+)
+
+elt = PropertyT.Adj(Δs)
+@assert elt == Δ^2 - PropertyT.Sq(Δs)
+unit = Δ
+
+@time model, varP = PropertyT.sos_problem_primal(
+    elt,
+    unit,
+    wd;
+    upper_bound = UPPER_BOUND,
+    augmented = true,
+    show_progress = true,
+)
+
+solve_in_loop(
+    model,
+    wd,
+    varP;
+    logdir = "./log/G2/r=$HALFRADIUS/Adj-$(UPPER_BOUND)Δ",
+    optimizer = scs_optimizer(;
+        linear_solver = SCS.MKLDirectSolver,
+        eps = 1e-9,
+        max_iters = 100_000,
+        accel = 50,
+        alpha = 1.95,
+    ),
+    data = (elt = elt, unit = unit, halfradius = HALFRADIUS),
+)
diff --git a/scripts/G₂_gens.jl b/scripts/G₂_gens.jl
new file mode 100644
index 0000000..48677d5
--- /dev/null
+++ b/scripts/G₂_gens.jl
@@ -0,0 +1,308 @@
+#= GAP code to generate matrices
+alg := SimpleLieAlgebra("G", 2, Rationals);
+root_sys := RootSystem(alg);
+pos_gens := PositiveRootVectors(root_sys);
+pos_rts := PositiveRoots(root_sys);
+
+neg_gens := NegativeRootVectors(root_sys);
+neg_rts := NegativeRoots(root_sys);
+
+alg_gens := ShallowCopy(pos_gens);;
+Append(alg_gens, neg_gens);
+grading := ShallowCopy(pos_rts);
+Append(grading, neg_rts);
+
+mats := List(alg_gens, x->AdjointMatrix(Basis(alg), x));
+
+W := WeylGroup(root_sys);
+PW := Action(W, grading, OnRight);
+=#
+
+using LinearAlgebra
+
+function matrix_exp(M::AbstractMatrix{<:Integer})
+    res = zeros(Rational{eltype(M)}, size(M))
+    res += I
+    k = 0
+    expM = one(M)
+    while !iszero(expM)
+        k += 1
+        expM *= M
+        @. res += 1 // factorial(k) * expM
+        if k == 20
+            @warn "matrix exponential did not converge" norm(expM - exp(M))
+            break
+        end
+    end
+    @debug "matrix_exp converged after $k iterations"
+    return res
+end
+
+const gap_adj_mats = [
+    [
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 1],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+    ],
+    [
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, -2],
+        [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
+    ],
+    [
+        [0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1],
+        [2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0],
+    ],
+    [
+        [0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0],
+        [3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0],
+    ],
+    [
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 1],
+        [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
+    ],
+    [
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0],
+    ],
+    [
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, -1],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+    ],
+    [
+        [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 2],
+        [0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+    ],
+    [
+        [0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1],
+        [0, 0, 0, 0, 0, 0, -2, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 0, 0, 0, 0],
+        [0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+    ],
+    [
+        [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
+        [0, 0, 0, 0, 0, 0, -3, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0],
+        [0, 0, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+    ],
+    [
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, -1],
+        [0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+    ],
+    [
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
+        [0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0],
+    ],
+]
+
+function G₂_matrices_roots()
+    adj_mats = map(gap_adj_mats) do m
+        return hcat(m...)
+    end
+    adj_mats = filter!(!isdiag, adj_mats) # remove the ones from center
+
+    gens_mats = [convert(Matrix{Int}, matrix_exp(m')) for m in adj_mats]
+
+    #=
+    The roots from
+
+    G₂roots_gap = [
+        [2, -1], # α = e₁ - e₂
+        [-3, 2], # A = -α + β = -e₁ + 2e₂ - e₃
+        [-1, 1], # β = e₂ - e₃
+        [1, 0], # α + β = e₁ - e₃
+        [3, -1], # B = 2α + β = 2e₁ - e₂ - e₃
+        [0, 1], # A + B = α + 2β = e₁ + e₂ - 2e₃
+        [-2, 1], # -α
+        [3, -2], # -A
+        [1, -1], # -β
+        [-1, 0], # -α - β
+        [-3, 1], # -B
+        [0, -1], # -A - B
+    ]
+
+    G₂roots_gap are the ones from cartan matrix. To obtain the standard
+    (hexagonal) picture map them by `T` defined as follows:
+    ```julia
+        cartan = hcat(G₂roots_gap[1:2]...)
+        rot(α) = [cos(α) -sin(α); sin(α) cos(α)]
+
+        c₁ = [√2, 0]
+        c₂ = rot(5π / 6) * [√2, 0] * √3 # (= 1/2[√6, 1])
+
+        T = hcat(c₁, c₂) * inv(cartan)
+    ```
+    By plotting one against the others (or by blind calculation) one can see
+    the following assignment. Here `⟨α, β⟩_ℤ = A₂` and `⟨A, B⟩_ℤ ≅ √3/√2 A₂`.
+    =#
+    e₁ = PropertyT.Roots.𝕖(3, 1)
+    e₂ = PropertyT.Roots.𝕖(3, 2)
+    e₃ = PropertyT.Roots.𝕖(3, 3)
+
+    α = e₁ - e₂
+    β = e₂ - e₃
+    A = -α + β
+    B = α + (α + β)
+
+    roots = [α, A, β, α + β, B, A + B, -α, -A, -β, -α - β, -B, -A - B]
+
+    return gens_mats, roots
+end
+
+function G₂_roots_weyl()
+    (mats, roots) = G₂_matrices_roots()
+    d = size(first(mats), 1)
+    G₂ = Groups.MatrixGroup{d}(mats)
+
+    m = Groups.gens(G₂)
+
+    σ = let w = m[1] * inv(m[7]) * m[1], m = union(m, inv.(m))
+        PermutationGroups.Perm([findfirst(==(inv(w) * x * w), m) for x in m])
+    end
+
+    τ = let w = m[2] * inv(m[8]) * m[2], m = union(m, inv.(m))
+        PermutationGroups.Perm([findfirst(==(inv(w) * x * w), m) for x in m])
+    end
+
+    W = PermGroup(σ, τ)
+
+    return G₂, roots, W
+end
diff --git a/scripts/G₂_has_T.jl b/scripts/G₂_has_T.jl
new file mode 100644
index 0000000..cde9048
--- /dev/null
+++ b/scripts/G₂_has_T.jl
@@ -0,0 +1,98 @@
+using LinearAlgebra
+BLAS.set_num_threads(8)
+
+ENV["OMP_NUM_THREADS"] = 4
+
+using Groups
+import Groups.MatrixGroups
+
+include(joinpath(@__DIR__, "../test/optimizers.jl"))
+using PropertyT
+
+using PropertyT.SymbolicWedderburn
+using PropertyT.PermutationGroups
+using PropertyT.StarAlgebras
+
+include(joinpath(@__DIR__, "argparse.jl"))
+include(joinpath(@__DIR__, "utils.jl"))
+
+# const N = parsed_args["N"]
+const HALFRADIUS = parsed_args["halfradius"]
+const UPPER_BOUND = parsed_args["upper_bound"]
+
+include(joinpath(@__DIR__, "./G₂_gens.jl"))
+
+G, roots, Weyl = G₂_roots_weyl()
+@info "Running Δ² - λ·Δ sum of squares decomposition for G₂"
+
+@info "computing group algebra structure"
+RG, S, sizes = @time PropertyT.group_algebra(G, halfradius = HALFRADIUS)
+
+@info "computing WedderburnDecomposition"
+wd = let Σ = Weyl, RG = RG
+    act = PropertyT.AlphabetPermutation{eltype(Σ),Int64}(
+        Dict(g => PermutationGroups.perm(g) for g in Σ),
+    )
+
+    @time SymbolicWedderburn.WedderburnDecomposition(
+        Float64,
+        Σ,
+        act,
+        basis(RG),
+        StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[HALFRADIUS]]),
+        semisimple = false,
+    )
+end
+@info wd
+
+Δ = RG(length(S)) - sum(RG(s) for s in S)
+elt = Δ^2
+unit = Δ
+
+@time model, varP = PropertyT.sos_problem_primal(
+    elt,
+    unit,
+    wd;
+    upper_bound = UPPER_BOUND,
+    augmented = true,
+    show_progress = false,
+)
+warm = nothing
+status = JuMP.OPTIMIZE_NOT_CALLED
+certified, λ = false, nothing
+
+while status ≠ JuMP.OPTIMAL
+    @time status, warm = PropertyT.solve(
+        model,
+        scs_optimizer(;
+            eps = 1e-10,
+            max_iters = 20_000,
+            accel = 50,
+            alpha = 1.95,
+        ),
+        warm,
+    )
+
+    @info "reconstructing the solution"
+    Q = @time let wd = wd, Ps = [JuMP.value.(P) for P in varP]
+        Qs = real.(sqrt.(Ps))
+        PropertyT.reconstruct(Qs, wd)
+    end
+
+    @info "certifying the solution"
+    @time certified, λ = PropertyT.certify_solution(
+        elt,
+        unit,
+        JuMP.objective_value(model),
+        Q;
+        halfradius = HALFRADIUS,
+        augmented = true,
+    )
+end
+
+if certified && λ > 0
+    Κ(λ, S) = round(sqrt(2λ / length(S)), Base.RoundDown; digits = 5)
+    @info "Certified result: G₂ has property (T):" N λ Κ(λ, S)
+else
+    @info "Could NOT certify the result:" certified λ
+end
diff --git a/scripts/SpN_Adj.jl b/scripts/SpN_Adj.jl
index c5922de..35b9a6b 100644
--- a/scripts/SpN_Adj.jl
+++ b/scripts/SpN_Adj.jl
@@ -24,8 +24,7 @@ const GENUS = 2N
 
 G = MatrixGroups.SymplecticGroup{GENUS}(Int8)
 
-RG, S, sizes =
-    @time PropertyT.group_algebra(G, halfradius=HALFRADIUS, twisted=true)
+RG, S, sizes = @time PropertyT.group_algebra(G, halfradius = HALFRADIUS)
 
 wd = let RG = RG, N = N
     G = StarAlgebras.object(RG)
@@ -42,6 +41,8 @@ wd = let RG = RG, N = N
         basis(RG),
         StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[HALFRADIUS]]),
     )
+    @info wdfl
+    wdfl
 end
 
 Δ = RG(length(S)) - sum(RG(s) for s in S)
@@ -58,23 +59,22 @@ unit = Δ
 @time model, varP = PropertyT.sos_problem_primal(
     elt,
     unit,
-    wd,
-    upper_bound=UPPER_BOUND,
-    augmented=true,
-    show_progress=true
+    wd;
+    upper_bound = UPPER_BOUND,
+    augmented = true,
+    show_progress = true,
 )
 
 solve_in_loop(
     model,
     wd,
-    varP,
-    logdir="./log/Sp($N,Z)/r=$HALFRADIUS/Adj_C₂-InfΔ",
-    optimizer=cosmo_optimizer(
-        eps=1e-10,
-        max_iters=20_000,
-        accel=50,
-        alpha=1.95,
+    varP;
+    logdir = "./log/Sp($N,Z)/r=$HALFRADIUS/Adj_C₂-$(UPPER_BOUND)Δ",
+    optimizer = cosmo_optimizer(;
+        eps = 1e-10,
+        max_iters = 20_000,
+        accel = 50,
+        alpha = 1.95,
     ),
-    data=(elt=elt, unit=unit, halfradius=HALFRADIUS)
+    data = (elt = elt, unit = unit, halfradius = HALFRADIUS),
 )
-
diff --git a/scripts/utils.jl b/scripts/utils.jl
index df891b3..4c4b394 100644
--- a/scripts/utils.jl
+++ b/scripts/utils.jl
@@ -11,12 +11,19 @@ function get_solution(model)
     return solution
 end
 
-function get_solution(model, wd, varP; logdir)
+function get_solution(model, wd, varP, eps = 1e-10)
     λ = JuMP.value(model[:λ])
 
-    Qs = [real.(sqrt(JuMP.value.(P))) for P in varP]
-    Q = PropertyT.reconstruct(Qs, wd)
+    @info "reconstructing the solution"
+    Q = @time let wd = wd, Ps = [JuMP.value.(P) for P in varP], eps = eps
+        PropertyT.__droptol!.(Ps, 100eps)
+        Qs = real.(sqrt.(Ps))
+        PropertyT.__droptol!.(Qs, eps)
+        PropertyT.reconstruct(Qs, wd)
+    end
+
     solution = Dict(:λ => λ, :Q => Q)
+
     return solution
 end
 
@@ -42,22 +49,23 @@ function solve_in_loop(model::JuMP.Model, args...; logdir, optimizer, data)
             # logstream = current_logger().logger.stream
             # v = @ccall setvbuf(logstream.handle::Ptr{Cvoid}, C_NULL::Ptr{Cvoid}, 1::Cint, 0::Cint)::Cint
             # @warn v
-            status, warm = @time PropertyT.solve(log_file, model, optimizer, warm)
-
-            solution = get_solution(model, args...; logdir=logdir)
-            solution[:warm] = warm
+            status, warm =
+                @time PropertyT.solve(log_file, model, optimizer, warm)
 
+            solution = get_solution(model, args...)
             serialize(joinpath(logdir, "solution_$date.sjl"), solution)
             serialize(joinpath(logdir, "solution.sjl"), solution)
 
-            flag, λ_cert = open(log_file, append=true) do io
+            solution[:warm] = warm
+
+            flag, λ_cert = open(log_file; append = true) do io
                 with_logger(SimpleLogger(io)) do
-                    PropertyT.certify_solution(
+                    return PropertyT.certify_solution(
                         data.elt,
                         data.unit,
                         solution[:λ],
-                        solution[:Q],
-                        halfradius=data.halfradius,
+                        solution[:Q];
+                        halfradius = data.halfradius,
                     )
                 end
             end
@@ -66,19 +74,20 @@ function solve_in_loop(model::JuMP.Model, args...; logdir, optimizer, data)
         end
 
         if flag == true && certified_λ ≥ 0
-            @info "Certification done with λ = $certified_λ"
+            @info "Certification done with λ = $certified_λ" certified_λ status
             return certified_λ
         else
-            rel_change = abs(certified_λ - old_lambda) / (abs(certified_λ) + abs(old_lambda))
-            @info "Certification failed with λ = $λ" certified_λ rel_change
+            rel_change =
+                abs(certified_λ - old_lambda) /
+                (abs(certified_λ) + abs(old_lambda))
+            @info "Certification failed with λ = $λ" certified_λ rel_change status
+            if rel_change < 1e-9
+                @info "No progress detected, breaking" certified_λ rel_change status
+                break
+            end
         end
 
         old_lambda = certified_λ
-
-        if rel_change < 1e-9
-            @info "No progress detected, breaking"
-            break
-        end
     end
 
     return status == JuMP.OPTIMAL ? old_lambda : NaN
diff --git a/src/PropertyT.jl b/src/PropertyT.jl
index 6b97d5d..60b0aef 100644
--- a/src/PropertyT.jl
+++ b/src/PropertyT.jl
@@ -3,7 +3,6 @@ module PropertyT
 using LinearAlgebra
 using SparseArrays
 
-using IntervalArithmetic
 using JuMP
 
 using Groups
@@ -26,17 +25,17 @@ include("gradings.jl")
 
 include("actions/actions.jl")
 
-function group_algebra(G::Groups.Group, S=gens(G); halfradius::Integer, twisted::Bool)
+function group_algebra(G::Groups.Group, S = gens(G); halfradius::Integer)
     S = union!(S, inv.(S))
     @info "generating wl-metric ball of radius $(2halfradius)"
-    @time E, sizes = Groups.wlmetric_ball(S, radius=2halfradius)
+    @time E, sizes = Groups.wlmetric_ball(S; radius = 2halfradius)
     @info "sizes = $(sizes)"
     @info "computing the *-algebra structure for G"
-    @time RG = StarAlgebras.StarAlgebra{twisted}(
+    @time RG = StarAlgebras.StarAlgebra(
         G,
         StarAlgebras.Basis{UInt32}(E),
-        (sizes[halfradius], sizes[halfradius]),
-        precompute=false,
+        (sizes[halfradius], sizes[halfradius]);
+        precompute = false,
     )
     return RG, S, sizes
 end
diff --git a/src/actions/actions.jl b/src/actions/actions.jl
index 9cc6872..8205e8b 100644
--- a/src/actions/actions.jl
+++ b/src/actions/actions.jl
@@ -1,5 +1,4 @@
 import SymbolicWedderburn.action
-StarAlgebras.star(g::Groups.GroupElement) = inv(g)
 
 include("alphabet_permutation.jl")
 
@@ -10,7 +9,7 @@ include("autfn_conjugation.jl")
 function SymbolicWedderburn.action(
     act::SymbolicWedderburn.ByPermutations,
     g::Groups.GroupElement,
-    α::StarAlgebras.AlgebraElement
+    α::StarAlgebras.AlgebraElement,
 )
     res = StarAlgebras.zero!(similar(α))
     B = basis(parent(α))
diff --git a/src/certify.jl b/src/certify.jl
index cdebcc9..b93c843 100644
--- a/src/certify.jl
+++ b/src/certify.jl
@@ -1,3 +1,6 @@
+import IntervalArithmetic
+import IntervalMatrices
+
 function augment_columns!(Q::AbstractMatrix)
     for c in eachcol(Q)
         c .-= sum(c) ./ length(c)
@@ -8,26 +11,25 @@ end
 function __sos_via_sqr!(
     res::StarAlgebras.AlgebraElement,
     P::AbstractMatrix;
-    augmented::Bool
+    augmented::Bool,
+    id = (b = basis(parent(res)); b[one(first(b))]),
 )
-    StarAlgebras.zero!(res)
     A = parent(res)
-    b = basis(A)
-    @assert size(A.mstructure) == size(P)
-    e = b[one(b[1])]
+    mstr = A.mstructure
+    @assert size(mstr) == size(P)
 
-    for i in axes(A.mstructure, 1)
-        x = StarAlgebras._istwisted(A.mstructure) ? StarAlgebras.star(b[i]) : b[i]
-        for j in axes(A.mstructure, 2)
+    StarAlgebras.zero!(res)
+    for j in axes(mstr, 2)
+        for i in axes(mstr, 1)
             p = P[i, j]
-            xy = b[A.mstructure[i, j]]
-            # either result += P[x,y]*(x*y)
-            res[xy] += p
+            x_star_y = mstr[-i, j]
+            res[x_star_y] += p
+            # either result += P[x,y]*(x'*y)
             if augmented
-                # or result += P[x,y]*(1-x)*(1-y) == P[x,y]*(2-x-y+xy)
-                y = b[j]
-                res[e] += p
-                res[x] -= p
+                # or result += P[x,y]*(1-x)'*(1-y) == P[x,y]*(1-x'-y+x'y)
+                res[id] += p
+                x_star, y = mstr[-i, id], j
+                res[x_star] -= p
                 res[y] -= p
             end
         end
@@ -36,7 +38,11 @@ function __sos_via_sqr!(
     return res
 end
 
-function __sos_via_cnstr!(res::StarAlgebras.AlgebraElement, Q²::AbstractMatrix, cnstrs)
+function __sos_via_cnstr!(
+    res::StarAlgebras.AlgebraElement,
+    Q²::AbstractMatrix,
+    cnstrs,
+)
     StarAlgebras.zero!(res)
     for (g, A_g) in cnstrs
         res[g] = dot(A_g, Q²)
@@ -44,10 +50,14 @@ function __sos_via_cnstr!(res::StarAlgebras.AlgebraElement, Q²::AbstractMatrix,
     return res
 end
 
-function compute_sos(A::StarAlgebras.StarAlgebra, Q::AbstractMatrix; augmented::Bool)
+function compute_sos(
+    A::StarAlgebras.StarAlgebra,
+    Q::AbstractMatrix;
+    augmented::Bool,
+)
     Q² = Q' * Q
     res = StarAlgebras.AlgebraElement(zeros(eltype(Q²), length(basis(A))), A)
-    res = __sos_via_sqr!(res, Q², augmented=augmented)
+    res = __sos_via_sqr!(res, Q²; augmented = augmented)
     return res
 end
 
@@ -56,7 +66,7 @@ function sufficient_λ(residual::StarAlgebras.AlgebraElement, λ; halfradius)
     suff_λ = λ - 2.0^(2ceil(log2(halfradius))) * L1_norm
 
     eq_sign = let T = eltype(residual)
-        if T <: Interval
+        if T <: IntervalArithmetic.Interval
             "∈"
         elseif T <: Union{Rational,Integer}
             "="
@@ -81,13 +91,12 @@ function sufficient_λ(
     order_unit::StarAlgebras.AlgebraElement,
     λ,
     sos::StarAlgebras.AlgebraElement;
-    halfradius
+    halfradius,
 )
-
     @assert parent(elt) === parent(order_unit) == parent(sos)
     residual = (elt - λ * order_unit) - sos
 
-    return sufficient_λ(residual, λ; halfradius=halfradius)
+    return sufficient_λ(residual, λ; halfradius = halfradius)
 end
 
 function certify_solution(
@@ -96,24 +105,27 @@ function certify_solution(
     λ,
     Q::AbstractMatrix{<:AbstractFloat};
     halfradius,
-    augmented=iszero(StarAlgebras.aug(elt)) && iszero(StarAlgebras.aug(orderunit))
+    augmented = iszero(StarAlgebras.aug(elt)) &&
+        iszero(StarAlgebras.aug(orderunit)),
 )
-
-    should_we_augment = !augmented && StarAlgebras.aug(elt) == StarAlgebras.aug(orderunit) == 0
+    should_we_augment =
+        !augmented && StarAlgebras.aug(elt) == StarAlgebras.aug(orderunit) == 0
 
     Q = should_we_augment ? augment_columns!(Q) : Q
-    @time sos = compute_sos(parent(elt), Q, augmented=augmented)
+    @time sos = compute_sos(parent(elt), Q; augmented = augmented)
 
     @info "Checking in $(eltype(sos)) arithmetic with" λ
 
-    λ_flpoint = sufficient_λ(elt, orderunit, λ, sos, halfradius=halfradius)
+    λ_flpoint = sufficient_λ(elt, orderunit, λ, sos; halfradius = halfradius)
 
     if λ_flpoint ≤ 0
         return false, λ_flpoint
     end
 
-    λ_int = @interval(λ)
-    Q_int = [@interval(q) for q in Q]
+    λ_int = IntervalArithmetic.@interval(λ)
+    Q_int = IntervalMatrices.IntervalMatrix([
+        IntervalArithmetic.@interval(q) for q in Q
+    ])
 
     check, sos_int = @time if should_we_augment
         @info("Projecting columns of Q to the augmentation ideal...")
@@ -123,16 +135,16 @@ function certify_solution(
         check_augmented || @error(
             "Augmentation failed! The following numbers are not certified!"
         )
-        sos_int = compute_sos(parent(elt), Q_int; augmented=augmented)
+        sos_int = compute_sos(parent(elt), Q_int; augmented = augmented)
         check_augmented, sos_int
     else
-        true, compute_sos(parent(elt), Q_int, augmented=augmented)
+        true, compute_sos(parent(elt), Q_int; augmented = augmented)
     end
 
     @info "Checking in $(eltype(sos_int)) arithmetic with" λ_int
 
     λ_certified =
-        sufficient_λ(elt, orderunit, λ_int, sos_int, halfradius=halfradius)
+        sufficient_λ(elt, orderunit, λ_int, sos_int; halfradius = halfradius)
 
-    return check && inf(λ_certified) > 0.0, inf(λ_certified)
+    return check && IntervalArithmetic.inf(λ_certified) > 0.0, λ_certified
 end
diff --git a/src/constraint_matrix.jl b/src/constraint_matrix.jl
index 8fbd36a..997a037 100644
--- a/src/constraint_matrix.jl
+++ b/src/constraint_matrix.jl
@@ -28,7 +28,12 @@ struct ConstraintMatrix{T,I} <: AbstractMatrix{T}
     size::Tuple{Int,Int}
     val::T
 
-    function ConstraintMatrix{T}(nzeros::AbstractArray{<:Integer}, n, m, val) where {T}
+    function ConstraintMatrix{T}(
+        nzeros::AbstractArray{<:Integer},
+        n,
+        m,
+        val,
+    ) where {T}
         @assert n ≥ 1
         @assert m ≥ 1
 
@@ -45,15 +50,26 @@ struct ConstraintMatrix{T,I} <: AbstractMatrix{T}
     end
 end
 
-ConstraintMatrix(nzeros::AbstractArray{<:Integer}, n, m, val::T) where {T} =
-    ConstraintMatrix{T}(nzeros, n, m, val)
+function ConstraintMatrix(
+    nzeros::AbstractArray{<:Integer},
+    n,
+    m,
+    val::T,
+) where {T}
+    return ConstraintMatrix{T}(nzeros, n, m, val)
+end
 
 Base.size(cm::ConstraintMatrix) = cm.size
 
-__get_positive(cm::ConstraintMatrix, idx::Integer) =
-    convert(eltype(cm), cm.val * length(searchsorted(cm.pos, idx)))
-__get_negative(cm::ConstraintMatrix, idx::Integer) =
-    convert(eltype(cm), cm.val * length(searchsorted(cm.neg, idx)))
+function __get_positive(cm::ConstraintMatrix, idx::Integer)
+    return convert(eltype(cm), cm.val * length(searchsorted(cm.pos, idx)))
+end
+function __get_negative(cm::ConstraintMatrix, idx::Integer)
+    return convert(
+        eltype(cm),
+        cm.val * length(searchsorted(cm.neg, idx; rev = true)),
+    )
+end
 
 Base.@propagate_inbounds function Base.getindex(
     cm::ConstraintMatrix,
@@ -67,26 +83,29 @@ Base.@propagate_inbounds function Base.getindex(
     return pos - neg
 end
 
-struct NZPairsIter{T}
-    m::ConstraintMatrix{T}
+struct NZPairsIter{T,I}
+    m::ConstraintMatrix{T,I}
 end
 
 Base.eltype(::Type{NZPairsIter{T}}) where {T} = Pair{Int,T}
 Base.IteratorSize(::Type{<:NZPairsIter}) = Base.SizeUnknown()
 
 # TODO: iterate over (idx=>val) pairs combining vals
-function Base.iterate(itr::NZPairsIter, state::Tuple{Int,Int}=(1, 1))
+function Base.iterate(
+    itr::NZPairsIter,
+    state::Tuple{Int,Nothing} = (1, nothing),
+)
     k = iterate(itr.m.pos, state[1])
-    isnothing(k) && return iterate(itr, state[2])
+    isnothing(k) && return iterate(itr, (nothing, 1))
     idx, st = k
-    return idx => itr.m.val, (st, 1)
+    return idx => itr.m.val, (st, nothing)
 end
 
-function Base.iterate(itr::NZPairsIter, state::Int)
-    k = iterate(itr.m.neg, state[1])
+function Base.iterate(itr::NZPairsIter, state::Tuple{Nothing,Int})
+    k = iterate(itr.m.neg, state[2])
     isnothing(k) && return nothing
     idx, st = k
-    return idx => -itr.m.val, st
+    return idx => -itr.m.val, (nothing, st)
 end
 
 """
@@ -127,3 +146,50 @@ function LinearAlgebra.dot(cm::ConstraintMatrix, m::AbstractMatrix{T}) where {T}
     neg = isempty(cm.neg) ? zero(first(m)) : sum(@view m[cm.neg])
     return convert(eltype(cm), cm.val) * (pos - neg)
 end
+
+function constraints(A::StarAlgebras.StarAlgebra; augmented::Bool)
+    return constraints(basis(A), A.mstructure; augmented = augmented)
+end
+
+function constraints(
+    basis::StarAlgebras.AbstractBasis,
+    mstr::StarAlgebras.MultiplicativeStructure;
+    augmented = false,
+)
+    cnstrs = _constraints(
+        mstr;
+        augmented = augmented,
+        num_constraints = length(basis),
+        id = basis[one(first(basis))],
+    )
+
+    return Dict(
+        basis[i] => ConstraintMatrix(c, size(mstr)..., 1) for
+        (i, c) in pairs(cnstrs)
+    )
+end
+
+function _constraints(
+    mstr::StarAlgebras.MultiplicativeStructure;
+    augmented::Bool = false,
+    num_constraints = maximum(mstr),
+    id,
+)
+    cnstrs = [signed(eltype(mstr))[] for _ in 1:num_constraints]
+    LI = LinearIndices(size(mstr))
+
+    for ci in CartesianIndices(size(mstr))
+        k = LI[ci]
+        i, j = Tuple(ci)
+        a_star_b = mstr[-i, j]
+        push!(cnstrs[a_star_b], k)
+        if augmented
+            # (1-a)'(1-b) = 1 - a' - b + a'b
+            push!(cnstrs[id], k)
+            a_star, b = mstr[-i, id], j
+            push!(cnstrs[a_star], -k)
+            push!(cnstrs[b], -k)
+        end
+    end
+    return cnstrs
+end
diff --git a/src/gradings.jl b/src/gradings.jl
index 8f34efa..b51de1a 100644
--- a/src/gradings.jl
+++ b/src/gradings.jl
@@ -1,7 +1,8 @@
 ## something about roots
 
-Roots.Root(e::MatrixGroups.ElementaryMatrix{N}) where {N} =
-    Roots.𝕖(N, e.i) - Roots.𝕖(N, e.j)
+function Roots.Root(e::MatrixGroups.ElementaryMatrix{N}) where {N}
+    return Roots.𝕖(N, e.i) - Roots.𝕖(N, e.j)
+end
 
 function Roots.Root(s::MatrixGroups.ElementarySymplectic{N}) where {N}
     if s.symbol === :A
@@ -51,20 +52,28 @@ function Adj(rootsystem::AbstractDict, subtype::Symbol)
         +,
         (
             Δα * Δβ for (α, Δα) in rootsystem for (β, Δβ) in rootsystem if
-            Roots.classify_sub_root_system(
-                roots,
-                first(α),
-                first(β),
-            ) == subtype
-        ),
-        init=zero(first(values(rootsystem))),
+            Roots.classify_sub_root_system(roots, first(α), first(β)) == subtype
+        );
+        init = zero(first(values(rootsystem))),
+    )
+end
+
+Adj(rootsystem::AbstractDict) = sum(values(rootsystem))^2 - Sq(rootsystem)
+
+function Sq(rootsystem::AbstractDict)
+    return reduce(
+        +,
+        Δα^2 for (_, Δα) in rootsystem;
+        init = zero(first(values(rootsystem))),
     )
 end
 
 function level(rootsystem, level::Integer)
     1 ≤ level ≤ 4 || throw("level is implemented only for i ∈{1,2,3,4}")
     level == 1 && return Adj(rootsystem, :C₁) # always positive
-    level == 2 && return Adj(rootsystem, :A₁) + Adj(rootsystem, Symbol("C₁×C₁")) + Adj(rootsystem, :C₂) # C₂ is not positive
+    level == 2 && return Adj(rootsystem, :A₁) +
+           Adj(rootsystem, Symbol("C₁×C₁")) +
+           Adj(rootsystem, :C₂) # C₂ is not positive
     level == 3 && return Adj(rootsystem, :A₂) + Adj(rootsystem, Symbol("A₁×C₁"))
     level == 4 && return Adj(rootsystem, Symbol("A₁×A₁")) # positive
 end
diff --git a/src/reconstruct.jl b/src/reconstruct.jl
index 75451ee..ce55619 100644
--- a/src/reconstruct.jl
+++ b/src/reconstruct.jl
@@ -7,33 +7,28 @@ end
 
 function reconstruct(
     Ms::AbstractVector{<:AbstractMatrix},
-    wbdec::WedderburnDecomposition;
-    atol=eps(real(eltype(wbdec))) * 10__outer_dim(wbdec)
+    wbdec::WedderburnDecomposition,
 )
     n = __outer_dim(wbdec)
     res = sum(zip(Ms, SymbolicWedderburn.direct_summands(wbdec))) do (M, ds)
         res = similar(M, n, n)
-        reconstruct!(res, M, ds, __group_of(wbdec), wbdec.hom, atol=atol)
+        res = _reconstruct!(res, M, ds)
+        return res
     end
+    res = average!(zero(res), res, __group_of(wbdec), wbdec.hom)
     return res
 end
 
-function reconstruct!(
+function _reconstruct!(
     res::AbstractMatrix,
     M::AbstractMatrix,
     ds::SymbolicWedderburn.DirectSummand,
-    G,
-    hom;
-    atol=eps(real(eltype(ds))) * 10max(size(res)...)
 )
     res .= zero(eltype(res))
-    U = SymbolicWedderburn.image_basis(ds)
-    d = SymbolicWedderburn.degree(ds)
-    tmp = (U' * M * U) .* d
-
-    res = average!(res, tmp, G, hom)
-    if eltype(res) <: AbstractFloat
-        __droptol!(res, atol) # TODO: is this really necessary?!
+    if !iszero(M)
+        U = SymbolicWedderburn.image_basis(ds)
+        d = SymbolicWedderburn.degree(ds)
+        res = (U' * M * U) .* d
     end
     return res
 end
@@ -52,18 +47,22 @@ function average!(
     res::AbstractMatrix,
     M::AbstractMatrix,
     G::Groups.Group,
-    hom::SymbolicWedderburn.InducedActionHomomorphism{<:SymbolicWedderburn.ByPermutations}
+    hom::SymbolicWedderburn.InducedActionHomomorphism{
+        <:SymbolicWedderburn.ByPermutations,
+    },
 )
+    res .= zero(eltype(res))
     @assert size(M) == size(res)
+    o = Groups.order(Int, G)
     for g in G
-        gext = SymbolicWedderburn.induce(hom, g)
+        p = SymbolicWedderburn.induce(hom, g)
         Threads.@threads for c in axes(res, 2)
             for r in axes(res, 1)
-                res[r, c] += M[r^gext, c^gext]
+                if !iszero(M[r, c])
+                    res[r^p, c^p] += M[r, c] / o
+                end
             end
         end
     end
-    o = Groups.order(Int, G)
-    res ./= o
     return res
 end
diff --git a/src/roots.jl b/src/roots.jl
index 3c5d6f3..92a915c 100644
--- a/src/roots.jl
+++ b/src/roots.jl
@@ -5,111 +5,104 @@ using LinearAlgebra
 
 export Root, isproportional, isorthogonal, ~, ⟂
 
-abstract type AbstractRoot{N,T} end
+abstract type AbstractRoot{N,T} end # <: AbstractVector{T} ?
 
-struct Root{N,T} <: AbstractRoot{N,T}
-    coord::SVector{N,T}
-end
-
-Root(a) = Root(SVector(a...))
-
-function Base.:(==)(r::Root{N}, s::Root{M}) where {M,N}
-    M == N || return false
-    r.coord == s.coord || return false
-    return true
-end
-
-Base.hash(r::Root, h::UInt) = hash(r.coord, hash(Root, h))
-
-Base.:+(r::Root{N,T}, s::Root{N,T}) where {N,T} = Root{N,T}(r.coord + s.coord)
-Base.:-(r::Root{N,T}, s::Root{N,T}) where {N,T} = Root{N,T}(r.coord - s.coord)
-Base.:-(r::Root{N}) where {N} = Root(-r.coord)
-
-Base.:*(a::Number, r::Root) = Root(a * r.coord)
-Base.:*(r::Root, a::Number) = a * r
-
-Base.length(r::AbstractRoot) = norm(r, 2)
-
-LinearAlgebra.norm(r::Root, p::Real=2) = norm(r.coord, p)
-LinearAlgebra.dot(r::Root, s::Root) = dot(r.coord, s.coord)
+ℓ₂length(r::AbstractRoot) = norm(r, 2)
+ambient_dim(r::AbstractRoot) = length(r)
+Base.:*(r::AbstractRoot, a::Number) = a * r
 
 cos_angle(a, b) = dot(a, b) / (norm(a) * norm(b))
 
-function isproportional(α::AbstractRoot{N}, β::AbstractRoot{M}) where {N,M}
-    N == M || return false
+function isproportional(α::AbstractRoot, β::AbstractRoot)
+    ambient_dim(α) == ambient_dim(β) || return false
     val = abs(cos_angle(α, β))
-    return isapprox(val, one(val), atol=eps(one(val)))
+    return isapprox(val, one(val); atol = eps(one(val)))
 end
 
-function isorthogonal(α::AbstractRoot{N}, β::AbstractRoot{M}) where {N,M}
-    N == M || return false
+function isorthogonal(α::AbstractRoot, β::AbstractRoot)
+    ambient_dim(α) == ambient_dim(β) || return false
     val = cos_angle(α, β)
-    return isapprox(val, zero(val), atol=eps(one(val)))
+    return isapprox(val, zero(val); atol = eps(one(val)))
 end
 
-function _positive_direction(α::Root{N}) where {N}
-    v = α.coord + 1 / (N * 100) * rand(N)
-    return Root{N,Float64}(v / norm(v, 2))
-end
-
-function positive(roots::AbstractVector{<:Root{N}}) where {N}
-    # return those roots for which dot(α, Root([½, ¼, …])) > 0.0
+function positive(roots::AbstractVector{<:AbstractRoot})
+    isempty(roots) && return empty(roots)
     pd = _positive_direction(first(roots))
     return filter(α -> dot(α, pd) > 0.0, roots)
 end
 
-function Base.show(io::IO, r::Root)
-    print(io, "Root$(r.coord)")
+function Base.show(io::IO, r::AbstractRoot)
+    return print(io, "Root $(r.coord)")
 end
 
-function Base.show(io::IO, ::MIME"text/plain", r::Root{N}) where {N}
-    lngth² = sum(x -> x^2, r.coord)
-    l = isinteger(sqrt(lngth²)) ? "$(sqrt(lngth²))" : "√$(lngth²)"
-    print(io, "Root in ℝ^$N of length $l\n", r.coord)
+function Base.show(io::IO, ::MIME"text/plain", r::AbstractRoot)
+    l₂l = ℓ₂length(r)
+    l = round(Int, l₂l) ≈ l₂l ? "$(round(Int, l₂l))" : "√$(round(Int, l₂l^2))"
+    return print(io, "Root in ℝ^$(length(r)) of length $l\n", r.coord)
 end
 
-𝕖(N, i) = Root(ntuple(k -> k == i ? 1 : 0, N))
-𝕆(N, ::Type{T}) where {T} = Root(ntuple(_ -> zero(T), N))
+function reflection(α::AbstractRoot, β::AbstractRoot)
+    return β - Int(2dot(α, β) // dot(α, α)) * α
+end
+function cartan(α::AbstractRoot, β::AbstractRoot)
+    ambient_dim(α) == ambient_dim(β) || throw("incompatible ambient dimensions")
+    return [
+        ℓ₂length(reflection(a, b) - b) / ℓ₂length(a) for a in (α, β),
+        b in (α, β)
+    ]
+end
 
 """
     classify_root_system(α, β)
 Return the symbol of smallest system generated by roots `α` and `β`.
 
-The classification is based only on roots length and
-proportionality/orthogonality.
+The classification is based only on roots length,
+proportionality/orthogonality and Cartan matrix.
 """
-function classify_root_system(α::AbstractRoot, β::AbstractRoot)
-    lα, lβ = length(α), length(β)
+function classify_root_system(
+    α::AbstractRoot,
+    β::AbstractRoot,
+    long::Tuple{Bool,Bool},
+)
     if isproportional(α, β)
-        if lα ≈ lβ ≈ √2
-            return :A₁
-        elseif lα ≈ lβ ≈ 2.0
+        if all(long)
             return :C₁
+        elseif all(.!long) # both short
+            return :A₁
         else
+            @error "Proportional roots of different length"
             error("Unknown root system ⟨α, β⟩:\n α = $α\n β = $β")
         end
     elseif isorthogonal(α, β)
-        if lα ≈ lβ ≈ √2
-            return Symbol("A₁×A₁")
-        elseif lα ≈ lβ ≈ 2.0
+        if all(long)
             return Symbol("C₁×C₁")
-        elseif (lα ≈ 2.0 && lβ ≈ √2) || (lα ≈ √2 && lβ ≈ 2)
+        elseif all(.!long) # both short
+            return Symbol("A₁×A₁")
+        elseif any(long)
             return Symbol("A₁×C₁")
-        else
-            error("Unknown root system ⟨α, β⟩:\n α = $α\n β = $β")
         end
     else # ⟨α, β⟩ is 2-dimensional, but they're not orthogonal
-        if lα ≈ lβ ≈ √2
+        a, b, c, d = abs.(cartan(α, β))
+        @assert a == d == 2
+        b, c = b < c ? (b, c) : (c, b)
+        if b == c == 1
             return :A₂
-        elseif (lα ≈ 2.0 && lβ ≈ √2) || (lα ≈ √2 && lβ ≈ 2)
+        elseif b == 1 && c == 2
             return :C₂
+        elseif b == 1 && c == 3
+            @warn ":G₂? really?"
+            return :G₂
         else
+            @error a, b, c, d
             error("Unknown root system ⟨α, β⟩:\n α = $α\n β = $β")
         end
     end
 end
 
-function proportional_root_from_system(Ω::AbstractVector{<:Root}, α::Root)
+function proportional_root_from_system(
+    Ω::AbstractVector{<:AbstractRoot},
+    α::AbstractRoot,
+)
     k = findfirst(v -> isproportional(α, v), Ω)
     if isnothing(k)
         error("Line L_α not contained in root system Ω:\n α = $α\n Ω = $Ω")
@@ -117,25 +110,31 @@ function proportional_root_from_system(Ω::AbstractVector{<:Root}, α::Root)
     return Ω[k]
 end
 
-struct Plane{R<:Root}
+struct Plane{R<:AbstractRoot}
     v1::R
     v2::R
     vectors::Vector{R}
 end
 
-Plane(α::R, β::R) where {R<:Root} =
-    Plane(α, β, [a * α + b * β for a in -3:3 for b in -3:3])
+function Plane(α::AbstractRoot, β::AbstractRoot)
+    return Plane(α, β, [a * α + b * β for a in -3:3 for b in -3:3])
+end
 
-function Base.in(r::R, plane::Plane{R}) where {R}
+function Base.in(r::AbstractRoot, plane::Plane)
     return any(isproportional(r, v) for v in plane.vectors)
 end
 
-function classify_sub_root_system(
-    Ω::AbstractVector{<:Root{N}},
-    α::Root{N},
-    β::Root{N},
-) where {N}
+function _islong(α::AbstractRoot, Ω)
+    lα = ℓ₂length(α)
+    return any(r -> lα - ℓ₂length(r) > eps(lα), Ω)
+end
 
+function classify_sub_root_system(
+    Ω::AbstractVector{<:AbstractRoot{N}},
+    α::AbstractRoot{N},
+    β::AbstractRoot{N},
+) where {N}
+    @assert 1 ≤ length(unique(ℓ₂length, Ω)) ≤ 2
     v = proportional_root_from_system(Ω, α)
     w = proportional_root_from_system(Ω, β)
 
@@ -146,32 +145,75 @@ function classify_sub_root_system(
     l = length(subsystem)
     if l == 1
         x = first(subsystem)
-        return classify_root_system(x, x)
+        long = _islong(x, Ω)
+        return classify_root_system(x, -x, (long, long))
     elseif l == 2
-        return classify_root_system(subsystem...)
+        x, y = subsystem
+        return classify_root_system(x, y, (_islong(x, Ω), _islong(y, Ω)))
     elseif l == 3
-        a = classify_root_system(subsystem[1], subsystem[2])
-        b = classify_root_system(subsystem[2], subsystem[3])
-        c = classify_root_system(subsystem[1], subsystem[3])
+        x, y, z = subsystem
+        l1, l2, l3 = _islong(x, Ω), _islong(y, Ω), _islong(z, Ω)
+        a = classify_root_system(x, y, (l1, l2))
+        b = classify_root_system(y, z, (l2, l3))
+        c = classify_root_system(x, z, (l1, l3))
 
-        if a == b == c # it's only A₂
+        if :A₂ == a == b == c # it's only A₂
             return a
         end
 
-        C = (:C₂, Symbol("C₁×C₁"))
-        if (a ∈ C && b ∈ C && c ∈ C) && (:C₂ ∈ (a, b, c))
-            return :C₂
-        end
+        throw("Unknown subroot system! $((x,y,z))")
     elseif l == 4
-        for i = 1:l
-            for j = (i+1):l
-                T = classify_root_system(subsystem[i], subsystem[j])
-                T == :C₂ && return :C₂
-            end
+        subtypes = [
+            classify_root_system(x, y, (_islong(x, Ω), _islong(y, Ω))) for
+            x in subsystem for y in subsystem if x ≠ y
+        ]
+        if :C₂ in subtypes
+            return :C₂
         end
     end
     @error "Unknown root subsystem generated by" α β
     throw("Unknown root system: $subsystem")
 end
 
+## concrete implementation:
+struct Root{N,T} <: AbstractRoot{N,T}
+    coord::SVector{N,T}
+end
+
+Root(a) = Root(SVector(a...))
+
+# convienience constructors
+𝕖(N, i) = Root(ntuple(k -> k == i ? 1 : 0, N))
+𝕆(N, ::Type{T}) where {T} = Root(ntuple(_ -> zero(T), N))
+
+function Base.:(==)(r::Root{N}, s::Root{M}) where {M,N}
+    M == N || return false
+    r.coord == s.coord || return false
+    return true
+end
+
+Base.hash(r::Root, h::UInt) = hash(r.coord, hash(Root, h))
+
+function Base.:+(r::Root, s::Root)
+    ambient_dim(r) == ambient_dim(s) || throw("incompatible ambient dimensions")
+    return Root(r.coord + s.coord)
+end
+
+function Base.:-(r::Root, s::Root)
+    ambient_dim(r) == ambient_dim(s) || throw("incompatible ambient dimensions")
+    return Root(r.coord - s.coord)
+end
+Base.:-(r::Root) = Root(-r.coord)
+
+Base.:*(a::Number, r::Root) = Root(a * r.coord)
+
+Base.length(r::Root) = length(r.coord)
+
+LinearAlgebra.norm(r::Root, p::Real = 2) = norm(r.coord, p)
+LinearAlgebra.dot(r::Root, s::Root) = dot(r.coord, s.coord)
+
+function _positive_direction(α::Root{N}) where {N}
+    v = α.coord + 1 / (N * 100) * rand(N)
+    return Root{N,Float64}(v / norm(v, 2))
+end
 end # of module Roots
diff --git a/src/solve.jl b/src/solve.jl
index e7df555..c872617 100644
--- a/src/solve.jl
+++ b/src/solve.jl
@@ -7,12 +7,15 @@ function setwarmstart!(model::JuMP.Model, warmstart)
         ct => JuMP.all_constraints(model, ct...) for
         ct in JuMP.list_of_constraint_types(model)
     )
+    try
+        JuMP.set_start_value.(JuMP.all_variables(model), warmstart.primal)
 
-    JuMP.set_start_value.(JuMP.all_variables(model), warmstart.primal)
-
-    for (ct, idx) in pairs(constraint_map)
-        JuMP.set_start_value.(idx, warmstart.slack[ct])
-        JuMP.set_dual_start_value.(idx, warmstart.dual[ct])
+        for (ct, idx) in pairs(constraint_map)
+            JuMP.set_start_value.(idx, warmstart.slack[ct])
+            JuMP.set_dual_start_value.(idx, warmstart.dual[ct])
+        end
+    catch e
+        @warn "Warmstarting was not succesfull!" e
     end
     return model
 end
@@ -28,11 +31,10 @@ function getwarmstart(model::JuMP.Model)
     slack = Dict(k => value.(v) for (k, v) in constraint_map)
     duals = Dict(k => JuMP.dual.(v) for (k, v) in constraint_map)
 
-    return (primal=primal, dual=duals, slack=slack)
+    return (primal = primal, dual = duals, slack = slack)
 end
 
-function solve(m::JuMP.Model, optimizer, warmstart=nothing)
-
+function solve(m::JuMP.Model, optimizer, warmstart = nothing)
     JuMP.set_optimizer(m, optimizer)
     MOIU.attach_optimizer(m)
 
@@ -45,8 +47,7 @@ function solve(m::JuMP.Model, optimizer, warmstart=nothing)
     return status, getwarmstart(m)
 end
 
-function solve(solverlog::String, m::JuMP.Model, optimizer, warmstart=nothing)
-
+function solve(solverlog::String, m::JuMP.Model, optimizer, warmstart = nothing)
     isdir(dirname(solverlog)) || mkpath(dirname(solverlog))
 
     Base.flush(Base.stdout)
@@ -55,7 +56,7 @@ function solve(solverlog::String, m::JuMP.Model, optimizer, warmstart=nothing)
         redirect_stdout(logfile) do
             status, warmstart = solve(m, optimizer, warmstart)
             Base.Libc.flush_cstdio()
-            status, warmstart
+            return status, warmstart
         end
     end
 
diff --git a/src/sos_sdps.jl b/src/sos_sdps.jl
index 2f6faea..7dfb306 100644
--- a/src/sos_sdps.jl
+++ b/src/sos_sdps.jl
@@ -6,15 +6,13 @@ See also [sos_problem_primal](@ref).
 """
 function sos_problem_dual(
     elt::StarAlgebras.AlgebraElement,
-    order_unit::StarAlgebras.AlgebraElement=zero(elt);
-    lower_bound=-Inf
+    order_unit::StarAlgebras.AlgebraElement = zero(elt);
+    lower_bound = -Inf,
 )
     @assert parent(elt) == parent(order_unit)
     algebra = parent(elt)
-    mstructure = if StarAlgebras._istwisted(algebra.mstructure)
-        algebra.mstructure
-    else
-        StarAlgebras.MTable{true}(basis(algebra), table_size=size(algebra.mstructure))
+    moment_matrix = let m = algebra.mstructure
+        [m[-i, j] for i in axes(m, 1), j in axes(m, 2)]
     end
 
     # 1 variable for every primal constraint
@@ -22,60 +20,21 @@ function sos_problem_dual(
     # Symmetrized:
     # 1 dual variable for every orbit of G acting on basis
     model = Model()
-    @variable model y[1:length(basis(algebra))]
-    @constraint model λ_dual dot(order_unit, y) == 1
-    @constraint(model, psd, y[mstructure] in PSDCone())
+    JuMP.@variable(model, y[1:length(basis(algebra))])
+    JuMP.@constraint(model, λ_dual, dot(order_unit, y) == 1)
+    JuMP.@constraint(model, psd, y[moment_matrix] in PSDCone())
 
     if !isinf(lower_bound)
         throw("Not Implemented yet")
-        @variable model λ_ub_dual
-        @objective model Min dot(elt, y) + lower_bound * λ_ub_dual
+        JuMP.@variable(model, λ_ub_dual)
+        JuMP.@objective(model, Min, dot(elt, y) + lower_bound * λ_ub_dual)
     else
-        @objective model Min dot(elt, y)
+        JuMP.@objective(model, Min, dot(elt, y))
     end
 
     return model
 end
 
-function constraints(
-    basis::StarAlgebras.AbstractBasis,
-    mstr::AbstractMatrix{<:Integer};
-    augmented::Bool=false,
-    table_size=size(mstr)
-)
-    cnstrs = [signed(eltype(mstr))[] for _ in basis]
-    LI = LinearIndices(table_size)
-
-    for ci in CartesianIndices(table_size)
-        k = LI[ci]
-        a_star_b = basis[mstr[k]]
-        push!(cnstrs[basis[a_star_b]], k)
-        if augmented
-            # (1-a_star)(1-b) = 1 - a_star - b + a_star_b
-
-            i, j = Tuple(ci)
-            a, b = basis[i], basis[j]
-
-            push!(cnstrs[basis[one(a)]], k)
-            push!(cnstrs[basis[StarAlgebras.star(a)]], -k)
-            push!(cnstrs[basis[b]], -k)
-        end
-    end
-
-    return Dict(
-        basis[i] => ConstraintMatrix(c, table_size..., 1) for (i, c) in pairs(cnstrs)
-    )
-end
-
-function constraints(A::StarAlgebras.StarAlgebra; augmented::Bool, twisted::Bool)
-    mstructure = if StarAlgebras._istwisted(A.mstructure) == twisted
-        A.mstructure
-    else
-        StarAlgebras.MTable{twisted}(basis(A), table_size=size(A.mstructure))
-    end
-    return constraints(basis(A), mstructure, augmented=augmented)
-end
-
 """
     sos_problem_primal(X, [u = zero(X); upper_bound=Inf])
 Formulate sum of squares decomposition problem for `X - λ·u`.
@@ -96,9 +55,10 @@ The default `u = zero(X)` formulates a simple feasibility problem.
 """
 function sos_problem_primal(
     elt::StarAlgebras.AlgebraElement,
-    order_unit::StarAlgebras.AlgebraElement=zero(elt);
-    upper_bound=Inf,
-    augmented::Bool=iszero(StarAlgebras.aug(elt)) && iszero(StarAlgebras.aug(order_unit))
+    order_unit::StarAlgebras.AlgebraElement = zero(elt);
+    upper_bound = Inf,
+    augmented::Bool = iszero(StarAlgebras.aug(elt)) &&
+                      iszero(StarAlgebras.aug(order_unit)),
 )
     @assert parent(elt) === parent(order_unit)
 
@@ -111,7 +71,7 @@ function sos_problem_primal(
         @warn "Setting `upper_bound` together with zero `order_unit` has no effect"
     end
 
-    A = constraints(parent(elt), augmented=augmented, twisted=true)
+    A = constraints(parent(elt); augmented = augmented)
 
     if !iszero(order_unit)
         λ = JuMP.@variable(model, λ)
@@ -135,11 +95,11 @@ end
 function invariant_constraint!(
     result::AbstractMatrix,
     basis::StarAlgebras.AbstractBasis,
-    cnstrs::AbstractDict{K,CM},
+    cnstrs::AbstractDict{K,<:ConstraintMatrix},
     invariant_vec::SparseVector,
-) where {K,CM<:ConstraintMatrix}
+) where {K}
     result .= zero(eltype(result))
-    for i in SparseArrays.nonzeroinds(invariant_vec)
+    @inbounds for i in SparseArrays.nonzeroinds(invariant_vec)
         g = basis[i]
         A = cnstrs[g]
         for (idx, v) in nzpairs(A)
@@ -149,25 +109,47 @@ function invariant_constraint!(
     return result
 end
 
+function invariant_constraint(basis, cnstrs, invariant_vec)
+    I = UInt32[]
+    J = UInt32[]
+    V = Float64[]
+    _M = first(values(cnstrs))
+    CI = CartesianIndices(_M)
+    @inbounds for i in SparseArrays.nonzeroinds(invariant_vec)
+        g = basis[i]
+        A = cnstrs[g]
+        for (idx, v) in nzpairs(A)
+            ci = CI[idx]
+            push!(I, ci[1])
+            push!(J, ci[2])
+            push!(V, invariant_vec[i] * v)
+        end
+    end
+    return sparse(I, J, V, size(_M)...)
+end
+
 function isorth_projection(ds::SymbolicWedderburn.DirectSummand)
     U = SymbolicWedderburn.image_basis(ds)
     return isapprox(U * U', I)
 end
 
-sos_problem_primal(
+function sos_problem_primal(
     elt::StarAlgebras.AlgebraElement,
     wedderburn::WedderburnDecomposition;
-    kwargs...
-) = sos_problem_primal(elt, zero(elt), wedderburn; kwargs...)
+    kwargs...,
+)
+    return sos_problem_primal(elt, zero(elt), wedderburn; kwargs...)
+end
 
 function __fast_recursive_dot!(
     res::JuMP.AffExpr,
     Ps::AbstractArray{<:AbstractMatrix{<:JuMP.VariableRef}},
-    Ms::AbstractArray{<:AbstractSparseMatrix};
+    Ms::AbstractArray{<:AbstractSparseMatrix},
+    weights,
 )
     @assert length(Ps) == length(Ms)
 
-    for (A, P) in zip(Ms, Ps)
+    for (w, A, P) in zip(weights, Ms, Ps)
         iszero(Ms) && continue
         rows = rowvals(A)
         vals = nonzeros(A)
@@ -175,13 +157,21 @@ function __fast_recursive_dot!(
             for i in nzrange(A, cidx)
                 ridx = rows[i]
                 val = vals[i]
-                JuMP.add_to_expression!(res, P[ridx, cidx], val)
+                JuMP.add_to_expression!(res, P[ridx, cidx], w * val)
             end
         end
     end
     return res
 end
 
+function _dot(
+    Ps::AbstractArray{<:AbstractMatrix{<:JuMP.VariableRef}},
+    Ms::AbstractArray{<:AbstractMatrix{T}},
+    weights = Iterators.repeated(one(T), length(Ms)),
+) where {T}
+    return __fast_recursive_dot!(JuMP.AffExpr(), Ps, Ms, weights)
+end
+
 import ProgressMeter
 __show_itrs(n, total) = () -> [(Symbol("constraint"), "$n/$total")]
 
@@ -189,21 +179,38 @@ function sos_problem_primal(
     elt::StarAlgebras.AlgebraElement,
     orderunit::StarAlgebras.AlgebraElement,
     wedderburn::WedderburnDecomposition;
-    upper_bound=Inf,
-    augmented=iszero(StarAlgebras.aug(elt)) && iszero(StarAlgebras.aug(orderunit)),
-    check_orthogonality=true,
-    show_progress=false
+    upper_bound = Inf,
+    augmented = iszero(StarAlgebras.aug(elt)) &&
+        iszero(StarAlgebras.aug(orderunit)),
+    check_orthogonality = true,
+    show_progress = false,
 )
-
     @assert parent(elt) === parent(orderunit)
     if check_orthogonality
         if any(!isorth_projection, direct_summands(wedderburn))
-            error("Wedderburn decomposition contains a non-orthogonal projection")
+            error(
+                "Wedderburn decomposition contains a non-orthogonal projection",
+            )
         end
     end
 
+    id_one = findfirst(invariant_vectors(wedderburn)) do v
+        b = basis(parent(elt))
+        return sparsevec([b[one(first(b))]], [1 // 1], length(v)) == v
+    end
+
+    prog = ProgressMeter.Progress(
+        length(invariant_vectors(wedderburn));
+        dt = 1,
+        desc = "Adding constraints: ",
+        enabled = show_progress,
+        barlen = 60,
+        showspeed = true,
+    )
+
     feasibility_problem = iszero(orderunit)
 
+    # problem creation starts here
     model = JuMP.Model()
     if !feasibility_problem # add λ or not?
         λ = JuMP.@variable(model, λ)
@@ -217,58 +224,52 @@ function sos_problem_primal(
         end
     end
 
-    P = map(direct_summands(wedderburn)) do ds
+    # semidefinite constraints as described by wedderburn
+    Ps = map(direct_summands(wedderburn)) do ds
         dim = size(ds, 1)
         P = JuMP.@variable(model, [1:dim, 1:dim], Symmetric)
         JuMP.@constraint(model, P in PSDCone())
-        P
+        return P
     end
 
-    begin # preallocating
+    begin # Ms are preallocated for the constraints loop
         T = eltype(wedderburn)
-        Ms = [spzeros.(T, size(p)...) for p in P]
-        M_orb = zeros(T, size(parent(elt).mstructure)...)
+        Ms = [spzeros.(T, size(p)...) for p in Ps]
+        _eps = 10 * eps(T) * max(size(parent(elt).mstructure)...)
     end
 
-    X = convert(Vector{T}, StarAlgebras.coeffs(elt))
-    U = convert(Vector{T}, StarAlgebras.coeffs(orderunit))
+    X = StarAlgebras.coeffs(elt)
+    U = StarAlgebras.coeffs(orderunit)
 
     # defining constraints based on the multiplicative structure
-    cnstrs = constraints(parent(elt), augmented=augmented, twisted=true)
-
-    prog = ProgressMeter.Progress(
-        length(invariant_vectors(wedderburn)),
-        dt=1,
-        desc="Adding constraints... ",
-        enabled=show_progress,
-        barlen=60,
-        showspeed=true
-    )
+    cnstrs = constraints(parent(elt); augmented = augmented)
 
+    # adding linear constraints: one per orbit
     for (i, iv) in enumerate(invariant_vectors(wedderburn))
-        ProgressMeter.next!(prog, showvalues=__show_itrs(i, prog.n))
+        ProgressMeter.next!(prog; showvalues = __show_itrs(i, prog.n))
+        augmented && i == id_one && continue
+        # i == 500 && break
 
         x = dot(X, iv)
         u = dot(U, iv)
 
-        M_orb = invariant_constraint!(M_orb, basis(parent(elt)), cnstrs, iv)
-        Ms = SymbolicWedderburn.diagonalize!(Ms, M_orb, wedderburn)
-        SparseArrays.droptol!.(Ms, 10 * eps(T) * max(size(M_orb)...))
-
-        # @info [nnz(m) / length(m) for m in Ms]
+        spM_orb = invariant_constraint(basis(parent(elt)), cnstrs, iv)
 
+        Ms = SymbolicWedderburn.diagonalize!(
+            Ms,
+            spM_orb,
+            wedderburn;
+            trace_preserving = true,
+        )
+        for M in Ms
+            SparseArrays.droptol!(M, _eps)
+        end
         if feasibility_problem
-            JuMP.@constraint(
-                model,
-                x == __fast_recursive_dot!(JuMP.AffExpr(), P, Ms)
-            )
+            JuMP.@constraint(model, x == _dot(Ps, Ms))
         else
-            JuMP.@constraint(
-                model,
-                x - λ * u == __fast_recursive_dot!(JuMP.AffExpr(), P, Ms)
-            )
+            JuMP.@constraint(model, x - λ * u == _dot(Ps, Ms))
         end
     end
     ProgressMeter.finish!(prog)
-    return model, P
+    return model, Ps
 end
diff --git a/test/1703.09680.jl b/test/1703.09680.jl
index 6ed0f79..e801ab9 100644
--- a/test/1703.09680.jl
+++ b/test/1703.09680.jl
@@ -1,9 +1,8 @@
 @testset "1703.09680 Examples" begin
-
     @testset "SL(2,Z)" begin
         N = 2
         G = MatrixGroups.SpecialLinearGroup{N}(Int8)
-        RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
+        RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2)
 
         Δ = let RG = RG, S = S
             RG(length(S)) - sum(RG(s) for s in S)
@@ -15,15 +14,15 @@
 
         status, certified, λ = check_positivity(
             elt,
-            unit,
-            upper_bound=ub,
-            halfradius=2,
-            optimizer=scs_optimizer(
-                eps=1e-10,
-                max_iters=5_000,
-                accel=50,
-                alpha=1.9,
-            )
+            unit;
+            upper_bound = ub,
+            halfradius = 2,
+            optimizer = scs_optimizer(;
+                eps = 1e-10,
+                max_iters = 5_000,
+                accel = 50,
+                alpha = 1.9,
+            ),
         )
 
         @test status == JuMP.ALMOST_OPTIMAL
@@ -33,8 +32,10 @@
 
     @testset "SL(3,F₅)" begin
         N = 3
-        G = MatrixGroups.SpecialLinearGroup{N}(SymbolicWedderburn.Characters.FiniteFields.GF{5})
-        RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
+        G = MatrixGroups.SpecialLinearGroup{N}(
+            SymbolicWedderburn.Characters.FiniteFields.GF{5},
+        )
+        RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2)
 
         Δ = let RG = RG, S = S
             RG(length(S)) - sum(RG(s) for s in S)
@@ -46,15 +47,15 @@
 
         status, certified, λ = check_positivity(
             elt,
-            unit,
-            upper_bound=ub,
-            halfradius=2,
-            optimizer=scs_optimizer(
-                eps=1e-10,
-                max_iters=5_000,
-                accel=50,
-                alpha=1.9,
-            )
+            unit;
+            upper_bound = ub,
+            halfradius = 2,
+            optimizer = scs_optimizer(;
+                eps = 1e-10,
+                max_iters = 5_000,
+                accel = 50,
+                alpha = 1.9,
+            ),
         )
 
         @test status == JuMP.OPTIMAL
@@ -62,12 +63,15 @@
         @test λ > 1
 
         m = PropertyT.sos_problem_dual(elt, unit)
-        PropertyT.solve(m, cosmo_optimizer(
-            eps=1e-6,
-            max_iters=5_000,
-            accel=50,
-            alpha=1.9,
-        ))
+        PropertyT.solve(
+            m,
+            cosmo_optimizer(;
+                eps = 1e-6,
+                max_iters = 5_000,
+                accel = 50,
+                alpha = 1.9,
+            ),
+        )
 
         @test JuMP.termination_status(m) in (JuMP.ALMOST_OPTIMAL, JuMP.OPTIMAL)
         @test JuMP.objective_value(m) ≈ 1.5 atol = 1e-2
@@ -76,7 +80,7 @@
     @testset "SAut(F₂)" begin
         N = 2
         G = SpecialAutomorphismGroup(FreeGroup(N))
-        RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
+        RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2)
 
         Δ = let RG = RG, S = S
             RG(length(S)) - sum(RG(s) for s in S)
@@ -88,53 +92,46 @@
 
         status, certified, λ = check_positivity(
             elt,
-            unit,
-            upper_bound=ub,
-            halfradius=2,
-            optimizer=scs_optimizer(
-                eps=1e-10,
-                max_iters=5_000,
-                accel=50,
-                alpha=1.9,
-            )
+            unit;
+            upper_bound = ub,
+            halfradius = 2,
+            optimizer = scs_optimizer(;
+                eps = 1e-10,
+                max_iters = 5_000,
+                accel = 50,
+                alpha = 1.9,
+            ),
         )
 
         @test status == JuMP.ALMOST_OPTIMAL
         @test λ < 0
         @test !certified
 
-        @time sos_problem =
-            PropertyT.sos_problem_primal(elt, upper_bound=ub)
+        @time sos_problem = PropertyT.sos_problem_primal(elt; upper_bound = ub)
 
         status, _ = PropertyT.solve(
             sos_problem,
-            cosmo_optimizer(
-                eps=1e-7,
-                max_iters=10_000,
-                accel=0,
-                alpha=1.9,
-            )
+            cosmo_optimizer(;
+                eps = 1e-7,
+                max_iters = 10_000,
+                accel = 0,
+                alpha = 1.9,
+            ),
         )
         @test status == JuMP.OPTIMAL
         P = JuMP.value.(sos_problem[:P])
         Q = real.(sqrt(P))
-        certified, λ_cert = PropertyT.certify_solution(
-            elt,
-            zero(elt),
-            0.0,
-            Q,
-            halfradius=2,
-        )
+        certified, λ_cert =
+            PropertyT.certify_solution(elt, zero(elt), 0.0, Q; halfradius = 2)
         @test !certified
         @test λ_cert < 0
-
     end
 
     @testset "SL(3,Z) has (T)" begin
         n = 3
 
         SL = MatrixGroups.SpecialLinearGroup{n}(Int8)
-        RSL, S, sizes = PropertyT.group_algebra(SL, halfradius=2, twisted=true)
+        RSL, S, sizes = PropertyT.group_algebra(SL; halfradius = 2)
 
         Δ = RSL(length(S)) - sum(RSL(s) for s in S)
 
@@ -145,18 +142,18 @@
 
             opt_problem = PropertyT.sos_problem_primal(
                 elt,
-                unit,
-                upper_bound=ub,
-                augmented=false,
+                unit;
+                upper_bound = ub,
+                augmented = false,
             )
 
             status, _ = PropertyT.solve(
                 opt_problem,
-                cosmo_optimizer(
-                    eps=1e-10,
-                    max_iters=10_000,
-                    accel=0,
-                    alpha=1.5,
+                cosmo_optimizer(;
+                    eps = 1e-10,
+                    max_iters = 10_000,
+                    accel = 0,
+                    alpha = 1.5,
                 ),
             )
 
@@ -170,13 +167,13 @@
                 elt,
                 unit,
                 λ,
-                Q,
-                halfradius=2,
-                augmented=false,
+                Q;
+                halfradius = 2,
+                augmented = false,
             )
 
             @test certified
-            @test isapprox(λ_cert, λ, rtol=1e-5)
+            @test isapprox(λ_cert, λ, rtol = 1e-5)
         end
 
         @testset "augmented formulation" begin
@@ -186,18 +183,18 @@
 
             opt_problem = PropertyT.sos_problem_primal(
                 elt,
-                unit,
-                upper_bound=ub,
-                augmented=true,
+                unit;
+                upper_bound = ub,
+                augmented = true,
             )
 
             status, _ = PropertyT.solve(
                 opt_problem,
-                scs_optimizer(
-                    eps=1e-10,
-                    max_iters=10_000,
-                    accel=-10,
-                    alpha=1.5,
+                scs_optimizer(;
+                    eps = 1e-10,
+                    max_iters = 10_000,
+                    accel = -10,
+                    alpha = 1.5,
                 ),
             )
 
@@ -210,13 +207,13 @@
                 elt,
                 unit,
                 λ,
-                Q,
-                halfradius=2,
-                augmented=true,
+                Q;
+                halfradius = 2,
+                augmented = true,
             )
 
             @test certified
-            @test isapprox(λ_cert, λ, rtol=1e-5)
+            @test isapprox(λ_cert, λ, rtol = 1e-5)
             @test λ_cert > 2 // 10
         end
     end
diff --git a/test/1712.07167.jl b/test/1712.07167.jl
index be479d9..abc19b8 100644
--- a/test/1712.07167.jl
+++ b/test/1712.07167.jl
@@ -1,9 +1,9 @@
 @testset "1712.07167 Examples" begin
-
     @testset "SAut(F₃)" begin
         N = 3
         G = SpecialAutomorphismGroup(FreeGroup(N))
-        RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
+        @info "running tests for" G
+        RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2)
 
         P = PermGroup(perm"(1,2)", Perm(circshift(1:N, -1)))
         Σ = PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P)
@@ -15,6 +15,7 @@
             basis(RG),
             StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]),
         )
+        @info wd
 
         Δ = let RG = RG, S = S
             RG(length(S)) - sum(RG(s) for s in S)
@@ -27,14 +28,14 @@
         status, certified, λ_cert = check_positivity(
             elt,
             unit,
-            wd,
-            upper_bound=ub,
-            halfradius=2,
-            optimizer=cosmo_optimizer(
-                eps=1e-7,
-                max_iters=10_000,
-                accel=50,
-                alpha=1.9,
+            wd;
+            upper_bound = ub,
+            halfradius = 2,
+            optimizer = cosmo_optimizer(;
+                eps = 1e-7,
+                max_iters = 10_000,
+                accel = 50,
+                alpha = 1.9,
             ),
         )
 
@@ -47,7 +48,8 @@
         n = 3
 
         SL = MatrixGroups.SpecialLinearGroup{n}(Int8)
-        RSL, S, sizes = PropertyT.group_algebra(SL, halfradius=2, twisted=true)
+        @info "running tests for" SL
+        RSL, S, sizes = PropertyT.group_algebra(SL; halfradius = 2)
 
         Δ = RSL(length(S)) - sum(RSL(s) for s in S)
 
@@ -62,6 +64,7 @@
                 basis(RSL),
                 StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]]),
             )
+            @info wd
 
             elt = Δ^2
             unit = Δ
@@ -71,8 +74,8 @@
                 elt,
                 unit,
                 wd,
-                upper_bound=ub,
-                augmented=false,
+                upper_bound = ub,
+                augmented = false,
             )
 
             wdfl = SymbolicWedderburn.WedderburnDecomposition(
@@ -86,18 +89,18 @@
             model, varP = PropertyT.sos_problem_primal(
                 elt,
                 unit,
-                wdfl,
-                upper_bound=ub,
-                augmented=false,
+                wdfl;
+                upper_bound = ub,
+                augmented = false,
             )
 
             status, warm = PropertyT.solve(
                 model,
-                cosmo_optimizer(
-                    eps=1e-10,
-                    max_iters=20_000,
-                    accel=50,
-                    alpha=1.9,
+                cosmo_optimizer(;
+                    eps = 1e-10,
+                    max_iters = 20_000,
+                    accel = 50,
+                    alpha = 1.9,
                 ),
             )
 
@@ -105,34 +108,34 @@
 
             status, _ = PropertyT.solve(
                 model,
-                scs_optimizer(
-                    eps=1e-10,
-                    max_iters=100,
-                    accel=-20,
-                    alpha=1.2,
+                scs_optimizer(;
+                    eps = 1e-10,
+                    max_iters = 100,
+                    accel = -20,
+                    alpha = 1.2,
                 ),
-                warm
+                warm,
             )
 
             @test status == JuMP.OPTIMAL
 
             Q = @time let varP = varP
                 Qs = map(varP) do P
-                    real.(sqrt(JuMP.value.(P)))
+                    return real.(sqrt(JuMP.value.(P)))
                 end
                 PropertyT.reconstruct(Qs, wdfl)
             end
             λ = JuMP.value(model[:λ])
 
-            sos = PropertyT.compute_sos(parent(elt), Q; augmented=false)
+            sos = PropertyT.compute_sos(parent(elt), Q; augmented = false)
 
             certified, λ_cert = PropertyT.certify_solution(
                 elt,
                 unit,
                 λ,
-                Q,
-                halfradius=2,
-                augmented=false,
+                Q;
+                halfradius = 2,
+                augmented = false,
             )
 
             @test certified
@@ -155,22 +158,23 @@
                 basis(RSL),
                 StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]]),
             )
+            @info wdfl
 
             opt_problem, varP = PropertyT.sos_problem_primal(
                 elt,
                 unit,
-                wdfl,
-                upper_bound=ub,
+                wdfl;
+                upper_bound = ub,
                 # augmented = true # since both elt and unit are augmented
             )
 
             status, _ = PropertyT.solve(
                 opt_problem,
-                scs_optimizer(
-                    eps=1e-8,
-                    max_iters=20_000,
-                    accel=0,
-                    alpha=1.9,
+                scs_optimizer(;
+                    eps = 1e-8,
+                    max_iters = 20_000,
+                    accel = 0,
+                    alpha = 1.9,
                 ),
             )
 
@@ -178,7 +182,7 @@
 
             Q = @time let varP = varP
                 Qs = map(varP) do P
-                    real.(sqrt(JuMP.value.(P)))
+                    return real.(sqrt(JuMP.value.(P)))
                 end
                 PropertyT.reconstruct(Qs, wdfl)
             end
@@ -187,8 +191,8 @@
                 elt,
                 unit,
                 JuMP.objective_value(opt_problem),
-                Q,
-                halfradius=2,
+                Q;
+                halfradius = 2,
                 # augmented = true # since both elt and unit are augmented
             )
 
diff --git a/test/1812.03456.jl b/test/1812.03456.jl
index b0d9482..383e01e 100644
--- a/test/1812.03456.jl
+++ b/test/1812.03456.jl
@@ -21,8 +21,9 @@ using SparseArrays
     @testset "Sq, Adj, Op in SL(4,Z)" begin
         N = 4
         G = MatrixGroups.SpecialLinearGroup{N}(Int8)
+        @info "running tests for" G
 
-        RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
+        RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2)
 
         Δ = let RG = RG, S = S
             RG(length(S)) - sum(RG(s) for s in S)
@@ -39,6 +40,7 @@ using SparseArrays
             basis(RG),
             StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]),
         )
+        @info wd
         ivs = SymbolicWedderburn.invariant_vectors(wd)
 
         sq, adj, op = PropertyT.SqAdjOp(RG, N)
@@ -82,7 +84,7 @@ using SparseArrays
     @testset "SAut(F₃)" begin
         n = 3
         G = SpecialAutomorphismGroup(FreeGroup(n))
-        RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
+        RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2)
         sq, adj, op = PropertyT.SqAdjOp(RG, n)
 
         @test sq(one(G)) == 216
@@ -98,7 +100,8 @@ end
         n = 3
 
         G = MatrixGroups.SpecialLinearGroup{n}(Int8)
-        RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
+        @info "running tests for" G
+        RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2)
 
         Δ = RG(length(S)) - sum(RG(s) for s in S)
 
@@ -113,6 +116,7 @@ end
             basis(RG),
             StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]),
         )
+        @info wd
 
         sq, adj, op = PropertyT.SqAdjOp(RG, n)
 
@@ -123,10 +127,10 @@ end
             status, certified, λ_cert = check_positivity(
                 elt,
                 Δ,
-                wd,
-                upper_bound=UB,
-                halfradius=2,
-                optimizer=cosmo_optimizer(accel=50, alpha=1.9)
+                wd;
+                upper_bound = UB,
+                halfradius = 2,
+                optimizer = cosmo_optimizer(; accel = 50, alpha = 1.9),
             )
             @test status == JuMP.OPTIMAL
             @test certified
@@ -140,10 +144,10 @@ end
             status, certified, λ_cert = check_positivity(
                 elt,
                 Δ,
-                wd,
-                upper_bound=UB,
-                halfradius=2,
-                optimizer=cosmo_optimizer(accel=50, alpha=1.9)
+                wd;
+                upper_bound = UB,
+                halfradius = 2,
+                optimizer = cosmo_optimizer(; accel = 50, alpha = 1.9),
             )
             @test status == JuMP.OPTIMAL
             @test certified
@@ -152,10 +156,11 @@ end
             m, _ = PropertyT.sos_problem_primal(elt, wd)
             PropertyT.solve(
                 m,
-                scs_optimizer(max_iters=5000, accel=50, alpha=1.9)
+                scs_optimizer(; max_iters = 5000, accel = 50, alpha = 1.9),
             )
 
-            @test JuMP.termination_status(m) in (JuMP.ALMOST_OPTIMAL, JuMP.OPTIMAL, JuMP.ITERATION_LIMIT)
+            @test JuMP.termination_status(m) in
+                  (JuMP.ALMOST_OPTIMAL, JuMP.OPTIMAL, JuMP.ITERATION_LIMIT)
             @test abs(JuMP.objective_value(m)) < 1e-3
         end
 
@@ -168,10 +173,10 @@ end
             status, certified, λ_cert = check_positivity(
                 elt,
                 Δ,
-                wd,
-                upper_bound=UB,
-                halfradius=2,
-                optimizer=cosmo_optimizer(accel=50, alpha=1.9)
+                wd;
+                upper_bound = UB,
+                halfradius = 2,
+                optimizer = scs_optimizer(; accel = 50, alpha = 1.9),
             )
             @test status == JuMP.OPTIMAL
             @test !certified
@@ -183,7 +188,8 @@ end
         n = 4
 
         G = MatrixGroups.SpecialLinearGroup{n}(Int8)
-        RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
+        @info "running tests for" G
+        RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2)
 
         Δ = RG(length(S)) - sum(RG(s) for s in S)
 
@@ -198,6 +204,7 @@ end
             basis(RG),
             StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]),
         )
+        @info wd
 
         sq, adj, op = PropertyT.SqAdjOp(RG, n)
 
@@ -208,10 +215,10 @@ end
             status, certified, λ_cert = check_positivity(
                 elt,
                 Δ,
-                wd,
-                upper_bound=UB,
-                halfradius=2,
-                optimizer=cosmo_optimizer(accel=50, alpha=1.9)
+                wd;
+                upper_bound = UB,
+                halfradius = 2,
+                optimizer = cosmo_optimizer(; accel = 50, alpha = 1.9),
             )
             @test status == JuMP.OPTIMAL
             @test certified
@@ -225,10 +232,10 @@ end
             status, certified, λ_cert = check_positivity(
                 elt,
                 Δ,
-                wd,
-                upper_bound=UB,
-                halfradius=2,
-                optimizer=cosmo_optimizer(accel=50, alpha=1.9)
+                wd;
+                upper_bound = UB,
+                halfradius = 2,
+                optimizer = cosmo_optimizer(; accel = 50, alpha = 1.9),
             )
             @test status == JuMP.OPTIMAL
             @test certified
@@ -242,10 +249,10 @@ end
             status, certified, λ_cert = check_positivity(
                 elt,
                 Δ,
-                wd,
-                upper_bound=UB,
-                halfradius=2,
-                optimizer=cosmo_optimizer(accel=50, alpha=1.9)
+                wd;
+                upper_bound = UB,
+                halfradius = 2,
+                optimizer = cosmo_optimizer(; accel = 50, alpha = 1.9),
             )
             @test status == JuMP.OPTIMAL
             @test !certified
diff --git a/test/Chevalley.jl b/test/Chevalley.jl
new file mode 100644
index 0000000..e4fa863
--- /dev/null
+++ b/test/Chevalley.jl
@@ -0,0 +1,168 @@
+countmap(v) = countmap(identity, v)
+function countmap(f, v)
+    counts = Dict{eltype(f(first(v))),Int}()
+    for x in v
+        fx = f(x)
+        counts[fx] = get!(counts, fx, 0) + 1
+    end
+    return counts
+end
+
+@testset "classify_root_system" begin
+    α = PropertyT.Roots.Root([1, -1, 0])
+    β = PropertyT.Roots.Root([0, 1, -1])
+    γ = PropertyT.Roots.Root([2, 0, 0])
+
+    @test PropertyT.Roots.classify_root_system(α, β, (false, false)) == :A₂
+    @test PropertyT.Roots.classify_root_system(α, γ, (false, true)) == :C₂
+    @test PropertyT.Roots.classify_root_system(β, γ, (false, true)) ==
+          Symbol("A₁×C₁")
+end
+
+@testset "Exceptional root systems" begin
+    @testset "F4" begin
+        F4 = let Σ = PermutationGroups.PermGroup(perm"(1,2,3,4)", perm"(1,2)")
+            long = let x = (1, 1, 0, 0) .// 1
+                PropertyT.Roots.Root.(
+                    union(
+                        (x^g for g in Σ),
+                        ((x .* (-1, 1, 1, 1))^g for g in Σ),
+                        ((-1 .* x)^g for g in Σ),
+                    ),
+                )
+            end
+
+            short = let x = (1, 0, 0, 0) .// 1
+                PropertyT.Roots.Root.(
+                    union((x^g for g in Σ), ((-1 .* x)^g for g in Σ))
+                )
+            end
+
+            signs = collect(Iterators.product(fill([-1, +1], 4)...))
+            halfs = let x = (1, 1, 1, 1) .// 2
+                PropertyT.Roots.Root.(union(x .* sgn for sgn in signs))
+            end
+
+            union(long, short, halfs)
+        end
+
+        @test length(F4) == 48
+
+        a = F4[1]
+        @test isapprox(PropertyT.Roots.ℓ₂length(a), sqrt(2))
+        b = F4[6]
+        @test isapprox(PropertyT.Roots.ℓ₂length(b), sqrt(2))
+        c = a + b
+        @test isapprox(PropertyT.Roots.ℓ₂length(c), 2.0)
+        @test PropertyT.Roots.classify_root_system(b, c, (false, true)) == :C₂
+
+        long = F4[findfirst(r -> PropertyT.Roots.ℓ₂length(r) == sqrt(2), F4)]
+        short = F4[findfirst(r -> PropertyT.Roots.ℓ₂length(r) == 1.0, F4)]
+
+        subtypes = Set([:C₂, :A₂, Symbol("A₁×C₁")])
+
+        let Ω = F4, α = long
+            counts = countmap([
+                PropertyT.Roots.classify_sub_root_system(Ω, α, γ) for
+                γ in Ω if !PropertyT.Roots.isproportional(α, γ)
+            ])
+            @test Set(keys(counts)) == subtypes
+            d, r = divrem(counts[:C₂], 6)
+            @test r == 0 && d == 3
+
+            d, r = divrem(counts[:A₂], 4)
+            @test r == 0 && d == 4
+        end
+
+        let Ω = F4, α = short
+            counts = countmap([
+                PropertyT.Roots.classify_sub_root_system(Ω, α, γ) for
+                γ in Ω if !PropertyT.Roots.isproportional(α, γ)
+            ])
+            @test Set(keys(counts)) == subtypes
+            d, r = divrem(counts[:C₂], 6)
+            @test r == 0 && d == 3
+
+            d, r = divrem(counts[:A₂], 4)
+            @test r == 0 && d == 4
+        end
+    end
+
+    @testset "E6-7-8 exceptional root systems" begin
+        E8 =
+            let Σ = PermutationGroups.PermGroup(
+                    perm"(1,2,3,4,5,6,7,8)",
+                    perm"(1,2)",
+                )
+                long = let x = (1, 1, 0, 0, 0, 0, 0, 0) .// 1
+                    PropertyT.Roots.Root.(
+                        union(
+                            (x^g for g in Σ),
+                            ((x .* (-1, 1, 1, 1, 1, 1, 1, 1))^g for g in Σ),
+                            ((-1 .* x)^g for g in Σ),
+                        ),
+                    )
+                end
+
+                signs = collect(
+                    p for p in Iterators.product(fill([-1, +1], 8)...) if
+                    iseven(count(==(-1), p))
+                )
+                halfs = let x = (1, 1, 1, 1, 1, 1, 1, 1) .// 2
+                    rts = unique(PropertyT.Roots.Root(x .* sgn) for sgn in signs)
+                end
+
+                union(long, halfs)
+            end
+
+        subtypes = Set([:A₂, Symbol("A₁×A₁")])
+
+        @testset "E8" begin
+            @test length(E8) == 240
+            @test all(r -> PropertyT.Roots.ℓ₂length(r) ≈ sqrt(2), E8)
+
+            let Ω = E8, α = first(Ω)
+                counts = countmap([
+                    PropertyT.Roots.classify_sub_root_system(Ω, α, γ) for
+                    γ in Ω if !PropertyT.Roots.isproportional(α, γ)
+                ])
+                @test Set(keys(counts)) == subtypes
+                d, r = divrem(counts[:A₂], 4)
+                @test r == 0 && d == 28
+            end
+        end
+        @testset "E7" begin
+            E7 = filter(r -> iszero(sum(r.coord)), E8)
+            @test length(E7) == 126
+
+            let Ω = E7, α = first(Ω)
+                counts = countmap([
+                    PropertyT.Roots.classify_sub_root_system(Ω, α, γ) for
+                    γ in Ω if !PropertyT.Roots.isproportional(α, γ)
+                ])
+                @test Set(keys(counts)) == subtypes
+                d, r = divrem(counts[:A₂], 4)
+                @test r == 0 && d == 16
+            end
+        end
+
+        @testset "E6" begin
+            E6 = filter(
+                r -> r.coord[end] == r.coord[end-1] == r.coord[end-2],
+                E8,
+            )
+            @test length(E6) == 72
+
+            let Ω = E6, α = first(Ω)
+                counts = countmap([
+                    PropertyT.Roots.classify_sub_root_system(Ω, α, γ) for
+                    γ in Ω if !PropertyT.Roots.isproportional(α, γ)
+                ])
+                @test Set(keys(counts)) == subtypes
+                d, r = divrem(counts[:A₂], 4)
+                @info d, r
+                @test r == 0 && d == 10
+            end
+        end
+    end
+end
diff --git a/test/actions.jl b/test/actions.jl
index af7d123..08e6a03 100644
--- a/test/actions.jl
+++ b/test/actions.jl
@@ -3,7 +3,7 @@ function test_action(basis, group, act)
     return @testset "action definition" begin
         @test all(basis) do b
             e = one(group)
-            action(act, e, b) == b
+            return action(act, e, b) == b
         end
 
         a = let a = rand(basis)
@@ -30,12 +30,12 @@ function test_action(basis, group, act)
         @test all([(g, h) for g in group for h in group]) do (g, h)
             x = action(act, h, action(act, g, a))
             y = action(act, g * h, a)
-            x == y
+            return x == y
         end
 
         if act isa SymbolicWedderburn.ByPermutations
             @test all(basis) do b
-                action(act, g, b) ∈ basis && action(act, h, b) ∈ basis
+                return action(act, g, b) ∈ basis && action(act, h, b) ∈ basis
             end
         end
     end
@@ -44,21 +44,18 @@ end
 ## Testing
 
 @testset "Actions on SL(3,ℤ)" begin
-
     n = 3
 
     SL = MatrixGroups.SpecialLinearGroup{n}(Int8)
-    RSL, S, sizes = PropertyT.group_algebra(SL, halfradius=2, twisted=true)
+    RSL, S, sizes = PropertyT.group_algebra(SL; halfradius = 2)
 
     @testset "Permutation action" begin
-
         Γ = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1)))
         ΓpA = PropertyT.action_by_conjugation(SL, Γ)
 
         test_action(basis(RSL), Γ, ΓpA)
 
         @testset "mps is successful" begin
-
             charsΓ =
                 SymbolicWedderburn.Character{
                     Rational{Int},
@@ -66,9 +63,9 @@ end
 
             RΓ = SymbolicWedderburn._group_algebra(Γ)
 
-            @time mps, simple =
+            @time mps, ranks =
                 SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ)
-            @test all(simple)
+            @test all(isone, ranks)
         end
 
         @testset "Wedderburn decomposition" begin
@@ -77,17 +74,17 @@ end
                 Γ,
                 ΓpA,
                 basis(RSL),
-                StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]])
+                StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]]),
             )
 
             @test length(invariant_vectors(wd)) == 918
-            @test SymbolicWedderburn.size.(direct_summands(wd), 1) == [40, 23, 18]
+            @test SymbolicWedderburn.size.(direct_summands(wd), 1) ==
+                  [40, 23, 18]
             @test all(issimple, direct_summands(wd))
         end
     end
 
     @testset "Wreath action" begin
-
         Γ = let P = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1)))
             PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P)
         end
@@ -97,7 +94,6 @@ end
         test_action(basis(RSL), Γ, ΓpA)
 
         @testset "mps is successful" begin
-
             charsΓ =
                 SymbolicWedderburn.Character{
                     Rational{Int},
@@ -105,9 +101,9 @@ end
 
             RΓ = SymbolicWedderburn._group_algebra(Γ)
 
-            @time mps, simple =
+            @time mps, ranks =
                 SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ)
-            @test all(simple)
+            @test all(isone, ranks)
         end
 
         @testset "Wedderburn decomposition" begin
@@ -116,32 +112,30 @@ end
                 Γ,
                 ΓpA,
                 basis(RSL),
-                StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]])
+                StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]]),
             )
 
             @test length(invariant_vectors(wd)) == 247
-            @test SymbolicWedderburn.size.(direct_summands(wd), 1) == [14, 9, 6, 14, 12]
+            @test SymbolicWedderburn.size.(direct_summands(wd), 1) ==
+                  [14, 9, 6, 14, 12]
             @test all(issimple, direct_summands(wd))
         end
     end
 end
 
 @testset "Actions on SAut(F4)" begin
-
     n = 4
 
     SAutFn = SpecialAutomorphismGroup(FreeGroup(n))
-    RSAutFn, S, sizes = PropertyT.group_algebra(SAutFn, halfradius=1, twisted=true)
+    RSAutFn, S, sizes = PropertyT.group_algebra(SAutFn; halfradius = 1)
 
     @testset "Permutation action" begin
-
         Γ = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1)))
         ΓpA = PropertyT.action_by_conjugation(SAutFn, Γ)
 
         test_action(basis(RSAutFn), Γ, ΓpA)
 
         @testset "mps is successful" begin
-
             charsΓ =
                 SymbolicWedderburn.Character{
                     Rational{Int},
@@ -149,9 +143,9 @@ end
 
             RΓ = SymbolicWedderburn._group_algebra(Γ)
 
-            @time mps, simple =
+            @time mps, ranks =
                 SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ)
-            @test all(simple)
+            @test all(isone, ranks)
         end
 
         @testset "Wedderburn decomposition" begin
@@ -160,17 +154,17 @@ end
                 Γ,
                 ΓpA,
                 basis(RSAutFn),
-                StarAlgebras.Basis{UInt16}(@view basis(RSAutFn)[1:sizes[1]])
+                StarAlgebras.Basis{UInt16}(@view basis(RSAutFn)[1:sizes[1]]),
             )
 
             @test length(invariant_vectors(wd)) == 93
-            @test SymbolicWedderburn.size.(direct_summands(wd), 1) == [4, 8, 5, 4]
+            @test SymbolicWedderburn.size.(direct_summands(wd), 1) ==
+                  [4, 8, 5, 4]
             @test all(issimple, direct_summands(wd))
         end
     end
 
     @testset "Wreath action" begin
-
         Γ = let P = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1)))
             PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P)
         end
@@ -180,7 +174,6 @@ end
         test_action(basis(RSAutFn), Γ, ΓpA)
 
         @testset "mps is successful" begin
-
             charsΓ =
                 SymbolicWedderburn.Character{
                     Rational{Int},
@@ -188,9 +181,9 @@ end
 
             RΓ = SymbolicWedderburn._group_algebra(Γ)
 
-            @time mps, simple =
+            @time mps, ranks =
                 SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ)
-            @test all(simple)
+            @test all(isone, ranks)
         end
 
         @testset "Wedderburn decomposition" begin
@@ -199,11 +192,12 @@ end
                 Γ,
                 ΓpA,
                 basis(RSAutFn),
-                StarAlgebras.Basis{UInt16}(@view basis(RSAutFn)[1:sizes[1]])
+                StarAlgebras.Basis{UInt16}(@view basis(RSAutFn)[1:sizes[1]]),
             )
 
             @test length(invariant_vectors(wd)) == 18
-            @test SymbolicWedderburn.size.(direct_summands(wd), 1) == [1, 1, 2, 2, 1, 2, 2, 1]
+            @test SymbolicWedderburn.size.(direct_summands(wd), 1) ==
+                  [1, 1, 2, 2, 1, 2, 2, 1]
             @test all(issimple, direct_summands(wd))
         end
     end
diff --git a/test/check_positivity.jl b/test/check_positivity.jl
index 974077b..d57e80f 100644
--- a/test/check_positivity.jl
+++ b/test/check_positivity.jl
@@ -1,6 +1,12 @@
-function check_positivity(elt, unit; upper_bound=Inf, halfradius=2, optimizer)
+function check_positivity(
+    elt,
+    unit;
+    upper_bound = Inf,
+    halfradius = 2,
+    optimizer,
+)
     @time sos_problem =
-        PropertyT.sos_problem_primal(elt, unit, upper_bound=upper_bound)
+        PropertyT.sos_problem_primal(elt, unit; upper_bound = upper_bound)
 
     status, _ = PropertyT.solve(sos_problem, optimizer)
     P = JuMP.value.(sos_problem[:P])
@@ -9,16 +15,23 @@ function check_positivity(elt, unit; upper_bound=Inf, halfradius=2, optimizer)
         elt,
         unit,
         JuMP.objective_value(sos_problem),
-        Q,
-        halfradius=halfradius,
+        Q;
+        halfradius = halfradius,
     )
     return status, certified, λ_cert
 end
 
-function check_positivity(elt, unit, wd; upper_bound=Inf, halfradius=2, optimizer)
+function check_positivity(
+    elt,
+    unit,
+    wd;
+    upper_bound = Inf,
+    halfradius = 2,
+    optimizer,
+)
     @assert aug(elt) == aug(unit) == 0
     @time sos_problem, Ps =
-        PropertyT.sos_problem_primal(elt, unit, wd, upper_bound=upper_bound)
+        PropertyT.sos_problem_primal(elt, unit, wd; upper_bound = upper_bound)
 
     @time status, _ = PropertyT.solve(sos_problem, optimizer)
 
@@ -29,13 +42,7 @@ function check_positivity(elt, unit, wd; upper_bound=Inf, halfradius=2, optimize
 
     λ = JuMP.value(sos_problem[:λ])
 
-    certified, λ_cert = PropertyT.certify_solution(
-        elt,
-        unit,
-        λ,
-        Q,
-        halfradius=halfradius
-    )
+    certified, λ_cert =
+        PropertyT.certify_solution(elt, unit, λ, Q; halfradius = halfradius)
     return status, certified, λ_cert
 end
-
diff --git a/test/constratint_matrices.jl b/test/constratint_matrices.jl
index afdee61..1d04876 100644
--- a/test/constratint_matrices.jl
+++ b/test/constratint_matrices.jl
@@ -1,12 +1,17 @@
 @testset "ConstraintMatrix" begin
-    @test PropertyT.ConstraintMatrix{Float64}([-1, 2, -1, 1, 4, 2, 6], 3, 2, π) isa AbstractMatrix
+    @test PropertyT.ConstraintMatrix{Float64}(
+        [-1, 2, -1, 1, 4, 2, 6],
+        3,
+        2,
+        π,
+    ) isa AbstractMatrix
 
     cm = PropertyT.ConstraintMatrix{Float64}([-1, 2, -1, 1, 4, 2, 6], 3, 2, π)
 
     @test cm == Float64[
         -π π
-        2π 0.0
-        0.0 π
+        2π 0
+        0 π
     ]
 
     @test collect(PropertyT.nzpairs(cm)) == [
@@ -18,4 +23,7 @@
         1 => -3.141592653589793
         1 => -3.141592653589793
     ]
+
+    @test PropertyT.ConstraintMatrix{Float64}([-9:-1; 1:9], 3, 3, 1.0) ==
+          zeros(3, 3)
 end
diff --git a/test/graded_adj.jl b/test/graded_adj.jl
index bf5a804..447b58c 100644
--- a/test/graded_adj.jl
+++ b/test/graded_adj.jl
@@ -1,10 +1,9 @@
 @testset "Adj via grading" begin
-
     @testset "SL(n,Z) & Aut(F₄)" begin
         n = 4
         halfradius = 1
         SL = MatrixGroups.SpecialLinearGroup{n}(Int8)
-        RSL, S, sizes = PropertyT.group_algebra(SL, halfradius=halfradius, twisted=true)
+        RSL, S, sizes = PropertyT.group_algebra(SL; halfradius = halfradius)
 
         Δ = RSL(length(S)) - sum(RSL(s) for s in S)
 
@@ -22,10 +21,9 @@
         @test PropertyT.Adj(Δs, :A₂) == adj
         @test PropertyT.Adj(Δs, Symbol("A₁×A₁")) == op
 
-
         halfradius = 1
         G = SpecialAutomorphismGroup(FreeGroup(n))
-        RG, S, sizes = PropertyT.group_algebra(G, halfradius=halfradius, twisted=true)
+        RG, S, sizes = PropertyT.group_algebra(G; halfradius = halfradius)
 
         Δ = RG(length(S)) - sum(RG(s) for s in S)
 
@@ -44,7 +42,6 @@
         @test PropertyT.Adj(Δs, Symbol("A₁×A₁")) == op
     end
 
-
     @testset "Symplectic group" begin
         @testset "Sp2(ℤ)" begin
             genus = 2
@@ -52,7 +49,8 @@
 
             SpN = MatrixGroups.SymplecticGroup{2genus}(Int8)
 
-            RSpN, S_sp, sizes_sp = PropertyT.group_algebra(SpN, halfradius=halfradius, twisted=true)
+            RSpN, S_sp, sizes_sp =
+                PropertyT.group_algebra(SpN; halfradius = halfradius)
 
             Δ, Δs = let RG = RSpN, S = S_sp, ψ = identity
                 Δ = RG(length(S)) - sum(RG(s) for s in S)
@@ -73,7 +71,8 @@
 
         SpN = MatrixGroups.SymplecticGroup{2genus}(Int8)
 
-        RSpN, S_sp, sizes_sp = PropertyT.group_algebra(SpN, halfradius=halfradius, twisted=true)
+        RSpN, S_sp, sizes_sp =
+            PropertyT.group_algebra(SpN; halfradius = halfradius)
 
         Δ, Δs = let RG = RSpN, S = S_sp, ψ = identity
             Δ = RG(length(S)) - sum(RG(s) for s in S)
@@ -86,9 +85,14 @@
         end
 
         @testset "Adj numerics for genus=$genus" begin
-
             all_subtypes = (
-                :A₁, :C₁, Symbol("A₁×A₁"), Symbol("C₁×C₁"), Symbol("A₁×C₁"), :A₂, :C₂
+                :A₁,
+                :C₁,
+                Symbol("A₁×A₁"),
+                Symbol("C₁×C₁"),
+                Symbol("A₁×C₁"),
+                :A₂,
+                :C₂,
             )
 
             @test PropertyT.Adj(Δs, :A₂)[one(SpN)] == 384
@@ -101,8 +105,8 @@
                     @test isinteger(PropertyT.Adj(Δs, subtype)[one(SpN)] / 16)
                 end
             end
-            @test sum(PropertyT.Adj(Δs, subtype) for subtype in all_subtypes) == Δ^2
+            @test sum(PropertyT.Adj(Δs, subtype) for subtype in all_subtypes) ==
+                  Δ^2
         end
     end
 end
-
diff --git a/test/optimizers.jl b/test/optimizers.jl
index f415de7..01b45db 100644
--- a/test/optimizers.jl
+++ b/test/optimizers.jl
@@ -4,12 +4,12 @@ import JuMP
 import SCS
 
 function scs_optimizer(;
-    accel=10,
-    alpha=1.5,
-    eps=1e-9,
-    max_iters=100_000,
-    verbose=true,
-    linear_solver=SCS.DirectSolver
+    accel = 10,
+    alpha = 1.5,
+    eps = 1e-9,
+    max_iters = 100_000,
+    verbose = true,
+    linear_solver = SCS.DirectSolver,
 )
     return JuMP.optimizer_with_attributes(
         SCS.Optimizer,
@@ -28,27 +28,27 @@ end
 import COSMO
 
 function cosmo_optimizer(;
-    accel=15,
-    alpha=1.6,
-    eps=1e-9,
-    max_iters=100_000,
-    verbose=true,
-    verbose_timing=verbose,
-    decompose=false
+    accel = 15,
+    alpha = 1.6,
+    eps = 1e-9,
+    max_iters = 100_000,
+    verbose = true,
+    verbose_timing = verbose,
+    decompose = false,
 )
     return JuMP.optimizer_with_attributes(
         COSMO.Optimizer,
-        "accelerator" => COSMO.with_options(
-            COSMO.AndersonAccelerator,
-            mem=max(accel, 2)
-        ),
+        "accelerator" =>
+            COSMO.with_options(COSMO.AndersonAccelerator; mem = max(accel, 2)),
         "alpha" => alpha,
         "decompose" => decompose,
         "eps_abs" => eps,
         "eps_rel" => eps,
+        "eps_prim_inf" => eps,
+        "eps_dual_inf" => eps,
         "max_iter" => max_iters,
         "verbose" => verbose,
         "verbose_timing" => verbose_timing,
-        "check_termination" => 200,
+        "check_termination" => 250,
     )
 end
diff --git a/test/quick_tests.jl b/test/quick_tests.jl
index 9a38b5c..9478286 100644
--- a/test/quick_tests.jl
+++ b/test/quick_tests.jl
@@ -1,11 +1,12 @@
 @testset "Quick tests" begin
-
     @testset "SL(2,F₇)" begin
         N = 2
         p = 7
         halfradius = 3
-        G = MatrixGroups.SpecialLinearGroup{N}(SymbolicWedderburn.Characters.FiniteFields.GF{p})
-        RG, S, sizes = PropertyT.group_algebra(G, halfradius=3, twisted=true)
+        G = MatrixGroups.SpecialLinearGroup{N}(
+            SymbolicWedderburn.Characters.FiniteFields.GF{p},
+        )
+        RG, S, sizes = PropertyT.group_algebra(G; halfradius = 3)
 
         Δ = let RG = RG, S = S
             RG(length(S)) - sum(RG(s) for s in S)
@@ -18,15 +19,15 @@
         @testset "standard formulation" begin
             status, certified, λ_cert = check_positivity(
                 elt,
-                unit,
-                upper_bound=ub,
-                halfradius=2,
-                optimizer=cosmo_optimizer(
-                    eps=1e-7,
-                    max_iters=5_000,
-                    accel=50,
-                    alpha=1.95,
-                )
+                unit;
+                upper_bound = ub,
+                halfradius = 2,
+                optimizer = cosmo_optimizer(;
+                    eps = 1e-7,
+                    max_iters = 5_000,
+                    accel = 50,
+                    alpha = 1.95,
+                ),
             )
 
             @test status == JuMP.OPTIMAL
@@ -34,14 +35,18 @@
             @test λ_cert > 5857 // 10000
 
             m = PropertyT.sos_problem_dual(elt, unit)
-            PropertyT.solve(m, cosmo_optimizer(
-                eps=1e-7,
-                max_iters=10_000,
-                accel=50,
-                alpha=1.95,
-            ))
+            PropertyT.solve(
+                m,
+                cosmo_optimizer(;
+                    eps = 1e-7,
+                    max_iters = 10_000,
+                    accel = 50,
+                    alpha = 1.95,
+                ),
+            )
 
-            @test JuMP.termination_status(m) in (JuMP.ALMOST_OPTIMAL, JuMP.OPTIMAL)
+            @test JuMP.termination_status(m) in
+                  (JuMP.ALMOST_OPTIMAL, JuMP.OPTIMAL)
             @test JuMP.objective_value(m) ≈ λ_cert atol = 1e-2
         end
 
@@ -55,20 +60,22 @@
                 Σ,
                 act,
                 basis(RG),
-                StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[halfradius]]),
+                StarAlgebras.Basis{UInt16}(
+                    @view basis(RG)[1:sizes[halfradius]]
+                ),
             )
 
             status, certified, λ_cert = check_positivity(
                 elt,
                 unit,
-                wd,
-                upper_bound=ub,
-                halfradius=2,
-                optimizer=cosmo_optimizer(
-                    eps=1e-7,
-                    max_iters=10_000,
-                    accel=50,
-                    alpha=1.9,
+                wd;
+                upper_bound = ub,
+                halfradius = 2,
+                optimizer = cosmo_optimizer(;
+                    eps = 1e-7,
+                    max_iters = 10_000,
+                    accel = 50,
+                    alpha = 1.9,
                 ),
             )
 
diff --git a/test/roots.jl b/test/roots.jl
new file mode 100644
index 0000000..7cb32d0
--- /dev/null
+++ b/test/roots.jl
@@ -0,0 +1,10 @@
+using PropertyT.Roots
+@testset "Roots" begin
+    @test Roots.Root{3,Int}([1, 2, 3]) isa Roots.AbstractRoot{}
+    @test Roots.Root([1, 2, 3]) isa Roots.AbstractRoot{3,Int}
+    # io
+    r = Roots.Root{3,Int}([1, 2, 3])
+    @test contains(sprint(show, MIME"text/plain"(), r), "of length √14\n")
+    r = Roots.Root{3,Int}([1, 2, 2])
+    @test contains(sprint(show, MIME"text/plain"(), r), "of length 3\n")
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 7003c1e..a306e55 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -24,6 +24,8 @@ if haskey(ENV, "FULL_TEST") || haskey(ENV, "CI")
         include("1712.07167.jl")
         include("1812.03456.jl")
 
+        include("roots.jl")
         include("graded_adj.jl")
+        include("Chevalley.jl")
     end
 end