Allow multiplication of characters

examples/ex_S4.jl

@@ -10,7 +10,7 @@ include(joinpath(@__DIR__, "solver.jl"))
 const N = 4
 @polyvar x[1:N]
 
-OPTIMIZER = csdp_optimizer(; max_iters = 2_000, accel = 10, eps = 1e-7)
+OPTIMIZER = scs_optimizer(; max_iters = 2_000, accel = 10, eps = 1e-7)
 
 f =
     1 +

src/Characters/Characters.jl

@@ -11,7 +11,11 @@ import PermutationGroups
 export AbstractClassFunction, Character, CharacterTable
 
 export conjugacy_classes,
-    constituents, degree, irreducible_characters, isirreducible, table
+    multiplicities,
+    degree,
+    irreducible_characters,
+    isirreducible,
+    table
 
 include("gf.jl")
 

src/Characters/character_tables.jl

@@ -24,11 +24,13 @@
 
 ## irreps
 function irreducible_characters(chtbl::CharacterTable)
-    return [Character(chtbl, i) for i in 1:size(chtbl, 1)]
+    return irreducible_characters(eltype(chtbl), chtbl)
 end
+
 function irreducible_characters(T::Type, chtbl::CharacterTable)
-    return Character{T}[Character{T}(chtbl, i) for i in 1:size(chtbl, 1)]
+    return [Character{T}(chtbl, i) for i in axes(chtbl, 1)]
 end
+
 function irreducible_characters(G::Group, cclasses = conjugacy_classes(G))
     return irreducible_characters(Rational{Int}, G, cclasses)
 end
@@ -41,23 +43,23 @@
     return irreducible_characters(CharacterTable(R, G, cclasses))
 end
 
-function trivial_character(chtbl::CharacterTable)
-    # return Character(chtbl, findfirst(r->all(isone, r), eachrow(chtbl)))
-    # can't use findfirst(f, eachrow(...)) on julia-1.6
-    for i in 1:size(chtbl, 1)
-        all(isone, @view(chtbl[i, :])) && return Character(chtbl, i)
-    end
-    # never hit, to keep compiler happy
-    return Character(chtbl, 0)
-end
+trivial_character(chtbl::CharacterTable) = Character(chtbl, 1)
 
 ## construcing tables
 
+function CharacterTable(G::Group, cclasses = conjugacy_classes(G))
+    return CharacterTable(Rational{Int}, G, cclasses)
+end
+
 function CharacterTable(
     Fp::Type{<:FiniteFields.GF},
     G::Group,
     cclasses = conjugacy_classes(G),
 )
+    # make sure that the first class contains the indentity
+    k = findfirst(cl -> one(G) in cl, cclasses)
+    cclasses[k], cclasses[1] = cclasses[1], cclasses[k]
+
     Ns = [CMMatrix(cclasses, i) for i in 1:length(cclasses)]
     esd = common_esd(Ns, Fp)
     @assert isdiag(esd)
@@ -71,6 +73,19 @@
     )
 
     tbl = normalize!(tbl)
+
+    let vals = tbl.values
+        # make order of characters deterministic
+        vals .= sortslices(vals; dims = 1)
+        # and that the trivial character is first
+        k = findfirst(i -> all(isone, @views vals[i, :]), axes(vals, 1))
+        # doesn't work on julia-1.6
+        # k = findfirst(r -> all(isone, r), eachrow(vals))
+        if k ≠ 1
+            _swap_rows!(vals, 1, k)
+        end
+    end
+
     return tbl
 end
 
@@ -131,7 +146,8 @@
 
 function normalize!(chtbl::CharacterTable{<:Group,<:FiniteFields.GF})
     id = one(parent(chtbl))
-    for (i, χ) in enumerate(irreducible_characters(chtbl))
+    Threads.@threads for i in axes(chtbl, 1)
+        χ = Character(chtbl, i)
         k = χ(id)
         if !isone(k)
             chtbl.values[i, :] .*= inv(k)

src/Characters/class_functions.jl

@@ -59,29 +59,30 @@
 Struct representing (possibly virtual) character of a group.
 
 Characters are backed by `table(χ)::CharacterTable` which actually stores the
