move unary/binary operations to arithmetic.jl

This commit is contained in:
kalmarek 2019-06-04 20:06:39 +02:00
parent 006c1fcfb2
commit a55d2e18f6
No known key found for this signature in database
GPG Key ID: 8BF1A3855328FC15
2 changed files with 265 additions and 204 deletions

View File

@ -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
###############################################################################
#

265
src/arithmetic.jl Normal file
View File

@ -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