Merge branch 'enh/pmfree_mul'

src/GroupRings.jl

@@ -1,7 +1,7 @@
 module GroupRings
 
 using Nemo
-import Nemo: Group, GroupElem, Ring, RingElem, parent, elem_type, parent_type
+import Nemo: Group, GroupElem, Ring, RingElem, parent, elem_type, parent_type, mul!, addeq!, divexact
 
 import Base: convert, show, hash, ==, +, -, *, //, /, length, norm, rationalize, deepcopy_internal, getindex, setindex!, eltype, one, zero
 
@@ -17,23 +17,40 @@ type GroupRing{Gr<:Group, T<:GroupElem} <: Ring
    basis_dict::Dict{T, Int}
    pm::Array{Int,2}
 
-   function GroupRing(G::Gr; initialise=true)
+   function GroupRing(G::Group, basis::Vector{T}; init::Bool=false)
-      A = new(G)
+      RG = new(G, basis, reverse_dict(basis))
-      if initialise
+      if init
-         complete(A)
+         RG.pm = try
+            create_pm(RG.basis, RG.basis_dict)
+         catch err
+            isa(err, KeyError) && throw("Product is not supported on basis")
+            throw(err)
+         end
+      else
+         RG.pm = zeros(Int, length(basis), length(basis))
       end
-      return A
+      return RG
    end
 
    function GroupRing(G::Gr, basis::Vector{T}, basis_dict::Dict{T,Int}, pm::Array{Int,2})
       return new(G, basis, basis_dict, pm)
    end
+
+   function GroupRing(G::Gr, pm::Array{Int,2})
+      RG = new(G)
+      RG.pm = pm
+      return RG
+   end
 end
 
-GroupRing{Gr<:Group}(G::Gr;initialise=true) = GroupRing{Gr, elem_type(G)}(G, initialise=initialise)
+GroupRing{Gr<:Group, T<:GroupElem}(G::Gr, basis::Vector{T}; init=false) =
+   GroupRing{Gr, T}(G, basis, init=init)
 
 GroupRing{Gr<:Group, T<:GroupElem}(G::Gr, b::Vector{T}, b_d::Dict{T,Int}, pm::Array{Int,2}) = GroupRing{Gr, T}(G, b, b_d, pm)
 
+GroupRing{Gr<:Group}(G::Gr, pm::Array{Int,2}) =
+   GroupRing{Gr, elem_type(G)}(G, pm)
+
 type GroupRingElem{T<:Number} <: RingElem
    coeffs::AbstractVector{T}
    parent::GroupRing
@@ -52,7 +69,7 @@ type GroupRingElem{T<:Number} <: RingElem
    end
 end
 
-export GroupRing, GroupRingElem, complete, create_pm
+export GroupRing, GroupRingElem, complete!, create_pm, star
 
 ###############################################################################
 #
@@ -65,8 +82,16 @@ elem_type(::GroupRing) = GroupRingElem
 parent_type(::GroupRingElem) = GroupRing
 parent_type(::Type{GroupRingElem}) = GroupRing
 
+eltype(X::GroupRingElem) = eltype(X.coeffs)
+
 parent{T}(g::GroupRingElem{T}) = g.parent
 
+Base.promote_rule{T<:Number,S<:Number}(::Type{GroupRingElem{T}}, ::Type{GroupRingElem{S}}) = GroupRingElem{promote_type(T,S)}
+
+function convert{T<:Number}(::Type{T}, X::GroupRingElem)
+   return GroupRingElem(convert(AbstractVector{T}, X.coeffs), parent(X))
+end
+
 ###############################################################################
 #
 #   GroupRing / GroupRingElem constructors
@@ -77,32 +102,13 @@ function GroupRingElem{T<:Number}(c::AbstractVector{T}, RG::GroupRing)
    return GroupRingElem{T}(c, RG)
 end
 
