222 lines
6.3 KiB
Julia
222 lines
6.3 KiB
Julia
module Roots
|
||
|
||
using StaticArrays
|
||
using LinearAlgebra
|
||
|
||
export Root, isproportional, isorthogonal, ~, ⟂
|
||
|
||
abstract type AbstractRoot{N,T} end # <: AbstractVector{T} ?
|
||
|
||
ℓ₂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, β::AbstractRoot)
|
||
ambient_dim(α) == ambient_dim(β) || return false
|
||
val = abs(cos_angle(α, β))
|
||
return isapprox(val, one(val); atol = eps(one(val)))
|
||
end
|
||
|
||
function isorthogonal(α::AbstractRoot, β::AbstractRoot)
|
||
ambient_dim(α) == ambient_dim(β) || return false
|
||
val = cos_angle(α, β)
|
||
return isapprox(val, zero(val); atol = eps(one(val)))
|
||
end
|
||
|
||
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::AbstractRoot)
|
||
return print(io, "Root $(r.coord)")
|
||
end
|
||
|
||
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
|
||
|
||
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,
|
||
proportionality/orthogonality and Cartan matrix.
|
||
"""
|
||
function classify_root_system(
|
||
α::AbstractRoot,
|
||
β::AbstractRoot,
|
||
long::Tuple{Bool,Bool},
|
||
)
|
||
if isproportional(α, β)
|
||
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 all(long)
|
||
return Symbol("C₁×C₁")
|
||
elseif all(.!long) # both short
|
||
return Symbol("A₁×A₁")
|
||
elseif any(long)
|
||
return Symbol("A₁×C₁")
|
||
end
|
||
else # ⟨α, β⟩ is 2-dimensional, but they're not orthogonal
|
||
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 b == 1 && c == 2
|
||
return :C₂
|
||
elseif b == 1 && c == 3
|
||
return :G₂
|
||
else
|
||
@error a, b, c, d
|
||
error("Unknown root system ⟨α, β⟩:\n α = $α\n β = $β")
|
||
end
|
||
end
|
||
end
|
||
|
||
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 Ω = $Ω")
|
||
end
|
||
return Ω[k]
|
||
end
|
||
|
||
struct Plane{R<:AbstractRoot}
|
||
v1::R
|
||
v2::R
|
||
vectors::Vector{R}
|
||
end
|
||
|
||
function Plane(α::AbstractRoot, β::AbstractRoot)
|
||
return Plane(α, β, [a * α + b * β for a in -3:3 for b in -3:3])
|
||
end
|
||
|
||
function Base.in(r::AbstractRoot, plane::Plane)
|
||
return any(isproportional(r, v) for v in plane.vectors)
|
||
end
|
||
|
||
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(Ω, β)
|
||
|
||
subsystem = filter(ω -> ω in Plane(v, w), Ω)
|
||
@assert length(subsystem) > 0
|
||
subsystem = positive(union(subsystem, -1 .* subsystem))
|
||
|
||
l = length(subsystem)
|
||
if l == 1
|
||
x = first(subsystem)
|
||
long = _islong(x, Ω)
|
||
return classify_root_system(x, -x, (long, long))
|
||
elseif l == 2
|
||
x, y = subsystem
|
||
return classify_root_system(x, y, (_islong(x, Ω), _islong(y, Ω)))
|
||
elseif l == 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₂ == a == b == c # it's only A₂
|
||
return a
|
||
end
|
||
|
||
throw("Unknown subroot system! $((x,y,z))")
|
||
elseif l == 4
|
||
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
|
||
@warn subtypes
|
||
elseif l == 6
|
||
return :G₂
|
||
end
|
||
@error "Unknown root subsystem generated by" α β l
|
||
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
|