Merge pull request #50 from kalmarek/mk/rand_via_pra

implement rand via Product Replacement Algorithm

Project.toml

@@ -1,7 +1,7 @@
 name = "GroupsCore"
 uuid = "d5909c97-4eac-4ecc-a3dc-fdd0858a4120"
-authors = ["Marek Kaluba <[email protected]> and contributors"]
-version = "0.4.2"
+authors = ["Marek Kaluba <[email protected]>"]
+version = "0.5.0"
 
 [deps]
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"

src/GroupsCore.jl

@@ -15,6 +15,7 @@ include("exceptions.jl")
 include("groups.jl")
 include("group_elements.jl")
 
+include("rand.jl")
 include("extensions.jl")
 
 end

src/group_elements.jl

@@ -27,17 +27,6 @@ Base.:(==)(g::GEl, h::GEl) where {GEl <: GroupElement} = throw(
     InterfaceNotImplemented(:Group, "Base.:(==)(::$GEl, ::$GEl)"),
 )
 
-Base.deepcopy_internal(g::GroupElement, stackdict::IdDict) = throw(
-    InterfaceNotImplemented(
-        :Group,
-        "Base.deepcopy_internal(::$(typeof(g)), ::IdDict)",
-    ),
-)
-# TODO: Technically, it is not necessary to implement `deepcopy_internal` method
-# if `parent(g)` can be reconstructed exactly from `g` (i.e. either it's cached,
-# or a singleton). However by defining this fallback we force everybody to
-# implement it, except isbits group elements.
-
 Base.copy(g::GroupElement) = deepcopy(g)
 
 @doc Markdown.doc"""

src/groups.jl

@@ -50,23 +50,6 @@ The result of this function is undefined unless `GroupsCore.hasgens(G)` return
 gens(G::Group) =
     throw(InterfaceNotImplemented(:Group, "GroupsCore.gens(::$(typeof(G)))"))
 
-@doc Markdown.doc"""
-    Base.rand(rng::Random.AbstractRNG, rs::Random.SamplerTrivial{Gr}) where {Gr <: Group}
-
-Return a random group element, treating the group as a collection.
-"""
-function Base.rand(
-    rng::Random.AbstractRNG,
-    rs::Random.SamplerTrivial{Gr}
-) where {Gr <: Group}
-    throw(
-        InterfaceNotImplemented(
-            :Random,
-            "Base.rand(::Random.AbstractRNG, ::Random.SamplerTrivial{$Gr}))",
-        ),
-    )
-end
-
 ################################################################################
 # Iterators
 ################################################################################
@@ -164,6 +147,3 @@ function ngens(G::Group)
         "Group does not seem to have generators. Did you alter `hasgens(::$(typeof(G)))`?",
     )
 end
-
-rand_pseudo(G::Group, args...) = rand_pseudo(Random.GLOBAL_RNG, G, args...)
-rand_pseudo(rng::Random.AbstractRNG, G::Group, args...) = rand(rng, G, args...)

src/rand.jl

@@ -0,0 +1,61 @@
+import Random
+
+"""
+    PRASampler
+Product Replacement Algorithm sampler for arbitrary group generated by
+an explicit finite set of generators.
+"""
+mutable struct PRASampler{T} <: Random.Sampler{T}
+    gentuple::Vector{T}
+    right::T
+    left::T
+end
+
+# constants taken from GAP
+function PRASampler(G, n::Integer=2ngens(G) + 10, scramble_time::Integer=10max(n, 10))
+    return PRASampler(Random.default_rng(), G, n, scramble_time)
+end
+
+function Random.Sampler(
+    RNG::Type{<:Random.AbstractRNG},
+    G::Group,
+    repetition::Random.Repetition=Val(Inf)
+)
+    return PRASampler(RNG(), G)
+end
+
+function PRASampler(
+    rng::Random.AbstractRNG,
+    G::Group,
+    n::Integer=2ngens(G) + 10,
+    scramble_time::Integer=10max(n, 10)
+)
+    if istrivial(G)
+        return PRASampler(fill(one(G), n), one(G), one(G))
+    end
+    @assert hasgens(G)
+    l = max(n, 2ngens(G), 2)
+    sampler = let S = gens(G)
+        S = union!(inv.(S), S)
+        append!(S, rand(rng, S, l - length(S)))
+        PRASampler(S, one(G), one(G))
+    end
+    for _ in 1:scramble_time
+        _ = rand(rng, sampler)
+    end
+    return sampler
+end
+
+function Random.rand(rng::Random.AbstractRNG, pra::PRASampler)
+    range = 1:length(pra.gentuple)
+
+    pra.right = pra.right * rand(rng, pra.gentuple)
+
+    i = rand(rng, range)
+    @inbounds pra.gentuple[i] = pra.gentuple[i] * pra.right
+
+    j = rand(rng, range)
+    @inbounds pra.left = pra.left * pra.gentuple[j]
+
+    return pra.left
+end

test/conformance_test.jl

@@ -77,12 +77,6 @@ function test_Group_interface(G::Group)
             @test rand(G, 2) isa AbstractVector{eltype(G)}
             g, h = rand(G, 2)
             @test parent(g) === parent(h) === G
-
-            @test GroupsCore.rand_pseudo(G) isa eltype(G)
-            @test GroupsCore.rand_pseudo(G, 2, 2) isa AbstractMatrix{eltype(G)}
-
-            g, h = GroupsCore.rand_pseudo(G, 2)
-            @test parent(g) === parent(h) === G
         end
     end
 end

test/cyclic.jl

@@ -21,14 +21,6 @@ Base.IteratorSize(::Type{CyclicGroup}) = Base.HasLength()
 GroupsCore.order(::Type{T}, C::CyclicGroup) where {T<:Integer} = T(C.order)
 GroupsCore.gens(C::CyclicGroup) = [CyclicGroupElement(1, C)]
 
-function Base.rand(
-    rng::Random.AbstractRNG,
-    rs::Random.SamplerTrivial{<:CyclicGroup},
-)
-    C = rs[]
-    return CyclicGroupElement(rand(rng, 0:C.order-1), C)
-end
-
 GroupsCore.parent(c::CyclicGroupElement) = c.parent
 Base.:(==)(g::CyclicGroupElement, h::CyclicGroupElement) =
     parent(g) === parent(h) && g.residual == h.residual

test/infinite_cyclic.jl

@@ -17,13 +17,6 @@ Base.IteratorSize(::Type{InfCyclicGroup}) = Base.IsInfinite()
 
 GroupsCore.gens(C::InfCyclicGroup) = [InfCyclicGroupElement(1)]
 
-function Base.rand(
-    rng::Random.AbstractRNG,
-    rs::Random.SamplerTrivial{<:InfCyclicGroup},
-)
-    return InfCyclicGroupElement(rand(rng, Int))
-end
-
 GroupsCore.parent(c::InfCyclicGroupElement) = InfCyclicGroup()
 Base.:(==)(g::InfCyclicGroupElement, h::InfCyclicGroupElement) = g.val == h.val
 

test/test_notsatisfied.jl

@@ -91,12 +91,9 @@ end
         @test_throws INI gens(G)
 
         Base.eltype(::Type{SomeGroup}) = SomeGroupElement
-        @test_throws INI rand(G, 2)
 
         @test_throws INI gens(G, 1)
         @test_throws INI ngens(G)
-
-        @test_throws INI GroupsCore.rand_pseudo(G, 2, 2)
     end
 
     @testset "GroupElem Interface" begin
@@ -116,7 +113,6 @@ end
 
         GroupsCore.isfiniteorder(::SomeGroupElement) = false
         @test_throws InfO order(g)
-        @test_throws INI deepcopy(g)
 
         @test_throws INI inv(g)
         @test_throws INI g * g