-function convert{T<:Number}(::Type{T}, X::GroupRingElem)
+function GroupRing(G::Group; init::Bool=false)
-   return GroupRingElem(convert(AbstractVector{T}, X.coeffs), parent(X))
+   return GroupRing(G, [elements(G)...], init=init)
-end
-
-function GroupRing(G::Group, pm::Array{Int,2})
-   size(pm,1) == size(pm,2) || throw("pm must be square, got $(size(pm))")
-   RG = GroupRing(G, initialise=false)
-   RG.pm = pm
-   return RG
-end
-
-function GroupRing(G::Group, basis::Vector)
-   basis_dict = reverse_dict(basis)
-   pm = try
-      create_pm(basis, basis_dict)
-   catch err
-      isa(err, KeyError) && throw("Products are not supported on basis")
-      throw(err)
-   end
-   return GroupRing(G, basis, basis_dict, pm)
 end
 
 function GroupRing(G::Group, basis::Vector, pm::Array{Int,2})
-   size(pm,1) == size(pm,2) || throw("pm must be of size (n,n), got
+   size(pm,1) == size(pm,2) || throw("pm must be square, got $(size(pm))")
-      $(size(pm))")
+   eltype(basis) == elem_type(G) || throw("Basis must consist of elements of $G")
-   eltype(basis) == elem_type(G) || throw("basis must consist of elements of $G")
    basis_dict = reverse_dict(basis)
    return GroupRing(G, basis, basis_dict, pm)
 end
@@ -113,23 +119,29 @@ end
 #
 ###############################################################################
 
+zero(RG::GroupRing, T::Type=Int) = RG(T)
+one(RG::GroupRing, T::Type=Int) = RG(RG.group(), T)
+
+function (RG::GroupRing)(i::Int, T::Type=Int)
+   elt = RG(T)
+   elt[RG.group()] = i
+   return elt
+end
+
 function (RG::GroupRing)(T::Type=Int)
-   isdefined(RG, :basis) || throw("Complete the definition of GroupRing first")
+   isdefined(RG, :basis) || throw("Can not coerce without basis of GroupRing")
    return GroupRingElem(spzeros(T,length(RG.basis)), RG)
 end
 
 function (RG::GroupRing)(g::GroupElem, T::Type=Int)
-   g = try
+   g = RG.group(g)
-      RG.group(g)
-   catch
-      throw("Can't coerce $g to the underlying group of $RG")
-   end
    result = RG(T)
    result[g] = one(T)
    return result
 end
 
 function (RG::GroupRing)(x::AbstractVector)
+   isdefined(RG, :basis) || throw("Can not coerce without basis of GroupRing")
    length(x) == length(RG.basis) || throw("Can not coerce to $RG: lengths differ")
    result = RG(eltype(x))
    result.coeffs = x
@@ -142,16 +154,19 @@ function (RG::GroupRing)(X::GroupRingElem)
 end
 
 function (RG::GroupRing)(X::GroupRingElem, emb::Function)
-    result = RG(eltype(X.coeffs))
+   isdefined(RG, :basis) || throw("Can not coerce without basis of GroupRing")
-    for g in parent(X).basis
+   result = RG(eltype(X.coeffs))
-        result[emb(g)] = X[g]
+   T = typeof(X.coeffs)
-    end
+   result.coeffs = T(result.coeffs)
-    return result
+   for g in parent(X).basis
+      result[emb(g)] = X[g]
+   end
+   return result
 end
 
 ###############################################################################
 #
-#   Basic manipulation
+#   Basic manipulation && Array protocol
 #
 ###############################################################################
 
@@ -160,7 +175,7 @@ function deepcopy_internal(X::GroupRingElem, dict::ObjectIdDict)
 end
 
 function hash(X::GroupRingElem, h::UInt)
-   return hash(X.coeffs, hash(parent(X), h))
+   return hash(full(X.coeffs), hash(parent(X), hash(GroupRingElem, h)))
 end
 
 function getindex(X::GroupRingElem, n::Int)
@@ -184,10 +199,8 @@ function setindex!(X::GroupRingElem, value, g::GroupElem)
    X.coeffs[RG.basis_dict[g]] = value
 end
 
-eltype(X::GroupRingElem) = eltype(X.coeffs)
+Base.size(X::GroupRingElem) = size(X.coeffs)
-
+Base.linearindexing{T<:GroupRingElem}(::Type{T}) = Base.LinearFast()
-one(RG::GroupRing) = RG(RG.group())
-zero(RG::GroupRing) = RG()
 
 ###############################################################################
 #
@@ -239,8 +252,8 @@ function (==)(A::GroupRing, B::GroupRing)
       A.basis == B.basis || return false
    else
       warn("Bases of GroupRings are not defined, comparing products mats.")
+      A.pm == B.pm || return false
    end
-   A.pm == B.pm || return false
    return true
 end
 
@@ -284,16 +297,23 @@ end
 #
 ###############################################################################
 
-function add{T<:Number}(X::GroupRingElem{T}, Y::GroupRingElem{T})
+function addeq!{T}(X::GroupRingElem{T}, Y::GroupRingElem{T})
-   parent(X) == parent(Y) || throw(ArgumentError(
+   X.coeffs .+= Y.coeffs
-   "Elements don't seem to belong to the same Group Ring!"))
+   return X
+end
+
+function add{T<:Number}(X::GroupRingElem{T}, Y::GroupRingElem{T}, check::Bool=true)
+   if check
+      parent(X) == parent(Y) || throw("Elements don't seem to belong to the same Group Ring!")
+   end
    return GroupRingElem(X.coeffs+Y.coeffs, parent(X))
 end
 
 function add{T<:Number, S<:Number}(X::GroupRingElem{T},
-   Y::GroupRingElem{S})
+   Y::GroupRingElem{S}, check::Bool=true)
-   parent(X) == parent(Y) || throw(ArgumentError(
+   if check
-   "Elements don't seem to belong to the same Group Ring!"))
+      parent(X) == parent(Y) || throw("Elements don't seem to belong to the same Group Ring!")
+   end
    warn("Adding elements with different base rings!")
    return GroupRingElem(+(promote(X.coeffs, Y.coeffs)...), parent(X))
 end
