rewrite the definitions of groups in symmetry.jl (#32)

* rewrite the definitions of groups in symmetry.jl

* add Random

Project.toml

@@ -10,6 +10,7 @@ GroupsCore = "d5909c97-4eac-4ecc-a3dc-fdd0858a4120"
 MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
 MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
 PermutationGroups = "8bc5a954-2dfc-11e9-10e6-cd969bffa420"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SumOfSquares = "4b9e565b-77fc-50a5-a571-1244f986bda1"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

src/symmetry.jl

@@ -1,13 +1,17 @@
 using SumOfSquares
+import Random
 using GroupsCore
 
 using PermutationGroups
-symmetric_group(n) = PermutationGroups.PermGroup(PermutationGroups.Perm([2,1]), PermutationGroups.Perm([2:n;1]))
-const G3 = symmetric_group(3)
-const P3 = collect(G3)
+function __symmetric_grp_gens(n::Integer)
+    return (Perm(UInt16[2, 1], false), Perm(UInt16[2:n; 1], false))
+end
+symmetric_group(n::Integer) = PermGroup(__symmetric_grp_gens(n)...)
 
 export Lattice1Group
 
+## Klein C₂⊕C₂ group
+struct KleinGroup <: Group end
 """
     struct KleinElement <: GroupElement
         id::Int
@@ -22,9 +26,17 @@ form a the Klein group.
 struct KleinElement <: GroupElement
     id::Int
 end
-Base.:(==)(a::KleinElement, b::KleinElement) = a.id == b.id
 
-PermutationGroups.order(el::KleinElement) = iszero(el.id) ? 1 : 2
+GroupsCore.gens(::KleinGroup) = [KleinElement(1), KleinElement(2)]
+GroupsCore.order(::Type{T}, G::KleinGroup) where {T} = convert(T, 4)
+Base.IteratorSize(::Type{KleinGroup}) = Base.HasLength()
+Base.eltype(::Type{KleinGroup}) = KleinElement
+Base.iterate(::KleinGroup, st=0) = ifelse(st ≥ 4, nothing, (KleinElement(st), st + 1))
+Base.one(::KleinGroup) = KleinElement(0)
+
+Base.parent(::KleinElement) = KleinGroup()
+GroupsCore.order(::Type{T}, el::KleinElement) where {T} = convert(T, ifelse(iszero(el.id), 1, 2))
+Base.:(==)(a::KleinElement, b::KleinElement) = a.id == b.id
 Base.inv(el::KleinElement) = el
 
 function Base.:*(a::KleinElement, b::KleinElement)
@@ -36,40 +48,52 @@ function Base.:*(a::KleinElement, b::KleinElement)
     return KleinElement(2)
 end
 
-Base.conj(a::KleinElement, b::KleinElement) = inv(b) * a * b
-Base.:^(a::KleinElement, b::KleinElement) = conj(a, b)
-
-struct KleinGroup <: Group end
-Base.one(::Union{KleinGroup, KleinElement}) = KleinElement(0)
-PermutationGroups.gens(::KleinGroup) = [KleinElement(1), KleinElement(2)]
-PermutationGroups.order(::Type{T}, G::KleinGroup) where {T} = convert(T, 4)
-function Base.iterate(::KleinGroup, prev::KleinElement = KleinElement(-1))
-    id = prev.id + 1
-    if id > 4
-        return nothing
-    else
-        next = KleinElement(id)
-        return next, next
-    end
+function Base.rand(rng::Random.AbstractRNG, st::Random.SamplerTrivial{KleinGroup})
+    klein = st[]
+    id = rand(rng, 1:order(Int, klein))
+    return KleinElement(id - 1)
 end
 
-#SymbolicWedderburn.conjugacy_classes_orbit(KleinGroup())
-
-function perm_klein(k::KleinElement, p::PermutationGroups.Permutation)
-    if k.id == 0
-        return k
-    else
-        return KleinElement(k.id^p)
-    end
-end
+## PermKleinGroup, aka the semi-direct product of Klein and S₃ over ϕ,
+# the homomorphism ϕ: S₃ → Aut(Klein) given by the permutation
+# of the three non-trivial elements in Klein.
+struct KleinPermGroup <: Group end
 
 """
 Group element `k * p = p * k^inv(p)`.
 """
 struct KleinPermElement <: GroupElement
-    p::eltype(P3)
+    p::PermutationGroups.Perm{UInt16}
     k::KleinElement
 end
+
+function GroupsCore.gens(::KleinPermGroup)
+    Kid = one(KleinGroup())
+    return [KleinPermElement(g, Kid) for g in __symmetric_grp_gens(3)]
+end
+GroupsCore.order(::Type{T}, G::KleinPermGroup) where {T} = convert(T, 6 * 4)
+Base.IteratorSize(::Type{KleinPermGroup}) = Base.HasShape{2}()
+Base.size(::KleinPermGroup) = (6, 4)
+Base.eltype(::Type{KleinPermGroup}) = KleinPermElement
+
+function Base.iterate(::KleinPermGroup)
+    P3 = map(PermutationGroups.Perms.perm, symmetric_group(3))
+    IT = Iterators.product(P3, KleinGroup())
+    (p, k), st = iterate(IT)
+    return KleinPermElement(p, k), (IT, st)
+end
+
+function Base.iterate(::KleinPermGroup, it_st)
+    IT, st = it_st
+    el_st = iterate(IT, st)
+    isnothing(el_st) && return nothing
+    (p, k), st = el_st
+    return KleinPermElement(p, k), (IT, st)
+end
+
+Base.one(::KleinPermGroup) = KleinPermElement(one(Perm{UInt16}), one(KleinGroup()))
+
+Base.parent(::KleinPermElement) = KleinPermGroup()
 Base.isone(el::KleinPermElement) = isone(el.k) && isone(el.p)
 
 function Base.hash(el::KleinPermElement, u::UInt64)
@@ -79,154 +103,157 @@ function Base.:(==)(a::KleinPermElement, b::KleinPermElement)
     return a.p == b.p && a.k == b.k
 end
 
+function __klein_ϕ(k::KleinElement, p::PermutationGroups.AbstractPermutation)
+    return ifelse(k.id == 0, k, KleinElement(k.id^p))
+end
+
 # k^(inv(p)) * inv(p) * k * p = k^(inv(p)) * k * inv(p) * p = 1
 function Base.inv(el::KleinPermElement)
     inv_p = inv(el.p)
-    KleinPermElement(inv_p, perm_klein(el.k, inv_p))
+    KleinPermElement(inv_p, __klein_ϕ(el.k, inv_p))
 end
 
 # p * k * q * k' = p * q * k^p * k'
 # k * p * k' * q = k * k'^p * p * q
 function Base.:*(a::KleinPermElement, b::KleinPermElement)
-    return KleinPermElement(a.p * b.p, perm_klein(a.k, b.p) * b.k)
+    return KleinPermElement(a.p * b.p, __klein_ϕ(a.k, b.p) * b.k)
 end
-function Base.:^(el::KleinPermElement, k::Integer)
-    return Base.power_by_squaring(el, k)
-end
-
-Base.conj(a::KleinPermElement, b::KleinPermElement) = inv(b) * a * b
-Base.:^(a::KleinPermElement, b::KleinPermElement) = conj(a, b)
 
-function PermutationGroups.order(el::KleinPermElement)
-    cur = el
-    i = 1
-    while !isone(cur)
-        i += 1
-        cur *= el
-    end
-    return i
+function Base.rand(rng::Random.AbstractRNG, ::Random.SamplerTrivial{KleinPermGroup})
+    p = Perm{UInt16}(Random.randperm(rng, 3), false)
+    k = rand(rng, KleinGroup())
+    return KleinPermElement(p, k)
 end
 
-struct KleinPermGroup <: Group end
-Base.one(::Union{KleinPermGroup, KleinPermElement}) = KleinPermElement(one(G3), one(KleinGroup()))
-# See https://github.com/blegat/CondensedMatterSOS.jl/pull/31#issuecomment-1717665659
-function PermutationGroups.gens(g::KleinPermGroup)
-    els = [one(g)]
-    for g in gens(G3)
-        push!(els, KleinPermElement(g, one(KleinGroup())))
-    end
-    return els
-end
-PermutationGroups.order(::Type{T}, G::KleinPermGroup) where {T} = convert(T, 6 * 4)
-const IT = Iterators.product(1:6, 0:3)
-function Base.iterate(::KleinPermGroup, args...)
-    el_st = iterate(IT, args...)
-    if el_st === nothing
-        return nothing
+function Base.deepcopy_internal(g::KleinPermElement, dict::IdDict)
+    p = if haskey(dict, objectid(g.p))
+        dict[objectid(g.p)]
     else
-        el, st = el_st
-        return KleinPermElement(P3[el[1]], KleinElement(el[2])), st
+        Base.deepcopy_internal(g.p, dict)
     end
+    return KleinPermElement(p, g.k)
 end
 
+# CyclicGroup
+struct CyclicGroup <: Group
+    n::Int
+end
 struct CyclicElem <: GroupElement
     n::Int
     id::Int
 end
+
+GroupsCore.gens(c::CyclicGroup) = [CyclicElem(c.n, 1)]
+GroupsCore.order(::Type{T}, c::CyclicGroup) where {T} = convert(T, c.n)
+Base.IteratorSize(::Type{CyclicGroup}) = Base.HasLength()
+Base.eltype(::Type{CyclicGroup}) = CyclicElem
+Base.iterate(c::CyclicGroup, st=0) = ifelse(st ≥ c.n, nothing, (CyclicElem(c.n, st), st + 1))
+Base.one(C::CyclicGroup) = CyclicElem(C.n, 0)
+
+Base.parent(c::CyclicElem) = CyclicGroup(c.n)
+GroupsCore.order(el::CyclicElem) = div(el.n, gcd(el.n, el.id))
 Base.:(==)(a::CyclicElem, b::CyclicElem) = a.n == b.n && a.id == b.id
 Base.inv(el::CyclicElem) = CyclicElem(el.n, (el.n - el.id) % el.n)
-
-function Base.:*(a::CyclicElem, b::CyclicElem)
-    return CyclicElem(a.n, (a.id + b.id) % a.n)
-end
+Base.:*(a::CyclicElem, b::CyclicElem) = CyclicElem(a.n, (a.id + b.id) % a.n)
 Base.:^(el::CyclicElem, k::Integer) = CyclicElem(el.n, (el.id * k) % el.n)
 
-Base.conj(a::CyclicElem, b::CyclicElem) = inv(b) * a * b
-Base.:^(a::CyclicElem, b::CyclicElem) = conj(a, b)
-
-function PermutationGroups.order(el::CyclicElem)
-    return div(el.n, gcd(el.n, el.id))
+function Base.rand(rng::Random.AbstractRNG, st::Random.SamplerTrivial{CyclicGroup})
+    C = st[]
+    id = rand(rng, 1:order(Int, C))
+    return CyclicElem(C.n, id - 1)
 end
 
-struct CyclicGroup <: Group
-    n::Int
+# General direct sum of two groups
+struct DirectSumGroup{G1<:Group,G2<:Group} <: Group
+    H::G1
+    K::G2
+end
+struct DirectSumElem{G1El<:GroupElement,G2El<:GroupElement} <: GroupElement
+    h::G1El
+    k::G2El
+end
+function GroupsCore.gens(G::DirectSumGroup)
+    elts = [DirectSumElem(h, one(G.K)) for h in gens(G.H)]
+    append!(elts, [DirectSumElem(one(G.H), k) for k in gens(G.K)])
+    return elts
 end
-Base.one(c::Union{CyclicGroup, CyclicElem}) = CyclicElem(c.n, 0)
-PermutationGroups.gens(c::CyclicGroup) = [CyclicElem(c.n, 1)]
-PermutationGroups.order(::Type{T}, c::CyclicGroup) where {T} = convert(T, c.n)
-function Base.iterate(c::CyclicGroup, prev::CyclicElem=CyclicElem(c.n, -1))
-    id = prev.id + 1
-    if id >= c.n
-        return nothing
+GroupsCore.order(::Type{T}, G::DirectSumGroup) where {T} =
+    convert(T, order(T, G.H) * order(T, G.K))
+
+function Base.IteratorSize(::Type{DirectSumGroup{H,K}}) where {H,K}
+    if Base.IteratorSize(H) isa Base.IsInfinite || Base.IteratorSize(K) isa Base.IsInfinite
+        return Base.IsInfinite()
+    elseif Base.IteratorSize(H) isa Base.SizeUnknown || Base.IteratorSize(K) isa Base.SizeUnknown
+        return Base.SizeUnknown()
     else
-        next = CyclicElem(c.n, id)
-        return next, next
+        return Base.HasShape{2}()
     end
 end
 
-#SymbolicWedderburn.conjugacy_classes_orbit(KleinPermGroup())
+Base.size(G::DirectSumGroup) = (order(Int, G.H), order(Int, G.K))
 
-struct DirectSum <: GroupElement
-    c::CyclicElem
-    k::KleinPermElement
-end
-function PermutationGroups.order(el::DirectSum)
-    return lcm(PermutationGroups.order(el.c), PermutationGroups.order(el.k))
+function Base.eltype(::Type{DirectSumGroup{H,K}}) where {H,K}
+    return DirectSumElem{eltype(H),eltype(K)}
 end
 
-function Base.hash(el::DirectSum, u::UInt64)
-    return hash(el.k, hash(el.c, u))
-end
-function Base.:(==)(a::DirectSum, b::DirectSum)
-    return a.c == b.c && a.k == b.k
+function Base.iterate(G::DirectSumGroup)
+    IT = Iterators.product(G.H, G.K)
+    el, st = iterate(IT)
+    return DirectSumElem(el...), (IT, st)
 end
 
-function Base.inv(el::DirectSum)
-    DirectSum(inv(el.c), inv(el.k))
+function Base.iterate(::DirectSumGroup, it_st)
+    el_st = iterate(it_st...)
+    isnothing(el_st) && return nothing
+    el, st = el_st
+    return DirectSumElem(el...), (first(it_st), st)
 end
 
-# k * p * k' * q = k * k'^p * p * q
-function Base.:*(a::DirectSum, b::DirectSum)
-    return DirectSum(a.c * b.c, a.k * b.k)
-end
-Base.:^(a::DirectSum, k::Integer) = DirectSum(a.c^k, a.k^k)
+Base.one(G::DirectSumGroup) = DirectSumElem(one(G.H), one(G.K))
 
-Base.conj(a::DirectSum, b::DirectSum) = inv(b) * a * b
-Base.:^(a::DirectSum, b::DirectSum) = conj(a, b)
+Base.parent(g::DirectSumElem) = DirectSumGroup(parent(g.h), parent(g.k))
+function GroupsCore.order(::Type{T}, g::DirectSumElem) where {T}
+    return convert(T, lcm(order(T, g.h), order(T, g.k)))
+end
+Base.:(==)(g1::DirectSumElem, g2::DirectSumElem) = g1.h == g2.h && g1.k == g2.k
+Base.hash(el::DirectSumElem, h::UInt) = hash(el.h, hash(el.k, hash(eltype(el), h)))
+Base.inv(g::DirectSumElem) = DirectSumElem(inv(g.h), inv(g.k))
+Base.:*(g1::DirectSumElem, g2::DirectSumElem) = DirectSumElem(g1.h * g2.h, g1.k * g2.k)
+Base.:^(g::DirectSumElem, n::Integer) = DirectSumElem(g.h^n, g.k^n)
 
-struct Lattice1Group <: Group
-    n::Int
+function Base.rand(rng::Random.AbstractRNG, st::Random.SamplerTrivial{<:DirectSumGroup})
+    G = st[]
+    h = rand(rng, G.H)
+    k = rand(rng, G.K)
+    return DirectSumElem(h, k)
 end
-Base.eltype(::Lattice1Group) = DirectSum
-Base.one(el::DirectSum) = DirectSum(one(el.c), one(el.k))
-Base.one(L::Lattice1Group) = DirectSum(CyclicElem(L.n, 0), one(KleinPermGroup()))
-function PermutationGroups.gens(L::Lattice1Group)
-    els = DirectSum[]
-    for g in gens(CyclicGroup(L.n))
-        push!(els, DirectSum(g, one(KleinPermGroup())))
-    end
-    for g in gens(KleinPermGroup())
-        push!(els, DirectSum(one(CyclicGroup(L.n)), g))
+
+function Base.deepcopy_internal(g::DirectSumElem, dict::IdDict)
+    h = if isbits(g.h)
+        g.h
+    elseif haskey(dict, objectid(g.h))
+        dict[objectid(g.h)]
+    else
+        Base.deepcopy_internal(g.h, dict)
     end
-    return els
-end
-PermutationGroups.order(L::Lattice1Group) = PermutationGroups.order(Int, L)
-function PermutationGroups.order(::Type{T}, L::Lattice1Group) where {T}
-    return order(T, CyclicGroup(L.n)) * order(T, KleinPermGroup())
-end
-Base.length(L::Lattice1Group) = PermutationGroups.order(L)
-function Base.iterate(L::Lattice1Group)
-    el_p, st_p = iterate(CyclicGroup(L.n))
-    el_k, st_k = iterate(KleinPermGroup())
-    return DirectSum(el_p, el_k), ((el_p, st_p), st_k)
-end
-function Base.iterate(L::Lattice1Group, st)
-    el_st_k = iterate(KleinPermGroup(), st[2])
-    if el_st_k === nothing
-        el_st_p = iterate(CyclicGroup(L.n), st[1][2])
-        el_st_p === nothing && return nothing
-        el_k, st_k = iterate(KleinPermGroup())
-        return DirectSum(el_st_p[1], el_k), (el_st_p, st_k)
+
+    k = if isbits(g.k)
+        g.k
+    elseif haskey(dict, objectid(g.k))
+        dict[objectid(g.k)]
+    else
+        Base.deepcopy_internal(g.k, dict)
     end
-    return DirectSum(st[1][1], el_st_k[1]), (st[1], el_st_k[2])
+    return DirectSumElem(h, k)
 end
+
+# two aliases
+const Lattice1Group = DirectSumGroup{CyclicGroup,KleinPermGroup}
+const DirectSum = DirectSumElem{CyclicElem,KleinPermElement}
+
+"""
+    Lattice1Group(n::Integer)
+
+The direct sum of the cyclic group of order `n` and `KleinPermGroup()`.
+"""
+(::Type{Lattice1Group})(n::Integer) = DirectSumGroup(CyclicGroup(n), KleinPermGroup())