From 3862e785c841e55586e5c906ef53d144efb2e47e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Lars=20G=C3=B6ttgens?= <lars.goettgens@rwth-aachen.de>
Date: Wed, 2 Oct 2024 12:10:36 +0200
Subject: [PATCH] LieAlgebras: reduce allocations (#4165)

---
 .gitignore                                    |   1 +
 .../LieAlgebras/src/AbstractLieAlgebra.jl     |  41 ++--
 experimental/LieAlgebras/src/LieAlgebra.jl    |  14 +-
 .../LieAlgebras/src/LieAlgebraIdeal.jl        |   2 +-
 .../LieAlgebras/src/LieAlgebraModule.jl       |   4 +-
 experimental/LieAlgebras/src/LieAlgebras.jl   |   1 +
 experimental/LieAlgebras/src/LieSubalgebra.jl |   2 +-
 .../LieAlgebras/src/LinearLieAlgebra.jl       |  14 +-
 experimental/LieAlgebras/src/RootSystem.jl    | 210 ++++++++++++++++--
 experimental/LieAlgebras/src/Types.jl         |  19 +-
 experimental/LieAlgebras/src/Util.jl          |  17 +-
 .../LieAlgebras/test/LinearLieAlgebra-test.jl |   2 +-
 .../LieAlgebras/test/RootSystem-test.jl       |  73 ++++++
 13 files changed, 331 insertions(+), 69 deletions(-)

diff --git a/.gitignore b/.gitignore
index ea513a4b2905..786f12eb31c3 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,6 +4,7 @@
 /docs/src/Nemo/
 /docs/src/Experimental/
 profile.pb.gz
+alloc-profile.pb.gz
 LocalPreferences.toml
 Manifest.toml
 .DS_Store
diff --git a/experimental/LieAlgebras/src/AbstractLieAlgebra.jl b/experimental/LieAlgebras/src/AbstractLieAlgebra.jl
index 4f5f16ffcd98..2cb4c4dbb848 100644
--- a/experimental/LieAlgebras/src/AbstractLieAlgebra.jl
+++ b/experimental/LieAlgebras/src/AbstractLieAlgebra.jl
@@ -265,9 +265,10 @@ function lie_algebra(
 
   npos = n_positive_roots(rs)
   # set Chevalley basis
-  r_plus = basis(L)[1:npos]
-  r_minus = basis(L)[(npos + 1):(2 * npos)]
-  h = basis(L)[(2 * npos + 1):end]
+  basis_L = basis(L)
+  r_plus = basis_L[1:npos]
+  r_minus = basis_L[(npos + 1):(2 * npos)]
+  h = basis_L[(2 * npos + 1):end]
   chev = (r_plus, r_minus, h)
   set_root_system_and_chevalley_basis!(L, rs, chev)
   return L
@@ -286,6 +287,7 @@ function _struct_consts(R::Field, rs::RootSystem, extraspecial_pair_signs)
   N = _N_matrix(rs, extraspecial_pair_signs)
 
   struct_consts = Matrix{sparse_row_type(R)}(undef, n, n)
+  beta_i_plus_beta_j = zero(RootSpaceElem, rs)
   for i in 1:nroots, j in i:nroots
     if i == j
       # [e_βi, e_βi] = 0
@@ -294,12 +296,14 @@ function _struct_consts(R::Field, rs::RootSystem, extraspecial_pair_signs)
     end
     beta_i = root(rs, i)
     beta_j = root(rs, j)
-    if iszero(beta_i + beta_j)
+    beta_i_plus_beta_j = add!(beta_i_plus_beta_j, beta_i, beta_j)
+    if iszero(beta_i_plus_beta_j)
       # [e_βi, e_-βi] = h_βi
       struct_consts[i, j] = sparse_row(
-        R, collect((nroots + 1):(nroots + nsimp)), _vec(coefficients(coroot(rs, i)))
+        R, collect((nroots + 1):(nroots + nsimp)), _vec(coefficients(coroot(rs, i)));
+        sort=false,
       )
-    elseif ((fl, k) = is_root_with_index(beta_i + beta_j); fl)
+    elseif ((fl, k) = is_root_with_index(beta_i_plus_beta_j); fl)
       # complicated case
       if i <= npos
         struct_consts[i, j] = sparse_row(R, [k], [N[i, j]])
@@ -313,7 +317,7 @@ function _struct_consts(R::Field, rs::RootSystem, extraspecial_pair_signs)
       struct_consts[i, j] = sparse_row(R)
     end
 
-    # # [e_βj, e_βi] = -[e_βi, e_βj]
+    # [e_βj, e_βi] = -[e_βi, e_βj]
     struct_consts[j, i] = -struct_consts[i, j]
   end
   for i in 1:nsimp, j in 1:nroots
@@ -321,7 +325,8 @@ function _struct_consts(R::Field, rs::RootSystem, extraspecial_pair_signs)
     struct_consts[nroots + i, j] = sparse_row(
       R,
       [j],
-      [dot(coefficients(root(rs, j)), view(cm, i, :))],
+      # [dot(coefficients(root(rs, j)), view(cm, i, :))], # currently the below is faster
+      [only(coefficients(root(rs, j)) * cm[i, :])],
     )
     # [e_βj, h_i] = -[h_i, e_βj]
     struct_consts[j, nroots + i] = -struct_consts[nroots + i, j]
@@ -352,25 +357,24 @@ function _N_matrix(rs::RootSystem, extraspecial_pair_signs::Vector{Bool})
         j -> !is_positive_root(alpha_i + beta_l - simple_root(rs, j)),
         1:(i - 1),
       ) || continue
-      p = 0
-      while is_root(beta_l - p * alpha_i)
-        p += 1
-      end
-      N[i, l] = (extraspecial_pair_signs[k - nsimp] ? 1 : -1) * p
+      p = _root_string_length_down(beta_l, alpha_i) + 1
+      N[i, l] = (extraspecial_pair_signs[k - nsimp] ? p : -p)
       N[l, i] = -N[i, l]
     end
   end
 
   # special pairs
+  alpha_i_plus_beta_j = zero(RootSpaceElem, rs)
   for (i, alpha_i) in enumerate(positive_roots(rs))
     for (j, beta_j) in enumerate(positive_roots(rs))
       i < j || continue
-      is_positive_root(alpha_i + beta_j) || continue
+      alpha_i_plus_beta_j = add!(alpha_i_plus_beta_j, alpha_i, beta_j)
+      is_positive_root(alpha_i_plus_beta_j) || continue
       l = findfirst(
-        l -> is_positive_root(alpha_i + beta_j - simple_root(rs, l)), 1:nsimp
+        l -> is_positive_root(alpha_i_plus_beta_j - simple_root(rs, l)), 1:nsimp
       )::Int
       l == i && continue # already extraspecial
-      fl, l_comp = is_positive_root_with_index(alpha_i + beta_j - simple_root(rs, l))
+      fl, l_comp = is_positive_root_with_index(alpha_i_plus_beta_j - simple_root(rs, l))
       @assert fl
       t1 = 0
       t2 = 0
@@ -383,10 +387,7 @@ function _N_matrix(rs::RootSystem, extraspecial_pair_signs::Vector{Bool})
         t2 = N[l, m] * N[j, m] * dot(root_m, root_m)//dot(alpha_i, alpha_i)
       end
       @assert t1 - t2 != 0
-      p = 0
-      while is_root(beta_j - p * alpha_i)
-        p += 1
-      end
+      p = _root_string_length_down(beta_j, alpha_i) + 1
       N[i, j] = Int(sign(t1 - t2) * sign(N[l, l_comp]) * p) # typo in CMT04
       N[j, i] = -N[i, j]
     end
diff --git a/experimental/LieAlgebras/src/LieAlgebra.jl b/experimental/LieAlgebras/src/LieAlgebra.jl
index e0855e1d006b..32cf8d0bbfd6 100644
--- a/experimental/LieAlgebras/src/LieAlgebra.jl
+++ b/experimental/LieAlgebras/src/LieAlgebra.jl
@@ -49,7 +49,9 @@ Return the `i`-th basis element of the Lie algebra `L`.
 function basis(L::LieAlgebra, i::Int)
   @req 1 <= i <= dim(L) "Index out of bounds."
   R = coefficient_ring(L)
-  return L([(j == i ? one(R) : zero(R)) for j in 1:dim(L)])
+  v = zero_matrix(R, 1, dim(L))
+  v[1, i] = one(R)
+  return L(v)
 end
 
 @doc raw"""
@@ -453,10 +455,12 @@ end
 
 function adjoint_matrix(x::LieAlgebraElem{C}) where {C<:FieldElem}
   L = parent(x)
-  return sum(
-    c * g for (c, g) in zip(coefficients(x), adjoint_matrices(L));
-    init=zero_matrix(coefficient_ring(L), dim(L), dim(L)),
-  )
+  mat = zero_matrix(coefficient_ring(L), dim(L), dim(L))
+  tmp = zero(mat)
+  for (c, g) in zip(coefficients(x), adjoint_matrices(L))
+    mat = addmul!(mat, g, c, tmp)
+  end
+  return mat
 end
 
 function _adjoint_matrix(S::LieSubalgebra{C}, x::LieAlgebraElem{C}) where {C<:FieldElem}
diff --git a/experimental/LieAlgebras/src/LieAlgebraIdeal.jl b/experimental/LieAlgebras/src/LieAlgebraIdeal.jl
index 73d8db05890b..35efee83d50a 100644
--- a/experimental/LieAlgebras/src/LieAlgebraIdeal.jl
+++ b/experimental/LieAlgebras/src/LieAlgebraIdeal.jl
@@ -270,7 +270,7 @@ where `L` is the Lie algebra where `I` lives in.
 function lie_algebra(I::LieAlgebraIdeal)
   LI = lie_algebra(basis(I))
   L = base_lie_algebra(I)
-  emb = hom(LI, L, basis(I))
+  emb = hom(LI, L, basis(I); check=false)
   return LI, emb
 end
 
diff --git a/experimental/LieAlgebras/src/LieAlgebraModule.jl b/experimental/LieAlgebras/src/LieAlgebraModule.jl
index d1aba5521814..500c2b90b0fa 100644
--- a/experimental/LieAlgebras/src/LieAlgebraModule.jl
+++ b/experimental/LieAlgebras/src/LieAlgebraModule.jl
@@ -53,7 +53,9 @@ Return the `i`-th basis element of the Lie algebra module `V`.
 function basis(V::LieAlgebraModule, i::Int)
   @req 1 <= i <= dim(V) "Index out of bounds."
   R = coefficient_ring(V)
-  return V([(j == i ? one(R) : zero(R)) for j in 1:dim(V)])
+  v = zero_matrix(R, 1, dim(V))
+  v[1, i] = one(R)
+  return V(v)
 end
 
 @doc raw"""
diff --git a/experimental/LieAlgebras/src/LieAlgebras.jl b/experimental/LieAlgebras/src/LieAlgebras.jl
index 9b264bf3ff43..d6cd094e36c6 100644
--- a/experimental/LieAlgebras/src/LieAlgebras.jl
+++ b/experimental/LieAlgebras/src/LieAlgebras.jl
@@ -84,6 +84,7 @@ import ..Oscar:
   tensor_product,
   weyl_vector,
   word,
+  zero!,
   zero_map,
   ⊕,
   ⊗
diff --git a/experimental/LieAlgebras/src/LieSubalgebra.jl b/experimental/LieAlgebras/src/LieSubalgebra.jl
index f83919ff27cb..dc6190154c5c 100644
--- a/experimental/LieAlgebras/src/LieSubalgebra.jl
+++ b/experimental/LieAlgebras/src/LieSubalgebra.jl
@@ -301,7 +301,7 @@ where `L` is the Lie algebra where `S` lives in.
 function lie_algebra(S::LieSubalgebra)
   LS = lie_algebra(basis(S))
   L = base_lie_algebra(S)
-  emb = hom(LS, L, basis(S))
+  emb = hom(LS, L, basis(S); check=false)
   return LS, emb
 end
 
diff --git a/experimental/LieAlgebras/src/LinearLieAlgebra.jl b/experimental/LieAlgebras/src/LinearLieAlgebra.jl
index 61d8183922ac..bab0fa3eee61 100644
--- a/experimental/LieAlgebras/src/LinearLieAlgebra.jl
+++ b/experimental/LieAlgebras/src/LinearLieAlgebra.jl
@@ -21,7 +21,7 @@ Return the basis `basis(L)` of the Lie algebra `L` in the underlying matrix
 representation.
 """
 function matrix_repr_basis(L::LinearLieAlgebra{C}) where {C<:FieldElem}
-  return Vector{dense_matrix_type(C)}(L.basis)
+  return L.basis::Vector{dense_matrix_type(C)}
 end
 
 @doc raw"""
@@ -31,7 +31,7 @@ Return the `i`-th element of the basis `basis(L)` of the Lie algebra `L` in the
 underlying matrix representation.
 """
 function matrix_repr_basis(L::LinearLieAlgebra{C}, i::Int) where {C<:FieldElem}
-  return (L.basis[i])::dense_matrix_type(C)
+  return matrix_repr_basis(L)[i]
 end
 
 ###############################################################################
@@ -130,10 +130,12 @@ Return the Lie algebra element `x` in the underlying matrix representation.
 """
 function Generic.matrix_repr(x::LinearLieAlgebraElem)
   L = parent(x)
-  return sum(
-    c * b for (c, b) in zip(_matrix(x), matrix_repr_basis(L));
-    init=zero_matrix(coefficient_ring(L), L.n, L.n),
-  )
+  mat = zero_matrix(coefficient_ring(L), L.n, L.n)
+  tmp = zero(mat)
+  for (c, b) in zip(coefficients(x), matrix_repr_basis(L))
+    mat = addmul!(mat, b, c, tmp)
+  end
+  return mat
 end
 
 function bracket(
diff --git a/experimental/LieAlgebras/src/RootSystem.jl b/experimental/LieAlgebras/src/RootSystem.jl
index 170a1b8f855c..bebe0380064e 100644
--- a/experimental/LieAlgebras/src/RootSystem.jl
+++ b/experimental/LieAlgebras/src/RootSystem.jl
@@ -115,6 +115,11 @@ Returns the `i`-th coroot of `R`, i.e. the `i`-th root of the dual root system o
 This is a more efficient version for `coroots(R)[i]`.
 
 Also see: `coroots`.
+
+!!! note
+    This function does not return a copy of the asked for object,
+    but the internal field of the root system.
+    Mutating the returned object will lead to undefined behavior.
 """
 function coroot(R::RootSystem, i::Int)
   if i <= n_positive_roots(R)
@@ -131,6 +136,11 @@ Returns the coroots of `R`, starting with the coroots of positive roots and then
 in the order of `positive_coroots` and `negative_coroots`.
 
 Also see: `coroot`.
+
+!!! note
+    This function does not return a copy of the asked for object,
+    but the internal field of the root system.
+    Mutating the returned object will lead to undefined behavior.
 """
 function coroots(R::RootSystem)
   return [[r for r in positive_coroots(R)]; [-r for r in positive_coroots(R)]]
@@ -172,6 +182,11 @@ Returns the `i`-th negative root of `R`.
 This is a more efficient version for `negative_roots(R)[i]`.
 
 Also see: `negative_roots`.
+
+!!! note
+    This function does not return a copy of the asked for object,
+    but the internal field of the root system.
+    Mutating the returned object will lead to undefined behavior.
 """
 function negative_root(R::RootSystem, i::Int)
   return -(R.positive_roots::Vector{RootSpaceElem})[i]
@@ -183,6 +198,11 @@ end
 Returns the negative roots of `R`. The $i$-th element of the returned vector is the negative root corresponding to the $i$-th positive root.
 
 Also see: `negative_root`.
+
+!!! note
+    This function does not return a copy of the asked for object,
+    but the internal field of the root system.
+    Mutating the returned object will lead to undefined behavior.
 """
 function negative_roots(R::RootSystem)
   return [-r for r in positive_roots(R)]
@@ -195,6 +215,11 @@ Returns the coroot corresponding to the `i`-th negative root of `R`
 This is a more efficient version for `negative_coroots(R)[i]`.
 
 Also see: `negative_coroots`.
+
+!!! note
+    This function does not return a copy of the asked for object,
+    but the internal field of the root system.
+    Mutating the returned object will lead to undefined behavior.
 """
 function negative_coroot(R::RootSystem, i::Int)
   return -(R.positive_coroots::Vector{DualRootSpaceElem})[i]
@@ -206,6 +231,11 @@ end
 Returns the coroots corresponding to the negative roots of `R`
 
 Also see: `negative_coroots`.
+
+!!! note
+    This function does not return a copy of the asked for object,
+    but the internal field of the root system.
+    Mutating the returned object will lead to undefined behavior.
 """
 function negative_coroots(R::RootSystem)
   return [-r for r in positive_coroots(R)]
@@ -251,6 +281,11 @@ Returns the `i`-th positive root of `R`.
 This is a more efficient version for `positive_roots(R)[i]`.
 
 Also see: `positive_roots`.
+
+!!! note
+    This function does not return a copy of the asked for object,
+    but the internal field of the root system.
+    Mutating the returned object will lead to undefined behavior.
 """
 function positive_root(R::RootSystem, i::Int)
   return (R.positive_roots::Vector{RootSpaceElem})[i]
@@ -263,6 +298,11 @@ Returns the positive roots of `R`, starting with the simple roots in the order o
 and then increasing in height.
 
 Also see: `positive_root`, `number_of_positive_roots`.
+
+!!! note
+    This function does not return a copy of the asked for object,
+    but the internal field of the root system.
+    Mutating the returned object will lead to undefined behavior.
 """
 function positive_roots(R::RootSystem)
   return R.positive_roots::Vector{RootSpaceElem}
@@ -275,6 +315,11 @@ Returns the coroot corresponding to the `i`-th positive root of `R`
 This is a more efficient version for `positive_coroots(R)[i]`.
 
 Also see: `positive_coroots`.
+
+!!! note
+    This function does not return a copy of the asked for object,
+    but the internal field of the root system.
+    Mutating the returned object will lead to undefined behavior.
 """
 function positive_coroot(R::RootSystem, i::Int)
   return (R.positive_coroots::Vector{DualRootSpaceElem})[i]
@@ -286,6 +331,11 @@ end
 Returns the coroots corresponding to the positive roots of `R`
 
 Also see: `positive_coroots`.
+
+!!! note
+    This function does not return a copy of the asked for object,
+    but the internal field of the root system.
+    Mutating the returned object will lead to undefined behavior.
 """
 function positive_coroots(R::RootSystem)
   return R.positive_coroots::Vector{DualRootSpaceElem}
@@ -336,6 +386,11 @@ Returns the `i`-th root of `R`.
 This is a more efficient version for `roots(R)[i]`.
 
 Also see: `roots`.
+
+!!! note
+    This function does not return a copy of the asked for object,
+    but the internal field of the root system.
+    Mutating the returned object will lead to undefined behavior.
 """
 function root(R::RootSystem, i::Int)
   if i <= n_positive_roots(R)
@@ -352,6 +407,11 @@ Returns the roots of `R`, starting with the positive roots and then the negative
 in the order of `positive_roots` and `negative_roots`.
 
 Also see: `root`.
+
+!!! note
+    This function does not return a copy of the asked for object,
+    but the internal field of the root system.
+    Mutating the returned object will lead to undefined behavior.
 """
 function roots(R::RootSystem)
   return [[r for r in positive_roots(R)]; [-r for r in positive_roots(R)]]
@@ -364,6 +424,11 @@ Returns the `i`-th simple root of `R`.
 This is a more efficient version for `simple_roots(R)[i]`.
 
 Also see: `simple_roots`.
+
+!!! note
+    This function does not return a copy of the asked for object,
+    but the internal field of the root system.
+    Mutating the returned object will lead to undefined behavior.
 """
 function simple_root(R::RootSystem, i::Int)
   @req 1 <= i <= rank(R) "Invalid index"
@@ -376,6 +441,11 @@ end
 Returns the simple roots of `R`.
 
 Also see: `simple_root`.
+
+!!! note
+    This function does not return a copy of the asked for object,
+    but the internal field of the root system.
+    Mutating the returned object will lead to undefined behavior.
 """
 function simple_roots(R::RootSystem)
   return positive_roots(R)[1:rank(R)]
@@ -388,6 +458,11 @@ Returns the coroot corresponding to the `i`-th simple root of `R`
 This is a more efficient version for `simple_coroots(R)[i]`.
 
 Also see: `simple_coroots`.
+
+!!! note
+    This function does not return a copy of the asked for object,
+    but the internal field of the root system.
+    Mutating the returned object will lead to undefined behavior.
 """
 function simple_coroot(R::RootSystem, i::Int)
   @req 1 <= i <= rank(R) "Invalid index"
@@ -400,6 +475,11 @@ end
 Returns the coroots corresponding to the simple roots of `R`
 
 Also see: `simple_coroots`.
+
+!!! note
+    This function does not return a copy of the asked for object,
+    but the internal field of the root system.
+    Mutating the returned object will lead to undefined behavior.
 """
 function simple_coroots(R::RootSystem)
   return positive_coroots(R)[1:rank(R)]
@@ -444,6 +524,14 @@ function RootSpaceElem(w::WeightLatticeElem)
   return RootSpaceElem(root_system(w), w)
 end
 
+function zero(::Type{RootSpaceElem}, R::RootSystem)
+  return RootSpaceElem(R, zero_matrix(QQ, 1, rank(R)))
+end
+
+function zero(r::RootSpaceElem)
+  return zero(RootSpaceElem, root_system(r))
+end
+
 function Base.:*(q::RationalUnion, r::RootSpaceElem)
   return RootSpaceElem(root_system(r), q * r.vec)
 end
@@ -464,6 +552,29 @@ function Base.:-(r::RootSpaceElem)
   return RootSpaceElem(root_system(r), -r.vec)
 end
 
+function zero!(r::RootSpaceElem)
+  r.vec = zero!(r.vec)
+  return r
+end
+
+function add!(rr::RootSpaceElem, r1::RootSpaceElem, r2::RootSpaceElem)
+  @req root_system(rr) === root_system(r1) === root_system(r2) "parent root system mismatch"
+  rr.vec = add!(rr.vec, r1.vec, r2.vec)
+  return rr
+end
+
+function neg!(rr::RootSpaceElem, r::RootSpaceElem)
+  @req root_system(rr) === root_system(r) "parent root system mismatch"
+  rr.vec = neg!(rr.vec, r.vec)
+  return rr
+end
+
+function sub!(rr::RootSpaceElem, r1::RootSpaceElem, r2::RootSpaceElem)
+  @req root_system(rr) === root_system(r1) === root_system(r2) "parent root system mismatch"
+  rr.vec = sub!(rr.vec, r1.vec, r2.vec)
+  return rr
+end
+
 function Base.:(==)(r::RootSpaceElem, r2::RootSpaceElem)
   return r.root_system === r2.root_system && r.vec == r2.vec
 end
@@ -505,7 +616,10 @@ end
 function dot(r1::RootSpaceElem, r2::RootSpaceElem)
   @req root_system(r1) === root_system(r2) "parent root system mismatch"
 
-  return dot(coefficients(r1) * _bilinear_form_QQ(root_system(r1)), coefficients(r2))
+  # return dot(coefficients(r1) * _bilinear_form_QQ(root_system(r1)), coefficients(r2)) # currently the below is faster
+  return only(
+    coefficients(r1) * _bilinear_form_QQ(root_system(r1)) * transpose(coefficients(r2))
+  )
 end
 
 @doc raw"""
@@ -602,6 +716,14 @@ function DualRootSpaceElem(root_system::RootSystem, vec::Vector{<:RationalUnion}
   return DualRootSpaceElem(root_system, matrix(QQ, 1, length(vec), vec))
 end
 
+function zero(::Type{DualRootSpaceElem}, R::RootSystem)
+  return DualRootSpaceElem(R, zero_matrix(QQ, 1, rank(R)))
+end
+
+function zero(r::DualRootSpaceElem)
+  return zero(DualRootSpaceElem, root_system(r))
+end
+
 function Base.:*(q::RationalUnion, r::DualRootSpaceElem)
   return DualRootSpaceElem(root_system(r), q * r.vec)
 end
@@ -622,6 +744,29 @@ function Base.:-(r::DualRootSpaceElem)
   return DualRootSpaceElem(root_system(r), -r.vec)
 end
 
+function zero!(r::DualRootSpaceElem)
+  r.vec = zero!(r.vec)
+  return r
+end
+
+function add!(rr::DualRootSpaceElem, r1::DualRootSpaceElem, r2::DualRootSpaceElem)
+  @req root_system(rr) === root_system(r1) === root_system(r2) "parent root system mismatch"
+  rr.vec = add!(rr.vec, r1.vec, r2.vec)
+  return rr
+end
+
+function neg!(rr::DualRootSpaceElem, r::DualRootSpaceElem)
+  @req root_system(rr) === root_system(r) "parent root system mismatch"
+  rr.vec = neg!(rr.vec, r.vec)
+  return rr
+end
+
+function sub!(rr::DualRootSpaceElem, r1::DualRootSpaceElem, r2::DualRootSpaceElem)
+  @req root_system(rr) === root_system(r1) === root_system(r2) "parent root system mismatch"
+  rr.vec = sub!(rr.vec, r1.vec, r2.vec)
+  return rr
+end
+
 function Base.:(==)(r::DualRootSpaceElem, r2::DualRootSpaceElem)
   return r.root_system === r2.root_system && r.vec == r2.vec
 end
@@ -759,18 +904,26 @@ function WeightLatticeElem(r::RootSpaceElem)
   return WeightLatticeElem(root_system(r), r)
 end
 
+function zero(::Type{WeightLatticeElem}, R::RootSystem)
+  return WeightLatticeElem(R, zero_matrix(ZZ, rank(R), 1))
+end
+
+function zero(r::WeightLatticeElem)
+  return zero(WeightLatticeElem, root_system(r))
+end
+
 function Base.:*(n::IntegerUnion, w::WeightLatticeElem)
   return WeightLatticeElem(root_system(w), n * w.vec)
 end
 
 function Base.:+(w::WeightLatticeElem, w2::WeightLatticeElem)
-  @req root_system(w) === root_system(w2) "parent weight lattics mismatch"
+  @req root_system(w) === root_system(w2) "parent root system mismatch"
 
   return WeightLatticeElem(root_system(w), w.vec + w2.vec)
 end
 
 function Base.:-(w::WeightLatticeElem, w2::WeightLatticeElem)
-  @req root_system(w) === root_system(w2) "parent weight lattics mismatch"
+  @req root_system(w) === root_system(w2) "parent root system mismatch"
 
   return WeightLatticeElem(root_system(w), w.vec - w2.vec)
 end
@@ -779,25 +932,29 @@ function Base.:-(w::WeightLatticeElem)
   return WeightLatticeElem(root_system(w), -w.vec)
 end
 
+function zero!(w::WeightLatticeElem)
+  w.vec = zero!(w.vec)
+  return w
+end
+
 function add!(wr::WeightLatticeElem, w1::WeightLatticeElem, w2::WeightLatticeElem)
-  add!(wr.vec, w1.vec, w2.vec)
+  @req root_system(wr) === root_system(w1) === root_system(w2) "parent root system mismatch"
+  wr.vec = add!(wr.vec, w1.vec, w2.vec)
   return wr
 end
 
-add!(w1::WeightLatticeElem, w2::WeightLatticeElem) = add!(w1, w1, w2)
-
-function neg!(w::WeightLatticeElem)
-  neg!(w.vec)
-  return w
+function neg!(wr::WeightLatticeElem, w::WeightLatticeElem)
+  @req root_system(wr) === root_system(w) "parent root system mismatch"
+  wr.vec = neg!(wr.vec, w.vec)
+  return wr
 end
 
 function sub!(wr::WeightLatticeElem, w1::WeightLatticeElem, w2::WeightLatticeElem)
-  sub!(wr.vec, w1.vec, w2.vec)
+  @req root_system(wr) === root_system(w1) === root_system(w2) "parent root system mismatch"
+  wr.vec = sub!(wr.vec, w1.vec, w2.vec)
   return wr
 end
 
-sub!(w1::WeightLatticeElem, w2::WeightLatticeElem) = sub!(w1, w1, w2)
-
 function Base.:(==)(w::WeightLatticeElem, w2::WeightLatticeElem)
   return w.root_system === w2.root_system && w.vec == w2.vec
 end
@@ -903,7 +1060,7 @@ function conjugate_dominant_weight_with_elem!(w::WeightLatticeElem)
 end
 
 function dot(w1::WeightLatticeElem, w2::WeightLatticeElem)
-  @req root_system(w1) === root_system(w2) "parent weight lattices mismatch"
+  @req root_system(w1) === root_system(w2) "parent root system mismatch"
   R = root_system(w1)
 
   return dot(
@@ -962,6 +1119,18 @@ function dot(w::WeightLatticeElem, r::RootSpaceElem)
   return dot(r, w)
 end
 
+# computes the maximum `p` such that `beta - p*alpha` is still a root
+# beta is assumed to be a root
+function _root_string_length_down(alpha::RootSpaceElem, beta::RootSpaceElem)
+  p = 0
+  beta_sub_p_alpha = beta - alpha
+  while is_root(beta_sub_p_alpha)
+    p += 1
+    beta_sub_p_alpha = sub!(beta_sub_p_alpha, alpha)
+  end
+  return p
+end
+
 @doc raw"""
     dim_of_simple_module([T = Int], R::RootSystem, hw::WeightLatticeElem) -> T
     dim_of_simple_module([T = Int], R::RootSystem, hw::Vector{<:IntegerUnion}) -> T
@@ -1300,8 +1469,8 @@ end
 function positive_roots_and_reflections(cartan_matrix::ZZMatrix)
   rank, _ = size(cartan_matrix)
 
-  roots = [[l == s ? 1 : 0 for l in 1:rank] for s in 1:rank]
-  coroots = [[l == s ? 1 : 0 for l in 1:rank] for s in 1:rank]
+  roots = [[l == s ? one(ZZ) : zero(ZZ) for l in 1:rank] for s in 1:rank]
+  coroots = [[l == s ? one(ZZ) : zero(ZZ) for l in 1:rank] for s in 1:rank]
   rootidx = Dict(roots[s] => s for s in 1:rank)
   refl = Dict((s, s) => 0 for s in 1:rank)
 
@@ -1312,12 +1481,11 @@ function positive_roots_and_reflections(cartan_matrix::ZZMatrix)
         continue
       end
 
-      pairing = let i = i
-        sum(roots[i][l] * cartan_matrix[s, l] for l in 1:rank; init=zero(ZZ))
-      end
-      copairing = let i = i
-        sum(coroots[i][l] * cartan_matrix[l, s] for l in 1:rank; init=zero(ZZ))
-      end
+      # pairing = dot(roots[i], view(cartan_matrix, s, :)) # currently the below is faster
+      pairing = only(view(cartan_matrix, s:s, :) * roots[i])
+      # copairing = dot(coroots[i], view(cartan_matrix, :, s)) # currently the below is faster
+      copairing = only(coroots[i] * view(cartan_matrix, :, s:s))
+
       if pairing * copairing >= 4
         refl[s, i] = 0
         continue
diff --git a/experimental/LieAlgebras/src/Types.jl b/experimental/LieAlgebras/src/Types.jl
index 3273a0bf4fa8..b32beb2f30aa 100644
--- a/experimental/LieAlgebras/src/Types.jl
+++ b/experimental/LieAlgebras/src/Types.jl
@@ -180,11 +180,19 @@ abstract type LieAlgebraElem{C<:FieldElem} <: AbstractAlgebra.SetElem end
       ) "Not anti-symmetric."
       @req all(
         iszero,
-        sum(
-          struct_consts[i, j][k] * struct_consts[k, l] +
-          struct_consts[j, l][k] * struct_consts[k, i] +
-          struct_consts[l, i][k] * struct_consts[k, j] for k in 1:dimL; init=sparse_row(R)
-        ) for i in 1:dimL, j in 1:dimL, l in 1:dimL
+        begin
+          row = sparse_row(R)
+          for (k, k_val) in struct_consts[i, j]
+            Hecke.add_scaled_row!(struct_consts[k, l], row, k_val)
+          end
+          for (k, k_val) in struct_consts[j, l]
+            Hecke.add_scaled_row!(struct_consts[k, i], row, k_val)
+          end
+          for (k, k_val) in struct_consts[l, i]
+            Hecke.add_scaled_row!(struct_consts[k, j], row, k_val)
+          end
+          row
+        end for i in 1:dimL, j in 1:dimL, l in 1:dimL
       ) "Jacobi identity does not hold."
     end
     return new{C}(R, dimL, struct_consts, s)
@@ -219,6 +227,7 @@ end
   ) where {C<:FieldElem}
     @req all(b -> size(b) == (n, n), basis) "Invalid basis element dimensions."
     @req length(s) == length(basis) "Invalid number of basis element names."
