use Cartan matrix to classify root-subsystems

src/roots.jl

@@ -71,39 +71,55 @@ end
 𝕖(N, i) = Root(ntuple(k -> k == i ? 1 : 0, N))
 𝕆(N, ::Type{T}) where {T} = Root(ntuple(_ -> zero(T), N))
 
+reflection(α::Root, β::Root) = β - Int(2dot(α, β) / dot(α, α)) * α
+function cartan(α, β)
+    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
@@ -130,12 +146,17 @@ function Base.in(r::R, plane::Plane{R}) where {R}
     return any(isproportional(r, v) for v in plane.vectors)
 end
 
+function _islong(α::Root, Ω)
+    lα = length(α)
+    return any(r -> lα - length(r) > eps(lα), Ω)
+end
+
 function classify_sub_root_system(
     Ω::AbstractVector{<:Root{N}},
     α::Root{N},
     β::Root{N},
 ) where {N}
-
+    @assert 1 ≤ length(unique(length, Ω)) ≤ 2
     v = proportional_root_from_system(Ω, α)
     w = proportional_root_from_system(Ω, β)
 
@@ -146,28 +167,30 @@ 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" α β