Enh/groups core

Project.toml

@@ -5,13 +5,14 @@ version = "0.1.0"
 
 [deps]
 Cyclotomics = "da8f5974-afbb-4dc8-91d8-516d5257c83b"
+GroupsCore = "d5909c97-4eac-4ecc-a3dc-fdd0858a4120"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 PermutationGroups = "8bc5a954-2dfc-11e9-10e6-cd969bffa420"
 Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"
 
 [compat]
-Cyclotomics = "0.2.1"
-PermutationGroups = "0.2.3"
+Cyclotomics = "0.2.2"
+PermutationGroups = "0.3"
 Primes = "0.4, 0.5"
 julia = "1"
 

src/SymbolicWedderburn.jl

@@ -3,8 +3,12 @@ module SymbolicWedderburn
 using LinearAlgebra
 using Primes
 
-using PermutationGroups
 using Cyclotomics
+using GroupsCore
+
+import PermutationGroups:
+    AbstractOrbit, AbstractPerm, AbstractPermutationGroup, Orbit, Perm, degree
+import PermutationGroups
 
 export conjugacy_classes, symmetry_adapted_basis
 

src/actions.jl

@@ -1,31 +1,32 @@
 abstract type AbstractActionHomomorphism{T} end
 
-struct ExtensionHomomorphism{T,V,Op} <: AbstractActionHomomorphism{T}
+struct ExtensionHomomorphism{T,I<:Integer,V,Op} <: AbstractActionHomomorphism{T}
     features::V
-    reversef::Dict{T,Int}
+    reversef::Dict{T,I}
     op::Op
 end
 
-ExtensionHomomorphism(features, op) = ExtensionHomomorphism(
-    features,
-    Dict(f => idx for (idx, f) in enumerate(features)),
-    op,
-)
+function ExtensionHomomorphism{I}(features, op) where I<:Integer
+    @assert typemax(I) >= length(features) "Use wider Integer type!"
+    reversef = Dict{eltype(features), I}(f => idx for (idx, f) in pairs(features))
+    return ExtensionHomomorphism(features, reversef, op)
+end
+
+ExtensionHomomorphism(features, op) = ExtensionHomomorphism{Int16}(features, op)
 
 Base.getindex(ehom::ExtensionHomomorphism, i::Integer) = ehom.features[i]
 Base.getindex(ehom::ExtensionHomomorphism{T}, f::T) where {T} = ehom.reversef[f]
-(ehom::ExtensionHomomorphism)(p::Perm{I}) where {I} =
-    Perm(vec(I[ehom[ehom.op(f, p)] for f in ehom.features]))
+function (ehom::ExtensionHomomorphism)(g::GroupsCore.GroupElement)
+    return Perm(vec([ehom[ehom.op(f, g)] for f in ehom.features]))
+end
 
-function (ehom::ExtensionHomomorphism)(
-    orb::O,
-) where {T,O<:AbstractOrbit{T,Nothing}}
+function (ehom::ExtensionHomomorphism)(orb::O) where {T,O<:AbstractOrbit{T,Nothing}}
     elts = ehom.(orb)
     dict = Dict(e => nothing for e in elts)
-    return O(elts, dict)
+    return Orbit(elts, dict)
 end
 
-function (ehom::ExtensionHomomorphism)(χ::CF) where {CF<:ClassFunction}
+function (ehom::ExtensionHomomorphism)(χ::Character)
     iccG = ehom.(conjugacy_classes(χ))
-    return CF(values(χ), χ.inv_of, iccG)
+    return Character(values(χ), χ.inv_of, iccG)
 end

src/characters.jl

@@ -29,7 +29,7 @@ end
 Base.:(==)(χ::AbstractClassFunction, ψ::AbstractClassFunction) =
     conjugacy_classes(χ) === conjugacy_classes(ψ) && values(χ) == values(ψ)
 
-Base.hash(χ::AbstractClassFunction, h::UInt = UInt(0)) =
+Base.hash(χ::AbstractClassFunction, h::UInt) =
     hash(conjugacy_classes(χ), hash(values(χ), hash(AbstractClassFunction, h)))
 
 ####################################
@@ -49,6 +49,23 @@ end
 
 const ClassFunction = Union{Character,VirtualCharacter}
 