+    @req eltype(basis) == dense_matrix_type(R) "Invalid basis element type."
     L = new{C}(R, n, length(basis), basis, s)
     if check
       @req all(b -> all(e -> parent(e) === R, b), basis) "Invalid matrices."
diff --git a/experimental/LieAlgebras/src/Util.jl b/experimental/LieAlgebras/src/Util.jl
index 3190fafa80d6..6c2ab263bba9 100644
--- a/experimental/LieAlgebras/src/Util.jl
+++ b/experimental/LieAlgebras/src/Util.jl
@@ -14,16 +14,17 @@ function coefficient_vector(M::MatElem{T}, basis::Vector{<:MatElem{T}}) where {T
   (nr, nc) = size(M)
   all(b -> size(b) == (nr, nc), basis) ||
     throw(DimensionMismatch("The matrices in B must have the same size as M."))
-  lgs = similar(M, nr * nc, length(basis))
+  @req all(b -> base_ring(b) == base_ring(M), basis) "Incompatible base rings"
+  lgs = zero(M, length(basis), nr * nc)
   for (k, b) in enumerate(basis)
-    for i in 1:nr, j in 1:nc
-      lgs[(i - 1) * nc + j, k] = b[i, j]
+    for i in 1:nr
+      lgs[k, (i - 1) * nc .+ (1:nc)] = view(b, i:i, :)
     end
   end
-  rhs = similar(M, nr * nc, 1)
-  for i in 1:nr, j in 1:nc
-    rhs[(i - 1) * nc + j, 1] = M[i, j]
+  rhs = similar(M, 1, nr * nc)
+  for i in 1:nr
+    rhs[1, (i - 1) * nc .+ (1:nc)] = view(M, i:i, :)
   end
-  sol = solve(lgs, rhs; side=:right)
-  return transpose(sol)
+  sol = solve(lgs, rhs; side=:left)
+  return sol
 end
diff --git a/experimental/LieAlgebras/test/LinearLieAlgebra-test.jl b/experimental/LieAlgebras/test/LinearLieAlgebra-test.jl
index 930f44737cdf..5a24fd851ce0 100644
--- a/experimental/LieAlgebras/test/LinearLieAlgebra-test.jl
+++ b/experimental/LieAlgebras/test/LinearLieAlgebra-test.jl
@@ -17,7 +17,7 @@
 
   @testset "conformance tests" begin
     @testset "0-dim Lie algebra /QQ" begin
-      L = lie_algebra(QQ, 0, MatElem{QQFieldElem}[], Symbol[])
+      L = lie_algebra(QQ, 0, QQMatrix[], Symbol[])
       lie_algebra_conformance_test(
         L, LinearLieAlgebra{QQFieldElem}, LinearLieAlgebraElem{QQFieldElem}
       )
diff --git a/experimental/LieAlgebras/test/RootSystem-test.jl b/experimental/LieAlgebras/test/RootSystem-test.jl
index 278000c23b73..71cd8001372f 100644
--- a/experimental/LieAlgebras/test/RootSystem-test.jl
+++ b/experimental/LieAlgebras/test/RootSystem-test.jl
@@ -132,6 +132,79 @@
         )
       end
 