@@ -301,7 +321,24 @@ end
 (+)(X::GroupRingElem, Y::GroupRingElem) = add(X,Y)
 (-)(X::GroupRingElem, Y::GroupRingElem) = add(X,-Y)
 
-function mul!{T}(result::AbstractVector{T}, X::AbstractVector, Y::AbstractVector, pm::Array{Int,2})
+doc"""
+    mul!{T}(result::AbstractArray{T},
+                 X::AbstractVector,
+                 Y::AbstractVector,
+                pm::Array{Int,2})
+> The most specialised multiplication for `X` and `Y` (`coeffs` of
+> `GroupRingElems`) using multiplication table `pm`.
+> Notes:
+> * this method will silently produce false results if `X[k]` is non-zero for
+> `k > size(pm,1)`.
+> * This method will fail if any zeros (i.e. uninitialised entries) are present
+> in `pm`.
+> * Use with extreme care!
+"""
+function mul!{T}(result::AbstractVector{T},
+                      X::AbstractVector,
+                      Y::AbstractVector,
+                     pm::Array{Int,2})
    z = zero(T)
    result .= z
    for j in eachindex(Y)
@@ -314,39 +351,95 @@ function mul!{T}(result::AbstractVector{T}, X::AbstractVector, Y::AbstractVector
          end
       end
    end
-end
-
-function mul{T<:Number}(X::AbstractVector{T}, Y::AbstractVector{T},
-   pm::Array{Int,2})
-   result = zeros(X)
-   mul!(X,Y,pm,result)
    return result
 end
 
-function mul(X::AbstractVector, Y::AbstractVector, pm::Array{Int,2})
+doc"""
-   T = promote_type(eltype(X), eltype(Y))
+    mul!{T}(result::GroupRingElem{T},
-   result = zeros(T, deepcopy(X))
+                 X::GroupRingElem,
-   mul!(X, Y, pm, result)
+                 Y::GroupRingElem)
+> In-place multiplication for `GroupRingElem`s `X` and `Y`.
+> `mul!` will make use the initialised entries of `pm` attribute of
+> `parent(X)::GroupRing` (if available), and will compute and store in `pm` the
+> remaining products.
+> The method will fail with `KeyError` if product `X*Y` is not supported on
+> `parent(X).basis`.
+"""
+function mul!{T}(result::GroupRingElem{T}, X::GroupRingElem, Y::GroupRingElem)
+   if result === X
+      result = deepcopy(result)
+   end
+
+   z = zero(T)
+   result.coeffs .= z
+
+   RG = parent(X)
+
+   for j::Int in eachindex(Y.coeffs)
+      if Y.coeffs[j] != z
+         for i::Int in eachindex(X.coeffs)
+            if X.coeffs[i] != z
+               if RG.pm[i,j] == 0
+                  g::elem_type(parent(X).group) = RG.basis[i]*RG.basis[j]
+                  RG.pm[i,j] = RG.basis_dict[g]
+               end
+               result.coeffs[RG.pm[i,j]] += X[i]*Y[j]
+            end
+         end
+      end
+   end
    return result
 end
 
-function *{T<:Number}(X::GroupRingElem{T}, Y::GroupRingElem{T})
+function *{T<:Number}(X::GroupRingElem{T}, Y::GroupRingElem{T}, check::Bool=true)
-   parent(X) == parent(Y) || throw(ArgumentError(
+   if check
-   "Elements don't seem to belong to the same Group Ring!"))
+      parent(X) == parent(Y) || throw("Elements don't seem to belong to the same Group Ring!")
-   RG = parent(X)
+   end
-   isdefined(RG, :pm) || complete(RG)
+   if isdefined(parent(X), :basis)
-   result = mul(X.coeffs, Y.coeffs, RG.pm)
+      result = parent(X)(similar(X.coeffs))
-   return GroupRingElem(result, RG)
+      result = mul!(result, X, Y)
+   else
+      result = mul!(similar(X.coeffs), X.coeffs, Y.coeffs, parent(X).pm)
+      result = GroupRingElem(result, parent(X))
+   end
+   return result
 end
 
-function *{T<:Number, S<:Number}(X::GroupRingElem{T}, Y::GroupRingElem{S})
+function *{T<:Number, S<:Number}(X::GroupRingElem{T}, Y::GroupRingElem{S}, check::Bool=true)
-   parent(X) == parent(Y) || throw("Elements don't seem to belong to the same
+   if true
-   Group Ring!")
+      parent(X) == parent(Y) || throw("Elements don't seem to belong to the same Group Ring!")
-   warn("Multiplying elements with different base rings!")
+   end
-   RG = parent(X)
+
-   isdefined(RG, :pm) || complete(RG)
+   TT = typeof(first(X.coeffs)*first(Y.coeffs))
-   result = mul(X.coeffs, Y.coeffs, RG.pm)
+   warn("Multiplying elements with different base rings! Promoting the result to $TT.")
-   return GroupRingElem(result, RG)
+
+   result = mul!(result, X, Y)
+   return result
+
+   if isdefined(parent(X), :basis)
+      result = parent(X)(similar(X.coeffs))
+      result = convert(TT, result)
+      result = mul!(result, X, Y)
+   else
+      result = convert(TT, similar(X.coeffs))
+      result = mul!(result, X.coeffs, Y.coeffs, parent(X).pm)
+      result = GroupRingElem(result, parent(X))
+   end
+   return result
+end
+
+
+
+function divexact{T}(X::GroupRingElem{T}, Y::GroupRingElem{T})
+   if length(Y) != 1
+      throw("Can not divide by a non-primitive element $(Y)!")
+   else
+      idx = findfirst(Y)
+      c = Y[idx]
+      c == 0 || throw("Can not invert")
+      g = parent(Y).basis[idx]
+      return X*1//c*parent(Y)(inv(g))
+   end
 end
 
 ###############################################################################
@@ -357,7 +450,7 @@ end
 
 function star{T}(X::GroupRingElem{T})
    RG = parent(X)
-   isdefined(RG, :basis) || complete(RG)
+   isdefined(RG, :basis) || throw("*-involution without basis is not possible")
    result = RG(T)
    for (i,c) in enumerate(X.coeffs)
       if c != zero(T)
@@ -413,22 +506,18 @@ end
 
 create_pm{T<:GroupElem}(b::Vector{T}) = create_pm(b, reverse_dict(b))
 
-function complete(A::GroupRing)
+function complete!(RG::GroupRing)
-   if !isdefined(A, :basis)
+   if !isdefined(RG, :basis)
-      A.basis = [elements(A.group)...]
+      RG.basis = [elements(RG.group)...]
    end
-   if !isdefined(A, :basis_dict)
+   if !isdefined(RG, :basis_dict)
-      A.basis_dict = reverse_dict(A.basis)
+      RG.basis_dict = reverse_dict(RG.basis)
    end
-   if !isdefined(A, :pm)
+   for linidx in find(RG.pm .== 0)
-      A.pm = try
+      i,j = ind2sub(size(RG.pm), linidx)
-         create_pm(A.basis, A.basis_dict)
+      RG.pm[i,j] = RG.basis_dict[RG.basis[i]*RG.basis[j]]
-      catch err
-         isa(err, KeyError) && throw("Product is not supported on basis")
-         throw(err)
-      end
    end
-   return A
+   return RG
 end
 
 end # of module GroupRings

test/runtests.jl

@@ -10,14 +10,16 @@ using Nemo
       @test isa(GroupRing(G), Nemo.Ring)
       @test isa(GroupRing(G), GroupRing)
 
-      RG = GroupRing(G, initialise=false)
+      RG = GroupRing(G, init=false)
-      @test isdefined(RG, :pm) == false
+      @test isdefined(RG, :basis) == true
-      @test isdefined(RG, :basis) == false
-      @test isdefined(RG, :basis_dict) == false
-
-      @test isa(complete(RG), GroupRing)
-      @test size(RG.pm) == (6,6)
       @test length(RG.basis) == 6
+      @test isdefined(RG, :basis_dict) == true
+      @test isdefined(RG, :pm) == true
+      @test RG.pm == zeros(Int, (6,6))
+
+      @test isa(complete!(RG), GroupRing)
+      @test all(RG.pm .> 0)
+      @test RG.pm == GroupRing(G, init=true).pm
 
       @test RG.basis_dict == GroupRings.reverse_dict(elements(G))
 
@@ -63,7 +65,7 @@ using Nemo
 
    @testset "GroupRingElems constructors/basic manipulation" begin
       G = PermutationGroup(3)
-      RG = GroupRing(G, initialise=true)
+      RG = GroupRing(G, init=true)
       a = rand(6)
       @test isa(GroupRingElem(a, RG), GroupRingElem)
       @test isa(RG(a), GroupRingElem)
@@ -120,6 +122,9 @@ using Nemo
          @test eltype(2.0*a) == typeof(2.0)
          @test (2.0*a).coeffs == 2.0.*(a.coeffs)
 
+         b = RG(1) + GroupRings.star(a)
+         @test a*b == mul!(a,a,b)
+
          @test isa(a/2, GroupRingElem)
          @test eltype(a/2) == typeof(1/2)
          @test (a/2).coeffs == 0.5*(a.coeffs)