From a55d2e18f6c3244786a201632deb4fdc6c5c14a3 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Tue, 4 Jun 2019 20:06:39 +0200 Subject: [PATCH] move unary/binary operations to arithmetic.jl --- src/GroupRings.jl | 204 ----------------------------------- src/arithmetic.jl | 265 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 265 insertions(+), 204 deletions(-) create mode 100644 src/arithmetic.jl diff --git a/src/GroupRings.jl b/src/GroupRings.jl index 531f360..f733533 100644 --- a/src/GroupRings.jl +++ b/src/GroupRings.jl @@ -125,210 +125,6 @@ function (==)(A::GroupRing, B::GroupRing) return true end -############################################################################### -# -# Scalar operators -# -############################################################################### - -(-)(X::GroupRingElem) = GroupRingElem(-X.coeffs, parent(X)) - -function mul!(a::T, X::GroupRingElem{T}) where {T<:Number} - X.coeffs .*= a - return X -end - -mul(a::T, X::GroupRingElem{T}) where {T<:Number} = GroupRingElem(a*X.coeffs, parent(X)) - -function mul(a::T, X::GroupRingElem{S}) where {T<:Number, S<:Number} - TT = promote_type(T,S) - TT == S || @warn("Scalar and coeffs are in different rings! Promoting result to $(TT)") - return GroupRingElem(a.*X.coeffs, parent(X)) -end - -(*)(a::Number, X::GroupRingElem) = mul(a, X) -(*)(X::GroupRingElem, a::Number) = mul(a, X) - -# disallow Rings to hijack *(::, ::GroupRingElem) -*(a::Union{AbstractFloat, Integer, RingElem, Rational}, X::GroupRingElem) = mul(a, X) - -(/)(X::GroupRingElem, a) = 1/a*X -(//)(X::GroupRingElem, a::Union{Integer, Rational}) = 1//a*X - -(^)(X::GroupRingElem, n::Integer) = Base.power_by_squaring(X, n) - -############################################################################### -# -# Binary operators -# -############################################################################### - -function addeq!(X::GroupRingElem, Y::GroupRingElem) - X.coeffs += Y.coeffs - return X -end - -function +(X::GroupRingElem{T}, Y::GroupRingElem{T}) where T - return GroupRingElem(X.coeffs+Y.coeffs, parent(X)) -end - -function +(X::GroupRingElem{S}, Y::GroupRingElem{T}) where {S, T} - @warn("Adding elements with different coefficient rings, Promoting result to $(promote_type(T,S))") - return GroupRingElem(X.coeffs+Y.coeffs, parent(X)) -end - --(X::GroupRingElem{T}, Y::GroupRingElem{T}) where T = addeq!((-Y), X) - -function -(X::GroupRingElem{S}, Y::GroupRingElem{T}) where {S, T} - @warn("Adding elements with different coefficient rings, Promoting result to $(promote_type(T,S))") - addeq!((-Y), X) -end - -@doc doc""" - fmac!(result::AbstractVector{T}, - X::AbstractVector, - Y::AbstractVector, - pm::Array{Int,2}) where {T<:Number} -> Fused multiply-add for group ring coeffs using multiplication table `pm`. -> The result of X*Y in GroupRing is added in-place to `result`. -> 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 fmac!(result::AbstractVector{T}, - X::AbstractVector, - Y::AbstractVector, - pm::Array{Int,2}) where {T<:Number} - z = zero(T) - s1 = size(pm,1) - s2 = size(pm,2) - - @inbounds for j in 1:s2 - if Y[j] != z - for i in 1:s1 - if X[i] != z - result[pm[i,j]] += X[i]*Y[j] - end - end - end - end - return result -end - -@doc doc""" - GRmul!(result::AbstractVector{T}, X::AbstractVector, Y::AbstractVector, - pm::Matrix{<:Integer}) where {T<:Number} -> The most specialised multiplication for `X` and `Y` (intended for `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 GRmul!(result::AbstractVector{T}, - X::AbstractVector, - Y::AbstractVector, - pm::AbstractMatrix{<:Integer}) where {T<:Number} - z = zero(T) - result .= z - - return fmac!(result, X, Y, pm) -end - -@doc doc""" - mul!(result::GroupRingElem, X::GroupRingElem, 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 necessary to perform the multiplication. -> The method will fail with `KeyError` if product `X*Y` is not supported on -> `parent(X).basis`. -""" -function mul!(result::GroupRingElem, X::GroupRingElem, Y::GroupRingElem) - if result === X - result = deepcopy(result) - end - T = eltype(result.coeffs) - z = zero(T) - result.coeffs .= z - - RG = parent(X) - - lX = length(X.coeffs) - lY = length(Y.coeffs) - - if isdefined(RG, :pm) - s = size(RG.pm) - k = findprev(!iszero, X.coeffs, lX) - (k == nothing ? 0 : k) <= s[1] || throw("Element in X outside of support of parents product") - k = findprev(!iszero, Y.coeffs, lY) - (k == nothing ? 0 : k) <= s[2] || throw("Element in Y outside of support of parents product") - - for j in 1:lY - if Y.coeffs[j] != z - for i in 1:lX - if X.coeffs[i] != z - if RG.pm[i,j] == 0 - RG.pm[i,j] = RG.basis_dict[RG.basis[i]*RG.basis[j]] - end - result.coeffs[RG.pm[i,j]] += X[i]*Y[j] - end - end - end - end - else - for j in 1:lY - if Y.coeffs[j] != z - for i in 1:lX - if X.coeffs[i] != z - result[RG.basis[i]*RG.basis[j]] += X[i]*Y[j] - end - end - end - end - end - return result -end - -function *(X::GroupRingElem{T}, Y::GroupRingElem{T}, check::Bool=true) where {T<:Number} - if check - parent(X) == parent(Y) || throw("Elements don't seem to belong to the same Group Ring!") - end - if isdefined(parent(X), :basis) - result = parent(X)(similar(X.coeffs)) - result = mul!(result, X, Y) - else - result = GRmul!(similar(X.coeffs), X.coeffs, Y.coeffs, parent(X).pm) - result = GroupRingElem(result, parent(X)) - end - return result -end - -function *(X::GroupRingElem{T}, Y::GroupRingElem{S}, check::Bool=true) where {T<:Number, S<:Number} - if check - parent(X) == parent(Y) || throw("Elements don't seem to belong to the same Group Ring!") - end - - TT = typeof(first(X.coeffs)*first(Y.coeffs)) - @warn("Multiplying elements with different base rings! Promoting the result to $TT.") - - 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 = RGmul!(result, X.coeffs, Y.coeffs, parent(X).pm) - result = GroupRingElem(result, parent(X)) - end - return result -end ############################################################################### # diff --git a/src/arithmetic.jl b/src/arithmetic.jl new file mode 100644 index 0000000..07c59f2 --- /dev/null +++ b/src/arithmetic.jl @@ -0,0 +1,265 @@ +############################################################################### +# +# Unsafe operators (NCRing interface) +# +############################################################################### + +@doc doc""" + zero!(X::GroupRingElem) +In-place set coefficients of `X` to `0`. The sparsity pattern of `X` will be preserved. +""" +zero!(X::GroupRingElem{T}) where T = ( fill!(X.coeffs, zero(T)); X ) + +@doc doc""" + mul!(result::GroupRingElem, X::GroupRingElem, Y::GroupRingElem) +Perform the multiplication `X * Y` and store the result in `result`. + +`mul!` will make use the initialised entries of `pm` attribute of `parent(X)` (if available), and will compute (and cache in `parent(X).pm`) the remaining products as necessary. + +# Notes: +* `result` is zeroed before use +* aliasing of `result` and `X` or `Y` is allowed (in the case a `similar` copy is created) +* `mul!` will throw `BoundsError` if `X * Y` is not supported on `parent(X).basis` +* no checks on arguments parents (i.e. mathematical correctns) are performed +""" +function mul!(result::GroupRingElem, X::GroupRingElem, Y::GroupRingElem) + + result = zero!(_dealias(result, X, Y)) + X_nzeros_idx = findall(!iszero, X.coeffs) + Y_nzeros_idx = findall(!iszero, Y.coeffs) + # X_nzeros_idx = [i for i in eachindex(X.coeffs) if X[i] != zero(eltype(X))] + # Y_nzeros_idx = [i for i in eachindex(Y.coeffs) if Y[i] != zero(eltype(Y))] + + RG = parent(X) + + if cachesmultiplication(RG) + complete!(RG, X_nzeros_idx, Y_nzeros_idx) + for j in Y_nzeros_idx + for i in X_nzeros_idx + result.coeffs[RG.pm[i,j]] += X[i]*Y[j] + end + end + else + for j in Y_nzeros_idx + for i in X_nzeros_idx + result[RG[i]*RG[j]] += X[i]*Y[j] + end + end + end + return result +end +@doc doc""" + add!(result::GroupRingElem, X::GroupRingElem, Y::GroupRingElem) +Perform te the addition `X + Y` and store the result in `result`. + +# Notes: +* `result` is zeroed before use +* aliasing of `result` and `X` or `Y` is allowed (in the case a `similar` copy is created) +* no checks on arguments parents (i.e. mathematical correctns) are performed +""" +function add!(result::GroupRingElem, X::GroupRingElem, Y::GroupRingElem) + result = _dealias(result, X, Y) + @inbounds for i in eachindex(result.coeffs) + result.coeffs[i] = X.coeffs[i] + Y.coeffs[i] + end + return result +end + +@doc doc""" + addeq!(X::GroupRingElem, Y::GroupRingElem) +Add `Y` to `X` in-place (the result is stored in `X`). + +# Notes: +* no checks on arguments parents (i.e. mathematical correctns) are performed +""" + +function addeq!(X::GroupRingElem, Y::GroupRingElem) + X.coeffs .+= Y.coeffs + return X +end + +function addmul!(result::GroupRingElem, X::GroupRingElem, Y::GroupRingElem, + tmp::GroupRingElem=similar(result)) + tmp = mul!(tmp, X, Y) + result = addeq!(result, tmp) + return result +end + +############################################################################### +# +# Arithmetic operations (NCRing interface) +# +############################################################################### + +-(X::GroupRingElem{T}) where T = GroupRingElem(-X.coeffs, parent(X)) +^(X::GroupRingElem, n::Integer) = Base.power_by_squaring(X, n) + +function Base.inv(X::GroupRingElem) + if isunit(X) + g = supp(X)[1] + i_g = eltype(X)(inv(X[g])) + return scalarmul!(base_ring(parent(X))(i_g), parent(X)(inv(g))) + end + throw(DivideError()) +end + +### Addition/Subtraction: + +function +(X::GroupRingElem{T, GR}, Y::GroupRingElem{T, GR}) where {T, GR<:GroupRing} + # @assert parent(X) == parent(Y) + return add!(X, X, Y) +end + +function +(X::GroupRingElem{S}, Y::GroupRingElem{T}) where {S,T} + # @assert parent(X) == parent(Y) + result = _promote(X, Y) + return add!(result, X, Y) +end + +# -Y creates a copy +-(X::GroupRingElem{T,GR}, Y::GroupRingElem{T,GR}) where {T,GR<:GroupRing} = addeq!(-Y, X) + +# X - Y => X + (-Y) for other parameters TODO: this allocates two elements instead of one + +### Multiplication: + +function _mul(result::GroupRingElem, X::GroupRingElem, Y::GroupRingElem) + RG = parent(X) + if hasbasis(RG) + result = mul!(result, X, Y) + else + result.coeffs = _mul!(result.coeffs, X.coeffs, Y.coeffs, RG.pm) + end + return result +end + +function *(X::GroupRingElem{T,GR}, Y::GroupRingElem{T,GR}) where {T, GR<:GroupRing} + # @assert parent(X) == parent(Y) + return _mul(X, X, Y) +end + +function *(X::GroupRingElem{S}, Y::GroupRingElem{T}) where {S,T} + # @assert parent(X) == parent(Y) + result = _promote(X, Y) + return _mul(result, X, Y) +end + +############################################################################### +# +# Scalar and Ad-hoc operators +# +############################################################################### + +### Addition/subtraction + +function addeq!(X::GroupRingElem{T}, a::T) where T + X[_identity_idx(parent(X))] += a + return X +end + +function addeq!(X::GroupRingElem{T, GroupRing{R,G,El}}, g::El, v=1) where {T,R,G,El} + @assert hasbasis(parent(X)) + X[g] += T(v) + return X +end + ++(X::GroupRingElem{T}, a::T) where T<:RingElement = addeq!(deepcopy(X), a) ++(a::T, X::GroupRingElem{T}) where T<:RingElement = addeq!(deepcopy(X), a) ++(X::GroupRingElem{T, GroupRing{R, G, El}}, g::El) where {T, R, G, El} = addeq!(deepcopy(X), g) ++(g::El, X::GroupRingElem{T, GroupRing{R, G, El}}) where {T, R, G, El} = addeq!(deepcopy(X), g) +# +(X::GroupRingElem{T}, a::S) where {S,T} = addeq!(deepcopy(X), base_ring(X)(a)) +# +(a::S, X::GroupRingElem{T}) where {S,T} = addeq!(deepcopy(X), base_ring(X)(a)) + +-(X::GroupRingElem{T}, a::T) where T<:RingElement = addeq!(deepcopy(X), -a) +-(a::T, X::GroupRingElem{T}) where T<:RingElement = addeq!(-X, a) +-(X::GroupRingElem{T, GroupRing{R, G, El}}, g::El) where {T, R, G, El} = addeq!(deepcopy(X), g, -1) +-(g::El, X::GroupRingElem{T, GroupRing{R, G, El}}) where {T, R, G, El} = addeq!(-X, g) +# -(X::GroupRingElem{T}, a::S) where {S,T} = addeq!(deepcopy(X), -base_ring(X)(a)) +# -(a::S, X::GroupRingElem{T}) where {S,T} = addeq!(-X, base_ring(X)(a)) + +### Scalar multiplication/scalar division + +scalarmul!(a::T, X::GroupRingElem{T}) where T<:RingElement = (X.coeffs .*= a; return X) + +function scalarmul(a::S, X::GroupRingElem{T}) where {S,T} + if promote_type(S, T) == T + return scalarmul!(base_ring(parent(X))(a), deepcopy(X)) + else + RG = change_base_ring(parent(X), parent(a)) + @warn "Coefficient ring does not contain scalar $a.\nThe result has coefficients in $(parent(a)) of type $(elem_type(parent(a)))." + return scalarmul!(a, GroupRingElem(base_ring(RG).(X.coeffs), RG)) + end +end + +*(a::T, X::GroupRingElem{T}) where T = scalarmul!(a, deepcopy(X)) +# assuming commutativity of base_ring +*(X::GroupRingElem{T}, a::T) where T = scalarmul!(a, deepcopy(X)) +*(a::T, X::GroupRingElem{S}) where {T, S} = scalarmul(a, X) +*(X::GroupRingElem{S}, a::T) where {T, S} = scalarmul(a, X) + +# deambiguations: +*(a::Union{AbstractFloat, Integer, RingElem, Rational}, X::GroupRingElem) = scalarmul(a, X) +*(X::GroupRingElem, a::Union{AbstractFloat, Integer, RingElem, Rational}) = scalarmul(a, X) + +# divisions +(/)(X::GroupRingElem, a) = inv(a)*X +(//)(X::GroupRingElem, a::Union{Integer, Rational}) = 1//a*X + +############################################################################### +# +# Unsafe operators on coefficients +# +############################################################################### + +@doc doc""" + _mul!(result::AbstractVector, X::AbstractVector, Y::AbstractVector, pm::AbstracctMatrix{<:Integer}) +Unsafe multiplication for group ring elements (represented by coefficient vectors `result`, `X` and `Y`) using the provided multiplication table `pm`. +Performs the multiplication `X * Y` and stores the result in `result`. + +# Notes: +* aliasing of `result` and `X` or `Y` is allowed (in the case a `similar` copy is created) +* `result` is zeroed before use +* `BoundsError` is thrown if `X * Y` is not determined by `pm` +* this method will silently produce false results if either `X[k]` or `Y[k]` +is non-zero for `k > size(pm,1)` +* use only if You know what you're doing! +""" +function _mul!(result::AbstractVector, X::AbstractVector, Y::AbstractVector, + pm::AbstractMatrix{<:Integer}) + # this mul! is used for base-less multiplication + if result === X || result === Y + result = zero(result) + else + result .= zero(eltype(result)) + end + return _addmul!(result, X, Y, pm) +end + +@doc doc""" + _addmul!(result::AbstractVector, X::AbstractVector, Y::AbstractVector, pm::AbstractMatrix{<:Integer}) +Unsafe fused multiply-add for group ring elements (represented by coefficient vectors `X` and `Y`) using the provided multiplication table `pm`. +Performs the multiplication `X * Y` and adds the result to `result` in-place. + +# Notes: +* aliasing of `result` and `X` or `Y` is NOT checked +* `BoundsError` is thrown if `X * Y` is not determined by `pm` +* this method will silently produce false results if either `X[k]` or `Y[k]` +is non-zero for `k > size(pm,1)` +* use only if You know what you're doing! +""" + +function _addmul!(result::AbstractVector, X::AbstractVector, Y::AbstractVector, + pm::AbstractMatrix{<:Integer}) + zX = zero(eltype(X)) + zY = zero(eltype(Y)) + @inbounds for j in 1:size(pm,2) + if Y[j] != zY + for i in 1:size(pm,1) + if X[i] != zX + result[pm[i,j]] += X[i]*Y[j] + end + end + end + end + return result +end