-character values. The constituents (decomposition into the irreducible summands)
-of a given character can be obtained by calling `constituents(χ)` which returns a
-vector of coefficients of `χ` in the basis of `irreducible_characters(table(χ))`.
+character values. The multiplicities (decomposition into the irreducible
+summands) of a given character can be obtained by calling `multiplicities(χ)`
+which returns a vector of coefficients of `χ` in the basis of
+`irreducible_characters(table(χ))`.
 
-It is assumed that equal class functions on the same group will have **identical**
-(ie. `===`) character tables.
+It is assumed that equal class functions on the same group will have
+**identical** (ie. `===`) character tables.
 """
 struct Character{T,S,ChT<:CharacterTable} <: AbstractClassFunction{T}
     table::ChT
-    constituents::Vector{S}
+    multips::Vector{S}
 end
 
 function Character{R}(
     chtbl::CharacterTable{Gr,T},
-    constituents::AbstractVector{S},
+    multips::AbstractVector{S},
 ) where {R,Gr,T,S}
-    return Character{R,S,typeof(chtbl)}(chtbl, constituents)
+    return Character{R,S,typeof(chtbl)}(chtbl, multips)
 end
 
-function Character(chtbl::CharacterTable, constituents::AbstractVector)
-    R = Base._return_type(*, Tuple{eltype(chtbl),eltype(constituents)})
+function Character(chtbl::CharacterTable, multips::AbstractVector)
+    R = Base._return_type(*, Tuple{eltype(chtbl),eltype(multips)})
     @assert R ≠ Any
-    return Character{R}(chtbl, constituents)
+    return Character{R}(chtbl, multips)
 end
 
 function Character(chtbl::CharacterTable, i::Integer)
@@ -95,14 +96,14 @@
 end
 
 function Character{T}(χ::Character) where {T}
-    S = eltype(constituents(χ))
+    S = eltype(multiplicities(χ))
     ChT = typeof(table(χ))
-    return Character{T,S,ChT}(table(χ), constituents(χ))
+    return Character{T,S,ChT}(table(χ), multiplicities(χ))
 end
 
 ## Accessors
 table(χ::Character) = χ.table
-constituents(χ::Character) = χ.constituents
+multiplicities(χ::Character) = χ.multips
 
 ## AbstractClassFunction api
 Base.parent(χ::Character) = parent(table(χ))
@@ -122,7 +123,7 @@
         T,
         sum(
             c * table(χ)[idx, i] for
-            (idx, c) in enumerate(constituents(χ)) if !iszero(c);
+            (idx, c) in enumerate(multiplicities(χ)) if !iszero(c);
             init = zero(T),
         ),
     )
@@ -131,12 +132,13 @@
 ## Basic functionality
 
 function Base.:(==)(χ::Character, ψ::Character)
-    return table(χ) === table(ψ) && constituents(χ) == constituents(ψ)
+    return table(χ) === table(ψ) && multiplicities(χ) == multiplicities(ψ)
 end
-Base.hash(χ::Character, h::UInt) = hash(table(χ), hash(constituents(χ), h))
+Base.hash(χ::Character, h::UInt) = hash(table(χ), hash(multiplicities(χ), h))
 
-function Base.deepcopy_internal(χ::Character, ::IdDict)
-    return Character(table(χ), copy(constituents(χ)))
+function Base.deepcopy_internal(χ::Character{T}, d::IdDict) where {T}
+    haskey(d, χ) && return d[χ]
+    return Character{T}(table(χ), copy(multiplicities(χ)))
 end
 
 ## Character arithmetic
@@ -145,17 +147,54 @@
     @eval begin
         function Base.$f(χ::Character, ψ::Character)
             @assert table(χ) === table(ψ)
-            return Character(table(χ), $f(constituents(χ), constituents(ψ)))
+            return Character(table(χ), $f(multiplicities(χ), multiplicities(ψ)))
         end
     end
 end
 
-Base.:*(χ::Character, c::Number) = Character(table(χ), c .* constituents(χ))
+Base.:*(χ::Character, c::Number) = Character(table(χ), c .* multiplicities(χ))
 Base.:*(c::Number, χ::Character) = χ * c
-Base.:/(χ::Character, c::Number) = Character(table(χ), constituents(χ) ./ c)
+Base.:/(χ::Character, c::Number) = Character(table(χ), multiplicities(χ) ./ c)
 
 Base.zero(χ::Character) = 0 * χ
 
+function __decompose(T::Type, values::AbstractVector, tbl::CharacterTable)
+    ψ = ClassFunction(values, conjugacy_classes(tbl), tbl.inv_of)
+    return decompose(T, ψ, tbl)
+end
+
+function decompose(cfun::AbstractClassFunction, tbl::CharacterTable)
+    return decompose(eltype(cfun), cfun, tbl)
+end
+
+function decompose(T::Type, cfun::AbstractClassFunction, tbl::CharacterTable)
+    vals = Vector{T}(undef, length(conjugacy_classes(tbl)))
+    return decompose!(vals, cfun, tbl)
+end
+
+function decompose!(
+    vals::AbstractVector,
+    cfun::AbstractClassFunction,
+    tbl::CharacterTable,
+)
+    @assert length(vals) == length(conjugacy_classes(tbl))
+    @assert conjugacy_classes(tbl) === conjugacy_classes(cfun)
+
+    for (i, idx) in enumerate(eachindex(vals))
+        χ = Character(tbl, i)
+        vals[idx] = dot(χ, cfun)
+    end
+    return vals
+end
+
+function Base.:*(χ::Character, ψ::Character)
+    @assert table(χ) == table(ψ)
+    values = Characters.values(χ) .* Characters.values(ψ)
+    return Character(table(χ), __decompose(Int, values, table(χ)))
+end
+
+Base.:^(χ::Character, n::Integer) = Base.power_by_squaring(χ, n)
+
 ## Group-theoretic functions:
 
 PermutationGroups.degree(χ::Character) = Int(χ(one(parent(χ))))
@@ -170,16 +209,16 @@
     all(isreal, vals) && return Character{T}(χ)
     tbl = table(χ)
     ψ = ClassFunction(vals[tbl.inv_of], conjugacy_classes(tbl), tbl.inv_of)
-    constituents = S[dot(ψ, χ) for χ in irreducible_characters(tbl)]
-    return Character{T,eltype(constituents),typeof(tbl)}(tbl, constituents)
+    multips = S[dot(ψ, χ) for χ in irreducible_characters(tbl)]
+    return Character{T,eltype(multips),typeof(tbl)}(tbl, multips)
 end
 
 function isvirtual(χ::Character)
-    return any(<(0), constituents(χ)) || any(!isinteger, constituents(χ))
+    return any(<(0), multiplicities(χ)) || any(!isinteger, multiplicities(χ))
 end
 
 function isirreducible(χ::Character)
-    C = constituents(χ)
+    C = multiplicities(χ)
     k = findfirst(!iszero, C)
     k !== nothing || return false # χ is zero
     isone(C[k]) || return false # muliplicity is ≠ 1
@@ -196,7 +235,7 @@
 function affordable_real!(χ::Character)
     ι = frobenius_schur(χ)
     if ι <= 0 # i.e. χ is complex or quaternionic
-        χ.constituents .+= constituents(conj(χ))
+        χ.multips .+= multiplicities(conj(χ))
     end
     return χ
 end
@@ -242,7 +281,7 @@
 
 function _print_char(io::IO, χ::Character)
     first = true
-    for (i, c) in enumerate(constituents(χ))
+    for (i, c) in enumerate(multiplicities(χ))
         iszero(c) && continue
         first || print(io, " ")
         print(io, ((c < 0 || first) ? "" : '+'))

src/action_characters.jl

@@ -63,8 +63,8 @@ function action_character(
     tbl::CharacterTable,
 ) where {T}
     ac_char = _action_class_fun(conjugacy_clss)
-    constituents = Int[dot(ac_char, χ) for χ in irreducible_characters(tbl)]
-    return Character{T}(tbl, constituents)
+    vals = Characters.decompose(Int, ac_char, tbl)
+    return Character{T}(tbl, vals)
 end
 
 function action_character(
@@ -73,14 +73,8 @@ function action_character(
     tbl::CharacterTable,
 ) where {T<:Union{AbstractFloat,ComplexF64}}
     ac_char = _action_class_fun(conjugacy_cls)
-
-    constituents = [dot(ac_char, χ) for χ in irreducible_characters(tbl)]
-    all(constituents) do c
-        ac = abs(c)
-        return abs(ac - round(ac)) < eps(real(T)) * length(conjugacy_cls)
-    end
-
-    return Character{T}(tbl, round.(Int, abs.(constituents)))
+    vals = Characters.decompose(T, ac_char, tbl)
+    return Character{T}(tbl, round.(Int, abs.(vals)))
 end
 
 """
