Merge pull request #9 from kalmarek/enh/dixon

Enh/dixon

.github/workflows/CompatHelper.yml

@@ -0,0 +1,16 @@
+name: CompatHelper
+on:
+  schedule:
+    - cron: '00 00 * * *'
+  workflow_dispatch:
+jobs:
+  CompatHelper:
+    runs-on: ubuntu-latest
+    steps:
+      - name: Pkg.add("CompatHelper")
+        run: julia -e 'using Pkg; Pkg.add("CompatHelper")'
+      - name: CompatHelper.main()
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+          COMPATHELPER_PRIV: ${{ secrets.COMPATHELPER_PRIV }}  # optional
+        run: julia -e 'using CompatHelper; CompatHelper.main()'

.github/workflows/TagBot.yml

@@ -0,0 +1,11 @@
+name: TagBot
+on:
+  schedule:
+    - cron: 0 * * * *
+jobs:
+  TagBot:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: JuliaRegistries/TagBot@v1
+        with:
+          token: ${{ secrets.GITHUB_TOKEN }}

.github/workflows/ci.yml

@@ -24,6 +24,16 @@ jobs:
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
+      - uses: actions/cache@v1
+        env:
+          cache-name: cache-artifacts
+        with:
+          path: ~/.julia/artifacts
+          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
+          restore-keys: |
+            ${{ runner.os }}-test-${{ env.cache-name }}-
+            ${{ runner.os }}-test-
+            ${{ runner.os }}-
       - uses: julia-actions/julia-buildpkg@latest
       - uses: julia-actions/julia-runtest@latest
       - uses: julia-actions/julia-processcoverage@v1

Project.toml

@@ -1,12 +1,16 @@
 name = "SymbolicWedderburn"
 uuid = "858aa9a9-4c7c-4c62-b466-2421203962a2"
-authors = ["tweisser <tillmann.weisser@web.de>"]
+authors = ["Marek Kaluba <kalmar@amu.edu.pl>", "tweisser <tillmann.weisser@web.de>"]
 version = "0.1.0"
 
 [deps]
-AbstractAlgebra = "c3fe647b-3220-5bb0-a1ea-a7954cac585d"
+Cyclotomics = "da8f5974-afbb-4dc8-91d8-516d5257c83b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+PermutationGroups = "8bc5a954-2dfc-11e9-10e6-cd969bffa420"
 Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"
-Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[compat]
+Cyclotomics = "0.1.2"
+PermutationGroups = "0.1.2"
+Primes = "0.4, 0.5"

src/SymbolicWedderburn.jl

@@ -1,16 +1,17 @@
 module SymbolicWedderburn
 
 using LinearAlgebra
-import AbstractAlgebra
 using Primes
 
-using Pkg
-Pkg.add(PackageSpec(url="https://github.com/kalmarek/PermutationGroups.jl"))
 using PermutationGroups
+using Cyclotomics
 
 include("gf.jl")
 include("eigenspacedecomposition.jl")
 include("ccmatrix.jl")
+
+include("characters.jl")
+include("powermap.jl")
 include("dixon.jl")
 
 end # module

src/ccmatrix.jl

@@ -21,7 +21,7 @@ function Base.getindex(M::CCMatrix, s::Integer, t::Integer)
 
         for g in M.cc[r]
             for h in M.cc[s]
-                out = AbstractAlgebra.mul!(out, g, h)
+                out = PermutationGroups.mul!(out, g, h)
                 for t in 1:size(M, 2)
                     if out == first(M.cc[t])
                         M.m[s, t] += 1

src/characters.jl

@@ -6,72 +6,86 @@ function LinearAlgebra.dot(
     χ::AbstractClassFunction{T},
     ψ::AbstractClassFunction{T},
 ) where {T}
-    # TODO: @assert v.cc == w.cc
 
-    R = parent(χ[1]) # TODO make something better here
-    val = zero(R)
-
-    for (i, cc) in enumerate(classes(χ))
-        val += R(length(cc)) * χ[i] * ψ[-i]
-    end
-
-    orderG = R(sum(length, classes(χ)))
-    val *= inv(orderG)
+    val = sum(length(cc) * χ[i] * ψ[-i] for (i, cc) in enumerate(conjugacy_classes(χ)))
+    orderG = sum(length, conjugacy_classes(χ))
+    val = div(val,orderG)
     return val
 end
 
-function (χ::AbstractClassFunction)(g::GroupElem)
-    for (i, cc) in enumerate(classes(χ))
-        g ∈ cc && return χ[i]
-    end
-    throw(DomainError(g, "element does not belong to conjugacy classes of χ"))
-end
-
 ####################################
 # Characters
 