+      @testset "Mutating arithmetics for $T" for T in (
+        RootSpaceElem, DualRootSpaceElem, WeightLatticeElem
+      )
+        rk = rank(R)
+
+        x1 = T(R, rand(-10:10, rk))
+        x2 = T(R, rand(-10:10, rk))
+        x3 = T(R, rand(-10:10, rk))
+        x1c = deepcopy(x1)
+        x2c = deepcopy(x2)
+        x3c = deepcopy(x3)
+
+        for (f, f!) in ((zero, zero!),)
+          x1 = f!(x1)
+          @test x1 == f(x1c)
+          x1 = deepcopy(x1c)
+        end
+
+        for (f, f!) in ((-, neg!),)
+          x1 = f!(x1, x2)
+          @test x1 == f(x2c)
+          @test x2 == x2c
+          x1 = deepcopy(x1c)
+          x2 = deepcopy(x2c)
+
+          x1 = f!(x1)
+          @test x1 == f(x1c)
+          x1 = deepcopy(x1c)
+        end
+
+        for (f, f!) in ((+, add!), (-, sub!))
+          x1 = f!(x1, x2, x3)
+          @test x1 == f(x2c, x3c)
+          @test x2 == x2c
+          @test x3 == x3c
+          x1 = deepcopy(x1c)
+          x2 = deepcopy(x2c)
+          x3 = deepcopy(x3c)
+
+          x1 = f!(x1, x1, x2)
+          @test x1 == f(x1c, x2c)
+          @test x2 == x2c
+          x1 = deepcopy(x1c)
+          x2 = deepcopy(x2c)
+
+          x1 = f!(x1, x2, x1)
+          @test x1 == f(x2c, x1c)
+          @test x2 == x2c
+          x1 = deepcopy(x1c)
+          x2 = deepcopy(x2c)
+
+          x1 = f!(x1, x2, x2)
+          @test x1 == f(x2c, x2c)
+          @test x2 == x2c
+          x1 = deepcopy(x1c)
+          x2 = deepcopy(x2c)
+
+          x1 = f!(x1, x1, x1)
+          @test x1 == f(x1c, x1c)
+          x1 = deepcopy(x1c)
+
+          x1 = f!(x1, x2)
+          @test x1 == f(x1c, x2c)
+          @test x2 == x2c
+          x1 = deepcopy(x1c)
+          x2 = deepcopy(x2c)
+
+          x1 = f!(x1, x1)
+          @test x1 == f(x1c, x1c)
+          x1 = deepcopy(x1c)
+        end
+      end
+
       @testset "Serialization" begin
         mktempdir() do path
           test_save_load_roundtrip(path, R) do loaded