@@ -100,8 +94,7 @@ function action_character(
     hom::InducedActionHomomorphism,
     tbl::CharacterTable,
 ) where {T}
-    ac_char = _action_class_fun(hom, conjugacy_classes(tbl))
-    vals = [dot(ac_char, χ) for χ in irreducible_characters(T, tbl)]
-    constituents = convert(Vector{Int}, vals)
-    return Character{T}(tbl, constituents)
+    act_character = _action_class_fun(hom, conjugacy_classes(tbl))
+    vals = Characters.decompose(Int, act_character, tbl)
+    return Character{T}(tbl, vals)
 end

src/matrix_projections.jl

@@ -110,7 +110,7 @@ _mproj_fitsT!(args...) = _mproj_outsT!(args...)
 function _mproj_outsT!(mproj::AbstractMatrix{T}, χ::Character) where {T}
     mproj .= sum(
         c .* matrix_projection_irr(Character(table(χ), i)) for
-        (i, c) in pairs(constituents(χ)) if !iszero(c)
+        (i, c) in pairs(multiplicities(χ)) if !iszero(c)
     )
     return mproj
 end
@@ -122,7 +122,7 @@ function _mproj_outsT!(
 ) where {T}
     mproj .= sum(
         c .* matrix_projection_irr(hom, Character(table(χ), i)) for
-        (i, c) in pairs(constituents(χ)) if !iszero(c)
+        (i, c) in pairs(multiplicities(χ)) if !iszero(c)
     )
     return mproj
 end