-mutable struct Character{T,CCl} <: AbstractClassFunction{T}
+mutable struct Character{T,CCl<:AbstractOrbit} <: AbstractClassFunction{T}
     vals::Vector{T}
     inv_of::Vector{Int}
     cc::Vector{CCl}
+end
 
-    function Character(
-        v::AbstractVector{T},
-        ccls::AbstractVector{CCl},
-        inv_of = _inv_of(ccls),
-    ) where {T,CCl}
-
-        χ = new{T,CCl}(v, inv_of, ccls)
-
-        # initial normalization
-        R = parent(first(v))
-        id = one(first(first(ccls)))
-        if !isone(χ(id))
-            χ.vals .*= inv(χ(id))
+if VERSION >= v"1.3.0"
+    function (χ::AbstractClassFunction)(g::PermutationGroups.AbstractPerm)
+        for (i, cc) in enumerate(conjugacy_classes(χ))
+            g ∈ cc && return χ[i]
         end
-
-        # computing the degree of χ:
-        deg_χ = sqrt(inv(dot(χ, χ)))
-        @debug χ.vals, deg_χ
-
-        # renormalizing χ
-        χ.vals .*= deg_χ
-
-        return χ
+        throw(DomainError(g, "element does not belong to conjugacy classes of χ"))
+    end
+else
+    function (χ::Character)(g::PermutationGroups.AbstractPerm)
+        for (i, cc) in enumerate(conjugacy_classes(χ))
+            g ∈ cc && return χ[i]
+        end
+        throw(DomainError(g, "element does not belong to conjugacy classes of χ"))
     end
 end
 
-function _inv_of(cc::AbstractVector)
+function Character(
+    v::AbstractVector{T},
+    ccls::AbstractVector{CCl},
+) where {T,CCl<:AbstractOrbit}
+    χ = Character{T,CCl}(v, _inv_of(ccls), ccls)
+    return χ
+end
+
+function _inv_of(cc::AbstractVector{<:AbstractOrbit})
     inv_of = zeros(Int, size(cc))
     for (i, c) in enumerate(cc)
         g = inv(first(c))
         inv_of[i] = something(findfirst(k -> g in k, cc), 0)
     end
-    @assert !any(iszero, inv_of) "Could not find the conjugacy class of $g."
+    any(iszero, inv_of) &&
+        throw(ArgumentError("Could not find the conjugacy class for inverse of $(first(cc[findfirst(iszero, inv_of)]))."))
     return inv_of
 end
 
-classes(χ::Character) = χ.cc
+function normalize!(χ::Character)
+
+    ccG = conjugacy_classes(χ)
+    id = one(first(first(ccG)))
+
+    k = χ(id)
+    if !isone(k)
+        k⁻¹ = inv(k)
+        for i in eachindex(χ.vals)
+            χ.vals *= k⁻¹
+        end
+    end
+
+    # ⟨χ, χ⟩ = 1/d²
+
+    deg = sqrt(inv(dot(χ,χ)))
+    # @debug "normalizing with" n dot(χ, χ) χ(id) χ
+
+    # normalizing χ
+    for i in eachindex(χ.vals)
+        χ.vals[i] *= deg
+    end
+    return χ
+end
+
+PermutationGroups.conjugacy_classes(χ::Character) = χ.cc
+
+PermutationGroups.degree(χ::Character) =
+    Int(χ(one(first(first(conjugacy_classes(χ))))))
 
 Base.@propagate_inbounds function Base.getindex(χ::Character, i::Integer)
     @boundscheck checkbounds(χ.vals, abs(i))
@@ -81,3 +95,29 @@ Base.@propagate_inbounds function Base.getindex(χ::Character, i::Integer)
         return @inbounds χ.vals[i]
     end
 end
+
+function Base.show(io::IO, ::MIME"text/plain", χ::Character{T}) where {T}
+    println(io, "Character over $T")
+    cc_reps = string.(first.(χ.cc))
+    k = maximum(length.(cc_reps))
+
+    for (c, v) in zip(cc_reps, χ.vals)
+        println(io, rpad("$c^G", k + 3), "→ \t", v)
+    end
+end
+
+function Base.show(io::IO, χ::Character{T}) where {T}
+    v = χ.vals
+    io = IOContext(io, :typeinfo => eltype(v))
+    limited = get(io, :limit, false)
+    opn, cls = '[', ']'
+
+    print(io, "Character: ")
+    if limited && length(v) > 20
+        Base.show_delim_array(io, v, opn, ",", "", false, 1, f+9)
+        print(io, "  …  ")
+        Base.show_delim_array(io, v, "", ",", cls, false, l-9, l)
+    else
+        Base.show_delim_array(io, v, opn, ",", cls, false)
+    end
+end