+"""
+    affordable_real!(χ::ClassFunction[, pmap::PowerMap])
+Return either `χ` or `2re(χ)` depending whether `χ` is afforded by a real representation, modifying `χ` in place.
+"""
+function affordable_real!(
+    χ::ClassFunction,
+    pmap = PowerMap(conjugacy_classes(χ)),
+)
+    ι = frobenius_schur_indicator(χ, pmap)
+    if ι <= 0 # i.e. χ is complex or quaternionic
+        for i in eachindex(χ.vals)
+            χ.vals[i] += conj(χ.vals[i])
+        end
+    end
+    return χ
+end
+
 """
     frobenius_schur_indicator(χ::AbstractClassFunction[, pmap::PowerMap])
 Return Frobenius-Schur indicator of `χ`, i.e. `Σχ(g²)` where sum is taken over
@@ -67,16 +84,21 @@ function frobenius_schur_indicator(
     χ::AbstractClassFunction,
     pmap::PowerMap = PowerMap(conjugacy_classes(χ)),
 )
-    res = sum(
+    ι = sum(
         length(c) * χ[pmap[i, 2]] for (i, c) in enumerate(conjugacy_classes(χ))
     )
-    return res / order(G)
+
+    ι_int = Int(ι)
+    ordG = sum(length, conjugacy_classes(χ))
+    d, r = divrem(ι_int, ordG)
+    @assert r == 0 "Non integral Frobenius Schur Indicator: $(ι_int) = $d * $ordG + $r"
+    return d
 end
 
 Base.isreal(χ::AbstractClassFunction) = frobenius_schur_indicator(χ) > 0
 
 if VERSION >= v"1.3.0"
-    function (χ::AbstractClassFunction)(g::PermutationGroups.GroupElem)
+    function (χ::AbstractClassFunction)(g::GroupsCore.GroupElement)
         for (i, cc) in enumerate(conjugacy_classes(χ))
             g ∈ cc && return χ[i]
         end
@@ -131,7 +153,7 @@ Base.@propagate_inbounds function Base.getindex(χ::ClassFunction, i::Integer)
 end
 
 if VERSION < v"1.3.0"
-    function (χ::Character)(g::PermutationGroups.GroupElem)
+    function (χ::Character)(g::GroupsCore.GroupElement)
         for (i, cc) in enumerate(conjugacy_classes(χ))
             g ∈ cc && return χ[i]
         end
@@ -142,7 +164,7 @@ if VERSION < v"1.3.0"
             ),
         )
     end
-    function (χ::VirtualCharacter)(g::PermutationGroups.GroupElem)
+    function (χ::VirtualCharacter)(g::GroupsCore.GroupElement)
         for (i, cc) in enumerate(conjugacy_classes(χ))
             g ∈ cc && return χ[i]
         end
@@ -199,38 +221,6 @@ function normalize!(χ::Character)
     return χ
 end
 
-"""
-    affordable_real!(χ::ClassFunction[, pmap::PowerMap])
-Return either `χ` or `2re(χ)` depending whether `χ` is afforded by a real representation, modifying `χ` in place.
-"""
-function affordable_real!(
-    χ::ClassFunction,
-    pmap = PowerMap(conjugacy_classes(χ)),
-)
-    ι = frobenius_schur_indicator(χ, pmap)
-    if !isone(ι) # χ is complex or quaternionic
-        for i in eachindex(χ.vals)
-            χ.vals[i] += conj(χ.vals[i])
-        end
-    end
-    return χ
-end
-
-function frobenius_schur_indicator(
-    χ::Character,
-    pmap::PowerMap = PowerMap(conjugacy_classes(χ)),
-)
-    ι = sum(
-        length(c) * χ[pmap[i, 2]] for (i, c) in enumerate(conjugacy_classes(χ))
-    )
-
-    ι_int = Int(ι)
-    ordG = sum(length, conjugacy_classes(χ))
-    d, r = divrem(ι_int, ordG)
-    @assert r == 0 "Non integral Frobenius Schur Indicator: $(ι_int) = $d * $ordG + $r"
-    return d
-end
-
 for C in (:Character, :VirtualCharacter)
     @eval begin
         function Base.show(io::IO, ::MIME"text/plain", χ::$C)

src/cmmatrix.jl

@@ -25,7 +25,7 @@ function Base.getindex(M::CMMatrix, s::Integer, t::Integer)
 
         for g in M.cc[r]
             for h in M.cc[s]
-                out = PermutationGroups.mul!(out, g, h)
+                out = GroupsCore.mul!(out, g, h)
                 for t = 1:size(M, 2)
                     if out == first(M.cc[t])
                         M.m[s, t] += 1
@@ -38,19 +38,19 @@ function Base.getindex(M::CMMatrix, s::Integer, t::Integer)
     return M.m[s, t]
 end
 
-conjugacy_classes(G::PermutationGroups.Group) = conjugacy_classes_orbit(G)
+conjugacy_classes(G::GroupsCore.Group) = conjugacy_classes_orbit(G)
 
-function conjugacy_classes_orbit(G::PermutationGroups.Group)
+function conjugacy_classes_orbit(G::GroupsCore.Group)
     id = one(G)
     S = gens(G)
     ordG = order(Int, G)
 
-    cclasses = [PermutationGroups.Orbit([id], Dict(id => nothing))]
+    cclasses = [Orbit([id], Dict(id => nothing))]
     elts_counted = 1
 
     for g in G
         any(ccl -> g ∈ ccl, cclasses) && continue
-        ccl_g = PermutationGroups.Orbit(S, g, ^)
+        ccl_g = Orbit(S, g, ^)
         elts_counted += length(ccl_g)
         push!(cclasses, ccl_g)
         elts_counted == ordG && break

src/dixon.jl

@@ -1,8 +1,8 @@
-Base.exponent(G::PermutationGroups.Group) = exponent(conjugacy_classes(G))
-Base.exponent(cclasses::AbstractVector) = lcm(order.(first.(cclasses)))
-dixon_prime(G::PermutationGroups.Group) = dixon_prime(order(G), exponent(G))
+Base.exponent(G::GroupsCore.Group) = exponent(conjugacy_classes(G))
+Base.exponent(cclasses::AbstractArray) = lcm(order.(first.(cclasses)))
+dixon_prime(G::GroupsCore.Group) = dixon_prime(order(G), exponent(G))
 
-function dixon_prime(cclasses::AbstractVector)
+function dixon_prime(cclasses::AbstractArray)
     ordG = sum(length, cclasses)
     m = exponent(cclasses)
     return dixon_prime(ordG, m)
@@ -32,9 +32,9 @@ function common_esd(Ns, F::Type{<:FiniteFields.GF})
     return esd
 end
 
-characters_dixon(G::PermutationGroups.Group) =
+characters_dixon(G::GroupsCore.Group) =
     characters_dixon(Rational{Int}, G)
-characters_dixon(::Type{R}, G::PermutationGroups.Group) where {R<:Real} =
+characters_dixon(::Type{R}, G::GroupsCore.Group) where {R<:Real} =
     characters_dixon(R, conjugacy_classes(G))
 
 function characters_dixon(::Type{R}, cclasses::AbstractVector) where {R}

src/eigenspacedecomposition.jl

@@ -1,27 +1,79 @@
+# Version over exact field:
 function row_echelon_form!(A::AbstractMatrix{T}) where {T<:Number}
-    c, d = size(A)
-    l = Int[]
+    pivots = Int[]
     pos = 0
-    i = 0
-    for i = 1:d
+    for i = 1:size(A, 2)
         j = findfirst(x -> !iszero(x), @view A[pos+1:end, i])
         j === nothing && continue
         j += pos
         pos += 1
-        push!(l, i)
+        push!(pivots, i)
 
         if (pos != j)
-            A[[pos, j], :] = A[[j, pos], :]
+            A[[pos, j], :] .= A[[j, pos], :]
         end
-        A[pos, :] ./= A[pos, i]
 
-        for j = 1:c
+        w = inv(A[pos, i])
+        A[pos, :] .*= w
+
+        for j = 1:size(A, 1)
             if j != pos
-                A[j, :] -= A[j, i] .* A[pos, :]
+                @. A[j, :] -= A[j, i] * A[pos, :]
             end
         end
     end
-    return A, l
+    return A, pivots
+end
+
+function row_echelon_form!(A::AbstractMatrix{C}) where {T<:AbstractFloat, C<:Cyclotomic{T}}
+    pivots = Int[]
+    pos = 0
+    for i = 1:size(A, 2)
+
+        j = findfirst(x -> abs(x) > sqrt(eps(T)), @view A[pos+1:end, i])
+        if j === nothing
+            # A[pos+1:end, :] .= Cyclotomics.zero!.(@view A[pos+1:end, :])
+            continue
+        end
+        j += pos
+        pos += 1
+        push!(pivots, i)
+
+        if (pos != j)
+            A[[pos, j], :] .= A[[j, pos], :]
+        end
+
+        @assert abs(A[pos, i]) >= sqrt(eps(T))
+        w = inv(A[pos, i])
+        # A[pos, :] .*= w
+
+        for idx in 1:size(A, 2)
+            A[pos, idx] = if abs(A[pos, idx]) <= sqrt(eps(T))
+                Cyclotomics.zero!(A[pos, idx])
+            elseif idx != i
+                A[pos, idx] * w
+            else
+                Cyclotomics.one!(A[pos, idx])
+            end
+        end
+
+        for j = 1:size(A, 1)
+            if j != pos
+                if abs(A[j, i]) < sqrt(eps(T))
+                    A[j, i] = Cyclotomics.zero!.(A[j, i])
+                else
+                    @. A[j, :] -= A[j, i] * A[pos, :]
+                end
+            end
+        end
+    end
+
+    for i in eachindex(A)
+        ndigs = floor(Int, -log10(eps(T)) - log10(max(size(A)...)) - 3)
+        A[i] = Cyclotomics.roundcoeffs!(A[i], digits=ndigs)
+    end
+
+    return A, pivots
 end
 
 function row_echelon_form(A::AbstractMatrix)
@@ -47,7 +99,7 @@ function right_nullspace(M::AbstractMatrix{T}) where {T}
 end
 
 function left_nullspace(M::AbstractMatrix)
-    return Matrix(transpose(right_nullspace(Matrix(transpose(M)))))
+    return transpose(right_nullspace(transpose(M)))
 end
 
 function left_eigen(M::AbstractMatrix{T}) where {T<:FiniteFields.GF}