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