src/sa_basis.jl

@@ -198,18 +198,18 @@ end
 function _constituents_decomposition(ψ::Character, tbl::CharacterTable)
     irr = irreducible_characters(tbl)
     degrees = degree.(irr)
-    multiplicities = constituents(ψ)
+    multips = multiplicities(ψ)
 
     @debug "Decomposition into character spaces:
     degrees:        $(join([lpad(d, 6) for d in degrees], ""))
-    multiplicities: $(join([lpad(m, 6) for m in multiplicities], ""))"
+    multiplicities: $(join([lpad(m, 6) for m in multips], ""))"
 
-    @assert dot(multiplicities, degrees) == degree(ψ)
+    @assert dot(multips, degrees) == degree(ψ)
     "Something went wrong: characters do not constitute a complete basis for action:
-    $(dot(multiplicities, degrees)) ≠ $(degree(ψ))"
+    $(dot(multips, degrees)) ≠ $(degree(ψ))"
 
-    present_irreps = [i for (i, m) in enumerate(multiplicities) if m ≠ 0]
-    return irr[present_irreps], multiplicities[present_irreps]
+    present_irreps = [i for (i, m) in enumerate(multips) if m ≠ 0]
+    return irr[present_irreps], multips[present_irreps]
 end
 
 function _symmetry_adapted_basis(
@@ -236,14 +236,14 @@ end
 function _symmetry_adapted_basis(
     T::Type,
     irr::AbstractVector{<:Character},
-    multiplicities::AbstractVector{<:Integer},
+    multips::AbstractVector{<:Integer},
     RG::StarAlgebra{<:Group},
     hom = nothing,
 )
     mps, ranks = minimal_projection_system(irr, RG)
     degrees = degree.(irr)
     @debug "ranks of projections obtained by mps:" degrees
-    res = map(zip(mps, multiplicities, degrees, ranks)) do (µ, m, deg, r)
+    res = map(zip(mps, multips, degrees, ranks)) do (µ, m, deg, r)
         Threads.@spawn begin
             µT = eltype(µ) == T ? µ : AlgebraElement{T}(µ)
             # here we use algebra to compute the dimension of image;
@@ -259,7 +259,7 @@ function _symmetry_adapted_basis(
     for (χ, ds) in zip(irr, direct_summands)
         if issimple(ds) &&
            (d = size(ds, 1)) !=
-           (e = multiplicity(ds) * sum(constituents(χ) .> 0))
+           (e = multiplicity(ds) * sum(multiplicities(χ) .> 0))
             throw(
                 "The dimension of the projection doesn't match with simple summand multiplicity: $d ≠ $e",
             )