reorganize Roots module

src/roots.jl

@@ -7,73 +7,48 @@ export Root, isproportional, isorthogonal, ~, ⟂
 
 abstract type AbstractRoot{N,T} end
 
-struct Root{N,T} <: AbstractRoot{N,T}
+ℓ₂length(r::AbstractRoot) = norm(r, 2)
-    coord::SVector{N,T}
+ambient_dim(r::AbstractRoot) = length(r)
-end
+Base.:*(r::AbstractRoot, a::Number) = a * r
-
-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)
 
 cos_angle(a, b) = dot(a, b) / (norm(a) * norm(b))
 
-function isproportional(α::AbstractRoot{N}, β::AbstractRoot{M}) where {N,M}
+function isproportional(α::AbstractRoot, β::AbstractRoot)
-    N == M || return false
+    ambient_dim(α) == ambient_dim(β) || return false
     val = abs(cos_angle(α, β))
     return isapprox(val, one(val); atol = eps(one(val)))
 end
 
-function isorthogonal(α::AbstractRoot{N}, β::AbstractRoot{M}) where {N,M}
+function isorthogonal(α::AbstractRoot, β::AbstractRoot)
-    N == M || return false
+    ambient_dim(α) == ambient_dim(β) || return false
     val = cos_angle(α, β)
     return isapprox(val, zero(val); atol = eps(one(val)))
 end
 
-function _positive_direction(α::Root{N}) where {N}
+function positive(roots::AbstractVector{<:AbstractRoot})
-    v = α.coord + 1 / (N * 100) * rand(N)
+    isempty(roots) && return empty(roots)
-    return Root{N,Float64}(v / norm(v, 2))
-end
-
-function positive(roots::AbstractVector{<:Root{N}}) where {N}
     pd = _positive_direction(first(roots))
     return filter(α -> dot(α, pd) > 0.0, roots)
 end
 
-function Base.show(io::IO, r::Root)
+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}
+function Base.show(io::IO, ::MIME"text/plain", r::AbstractRoot)
-    lngth² = sum(x -> x^2, r.coord)
+    l₂l = ℓ₂length(r)
-    l = isinteger(sqrt(lngth²)) ? "$(sqrt(lngth²))" : "√$(lngth²)"
+    l = isinteger(l₂l) ? "$(l₂l)" : "√$(l₂l^2)"
     return print(io, "Root in ℝ^$N of length $l\n", r.coord)
 end
 