src/dixon.jl

@@ -1,21 +1,107 @@
-AbstractAlgebra.exponent(G::AbstractAlgebra.AbstractPermutationGroup) =
-    AbstractAlgebra.exponent(conjugacy_classes(G))
-AbstractAlgebra.exponent(cclasses::AbstractVector) =
-    lcm(AbstractAlgebra.order.(first.(cclasses)))
-dixon_prime(G::AbstractAlgebra.Group) =
-    dixon_prime(AbstractAlgebra.order(G), exponent(G))
+Base.exponent(G::AbstractPermutationGroup) = exponent(conjugacy_classes(G))
+Base.exponent(cclasses::AbstractVector) = lcm(order.(first.(cclasses)))
+dixon_prime(G::AbstractPermutationGroup) = dixon_prime(order(G), exponent(G))
 
 function dixon_prime(cclasses::AbstractVector)
     ordG = sum(length, cclasses)
-    m = AbstractAlgebra.exponent(cclasses)
+    m = exponent(cclasses)
     return dixon_prime(ordG, m)
 end
 
 function dixon_prime(ordG::Integer, exponent::Integer)
-    p = 2floor(Int, sqrt(ordG))
+    p = 2 * floor(Int, sqrt(ordG))
     while true
-        p = nextprime(p+1)
+        p = nextprime(p + 1)
         isone(p % exponent) && break # we need -1 to be in the field
     end
     return p
 end
+
+function common_esd(Ns, F::Type{<:FiniteFields.GF})
+    @assert !isempty(Ns)
+    esd = EigenSpaceDecomposition(F.(first(Ns)))
+    for N in Iterators.rest(Ns, 2)
+        esd = refine(esd, F.(N))
+        @debug N esd.eigspace_ptrs
+        isdiag(esd) && return esd
+    end
+    return esd
+end
+
+characters_dixon(G::AbstractPermutationGroup) =
+    characters_dixon(conjugacy_classes(G))
+
+function characters_dixon(cclasses::AbstractVector{<:AbstractOrbit})
+    p = dixon_prime(cclasses)
+    chars_𝔽p = characters_dixon(cclasses, FiniteFields.GF{p})
+    return complex_characters(chars_𝔽p)
+end
+
+function characters_dixon(
+    cclasses::AbstractVector{<:AbstractOrbit},
+    F::Type{<:FiniteFields.GF},
+)
+    Ns = [CCMatrix(cclasses, i) for i = 1:length(cclasses)]
+    esd = common_esd(Ns, F)
+    @assert isdiag(esd) "Class Matricies failed to diagonalize! $esd"
+    inv_ccls = _inv_of(cclasses)
+    return [
+        normalize!(Character(vec(eigensubspace), inv_ccls, cclasses))
+        for eigensubspace in esd
+    ]
+end
+
+function _multiplicities(
+    chars::AbstractVector{<:Character{F}},
+    cclasses = conjugacy_classes(first(chars)),
+) where {F<:FiniteFields.GF}
+
+    e = Int(exponent(cclasses))
+    ie = inv(F(e))
+    ω = FiniteFields.rootofunity(F, e)
+
+
+    multiplicities = zeros(Int, length(chars), length(cclasses), e)
+    powermap = PowerMap(cclasses)
+    for (i, χ) in enumerate(chars)
+        for j = 1:length(cclasses), k = 0:e-1
+            multiplicities[i, j, k+1] =
+                Int(ie * sum(χ[powermap[j, l]] * ω^-(k * l) for l = 0:e-1))
+        end
+    end
+
+    return multiplicities
+end
+
+
+function complex_characters(
+    chars::AbstractVector{<:Character{F}},
+) where {F<:FiniteFields.GF}
+
+    cclasses = conjugacy_classes(first(chars))
+    lccl = length(cclasses)
+    mult_c = _multiplicities(chars, cclasses)
+    e = size(mult_c, 3) # the exponent
+
+    inv_of_cls = first(chars).inv_of
+
+    C = Cyclotomics.Cyclotomic{Int,Cyclotomics.SparseVector{Int,Int}}
+    # C = typeof(Cyclotomics.E(5))
+
+    complex_chars = Vector{Character{C,eltype(cclasses)}}(undef, length(chars))
+
+
+    for i = 1:length(complex_chars)
+        complex_chars[i] = Character(
+            [
+                Cyclotomics.reduced_embedding(sum(
+                    mult_c[i, j, k+1] * E(e, k) for k = 0:e-1
+                )) for j = 1:lccl
+            ],
+            inv_of_cls,
+            cclasses,
+        )
+    end
+
+    return complex_chars
+end