-𝕖(N, i) = Root(ntuple(k -> k == i ? 1 : 0, N))
+function reflection(α::AbstractRoot, β::AbstractRoot)
-𝕆(N, ::Type{T}) where {T} = Root(ntuple(_ -> zero(T), N))
+    return β - Int(2dot(α, β) // dot(α, α)) * α
-
+end
-reflection(α::Root, β::Root) = β - Int(2dot(α, β) / dot(α, α)) * α
+function cartan(α::AbstractRoot, β::AbstractRoot)
-function cartan(α, β)
+    ambient_dim(α) == ambient_dim(β) || throw("incompatible ambient dimensions")
     return [
-        length(reflection(a, b) - b) / length(a) for a in (α, β), b in (α, β)
+        ℓ₂length(reflection(a, b) - b) / ℓ₂length(a) for a in (α, β),
+        b in (α, β)
     ]
 end
 
@@ -124,7 +99,10 @@ function classify_root_system(
     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 Ω = $Ω")
@@ -132,31 +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
 
-function Plane(α::Root, β::Root)
+function Plane(α::AbstractRoot, β::AbstractRoot)
     return Plane(α, β, [a * α + b * β for a in -3:3 for b in -3:3])
 end
 
-function Base.in(r::Root, plane::Plane)
+function Base.in(r::AbstractRoot, plane::Plane)
     return any(isproportional(r, v) for v in plane.vectors)
 end
 
-function _islong(α::Root, Ω)
+function _islong(α::AbstractRoot, Ω)
-    lα = length(α)
+    lα = ℓ₂length(α)
-    return any(r -> lα - length(r) > eps(lα), Ω)
+    return any(r -> lα - ℓ₂length(r) > eps(lα), Ω)
 end
 
 function classify_sub_root_system(
-    Ω::AbstractVector{<:Root{N}},
+    Ω::AbstractVector{<:AbstractRoot{N}},
-    α::Root{N},
+    α::AbstractRoot{N},
-    β::Root{N},
+    β::AbstractRoot{N},
 ) where {N}
-    @assert 1 ≤ length(unique(length, Ω)) ≤ 2
+    @assert 1 ≤ length(unique(ℓ₂length, Ω)) ≤ 2
     v = proportional_root_from_system(Ω, α)
     w = proportional_root_from_system(Ω, β)
 
@@ -197,4 +175,45 @@ function classify_sub_root_system(
     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

test/Chevalley.jl

@@ -22,7 +22,7 @@ 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.0, 1.0, 0.0, 0.0)
+            long = let x = (1, 1, 0, 0) .// 1
                 PropertyT.Roots.Root.(
                     union(
                         (x^g for g in Σ),
@@ -32,14 +32,14 @@ end
                 )
             end
 
-            short = let x = (1.0, 0.0, 0.0, 0.0)
+            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 / 2 .* (1.0, 1.0, 1.0, 1.0)
+            halfs = let x = (1, 1, 1, 1) .// 2
                 PropertyT.Roots.Root.(union(x .* sgn for sgn in signs))
             end
 
@@ -49,15 +49,15 @@ end
         @test length(F4) == 48
 
         a = F4[1]
-        @test isapprox(length(a), sqrt(2))
+        @test isapprox(PropertyT.Roots.ℓ₂length(a), sqrt(2))
         b = F4[6]
-        @test isapprox(length(b), sqrt(2))
+        @test isapprox(PropertyT.Roots.ℓ₂length(b), sqrt(2))
         c = a + b
-        @test isapprox(length(c), 2.0)
+        @test isapprox(PropertyT.Roots.ℓ₂length(c), 2.0)
         @test PropertyT.Roots.classify_root_system(b, c, (false, true)) == :C₂
 
-        long = F4[findfirst(r -> length(r) == sqrt(2), F4)]
+        long = F4[findfirst(r -> PropertyT.Roots.ℓ₂length(r) == sqrt(2), F4)]
-        short = F4[findfirst(r -> length(r) == 1.0, F4)]
+        short = F4[findfirst(r -> PropertyT.Roots.ℓ₂length(r) == 1.0, F4)]
 
         subtypes = Set([:C₂, :A₂, Symbol("A₁×C₁")])
 
@@ -94,7 +94,7 @@ end
                     perm"(1,2,3,4,5,6,7,8)",
                     perm"(1,2)",
                 )
-                long = let x = (1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
+                long = let x = (1, 1, 0, 0, 0, 0, 0, 0) .// 1
                     PropertyT.Roots.Root.(
                         union(
                             (x^g for g in Σ),
@@ -108,7 +108,7 @@ end
                     p for p in Iterators.product(fill([-1, +1], 8)...) if
                     iseven(count(==(-1), p))
                 )
-                halfs = let x = 1 / 2 .* ntuple(i -> 1.0, 8)
+                halfs = let x = (1, 1, 1, 1, 1, 1, 1, 1) .// 2
                     rts = unique(PropertyT.Roots.Root(x .* sgn) for sgn in signs)
                 end
 
@@ -119,7 +119,7 @@ end
 
         @testset "E8" begin
             @test length(E8) == 240
-            @test all(r -> length(r) ≈ sqrt(2), E8)
+            @test all(r -> PropertyT.Roots.ℓ₂length(r) ≈ sqrt(2), E8)
 
             let Ω = E8, α = first(Ω)
                 counts = countmap([