diff --git a/experimental/BasisLieHighestWeight/docs/doc.main b/experimental/BasisLieHighestWeight/docs/doc.main
new file mode 100644
index 000000000000..670207105f72
--- /dev/null
+++ b/experimental/BasisLieHighestWeight/docs/doc.main
@@ -0,0 +1,6 @@
+[
+   "Bases for highest weight modules" => [
+      "introduction.md",
+      "user_functions.md"
+   ]
+]
diff --git a/experimental/BasisLieHighestWeight/docs/src/introduction.md b/experimental/BasisLieHighestWeight/docs/src/introduction.md
new file mode 100644
index 000000000000..48ff61e54637
--- /dev/null
+++ b/experimental/BasisLieHighestWeight/docs/src/introduction.md
@@ -0,0 +1,25 @@
+```@meta
+CurrentModule = Oscar
+DocTestSetup = quote
+  using Oscar
+end
+```
+
+# Introduction
+
+This subproject contains the code for the OSCAR book chapter by Ghislain Fourier and Xin Fang on the bases of highest weight modules.
+
+## Status
+
+This part of OSCAR is in an experimental state; please see [Adding new projects to experimental](@ref) for what this means.
+More documentation is to come in the future.
+
+## Contact
+
+Please direct questions about this part of OSCAR to the following people:
+* [Ghislain Fourier](https://www.art.rwth-aachen.de/cms/~rnko/)
+* [Lars Göttgens](https://lgoe.li/)
+
+You can ask questions in the [OSCAR Slack](https://www.oscar-system.org/community/#slack).
+
+Alternatively, you can [raise an issue on github](https://www.oscar-system.org/community/#how-to-report-issues).
diff --git a/experimental/BasisLieHighestWeight/docs/src/user_functions.md b/experimental/BasisLieHighestWeight/docs/src/user_functions.md
new file mode 100644
index 000000000000..c450d9c45f26
--- /dev/null
+++ b/experimental/BasisLieHighestWeight/docs/src/user_functions.md
@@ -0,0 +1,16 @@
+```@meta
+CurrentModule = Oscar
+DocTestSetup = quote
+  using Oscar
+end
+```
+
+# User-facing functions
+```@docs
+basis_lie_highest_weight_operators
+basis_lie_highest_weight
+basis_lie_highest_weight_lusztig
+basis_lie_highest_weight_nz
+basis_lie_highest_weight_pbw
+basis_lie_highest_weight_string
+```
diff --git a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl
new file mode 100644
index 000000000000..e4429e9ae401
--- /dev/null
+++ b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl
@@ -0,0 +1,59 @@
+module BasisLieHighestWeight
+
+using ..Oscar
+using ..Oscar: GAPWrap, IntegerUnion, _is_weighted
+
+using AbstractAlgebra.PrettyPrinting
+
+import Oscar: dim, monomial_ordering, monomials
+
+import Base: length
+
+# TODO: Test if ZZx should be a graded_polynomial_ring with weights_w as weights
+
+# TODO (?) Maybe export and docstring: 
+# get_dim_weightspace
+# orbit_weylgroup
+# get_lattice_points_of_weightspace
+# convert_lattice_points_to_monomials
+# convert_monomials_to_lattice_points
+# action_matrices_of_operators
+# weights_for_operators
+
+# TODO Use Oscar-Lie-Algebra type instead of LieAlgebra
+# TODO Data-Type for weights of Lie-Algebras? Two types, in alpha_i and w_i, conversion is defined in RootConversion
+# w_to_aplha
+# alpha_to_w
+
+# TODO GAPWrap-wrappers are missing
+
+include("LieAlgebras.jl")
+include("BirationalSequence.jl")
+include("MonomialBasis.jl")
+include("NewMonomial.jl")
+include("TensorModels.jl")
+include("MonomialOrder.jl")
+include("RootConversion.jl")
+include("WeylPolytope.jl")
+include("MainAlgorithm.jl")
+include("UserFunctions.jl")
+
+export basis_lie_highest_weight_operators
+export basis_lie_highest_weight
+export basis_lie_highest_weight_lusztig
+export basis_lie_highest_weight_nz
+export basis_lie_highest_weight_pbw
+export basis_lie_highest_weight_string
+
+end
+
+using .BasisLieHighestWeight
+
+export BasisLieHighestWeight
+
+export basis_lie_highest_weight_operators
+export basis_lie_highest_weight
+export basis_lie_highest_weight_lusztig
+export basis_lie_highest_weight_nz
+export basis_lie_highest_weight_pbw
+export basis_lie_highest_weight_string
diff --git a/experimental/BasisLieHighestWeight/src/BirationalSequence.jl b/experimental/BasisLieHighestWeight/src/BirationalSequence.jl
new file mode 100644
index 000000000000..0d0de239c830
--- /dev/null
+++ b/experimental/BasisLieHighestWeight/src/BirationalSequence.jl
@@ -0,0 +1,12 @@
+struct BirationalSequence
+  operators::Vector{GAP.Obj}
+  operators_vectors::Vector{Vector{Any}}
+  weights_w::Vector{Vector{ZZRingElem}}
+  weights_alpha::Vector{Vector{QQFieldElem}}
+end
+
+function Base.show(io::IO, birational_sequence::BirationalSequence)
+  println(io, "BirationalSequence")
+  println(io, "Operators: ", birational_sequence.operators)
+  print(io, "Weights in alpha_i:", birational_sequence.weights_alpha)
+end
diff --git a/experimental/BasisLieHighestWeight/src/LieAlgebras.jl b/experimental/BasisLieHighestWeight/src/LieAlgebras.jl
new file mode 100644
index 000000000000..3ef6b6c751d6
--- /dev/null
+++ b/experimental/BasisLieHighestWeight/src/LieAlgebras.jl
@@ -0,0 +1,83 @@
+@attributes mutable struct LieAlgebraStructure
+  lie_type::Symbol
+  rank::Int
+  lie_algebra_gap::GAP.Obj
+  chevalley_basis::NTuple{3,Vector{GAP.Obj}}
+
+  function LieAlgebraStructure(lie_type::Symbol, rank::Int)
+    lie_algebra_gap = GAP.Globals.SimpleLieAlgebra(
+      GAP.Obj(lie_type), rank, GAP.Globals.Rationals
+    )
+    chevalley_basis = NTuple{3,Vector{GAP.Obj}}(GAP.Globals.ChevalleyBasis(lie_algebra_gap))
+    return new(lie_type, rank, lie_algebra_gap, chevalley_basis)
+  end
+end
+
+rank(L::LieAlgebraStructure) = L.rank
+
+@attr QQMatrix function cartan_matrix(L::LieAlgebraStructure)
+  R = GAP.Globals.RootSystem(L.lie_algebra_gap)
+  C = matrix(QQ, GAP.Globals.CartanMatrix(R))
+  return C
+end
+
+@attr QQMatrix function inv_cartan_matrix(L::LieAlgebraStructure)
+  return inv(cartan_matrix(L))
+end
+
+function Base.show(io::IO, L::LieAlgebraStructure)
+  io = pretty(io)
+  print(io, LowercaseOff(), "Lie algebra of type $(L.lie_type)$(L.rank)")
+end
+
+function lie_algebra(type::Symbol, rk::Int)
+  return LieAlgebraStructure(type, rk)
+end
+
+function chevalley_basis_gap(L::LieAlgebraStructure)
+  return L.chevalley_basis
+end
+
+function cartan_sub_basis(L::LieAlgebraStructure)
+  return L.chevalley_basis[3]
+end
+
+function root_system_gap(L::LieAlgebraStructure)
+  return GAP.Globals.RootSystem(L.lie_algebra_gap)
+end
+
+function num_positive_roots(L::LieAlgebraStructure)
+  return length(GAP.Globals.PositiveRoots(root_system_gap(L)))
+end
+
+function matrices_of_operators_gap(
+  L::LieAlgebraStructure, highest_weight::Vector{ZZRingElem}, operators::Vector{GAP.Obj}
+)
+  """
+  used to create action_matrices_of_operators
+  """
+  M = GAP.Globals.HighestWeightModule(L.lie_algebra_gap, GAP.Obj(Int.(highest_weight)))
+  matrices_of_operators = [
+    sparse_matrix(transpose(matrix(QQ, GAP.Globals.MatrixOfAction(GAPWrap.Basis(M), o)))) # TODO: remove transpose?
+    for o in operators
+  ]
+  denominators = map(y -> denominator(y[2]), union(union(matrices_of_operators...)...))
+  common_denominator = lcm(denominators)# // 1
+  matrices_of_operators =
+    (A -> change_base_ring(ZZ, common_denominator * A)).(matrices_of_operators)
+  return matrices_of_operators
+end
+
+function weight(L::LieAlgebraStructure, operator::GAP.Obj)
+  """
+  Calculates the weight in w_i for operator
+  """
+  @req !iszero(operator) "Operators should be non-zero"
+  basis = GAP.Globals.Basis(L.lie_algebra_gap)
+  basis_ind = GAP.Globals.Position(basis, operator)
+  denom = GAP.Globals.Coefficients(basis, operator)[basis_ind]
+  return [
+    ZZ(GAP.Globals.Coefficients(basis, h * operator)[basis_ind]//denom) for
+    h in cartan_sub_basis(L)
+  ]
+end
diff --git a/experimental/BasisLieHighestWeight/src/MainAlgorithm.jl b/experimental/BasisLieHighestWeight/src/MainAlgorithm.jl
new file mode 100644
index 000000000000..08d09091d535
--- /dev/null
+++ b/experimental/BasisLieHighestWeight/src/MainAlgorithm.jl
@@ -0,0 +1,453 @@
+
+function basis_lie_highest_weight_compute(
+  L::LieAlgebraStructure,
+  chevalley_basis::NTuple{3,Vector{GAP.Obj}},
+  highest_weight::Vector{Int},
+  operators::Vector{<:GAP.Obj},     # operators are represented by our monomials. x_i is connected to operators[i]
+  monomial_ordering_symb::Symbol,
+)
+  """
+  Pseudocode:
+
+  basis_lie_highest_weight(highest_weight)
+      return compute_monomials(highest_weight)
+
+  compute_monomials(highest_weight)
+      if highest_weight was already computed 
+          return old results
+      if highest_weight = [0, ..., 0] or [0, ..., 1, ..., 0]
+          return add_by_hand(highest_weight, {})
+      else
+          set_mon = {}
+          go through all partitions lambda_1 + lambda_2 = highest_weight
+              add compute_monomials(lambda_1) (+) compute_monomials(lambda_1) to set_mon 
+          if set_mon too small
+              add_by_hand(highest_weight, set_mon)
+          return set_mon
+
+  add_by_hand(highest_weight, set_mon)
+      add_known_monomials(set_mon)
+      go through all weightspaces that are not full
+          add_new_monomials(weightspace, set_mon)
+      return set_mon
+    
+  add_known_monomials(set_mon)
+      add all monomials from set_mon to basis
+
+  add_new_monomials(weightspace, set_mon)
+      calculate monomials with weight in weightspace
+      go through them one by one in monomial_ordering until basis is full
+      return set_mon
+  """
+  highest_weight = ZZ.(highest_weight)
+  # The function precomputes objects that are independent of the highest weight and that can be used in all recursion 
+  # steps. Then it starts the recursion and returns the result.
+
+  weights_w = [weight(L, op) for op in operators] # weights of the operators
+  weights_alpha = [w_to_alpha(L, weight_w) for weight_w in weights_w] # other root system
+
+  asVec(v) = GAP.gap_to_julia(GAPWrap.ExtRepOfObj(v)) # TODO
+  birational_sequence = BirationalSequence(
+    operators, [asVec(v) for v in operators], weights_w, weights_alpha
+  )
+
+  ZZx, _ = PolynomialRing(ZZ, length(operators)) # for our monomials
+  monomial_ordering = get_monomial_ordering(monomial_ordering_symb, ZZx, weights_alpha)
+
+  # save computations from recursions
+  calc_highest_weight = Dict{Vector{ZZRingElem},Set{ZZMPolyRingElem}}(
+    [ZZ(0) for i in 1:rank(L)] => Set([ZZx(1)])
+  )
+  # save all highest weights, for which the Minkowski-sum did not suffice to gain all monomials
+  no_minkowski = Set{Vector{ZZRingElem}}()
+
+  # start recursion over highest_weight
+  set_mon = compute_monomials(
+    L,
+    birational_sequence,
+    ZZx,
+    highest_weight,
+    monomial_ordering,
+    calc_highest_weight,
+    no_minkowski,
+  )
+
+  # monomials = sort(collect(set_mon); lt=((m1, m2) -> cmp(monomial_ordering, m1, m2) < 0))
+  minkowski_gens = sort(collect(no_minkowski); by=(gen -> (sum(gen), reverse(gen))))
+
+  # output
+  return MonomialBasis(
+    L, highest_weight, monomial_ordering, set_mon, minkowski_gens, birational_sequence
+  )
+end
+
+function compute_monomials(
+  L::LieAlgebraStructure,
+  birational_sequence::BirationalSequence,
+  ZZx::ZZMPolyRing,
+  highest_weight::Vector{ZZRingElem},
+  monomial_ordering::MonomialOrdering,
+  calc_highest_weight::Dict{Vector{ZZRingElem},Set{ZZMPolyRingElem}},
+  no_minkowski::Set{Vector{ZZRingElem}},
+)
+  """
+  This function calculates the monomial basis M_{highest_weight} recursively. The recursion saves all computed 
+  results in calc_highest_weight and we first check, if we already encountered this highest weight in a prior step. 
+  If this is not the case, we need to perform computations. The recursion works by using the Minkowski-sum. 
+  If M_{highest_weight} is the desired set of monomials (identified by the exponents as lattice points), it is know 
+  that for lambda_1 + lambda_2 = highest_weight we have M_{lambda_1} + M_{lambda_2} subseteq M_{highest_weight}. 
+  The complexity grows exponentially in the size of highest_weight. Therefore, it is very helpful to obtain a part of
+  M_{highest_weight} by going through all partitions of highest_weight and using the Minkowski-property. The base 
+  cases of the recursion are the fundamental weights highest_weight = [0, ..., 1, ..., 0]. In this case, or if the 
+  Minkowski-property did not find enough monomials, we need to perform the computations "by hand".
+  """
+  # simple cases
+  # we already computed the highest_weight result in a prior recursion step
+  if haskey(calc_highest_weight, highest_weight)
+    return calc_highest_weight[highest_weight]
+  elseif highest_weight == [ZZ(0) for i in 1:(L.rank)] # we mathematically know the solution
+    return Set(ZZx(1))
+  end
+
+  # calculation required
+  # gap_dim is number of monomials that we need to find, i.e. |M_{highest_weight}|.
+  # if highest_weight is a fundamental weight, partition into smaller summands is possible. This is the basecase of 
+  # the recursion.
+  gap_dim = GAP.Globals.DimensionOfHighestWeightModule(
+    L.lie_algebra_gap, GAP.Obj(Int.(highest_weight))
+  ) # fundamental weights
+  if is_fundamental(highest_weight) || sum(abs.(highest_weight)) == 0
+    push!(no_minkowski, highest_weight)
+    set_mon = add_by_hand(
+      L, birational_sequence, ZZx, highest_weight, monomial_ordering, Set{ZZMPolyRingElem}()
+    )
+    push!(calc_highest_weight, highest_weight => set_mon)
+    return set_mon
+  else
+    # use Minkowski-Sum for recursion
+    set_mon = Set{ZZMPolyRingElem}()
+    i = 0
+    sub_weights_w = compute_sub_weights(highest_weight)
+    l = length(sub_weights_w)
+    # go through all partitions lambda_1 + lambda_2 = highest_weight until we have enough monomials or used all 
+    # partitions
+    while length(set_mon) < gap_dim && i < l
+      i += 1
+      lambda_1 = sub_weights_w[i]
+      lambda_2 = highest_weight .- lambda_1
+      mon_lambda_1 = compute_monomials(
+        L,
+        birational_sequence,
+        ZZx,
+        lambda_1,
+        monomial_ordering,
+        calc_highest_weight,
+        no_minkowski,
+      )
+      mon_lambda_2 = compute_monomials(
+        L,
+        birational_sequence,
+        ZZx,
+        lambda_2,
+        monomial_ordering,
+        calc_highest_weight,
+        no_minkowski,
+      )
+      # Minkowski-sum: M_{lambda_1} + M_{lambda_2} \subseteq M_{highest_weight}, if monomials get identified with 
+      # points in ZZ^n
+      mon_sum = Set([p * q for p in mon_lambda_1 for q in mon_lambda_2])
+      union!(set_mon, mon_sum)
+    end
+    # check if we found enough monomials
+
+    if length(set_mon) < gap_dim
+      push!(no_minkowski, highest_weight)
+      set_mon = add_by_hand(
+        L, birational_sequence, ZZx, highest_weight, monomial_ordering, set_mon
+      )
+    end
+    push!(calc_highest_weight, highest_weight => set_mon)
+    return set_mon
+  end
+end
+
+function add_known_monomials!(
+  weight_w::Vector{ZZRingElem},
+  set_mon_in_weightspace::Dict{Vector{ZZRingElem},Set{ZZMPolyRingElem}},
+  matrices_of_operators::Vector{<:SMat{ZZRingElem}},
+  space::Dict{Vector{ZZRingElem},<:SMat{QQFieldElem}},
+  v0::SRow{ZZRingElem},
+)
+  """
+  By using the Minkowski-sum, we know that all monomials in set_mon_in_weightspace are in our basis. Since we want to
+  extend the weightspace with missing monomials, we need to calculate and add the vector of each monomial to our 
+  basis.
+  """
+  for mon in set_mon_in_weightspace[weight_w]
+    # calculate the vector vec associated with mon
+    vec = calc_vec(v0, mon, matrices_of_operators)
+
+    # check if vec extends the basis
+    if !haskey(space, weight_w)
+      space[weight_w] = sparse_matrix(QQ)
+    end
+    Hecke._add_row_to_rref!(space[weight_w], change_base_ring(QQ, vec))
+  end
+end
+
+function add_new_monomials!(
+  L::LieAlgebraStructure,
+  birational_sequence::BirationalSequence,
+  ZZx::ZZMPolyRing,
+  matrices_of_operators::Vector{<:SMat{ZZRingElem}},
+  monomial_ordering::MonomialOrdering,
+  dim_weightspace::Int,
+  weight_w::Vector{ZZRingElem},
+  set_mon_in_weightspace::Dict{Vector{ZZRingElem},Set{ZZMPolyRingElem}},
+  space::Dict{Vector{ZZRingElem},<:SMat{QQFieldElem}},
+  v0::SRow{ZZRingElem},
+  set_mon::Set{ZZMPolyRingElem},
+)
+  """
+  If a weightspace is missing monomials, we need to calculate them by trial and error. We would like to go through all
+  monomials in the order monomial_ordering and calculate the corresponding vector. If it extends the basis, we add it 
+  to the result and else we try the next one. We know, that all monomials that work lay in the weyl-polytope. 
+  Therefore, we only inspect the monomials that lie both in the weyl-polytope and the weightspace. Since the weyl-
+  polytope is bounded these are finitely many and we can sort them and then go trough them, until we found enough. 
+  """
+  # get monomials that are in the weightspace, sorted by monomial_ordering
+  poss_mon_in_weightspace = convert_lattice_points_to_monomials(
+    ZZx,
+    get_lattice_points_of_weightspace(
+      birational_sequence.weights_alpha, w_to_alpha(L, weight_w)
+    ),
+  )
+  isempty(poss_mon_in_weightspace) && error("The input seems to be invalid.")
+  poss_mon_in_weightspace = sort(
+    poss_mon_in_weightspace; lt=((m1, m2) -> cmp(monomial_ordering, m1, m2) < 0)
+  )
+
+  # check which monomials should get added to the basis
+  i = 0
+  if weight_w == 0 # check if [0 0 ... 0] already in basis
+    i += 1
+  end
+  number_mon_in_weightspace = length(set_mon_in_weightspace[weight_w])
+  # go through possible monomials one by one and check if it extends the basis
+  while number_mon_in_weightspace < dim_weightspace
+    i += 1
+
+    mon = poss_mon_in_weightspace[i]
+    if mon in set_mon
+      continue
+    end
+
+    # calculate the vector vec associated with mon
+    vec = calc_vec(v0, mon, matrices_of_operators)
+
+    # check if vec extends the basis
+    if !haskey(space, weight_w)
+      space[weight_w] = sparse_matrix(QQ)
+    end
+    fl = Hecke._add_row_to_rref!(space[weight_w], change_base_ring(QQ, vec))
+    if !fl
+      continue
+    end
+
+    # save monom
+    number_mon_in_weightspace += 1
+    push!(set_mon, mon)
+  end
+end
+
+function add_by_hand(
+  L::LieAlgebraStructure,
+  birational_sequence::BirationalSequence,
+  ZZx::ZZMPolyRing,
+  highest_weight::Vector{ZZRingElem},
+  monomial_ordering::MonomialOrdering,
+  set_mon::Set{ZZMPolyRingElem},
+)
+  """
+  This function calculates the missing monomials by going through each non full weightspace and adding possible 
+  monomials manually by computing their corresponding vectors and checking if they enlargen the basis.
+  """
+  # initialization
+  # matrices g_i for (g_1^a_1 * ... * g_k^a_k)*v
+  matrices_of_operators = tensor_matrices_of_operators(
+    L, highest_weight, birational_sequence.operators
+  )
+  space = Dict(ZZ(0) * birational_sequence.weights_w[1] => sparse_matrix(QQ)) # span of basis vectors to keep track of the basis
+  v0 = sparse_row(ZZ, [(1, 1)])  # starting vector v
+
+  push!(set_mon, ZZx(1))
+  # required monomials of each weightspace
+  weightspaces = get_dim_weightspace(L, highest_weight)
+
+  # sort the monomials from the minkowski-sum by their weightspaces
+  set_mon_in_weightspace = Dict{Vector{ZZRingElem},Set{ZZMPolyRingElem}}()
+  for (weight_w, _) in weightspaces
+    set_mon_in_weightspace[weight_w] = Set{ZZMPolyRingElem}()
+  end
+  for mon in set_mon
+    weight_w = weight(mon, birational_sequence.weights_w)
+    push!(set_mon_in_weightspace[weight_w], mon)
+  end
+
+  # only inspect weightspaces with missing monomials
+  weights_with_full_weightspace = Set{Vector{ZZRingElem}}()
+  for (weight_w, dim_weightspace) in weightspaces
+    if (length(set_mon_in_weightspace[weight_w]) == dim_weightspace)
+      push!(weights_with_full_weightspace, weight_w)
+    end
+  end
+  delete!(weightspaces, weights_with_full_weightspace)
+
+  # The weightspaces could be calculated completely indepent (except for
+  # the caching). This is not implemented, since I used the package Distributed.jl for this, which is not in the 
+  # Oscar dependencies. But I plan to reimplement this. 
+  # insert known monomials into basis
+  for (weight_w, _) in weightspaces
+    add_known_monomials!(weight_w, set_mon_in_weightspace, matrices_of_operators, space, v0)
+  end
+
+  # calculate new monomials
+  for (weight_w, dim_weightspace) in weightspaces
+    add_new_monomials!(
+      L,
+      birational_sequence,
+      ZZx,
+      matrices_of_operators,
+      monomial_ordering,
+      dim_weightspace,
+      weight_w,
+      set_mon_in_weightspace,
+      space,
+      v0,
+      set_mon,
+    )
+  end
+
+  return set_mon
+end
+
+function operators_by_index(
+  L::LieAlgebraStructure,
+  chevalley_basis::NTuple{3,Vector{GAP.Obj}},
+  birational_sequence::Vector{Int},
+)
+  @req all(i -> 1 <= i <= num_positive_roots(L), birational_sequence) "Entry of birational_sequence out of bounds"
+
+  return [chevalley_basis[1][i] for i in birational_sequence] # TODO: change to [2]
+end
+
+function operators_by_simple_roots(
+  L::LieAlgebraStructure,
+  chevalley_basis::NTuple{3,Vector{GAP.Obj}},
+  birational_sequence::Vector{Vector{Int}},
+)
+  rs = root_system_gap(L)
+  simple_roots = Vector{Vector{Int}}(GAP.Globals.SimpleSystem(rs))
+  positive_roots = Vector{Vector{Int}}(GAP.Globals.PositiveRoots(rs))
+
+  root_inds = Int[]
+  for whgt_alpha in birational_sequence
+    @req length(whgt_alpha) == rank(L) "Length mismatch"
+    @req all(>=(0), whgt_alpha) "Only positive roots are allowed as input"
+    root = sum(whgt_alpha .* simple_roots)
+    root_ind = findfirst(==(root), positive_roots)
+    @req !isnothing(root_ind) "$whgt_alpha is not a positive root"
+    push!(root_inds, root_ind)
+  end
+
+  return operators_by_index(L, chevalley_basis, root_inds)
+end
+
+function operators_lusztig(
+  L::LieAlgebraStructure,
+  chevalley_basis::NTuple{3,Vector{GAP.Obj}},
+  reduced_expression::Vector{Int},
+)
+  root_inds = operators_lusztig_indices(L, reduced_expression)
+  return operators_by_index(L, chevalley_basis, root_inds)
+end
+
+function operators_lusztig_indices(L::LieAlgebraStructure, word::Vector{Int})
+  """
+  Computes the operators for the lusztig polytopes for a longest weyl-word 
+  reduced_expression.
+
+  \beta_k := s_{i_1} … s_{i_{k-1}} (\alpha_{i_k})
+
+  F.e. for A, 2, [1, 2, 1], we get
+  \beta_1 = \alpha_1
+  \beta_2 = \alpha_1 + \alpha_2
+  \beta_3 = \alpha_2
+  """
+  rs = root_system_gap(L)
+
+  simple_roots = GAP.Globals.SimpleSystem(rs)
+  positive_roots = Vector{Vector{Int}}(GAP.Globals.PositiveRoots(rs))
+  sparse_cartan_matrix = GAP.Globals.SparseCartanMatrix(GAP.Globals.WeylGroup(rs))
+
+  root_inds = Int[]
+
+  for k in 1:length(word)
+    # Calculate betas by applying simple reflections step-by-step.
+    root = copy(simple_roots[word[k]])
+    for j in (k - 1):-1:1
+      GAP.Globals.ApplySimpleReflection(sparse_cartan_matrix, word[j], root)
+    end
+    root_ind = findfirst(==(Vector{Int}(root)), positive_roots)
+    @req !isnothing(root_ind) "$root is not a positive root"
+    push!(root_inds, root_ind)
+  end
+  return root_inds
+end
+
+@doc """
+    is_fundamental(highest_weight::Vector{IntegerUnion}) -> Bool
+
+    returns true if ``highest_weight`` is fundamental, i.e. [0, ..., 1, ..., 0]
+
+# Examples
+```jldoctest
+julia> BasisLieHighestWeight.is_fundamental([0, 1, 0])
+true
+
+julia> BasisLieHighestWeight.is_fundamental([0, 1, 1])
+false
+```
+"""
+function is_fundamental(highest_weight::Vector{<:IntegerUnion})
+  hasone = false
+  for i in highest_weight
+    if iszero(i)
+      continue
+    elseif isone(i)
+      hasone && return false
+      hasone = true
+    else
+      return false
+    end
+  end
+  return hasone
+end
+
+function compute_sub_weights(highest_weight::Vector{ZZRingElem})
+  """
+  returns list of weights w != 0, highest_weight with 0 <= w <= highest_weight elementwise, ordered by l_2-norm
+  """
+  sub_weights_w = []
+  foreach(Iterators.product((0:x for x in highest_weight)...)) do i
+    push!(sub_weights_w, [i...])
+  end
+  if isempty(sub_weights_w) || length(sub_weights_w) == 1 # case [] or [[0, ..., 0]]
+    return []
+  else
+    popfirst!(sub_weights_w) # [0, ..., 0]
+    pop!(sub_weights_w) # highest_weight
+    sort!(sub_weights_w; by=x -> sum((x) .^ 2))
+    return sub_weights_w
+  end
+end
diff --git a/experimental/BasisLieHighestWeight/src/MonomialBasis.jl b/experimental/BasisLieHighestWeight/src/MonomialBasis.jl
new file mode 100644
index 000000000000..ec253fac0f57
--- /dev/null
+++ b/experimental/BasisLieHighestWeight/src/MonomialBasis.jl
@@ -0,0 +1,83 @@
+struct MonomialBasis
+  lie_algebra::LieAlgebraStructure
+  # birational_sequence::BirationalSequence
+  highest_weight::Vector{Int}
+  monomial_ordering::MonomialOrdering
+  dimension::Int
+  monomials::Set{ZZMPolyRingElem}
+  monomials_parent::ZZMPolyRing
+  minkowski_gens::Vector{Vector{ZZRingElem}} # TODO: put in attribute storage
+  birational_sequence::BirationalSequence # TODO: put in attribute storage
+
+  function MonomialBasis(
+    lie_algebra::LieAlgebraStructure,
+    highest_weight::Vector{<:IntegerUnion},
+    monomial_ordering::MonomialOrdering,
+    monomials::Set{ZZMPolyRingElem},
+    minkowski_gens::Vector{Vector{ZZRingElem}},
+    birational_sequence::BirationalSequence,
+  )
+    return new(
+      lie_algebra,
+      Int.(highest_weight),
+      monomial_ordering,
+      length(monomials),
+      monomials,
+      parent(first(monomials)),
+      minkowski_gens,
+      birational_sequence,
+    )
+  end
+end
+
+base_lie_algebra(basis::MonomialBasis) = basis.lie_algebra
+
+highest_weight(basis::MonomialBasis) = basis.highest_weight
+
+dim(basis::MonomialBasis) = basis.dimension
+length(basis::MonomialBasis) = dim(basis)
+
+monomials(basis::MonomialBasis) = basis.monomials
+
+monomial_ordering(basis::MonomialBasis) = basis.monomial_ordering
+
+function Base.show(io::IO, ::MIME"text/plain", basis::MonomialBasis)
+  io = pretty(io)
+  println(io, "Monomial basis of a highest weight module")
+  println(io, Indent(), "of highest weight $(highest_weight(basis))", Dedent())
+  println(io, Indent(), "of dimension $(dim(basis))", Dedent())
+  println(io, Indent(), "with monomial ordering $(monomial_ordering(basis))", Dedent())
+  println(io, "over ", Lowercase(), base_lie_algebra(basis))
+  print(
+    io,
+    Indent(),
+    "where the used birational sequence consists of the following roots (given as coefficients w.r.t. alpha_i):",
+    Indent(),
+  )
+  for weight in basis.birational_sequence.weights_alpha
+    print(io, '\n', Int.(weight))
+  end
+  println(io, Dedent(), Dedent())
+  print(
+    io,
+    Indent(),
+    "and the basis was generated by Minkowski sums of the bases of the following highest weight modules:",
+    Indent(),
+  )
+  for gen in basis.minkowski_gens
+    print(io, '\n', Int.(gen))
+  end
+  print(io, Dedent(), Dedent())
+end
+
+function Base.show(io::IO, basis::MonomialBasis)
+  if get(io, :supercompact, false)
+    print(io, "Monomial basis of a highest weight module")
+  else
+    print(
+      io,
+      "Monomial basis of a highest weight module with highest weight $(highest_weight(basis)) over ",
+    )
+    print(IOContext(io, :supercompact => true), base_lie_algebra(basis))
+  end
+end
diff --git a/experimental/BasisLieHighestWeight/src/MonomialOrder.jl b/experimental/BasisLieHighestWeight/src/MonomialOrder.jl
new file mode 100644
index 000000000000..0b4ecb72c63b
--- /dev/null
+++ b/experimental/BasisLieHighestWeight/src/MonomialOrder.jl
@@ -0,0 +1,17 @@
+function get_monomial_ordering(
+  ordering_input::Union{Symbol,Function},
+  ZZx::ZZMPolyRing,
+  weights_alpha::Vector{Vector{QQFieldElem}},
+)
+  """
+  Returns the desired monomial_ordering
+  """
+  if _is_weighted(ordering_input)
+    choosen_monomial_order = monomial_ordering(
+      ZZx, ordering_input, Int[Int(sum(w)) for w in weights_alpha]
+    )
+  else
+    choosen_monomial_order = monomial_ordering(ZZx, ordering_input)
+  end
+  return choosen_monomial_order
+end
diff --git a/experimental/BasisLieHighestWeight/src/NewMonomial.jl b/experimental/BasisLieHighestWeight/src/NewMonomial.jl
new file mode 100644
index 000000000000..e0c81b7344ed
--- /dev/null
+++ b/experimental/BasisLieHighestWeight/src/NewMonomial.jl
@@ -0,0 +1,21 @@
+function weight(mon::ZZMPolyRingElem, weights_w::Vector{Vector{ZZRingElem}})
+  @assert length(weights_w) == length(degrees(mon))
+  return sum(exp * weight for (exp, weight) in zip(degrees(mon), weights_w))
+end
+
+function calc_vec(
+  v0::SRow{ZZRingElem},
+  mon::ZZMPolyRingElem,
+  matrices_of_operators::Vector{<:SMat{ZZRingElem}},
+)
+  v = v0
+  degree_mon = degrees(mon)
+  for i in length(degree_mon):-1:1
+    for _ in 1:degree_mon[i]
+      # currently there is no sparse matrix * vector mult
+      # this is also the line that takes up almost all the computation time for big examples
+      v = v * transpose(matrices_of_operators[i]) # TODO: remove transpose?
+    end
+  end
+  return v
+end
diff --git a/experimental/BasisLieHighestWeight/src/RootConversion.jl b/experimental/BasisLieHighestWeight/src/RootConversion.jl
new file mode 100644
index 000000000000..c1f83d9a0964
--- /dev/null
+++ b/experimental/BasisLieHighestWeight/src/RootConversion.jl
@@ -0,0 +1,9 @@
+function w_to_alpha(
+  L::LieAlgebraStructure, weight_w::Union{Vector{ZZRingElem},Vector{QQFieldElem}}
+)
+  return weight_w * inv_cartan_matrix(L)
+end
+
+function alpha_to_w(L::LieAlgebraStructure, weight_alpha::Vector{QQFieldElem})
+  return weight_alpha * cartan_matrix(L)
+end
diff --git a/experimental/BasisLieHighestWeight/src/TensorModels.jl b/experimental/BasisLieHighestWeight/src/TensorModels.jl
new file mode 100644
index 000000000000..c60fa29e39e6
--- /dev/null
+++ b/experimental/BasisLieHighestWeight/src/TensorModels.jl
@@ -0,0 +1,39 @@
+
+function _tensor_product(A, B)
+  return kronecker_product(A, identity_matrix(SMat, ZZ, nrows(B))) +
+         kronecker_product(identity_matrix(SMat, ZZ, nrows(A)), B)
+end
+
+function _tensor_power(A, k)
+  return sum(
+    kronecker_product(
+      kronecker_product(identity_matrix(SMat, ZZ, nrows(A)^(j - 1)), A),
+      identity_matrix(SMat, ZZ, nrows(A)^(k - j)),
+    ) for j in 1:k
+  )
+end
+
+@doc raw"""
+    tensor_matrices_of_operators(L::LieAlgebraStructure, highest_weight::Vector{ZZRingElem}, operators::Vector{GAP.Obj}) -> Vector{SMat{ZZRingElem}}
+
+Calculates the action matrices of the operators in `operators` on
+the tensor product of the fundamental modules (with multiplicities in `highest_weight`).
+Note that the highest weight module with highest weight `highest_weight` is a submodule of this tensor product.
+"""
+function tensor_matrices_of_operators(
+  L::LieAlgebraStructure, highest_weight::Vector{ZZRingElem}, operators::Vector{GAP.Obj}
+)
+  matrices_of_operators = [zero_matrix(SMat, ZZ, 1) for _ in operators]
+  for (i, highest_weight_i) in enumerate(Int.(highest_weight))
+    if highest_weight_i <= 0
+      continue
+    end
+    wi = ZZ.(1:length(highest_weight) .== i) # i-th fundamental weight
+    matrices_of_operators = [
+      _tensor_product(mat_temp, _tensor_power(mat_wi, highest_weight_i)) for
+      (mat_temp, mat_wi) in
+      zip(matrices_of_operators, matrices_of_operators_gap(L, wi, operators))
+    ]
+  end
+  return matrices_of_operators
+end
diff --git a/experimental/BasisLieHighestWeight/src/UserFunctions.jl b/experimental/BasisLieHighestWeight/src/UserFunctions.jl
new file mode 100644
index 000000000000..a9f097e53cfd
--- /dev/null
+++ b/experimental/BasisLieHighestWeight/src/UserFunctions.jl
@@ -0,0 +1,408 @@
+@doc raw"""
+    basis_lie_highest_weight_operators(type::Symbol, rank::Int)
+
+Lists the operators available for a given simple Lie algebra of type `type_rank`,
+together with their index.
+Operators $f_\alpha$ of negative roots are shown as the coefficients of the corresponding positive root.
+w.r.t. the simple roots $\alpha_i$.
+
+# Example
+```jldoctest
+julia> basis_lie_highest_weight_operators(:B, 2)
+4-element Vector{Tuple{Int64, Vector{QQFieldElem}}}:
+ (1, [1, 0])
+ (2, [0, 1])
+ (3, [1, 1])
+ (4, [1, 2])
+```
+"""
+function basis_lie_highest_weight_operators(type::Symbol, rank::Int)
+  L = lie_algebra(type, rank)
+  chevalley_basis = chevalley_basis_gap(L)
+  operators = chevalley_basis[1]  # TODO: change to [2]
+  weights_alpha = [w_to_alpha(L, weight(L, op)) for op in operators]
+  return collect(enumerate(weights_alpha))
+end
+
+@doc raw"""
+    basis_lie_highest_weight(type::Symbol, rank::Int, highest_weight::Vector{Int}; monomial_ordering::Symbol=:degrevlex)
+    basis_lie_highest_weight(type::Symbol, rank::Int, highest_weight::Vector{Int}, birational_sequence::Vector{Int}; monomial_ordering::Symbol=:degrevlex)
+    basis_lie_highest_weight(type::Symbol, rank::Int, highest_weight::Vector{Int}, birational_sequence::Vector{Vector{Int}}; monomial_ordering::Symbol=:degrevlex)
+
+Computes a monomial basis for the highest weight module with highest weight
+`highest_weight` (in terms of the fundamental weights $\omega_i$),
+for a simple Lie algebra of type `type_rank`.
+
+If no birational sequence is specified, all operators in the order of `basis_lie_highest_weight_operators` are used.
+A birational sequence of type `Vector{Int}` is a sequence of indices of operators in `basis_lie_highest_weight_operators`.
+A birational sequence of type `Vector{Vector{Int}}` is a sequence of weights in terms of the simple roots $\alpha_i$.
+
+`monomial_ordering` describes the monomial ordering used for the basis.
+If this is a weighted ordering, the height of the corresponding root is used as weight.
+
+# Examples
+```jldoctest
+julia> base = basis_lie_highest_weight(:A, 2, [1, 1])
+Monomial basis of a highest weight module
+  of highest weight [1, 1]
+  of dimension 8
+  with monomial ordering degrevlex([x1, x2, x3])
+over Lie algebra of type A2
+  where the used birational sequence consists of the following roots (given as coefficients w.r.t. alpha_i):
+    [1, 0]
+    [0, 1]
+    [1, 1]
+  and the basis was generated by Minkowski sums of the bases of the following highest weight modules:
+    [1, 0]
+    [0, 1]
+
+julia> base = basis_lie_highest_weight(:A, 3, [2, 2, 3]; monomial_ordering = :lex)
+Monomial basis of a highest weight module
+  of highest weight [2, 2, 3]
+  of dimension 1260
+  with monomial ordering lex([x1, x2, x3, x4, x5, x6])
+over Lie algebra of type A3
+  where the used birational sequence consists of the following roots (given as coefficients w.r.t. alpha_i):
+    [1, 0, 0]
+    [0, 1, 0]
+    [0, 0, 1]
+    [1, 1, 0]
+    [0, 1, 1]
+    [1, 1, 1]
+  and the basis was generated by Minkowski sums of the bases of the following highest weight modules:
+    [1, 0, 0]
+    [0, 1, 0]
+    [0, 0, 1]
+
+julia> base = basis_lie_highest_weight(:A, 2, [1, 0], [1,2,1])
+Monomial basis of a highest weight module
+  of highest weight [1, 0]
+  of dimension 3
+  with monomial ordering degrevlex([x1, x2, x3])
+over Lie algebra of type A2
+  where the used birational sequence consists of the following roots (given as coefficients w.r.t. alpha_i):
+    [1, 0]
+    [0, 1]
+    [1, 0]
+  and the basis was generated by Minkowski sums of the bases of the following highest weight modules:
+    [1, 0]
+
+julia> base = basis_lie_highest_weight(:A, 2, [1, 0], [[1,0], [0,1], [1,0]])
+Monomial basis of a highest weight module
+  of highest weight [1, 0]
+  of dimension 3
+  with monomial ordering degrevlex([x1, x2, x3])
+over Lie algebra of type A2
+  where the used birational sequence consists of the following roots (given as coefficients w.r.t. alpha_i):
+    [1, 0]
+    [0, 1]
+    [1, 0]
+  and the basis was generated by Minkowski sums of the bases of the following highest weight modules:
+    [1, 0]
+
+julia> base = basis_lie_highest_weight(:C, 3, [1, 1, 1]; monomial_ordering = :lex)
+Monomial basis of a highest weight module
+  of highest weight [1, 1, 1]
+  of dimension 512
+  with monomial ordering lex([x1, x2, x3, x4, x5, x6, x7, x8, x9])
+over Lie algebra of type C3
+  where the used birational sequence consists of the following roots (given as coefficients w.r.t. alpha_i):
+    [1, 0, 0]
+    [0, 1, 0]
+    [0, 0, 1]
+    [1, 1, 0]
+    [0, 1, 1]
+    [1, 1, 1]
+    [0, 2, 1]
+    [1, 2, 1]
+    [2, 2, 1]
+  and the basis was generated by Minkowski sums of the bases of the following highest weight modules:
+    [1, 0, 0]
+    [0, 1, 0]
+    [0, 0, 1]
+    [0, 1, 1]
+    [1, 1, 1]
+```
+"""
+function basis_lie_highest_weight(
+  type::Symbol, rank::Int, highest_weight::Vector{Int}; monomial_ordering::Symbol=:degrevlex
+)
+  L = lie_algebra(type, rank)
+  chevalley_basis = chevalley_basis_gap(L)
+  operators = chevalley_basis[1] # TODO: change to [2]
+  return basis_lie_highest_weight_compute(
+    L, chevalley_basis, highest_weight, operators, monomial_ordering
+  )
+end
+
+function basis_lie_highest_weight(
+  type::Symbol,
+  rank::Int,
+  highest_weight::Vector{Int},
+  birational_sequence::Vector{Int};
+  monomial_ordering::Symbol=:degrevlex,
+)
+  L = lie_algebra(type, rank)
+  chevalley_basis = chevalley_basis_gap(L)
+  operators = operators_by_index(L, chevalley_basis, birational_sequence)
+  return basis_lie_highest_weight_compute(
+    L, chevalley_basis, highest_weight, operators, monomial_ordering
+  )
+end
+
+function basis_lie_highest_weight(
+  type::Symbol,
+  rank::Int,
+  highest_weight::Vector{Int},
+  birational_sequence::Vector{Vector{Int}};
+  monomial_ordering::Symbol=:degrevlex,
+)
+  L = lie_algebra(type, rank)
+  chevalley_basis = chevalley_basis_gap(L)
+  operators = operators_by_simple_roots(L, chevalley_basis, birational_sequence)
+  return basis_lie_highest_weight_compute(
+    L, chevalley_basis, highest_weight, operators, monomial_ordering
+  )
+end
+
+@doc raw"""
+    basis_lie_highest_weight_lusztig(type::Symbol, rank::Int, highest_weight::Vector{Int}, reduced_expression::Vector{Int})
+
+Computes a monomial basis for the highest weight module with highest weight
+`highest_weight` (in terms of the fundamental weights $\omega_i$),
+for a simple Lie algebra $L$ of type `type_rank`.
+
+Let $\omega_0 = s_{i_1} \cdots s_{i_N}$ be a reduced expression of the longest element in the Weyl group of $L$
+given as indices $[i_1, \dots, i_N]$ in `reduced_expression`.
+Then the birational sequence used consists of $\beta_1, \dots, \beta_N$ where $\beta_1 := \alpha_{i_1}$ and \beta_k := s_{i_1} \cdots s_{i_{k-1}} \alpha_{i_k}$ for $k = 2, \dots, N$.
+
+The monomial ordering is fixed to `wdegrevlex` (weighted degree reverse lexicographic order).
+
+# Examples
+```jldoctest
+julia> base = basis_lie_highest_weight_lusztig(:D, 4, [1,1,1,1], [4,3,2,4,3,2,1,2,4,3,2,1])
+Monomial basis of a highest weight module
+  of highest weight [1, 1, 1, 1]
+  of dimension 4096
+  with monomial ordering wdegrevlex([x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12], [1, 1, 3, 2, 2, 1, 5, 4, 3, 3, 2, 1])
+over Lie algebra of type D4
+  where the used birational sequence consists of the following roots (given as coefficients w.r.t. alpha_i):
+    [0, 0, 0, 1]
+    [0, 0, 1, 0]
+    [0, 1, 1, 1]
+    [0, 1, 1, 0]
+    [0, 1, 0, 1]
+    [0, 1, 0, 0]
+    [1, 2, 1, 1]
+    [1, 1, 1, 1]
+    [1, 1, 0, 1]
+    [1, 1, 1, 0]
+    [1, 1, 0, 0]
+    [1, 0, 0, 0]
+  and the basis was generated by Minkowski sums of the bases of the following highest weight modules:
+    [1, 0, 0, 0]
+    [0, 1, 0, 0]
+    [0, 0, 1, 0]
+    [0, 0, 0, 1]
+    [0, 0, 1, 1]
+```
+"""
+function basis_lie_highest_weight_lusztig(
+  type::Symbol, rank::Int, highest_weight::Vector{Int}, reduced_expression::Vector{Int}
+)
+  monomial_ordering = :wdegrevlex
+  L = lie_algebra(type, rank)
+  chevalley_basis = chevalley_basis_gap(L)
+  operators = operators_lusztig(L, chevalley_basis, reduced_expression)
+  return basis_lie_highest_weight_compute(
+    L, chevalley_basis, highest_weight, operators, monomial_ordering
+  )
+end
+
+@doc raw"""
+    basis_lie_highest_weight_string(type::Symbol, rank::Int, highest_weight::Vector{Int}, reduced_expression::Vector{Int})
+
+Computes a monomial basis for the highest weight module with highest weight
+`highest_weight` (in terms of the fundamental weights $\omega_i$),
+for a simple Lie algebra $L$ of type `type_rank`.
+
+Let $\omega_0 = s_{i_1} \cdots s_{i_N}$ be a reduced expression of the longest element in the Weyl group of $L$
+given as indices $[i_1, \dots, i_N]$ in `reduced_expression`.
+Then the birational sequence used consists of $\alpha_{i_1}, \dots, \alpha_{i_N}$.
+
+The monomial ordering is fixed to `neglex` (negative lexicographic order).      
+
+# Examples
+```jldoctest
+julia> basis_lie_highest_weight_string(:B, 3, [1,1,1], [3,2,3,2,1,2,3,2,1])
+Monomial basis of a highest weight module
+  of highest weight [1, 1, 1]
+  of dimension 512
+  with monomial ordering neglex([x1, x2, x3, x4, x5, x6, x7, x8, x9])
+over Lie algebra of type B3
+  where the used birational sequence consists of the following roots (given as coefficients w.r.t. alpha_i):
+    [0, 0, 1]
+    [0, 1, 0]
+    [0, 0, 1]
+    [0, 1, 0]
+    [1, 0, 0]
+    [0, 1, 0]
+    [0, 0, 1]
+    [0, 1, 0]
+    [1, 0, 0]
+  and the basis was generated by Minkowski sums of the bases of the following highest weight modules:
+    [1, 0, 0]
+    [0, 1, 0]
+    [0, 0, 1]
+
+julia> basis_lie_highest_weight_string(:A, 4, [1,1,1,1], [4,3,2,1,2,3,4,3,2,3])
+Monomial basis of a highest weight module
+  of highest weight [1, 1, 1, 1]
+  of dimension 1024
+  with monomial ordering neglex([x1, x2, x3, x4, x5, x6, x7, x8, x9, x10])
+over Lie algebra of type A4
+  where the used birational sequence consists of the following roots (given as coefficients w.r.t. alpha_i):
+    [0, 0, 0, 1]
+    [0, 0, 1, 0]
+    [0, 1, 0, 0]
+    [1, 0, 0, 0]
+    [0, 1, 0, 0]
+    [0, 0, 1, 0]
+    [0, 0, 0, 1]
+    [0, 0, 1, 0]
+    [0, 1, 0, 0]
+    [0, 0, 1, 0]
+  and the basis was generated by Minkowski sums of the bases of the following highest weight modules:
+    [1, 0, 0, 0]
+    [0, 1, 0, 0]
+    [0, 0, 1, 0]
+    [0, 0, 0, 1]
+    [0, 1, 0, 1]
+```
+"""
+function basis_lie_highest_weight_string(
+  type::Symbol, rank::Int, highest_weight::Vector{Int}, reduced_expression::Vector{Int}
+)
+  monomial_ordering = :neglex
+  L = lie_algebra(type, rank)
+  chevalley_basis = chevalley_basis_gap(L)
+  operators = operators_by_index(L, chevalley_basis, reduced_expression)
+  return basis_lie_highest_weight_compute(
+    L, chevalley_basis, highest_weight, operators, monomial_ordering
+  )
+end
+
+@doc raw"""
+    basis_lie_highest_weight_pbw(type::Symbol, rank::Int, highest_weight::Vector{Int})
+
+Computes a monomial basis for the highest weight module with highest weight
+`highest_weight` (in terms of the fundamental weights $\omega_i$),
+for a simple Lie algebra $L$ of type `type_rank`.
+
+Then the birational sequence used consists of all operators in descening height of the corresponding root.
+
+The monomial ordering is fixed to `neglex` (negative lexicographic order).      
+      
+# Examples
+```jldoctest
+julia> basis_lie_highest_weight_pbw(:A, 3, [1,1,1])
+Monomial basis of a highest weight module
+  of highest weight [1, 1, 1]
+  of dimension 64
+  with monomial ordering neglex([x1, x2, x3, x4, x5, x6])
+over Lie algebra of type A3
+  where the used birational sequence consists of the following roots (given as coefficients w.r.t. alpha_i):
+    [1, 1, 1]
+    [0, 1, 1]
+    [1, 1, 0]
+    [0, 0, 1]
+    [0, 1, 0]
+    [1, 0, 0]
+  and the basis was generated by Minkowski sums of the bases of the following highest weight modules:
+    [1, 0, 0]
+    [0, 1, 0]
+    [0, 0, 1]
+```
+"""
+function basis_lie_highest_weight_pbw(type::Symbol, rank::Int, highest_weight::Vector{Int})
+  monomial_ordering = :neglex
+  L = lie_algebra(type, rank)
+  chevalley_basis = chevalley_basis_gap(L)
+  operators = reverse(chevalley_basis[1]) # TODO: change to [2]
+  return basis_lie_highest_weight_compute(
+    L, chevalley_basis, highest_weight, operators, monomial_ordering
+  )
+end
+
+@doc raw"""
+    basis_lie_highest_weight_nz(type::Symbol, rank::Int, highest_weight::Vector{Int}, reduced_expression::Vector{Int})
+
+Computes a monomial basis for the highest weight module with highest weight
+`highest_weight` (in terms of the fundamental weights $\omega_i$),
+for a simple Lie algebra $L$ of type `type_rank`.
+
+Let $\omega_0 = s_{i_1} \cdots s_{i_N}$ be a reduced expression of the longest element in the Weyl group of $L$
+given as indices $[i_1, \dots, i_N]$ in `reduced_expression`.
+Then the birational sequence used consists of $\alpha_{i_1}, \dots, \alpha_{i_N}$.
+
+The monomial ordering is fixed to `degrevlex` (degree reverse lexicographic order).      
+
+# Examples
+```jldoctest
+julia> basis_lie_highest_weight_nz(:C, 3, [1,1,1], [3,2,3,2,1,2,3,2,1])
+Monomial basis of a highest weight module
+  of highest weight [1, 1, 1]
+  of dimension 512
+  with monomial ordering degrevlex([x1, x2, x3, x4, x5, x6, x7, x8, x9])
+over Lie algebra of type C3
+  where the used birational sequence consists of the following roots (given as coefficients w.r.t. alpha_i):
+    [0, 0, 1]
+    [0, 1, 0]
+    [0, 0, 1]
+    [0, 1, 0]
+    [1, 0, 0]
+    [0, 1, 0]
+    [0, 0, 1]
+    [0, 1, 0]
+    [1, 0, 0]
+  and the basis was generated by Minkowski sums of the bases of the following highest weight modules:
+    [1, 0, 0]
+    [0, 1, 0]
+    [0, 0, 1]
+
+julia> basis_lie_highest_weight_nz(:A, 4, [1,1,1,1], [4,3,2,1,2,3,4,3,2,3])
+Monomial basis of a highest weight module
+  of highest weight [1, 1, 1, 1]
+  of dimension 1024
+  with monomial ordering degrevlex([x1, x2, x3, x4, x5, x6, x7, x8, x9, x10])
+over Lie algebra of type A4
+  where the used birational sequence consists of the following roots (given as coefficients w.r.t. alpha_i):
+    [0, 0, 0, 1]
+    [0, 0, 1, 0]
+    [0, 1, 0, 0]
+    [1, 0, 0, 0]
+    [0, 1, 0, 0]
+    [0, 0, 1, 0]
+    [0, 0, 0, 1]
+    [0, 0, 1, 0]
+    [0, 1, 0, 0]
+    [0, 0, 1, 0]
+  and the basis was generated by Minkowski sums of the bases of the following highest weight modules:
+    [1, 0, 0, 0]
+    [0, 1, 0, 0]
+    [0, 0, 1, 0]
+    [0, 0, 0, 1]
+    [0, 1, 0, 1] 
+```
+"""
+function basis_lie_highest_weight_nz(
+  type::Symbol, rank::Int, highest_weight::Vector{Int}, reduced_expression::Vector{Int}
+)
+  monomial_ordering = :degrevlex
+  L = lie_algebra(type, rank)
+  chevalley_basis = chevalley_basis_gap(L)
+  operators = operators_by_index(L, chevalley_basis, reduced_expression)
+  return basis_lie_highest_weight_compute(
+    L, chevalley_basis, highest_weight, operators, monomial_ordering
+  )
+end
diff --git a/experimental/BasisLieHighestWeight/src/WeylPolytope.jl b/experimental/BasisLieHighestWeight/src/WeylPolytope.jl
new file mode 100644
index 000000000000..8650b8a3e43e
--- /dev/null
+++ b/experimental/BasisLieHighestWeight/src/WeylPolytope.jl
@@ -0,0 +1,90 @@
+
+@doc raw"""
+    orbit_weylgroup(L::LieAlgebraStructure, weight_w::Vector{ZZRingElem}) -> Vector{Vector{ZZRingElem}}
+
+Computes the orbit of the weight `weight_w` given as coefficients to the fundamental weights
+$\omega_i$ under the action of the Weyl group of the Lie algebra `L`.
+"""
+function orbit_weylgroup(L::LieAlgebraStructure, weight_w::Vector{ZZRingElem})
+  # initialization
+  weyl_group = GAP.Globals.WeylGroup(GAP.Globals.RootSystem(L.lie_algebra_gap))
+  orbit_iterator = GAP.Globals.WeylOrbitIterator(weyl_group, GAP.Obj(Int.(weight_w)))
+  vertices = Vector{ZZRingElem}[]
+
+  # operate with the weylgroup on weight_vector
+  while !(GAPWrap.IsDoneIterator(orbit_iterator))
+    w = GAPWrap.NextIterator(orbit_iterator)
+    push!(vertices, Vector{ZZRingElem}(w))
+  end
+
+  return vertices
+end
+
+@doc raw"""
+    get_dim_weightspace(L::LieAlgebraStructure, highest_weight::Vector{ZZRingElem}) -> Dict{Vector{ZZRingElem},Int}
+
+Computes the dimension of the weight spaces of the Lie algebra `L` module with highest weight `highest_weight`.
+For all dominant weights, the dimension is computed with GAP. For the remaining weights, the dimension is
+calculated by checking which dominant weight lies in the orbit of the weight under the action of the Weyl group.
+
+The weights are given as coefficients to the fundamental weights $\omega_i$.
+"""
+function get_dim_weightspace(L::LieAlgebraStructure, highest_weight::Vector{ZZRingElem})
+  # calculate dimension for dominant weights with GAP
+  root_system = root_system_gap(L)
+  dominant_char = GAP.Globals.DominantCharacter(root_system, GAP.Obj(Int.(highest_weight)))
+  dominant_weights = map(weight -> ZZ.(weight), dominant_char[1])
+  dominant_weights_dim = Int.(dominant_char[2])
+  weightspaces = Dict{Vector{ZZRingElem},Int}()
+
+  # calculate dimension for the rest by checking which dominant weight lies in the orbit
+  for (dominant_weight, dim) in zip(dominant_weights, dominant_weights_dim)
+    for weight in orbit_weylgroup(L, dominant_weight)
+      weightspaces[highest_weight - weight] = dim
+    end
+  end
+  return weightspaces
+end
+
+function convert_lattice_points_to_monomials(
+  ZZx::ZZMPolyRing, lattice_points_weightspace::Vector{Vector{ZZRingElem}}
+)
+  return [ZZx([ZZ(1)], [lattice_point]) for lattice_point in lattice_points_weightspace]
+end
+
+@doc raw"""
+    get_lattice_points_of_weightspace(weights::Vector{Vector{QQFieldElem}}, wght::Vector{QQFieldElem}) -> Vector{Vector{ZZRingElem}}
+
+Calculates all lattice points in a given weightspace for a Lie algebra highest weight module.
+This is equivalent to finding $\mathbb{Z}$-linear combinations of `weights` that equal `wght`.
+All weights are given as coefficients to the simple roots $\alpha_i$.
+"""
+function get_lattice_points_of_weightspace(
+  weights::Vector{Vector{QQFieldElem}}, wght::Vector{QQFieldElem}
+)
+  # calculate all integer solutions to the following linear program:
+  # [   |              |    ]       [   x   ]      
+  # [weights[1]...weights[k]]   *   [   |   ]   =   weight 
+  # [   |              |    ]       [  res  ] 
+  # [   |              |    ]       [   |   ]
+  # where res[i] >= 0 for all i
+
+  n = length(weights)
+  m = length(wght)
+  A = zero_matrix(QQ, 2m + n, n)
+  b = [zero(QQ) for _ in 1:(2m + n)]
+  # equalities
+  for i in 1:n
+    w = matrix(QQ, m, 1, weights[i])
+    A[1:m, i] = w
+    A[(m + 1):(2m), i] = -w
+  end
+  b[1:m] = wght
+  b[(m + 1):(2m)] = -wght
+  # non-negativity
+  for i in 1:n
+    A[2m + i, i] = -1
+  end
+
+  return Vector{ZZRingElem}.(lattice_points(polyhedron(A, b)))
+end
diff --git a/experimental/BasisLieHighestWeight/test/MBOld.jl b/experimental/BasisLieHighestWeight/test/MBOld.jl
new file mode 100644
index 000000000000..ebb31573be7d
--- /dev/null
+++ b/experimental/BasisLieHighestWeight/test/MBOld.jl
@@ -0,0 +1,261 @@
+# This file is an old version of the algorithm that can compute (not all cases) of 
+# basis_lie_highest_weight and is used only in runtests.jl to check that the newer algorithm matches
+# There is code doubling, but I am not sure how the src part is going to change when its integrated with the other
+# lie algebra work and would prefer to change this file after doing the integration to make sure that everything stays
+# correct.
+
+module MBOld
+
+using Oscar
+
+struct SparseVectorSpaceBasis
+  A::Vector{SRow{ZZRingElem}}
+  pivot::Vector{Int}
+end
+
+function normalize(v::SRow{ZZRingElem})
+  """
+  divides vector by gcd of nonzero entries, returns vector and first nonzero index
+  used: addAndReduce!
+  """
+  if isempty(v)
+    return (0, 0)
+  end
+
+  pivot = first(v)[1]
+
+  return divexact(v, gcd(map(y -> y[2], union(v)))), pivot
+end
+
+reduceCol(a, b, i::Int) = b[i] * a - a[i] * b
+
+function addAndReduce!(sp::SparseVectorSpaceBasis, v::SRow{ZZRingElem})
+  """
+  for each pivot of sp.A we make entry of v zero and return the result
+  0 => linear dependent
+  * => linear independent, new column element of sp.A since it increases basis
+  invariants: the row of a pivotelement in any column in A is 0 (except the pivotelement)
+             elements of A are integers, gcd of each column is 1
+  """
+  A = sp.A
+  pivot = sp.pivot
+  v, newPivot = normalize(v)
+  if newPivot == 0
+    return 0
+  end
+
+  for j in 1:length(A)
+    i = pivot[j]
+    if i != newPivot
+      continue
+    end
+    v = reduceCol(v, A[j], i)
+    v, newPivot = normalize(v)
+    if newPivot == 0
+      return 0
+    end
+  end
+
+  pos = findfirst(pivot .> newPivot)
+  if (pos === nothing)
+    pos = length(pivot) + 1
+  end
+
+  insert!(A, pos, v)
+  insert!(pivot, pos, newPivot)
+
+  return v
+end
+
+#### Lie algebras
+
+function lieAlgebra(t::String, n::Int)
+  L = GAP.Globals.SimpleLieAlgebra(GAP.Obj(t), n, GAP.Globals.Rationals)
+  return L, NTuple{3,Vector{GAP.Obj}}(GAP.Globals.ChevalleyBasis(L))
+end
+
+function matricesForOperators(L, hw, ops)
+  """
+  used to create tensorMatricesForOperators
+  """
+  M = GAP.Globals.HighestWeightModule(L, GAP.Obj(hw))
+  mats = map(
+    o -> sparse_matrix(
+      transpose(matrix(QQ, GAP.Globals.MatrixOfAction(GAP.Globals.Basis(M), o)))
+    ),
+    ops,
+  )
+  denominators = map(y -> denominator(y[2]), union(union(mats...)...))
+  d = lcm(denominators)# // 1
+  mats = (A -> change_base_ring(ZZ, d * A)).(mats)
+  return mats
+end
+
+function weightsForOperators(L, cartan, ops)
+  asVec(v) = Vector{Int}(GAP.Globals.Coefficients(GAP.Globals.Basis(L), v))
+  if any(iszero.(asVec.(ops)))
+    error("ops should be non-zero")
+  end
+  nzi(v) = findfirst(asVec(v) .!= 0)
+  return [[asVec(h * v)[nzi(v)] / asVec(v)[nzi(v)] for h in cartan] for v in ops]
+end
+
+#### tensor model
+
+function kron(A, B)
+  res = sparse_matrix(ZZ, nrows(A) * nrows(B), ncols(A) * ncols(B))
+  for i in 1:nrows(B)
+    for j in 1:nrows(A)
+      new_row_tuples = Vector{Tuple{Int,ZZRingElem}}([(1, ZZ(0))])
+      for (index_A, element_A) in union(getindex(A, j))
+        for (index_B, element_B) in union(getindex(B, i))
+          push!(new_row_tuples, ((index_A - 1) * ncols(B) + index_B, element_A * element_B))
+        end
+      end
+      new_row = sparse_row(ZZ, new_row_tuples)
+      setindex!(res, new_row, (j - 1) * nrows(B) + i)
+    end
+  end
+  return res
+end
+
+# temprary fix sparse in Oscar does not work
+function tensorProduct(A, B)
+  temp_mat = kron(A, spid(nrows(B))) + kron(spid(nrows(A)), B)
+  res = sparse_matrix(ZZ, nrows(A) * nrows(B), ncols(A) * ncols(B))
+  for i in 1:nrows(temp_mat)
+    setindex!(res, getindex(temp_mat, i), i)
+  end
+  return res
+end
+
+spid(n) = identity_matrix(SMat, ZZ, n)
+tensorProducts(As, Bs) = (AB -> tensorProduct(AB[1], AB[2])).(zip(As, Bs))
+tensorPower(A, n) = (n == 1) ? A : tensorProduct(tensorPower(A, n - 1), A)
+tensorPowers(As, n) = (A -> tensorPower(A, n)).(As)
+
+function tensorMatricesForOperators(L, hw, ops)
+  """
+  Calculates the matrices g_i corresponding to the operator ops[i].
+  """
+  mats = []
+
+  for i in 1:length(hw)
+    if hw[i] <= 0
+      continue
+    end
+    wi = Int.(1:length(hw) .== i) # i-th fundamental weight
+    _mats = matricesForOperators(L, wi, ops)
+    _mats = tensorPowers(_mats, hw[i])
+    if size(mats)[1] > 0
+      A = _mats[1]
+      B = mats[1]
+    end
+    mats = mats == [] ? _mats : tensorProducts(mats, _mats)
+  end
+  return mats
+end
+
+"""
+    basisLieHighestWeight(t::String, n::Int, hw::Vector{Int}; parallel::Bool = true) :: Tuple{Vector{Vector{UInt8}},Vector{SRow{ZZRingElem}}}
+
+Compute a monomial basis for the highest weight module with highest weight ``hw`` (in terms of the fundamental weights), for a simple Lie algebra of type ``t`` and rank ``n``.
+
+Example
+======
+```jldoctest
+julia> dim, monomials, vectors = PolyBases.MB.basisLieHighestWeight(:A, 2, [1,0])
+(3, Vector{UInt8}[[0x00, 0x00, 0x00], [0x01, 0x00, 0x00], [0x00, 0x00, 0x01]], SparseArrays.SparseVector{Int64, Int64}[  [1]  =  1,   [2]  =  1,   [3]  =  1])
+```
+"""
+
+function basisLieHighestWeight(t::String, n::Int, hw::Vector{Int}; roots=[]) #--- :: Tuple{Int64,Vector{Vector{UInt8}},Vector{SRow{ZZRingElem}}}
+  L, CH = lieAlgebra(t, n)
+  ops = CH[1] # positive root vectors
+  # .. reorder..
+  wts = weightsForOperators(L, CH[3], ops)
+  wts = (v -> Int.(v)).(wts)
+
+  mats = tensorMatricesForOperators(L, hw, ops)
+  hwv = sparse_row(ZZ, [(1, 1)])
+
+  monomials = compute(hwv, mats, wts)
+  ZZx, x = PolynomialRing(ZZ, length(monomials[1]))
+  monomials = [
+    finish(push_term!(MPolyBuildCtx(ZZx), ZZ(1), Int.(mon))) for mon in monomials
+  ]
+  monomials = Set(monomials)
+  return monomials
+end
+
+nullMon(m) = zeros(UInt8, m)
+
+function compute(v0, mats, wts::Vector{Vector{Int}})
+  m = length(mats)
+  monomials = [nullMon(m)]
+  lastPos = 0
+  id(mon) = sum((1 << (sum(mon[1:i]) + i - 1) for i in 1:(m - 1)); init=1)
+  e = [UInt8.(1:m .== i) for i in 1:m]
+  maxid(deg) = id(deg .* e[1])
+
+  blacklists = [falses(maxid(0))]
+  blacklists[end][id(monomials[1])] = 1
+  newMons(deg) = (push!(blacklists, falses(maxid(deg))))
+  checkMon(mon) = blacklists[end - 1][id(mon)]
+  setMon(mon) = (blacklists[end][id(mon)] = true)
+
+  vectors = [v0]
+  weights = [0 * wts[1]]
+  space = Dict(weights[1] => SparseVectorSpaceBasis([], []))
+  addAndReduce!(space[weights[1]], v0)
+
+  deg = 0
+  while length(monomials) > lastPos
+    startPos = lastPos + 1
+    newPos = length(monomials)
+    deg = deg + 1
+    newMons(deg)
+    for i in 1:m, di in deg:-1:1
+      for p in startPos:newPos
+        if !all(monomials[p][1:(i - 1)] .== 0)
+          continue
+        end
+
+        if monomials[p][i] + 1 > di
+          startPos = p + 1
+          continue
+        end
+        if monomials[p][i] + 1 < di
+          break
+        end
+
+        mon = monomials[p] + e[i]
+
+        if any(i != j && mon[j] > 0 && !checkMon(mon - e[j]) for j in 1:m)
+          continue
+        end
+
+        wt = weights[p] + wts[i]
+        if !haskey(space, wt)
+          space[wt] = SparseVectorSpaceBasis([], [])
+        end
+
+        vec = vectors[p] * transpose(mats[i])
+        vec = addAndReduce!(space[wt], vec)
+        if vec == 0
+          continue
+        end
+
+        setMon(mon)
+        push!(monomials, mon)
+        push!(weights, wt)
+        push!(vectors, vec)
+      end
+    end
+    lastPos = newPos
+  end
+
+  return monomials
+end
+
+end
diff --git a/experimental/BasisLieHighestWeight/test/MainAlgorithm-test.jl b/experimental/BasisLieHighestWeight/test/MainAlgorithm-test.jl
new file mode 100644
index 000000000000..a3a8959bb35a
--- /dev/null
+++ b/experimental/BasisLieHighestWeight/test/MainAlgorithm-test.jl
@@ -0,0 +1,566 @@
+include("MBOld.jl")
+
+"""
+We are testing our code in multiple ways. First, we calculated two small examples per hand and compare those. Then we 
+check basic properties of the result. For example we know the size of our monomial basis. These properties get partially
+used in the algorithm and could therefore be true for false results. We have another basic algorithm that solves the 
+problem without the recursion, weightspaces and saving of computations. The third test compares the results we can 
+compute with the weaker version.
+"""
+
+function compare_algorithms(dynkin::Symbol, n::Int64, lambda::Vector{Int64})
+  # old algorithm
+  mons_old = MBOld.basisLieHighestWeight(string(dynkin), n, lambda) # basic algorithm
+
+  # new algorithm
+  basis = basis_lie_highest_weight(dynkin, n, lambda)
+  mons_new = monomials(basis)
+  L = GAP.Globals.SimpleLieAlgebra(GAP.Obj(dynkin), n, GAP.Globals.Rationals)
+  gap_dim = GAP.Globals.DimensionOfHighestWeightModule(L, GAP.Obj(lambda)) # dimension
+
+  # comparison
+  # convert set of monomials over different ring objects to string representation to compare for equality
+  @test issetequal(string.(mons_old), string.(mons_new)) # compare if result of old and new algorithm match
+  @test gap_dim == length(mons_new) # check if dimension is correct
+end
+
+function check_dimension(
+  dynkin::Symbol, n::Int64, lambda::Vector{Int64}, monomial_ordering::Symbol
+)
+  basis = basis_lie_highest_weight(dynkin, n, lambda; monomial_ordering)
+  L = GAP.Globals.SimpleLieAlgebra(GAP.Obj(dynkin), n, GAP.Globals.Rationals)
+  gap_dim = GAP.Globals.DimensionOfHighestWeightModule(L, GAP.Obj(lambda)) # dimension
+  @test gap_dim == dim(basis) == length(monomials(basis)) # check if dimension is correct
+end
+
+@testset "Test BasisLieHighestWeight" begin
+  @testset "is_fundamental" begin
+    @test BasisLieHighestWeight.is_fundamental([ZZ(0), ZZ(1), ZZ(0)])
+    @test !BasisLieHighestWeight.is_fundamental([ZZ(0), ZZ(1), ZZ(1)])
+  end
+
+  @testset "compute_sub_weights" begin
+    @test isequal(BasisLieHighestWeight.compute_sub_weights([ZZ(0), ZZ(0), ZZ(0)]), [])
+    sub_weights = Vector{Vector{ZZRingElem}}([
+      [1, 0, 0],
+      [0, 1, 0],
+      [0, 0, 1],
+      [1, 1, 0],
+      [1, 0, 1],
+      [0, 1, 1],
+      [1, 1, 1],
+      [2, 0, 0],
+      [0, 2, 0],
+      [2, 1, 0],
+      [1, 2, 0],
+      [2, 0, 1],
+      [0, 2, 1],
+      [2, 1, 1],
+      [1, 2, 1],
+      [2, 2, 0],
+      [0, 3, 0],
+      [2, 2, 1],
+      [1, 3, 0],
+      [0, 3, 1],
+      [1, 3, 1],
+      [2, 3, 0],
+    ])
+    @test isequal(
+      BasisLieHighestWeight.compute_sub_weights([ZZ(2), ZZ(3), ZZ(1)]), sub_weights
+    )
+  end
+
+  @testset "Known examples basis_lie_highest_weight" begin
+    base = basis_lie_highest_weight(:A, 2, [1, 0])
+    mons = monomials(base)
+    @test issetequal(string.(mons), Set(["1", "x3", "x1"]))
+    base = basis_lie_highest_weight(:A, 2, [1, 0], [1, 2, 1])
+    mons = monomials(base)
+    @test issetequal(string.(mons), Set(["1", "x2*x3", "x3"]))
+  end
+
+  @testset "Compare basis_lie_highest_weight with algorithm of Johannes and check dimension" begin
+    @testset "Dynkin type $dynkin" for dynkin in (:A, :B, :C, :D)
+      @testset "n = $n" for n in 1:4
+        if (
+          !(dynkin == :B && n < 2) && !(dynkin == :C && n < 2) && !(dynkin == :D && n < 4)
+        )
+          for i in 1:n                                # w_i
+            lambda = zeros(Int64, n)
+            lambda[i] = 1
+            compare_algorithms(dynkin, n, lambda)
+          end
+
+          if (n > 1)
+            lambda = [1, (0 for i in 1:(n - 2))..., 1]  # w_1 + w_n
+            compare_algorithms(dynkin, n, lambda)
+          end
+
+          if (n < 4)
+            lambda = ones(Int64, n)                  # w_1 + ... + w_n
+            compare_algorithms(dynkin, n, lambda)
+          end
+        end
+      end
+    end
+  end
+
+  @testset "Compare against GAP algorithm of Xin on some examples" begin
+    basis_lusztig = basis_lie_highest_weight_lusztig(:A, 3, [2, 1, 1], [2, 3, 1, 2, 3, 1])
+
+    @test issetequal(
+      [only(exponents(m)) for m in monomials(basis_lusztig)],
+      [
+        [0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 1, 0],
+        [1, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 1],
+        [1, 0, 0, 0, 1, 0],
+        [0, 0, 1, 0, 0, 0],
+        [1, 0, 0, 0, 0, 1],
+        [0, 1, 0, 0, 0, 0],
+        [0, 0, 0, 0, 2, 0],
+        [0, 0, 0, 0, 1, 1],
+        [1, 0, 0, 0, 1, 1],
+        [0, 0, 1, 0, 0, 1],
+        [0, 1, 0, 0, 1, 0],
+        [0, 0, 0, 1, 0, 0],
+        [1, 0, 0, 0, 2, 0],
+        [0, 0, 1, 0, 1, 0],
+        [2, 0, 0, 0, 1, 0],
+        [2, 0, 0, 0, 0, 1],
+        [0, 1, 0, 0, 0, 1],
+        [0, 0, 0, 0, 2, 1],
+        [1, 0, 0, 0, 2, 1],
+        [0, 0, 1, 0, 1, 1],
+        [0, 1, 0, 0, 2, 0],
+        [0, 0, 0, 1, 1, 0],
+        [2, 0, 0, 0, 1, 1],
+        [1, 0, 1, 0, 0, 1],
+        [1, 1, 0, 0, 1, 0],
+        [1, 0, 0, 1, 0, 0],
+        [0, 1, 0, 0, 1, 1],
+        [0, 0, 0, 1, 0, 1],
+        [2, 0, 0, 0, 2, 0],
+        [1, 0, 1, 0, 1, 0],
+        [1, 1, 0, 0, 0, 1],
+        [0, 0, 1, 0, 2, 0],
+        [2, 0, 0, 0, 2, 1],
+        [1, 0, 1, 0, 1, 1],
+        [0, 0, 2, 0, 0, 1],
+        [1, 1, 0, 0, 2, 0],
+        [1, 0, 0, 1, 1, 0],
+        [0, 1, 1, 0, 1, 0],
+        [1, 1, 0, 0, 1, 1],
+        [1, 0, 0, 1, 0, 1],
+        [0, 1, 1, 0, 0, 1],
+        [0, 2, 0, 0, 1, 0],
+        [0, 0, 1, 0, 2, 1],
+        [0, 0, 0, 1, 2, 0],
+        [0, 1, 0, 0, 2, 1],
+        [0, 0, 0, 1, 1, 1],
+        [1, 0, 1, 0, 2, 0],
+        [3, 0, 0, 0, 2, 0],
+        [3, 0, 0, 0, 1, 1],
+        [1, 1, 0, 0, 2, 1],
+        [1, 0, 0, 1, 1, 1],
+        [0, 1, 1, 0, 1, 1],
+        [0, 0, 1, 1, 0, 1],
+        [0, 2, 0, 0, 2, 0],
+        [0, 1, 0, 1, 1, 0],
+        [1, 0, 1, 0, 2, 1],
+        [0, 0, 2, 0, 1, 1],
+        [1, 0, 0, 1, 2, 0],
+        [0, 1, 1, 0, 2, 0],
+        [3, 0, 0, 0, 2, 1],
+        [2, 0, 1, 0, 1, 1],
+        [2, 1, 0, 0, 2, 0],
+        [2, 0, 0, 1, 1, 0],
+        [2, 1, 0, 0, 1, 1],
+        [2, 0, 0, 1, 0, 1],
+        [0, 2, 0, 0, 1, 1],
+        [0, 0, 0, 1, 2, 1],
+        [2, 0, 1, 0, 2, 0],
+        [1, 0, 0, 1, 2, 1],
+        [0, 1, 1, 0, 2, 1],
+        [0, 0, 1, 1, 1, 1],
+        [0, 1, 0, 1, 2, 0],
+        [2, 1, 0, 0, 2, 1],
+        [2, 0, 0, 1, 1, 1],
+        [1, 1, 1, 0, 1, 1],
+        [1, 0, 1, 1, 0, 1],
+        [1, 2, 0, 0, 2, 0],
+        [1, 1, 0, 1, 1, 0],
+        [0, 2, 0, 0, 2, 1],
+        [0, 1, 0, 1, 1, 1],
+        [2, 0, 1, 0, 2, 1],
+        [1, 0, 2, 0, 1, 1],
+        [2, 0, 0, 1, 2, 0],
+        [1, 1, 1, 0, 2, 0],
+        [1, 2, 0, 0, 1, 1],
+        [0, 0, 2, 0, 2, 1],
+        [4, 0, 0, 0, 2, 1],
+        [2, 0, 0, 1, 2, 1],
+        [1, 1, 1, 0, 2, 1],
+        [1, 0, 1, 1, 1, 1],
+        [0, 1, 2, 0, 1, 1],
+        [1, 1, 0, 1, 2, 0],
+        [0, 2, 1, 0, 2, 0],
+        [1, 2, 0, 0, 2, 1],
+        [1, 1, 0, 1, 1, 1],
+        [0, 2, 1, 0, 1, 1],
+        [0, 3, 0, 0, 2, 0],
+        [0, 0, 1, 1, 2, 1],
+        [0, 1, 0, 1, 2, 1],
+        [1, 0, 2, 0, 2, 1],
+        [3, 0, 1, 0, 2, 1],
+        [3, 0, 0, 1, 2, 0],
+        [3, 1, 0, 0, 2, 1],
+        [3, 0, 0, 1, 1, 1],
+        [1, 1, 0, 1, 2, 1],
+        [0, 2, 1, 0, 2, 1],
+        [0, 1, 1, 1, 1, 1],
+        [0, 2, 0, 1, 2, 0],
+        [1, 0, 1, 1, 2, 1],
+        [0, 1, 2, 0, 2, 1],
+        [3, 0, 0, 1, 2, 1],
+        [2, 1, 1, 0, 2, 1],
+        [2, 0, 1, 1, 1, 1],
+        [2, 1, 0, 1, 2, 0],
+        [2, 2, 0, 0, 2, 1],
+        [2, 1, 0, 1, 1, 1],
+        [0, 3, 0, 0, 2, 1],
+        [2, 0, 2, 0, 2, 1],
+        [0, 1, 1, 1, 2, 1],
+        [2, 1, 0, 1, 2, 1],
+        [1, 2, 1, 0, 2, 1],
+        [1, 1, 1, 1, 1, 1],
+        [1, 2, 0, 1, 2, 0],
+        [0, 2, 0, 1, 2, 1],
+        [2, 0, 1, 1, 2, 1],
+        [1, 1, 2, 0, 2, 1],
+        [1, 3, 0, 0, 2, 1],
+        [4, 0, 0, 1, 2, 1],
+        [1, 1, 1, 1, 2, 1],
+        [0, 2, 2, 0, 2, 1],
+        [1, 2, 0, 1, 2, 1],
+        [0, 3, 1, 0, 2, 1],
+        [3, 0, 1, 1, 2, 1],
+        [3, 1, 0, 1, 2, 1],
+        [0, 2, 1, 1, 2, 1],
+        [2, 1, 1, 1, 2, 1],
+        [2, 2, 0, 1, 2, 1],
+        [1, 2, 1, 1, 2, 1],
+      ],
+    )
+
+    basis_string = basis_lie_highest_weight_string(:A, 3, [2, 1, 1], [2, 3, 1, 2, 3, 1])
+
+    @test issetequal(
+      [only(exponents(m)) for m in monomials(basis_string)],
+      [
+        [0, 0, 0, 0, 0, 0],
+        [0, 0, 1, 0, 0, 0],
+        [1, 0, 0, 0, 0, 0],
+        [0, 1, 0, 0, 0, 0],
+        [1, 0, 1, 0, 0, 0],
+        [0, 0, 1, 1, 0, 0],
+        [1, 1, 0, 0, 0, 0],
+        [0, 1, 0, 1, 0, 0],
+        [0, 0, 2, 0, 0, 0],
+        [0, 1, 1, 0, 0, 0],
+        [1, 1, 1, 0, 0, 0],
+        [0, 1, 1, 1, 0, 0],
+        [0, 1, 0, 1, 0, 1],
+        [0, 0, 1, 1, 1, 0],
+        [1, 0, 2, 0, 0, 0],
+        [0, 0, 2, 1, 0, 0],
+        [2, 0, 1, 0, 0, 0],
+        [2, 1, 0, 0, 0, 0],
+        [0, 2, 0, 1, 0, 0],
+        [0, 1, 2, 0, 0, 0],
+        [1, 1, 2, 0, 0, 0],
+        [0, 1, 2, 1, 0, 0],
+        [0, 1, 1, 1, 0, 1],
+        [0, 0, 2, 1, 1, 0],
+        [2, 1, 1, 0, 0, 0],
+        [1, 1, 1, 1, 0, 0],
+        [1, 1, 0, 1, 0, 1],
+        [1, 0, 1, 1, 1, 0],
+        [0, 2, 1, 1, 0, 0],
+        [0, 2, 0, 1, 0, 1],
+        [2, 0, 2, 0, 0, 0],
+        [1, 0, 2, 1, 0, 0],
+        [1, 2, 0, 1, 0, 0],
+        [0, 0, 3, 1, 0, 0],
+        [2, 1, 2, 0, 0, 0],
+        [1, 1, 2, 1, 0, 0],
+        [1, 1, 1, 1, 0, 1],
+        [1, 0, 2, 1, 1, 0],
+        [0, 1, 1, 2, 0, 1],
+        [0, 0, 2, 2, 1, 0],
+        [1, 2, 1, 1, 0, 0],
+        [1, 2, 0, 1, 0, 1],
+        [0, 2, 0, 2, 0, 1],
+        [0, 1, 1, 2, 1, 0],
+        [0, 1, 3, 1, 0, 0],
+        [0, 0, 3, 1, 1, 0],
+        [0, 2, 2, 1, 0, 0],
+        [0, 2, 1, 1, 0, 1],
+        [1, 0, 3, 1, 0, 0],
+        [3, 0, 2, 0, 0, 0],
+        [3, 1, 1, 0, 0, 0],
+        [1, 2, 2, 1, 0, 0],
+        [1, 2, 1, 1, 0, 1],
+        [0, 2, 1, 2, 0, 1],
+        [0, 2, 0, 2, 0, 2],
+        [0, 1, 2, 2, 1, 0],
+        [0, 1, 1, 2, 1, 1],
+        [1, 1, 3, 1, 0, 0],
+        [1, 0, 3, 1, 1, 0],
+        [0, 1, 2, 2, 0, 1],
+        [0, 0, 3, 2, 1, 0],
+        [3, 1, 2, 0, 0, 0],
+        [2, 1, 2, 1, 0, 0],
+        [2, 1, 1, 1, 0, 1],
+        [2, 0, 2, 1, 1, 0],
+        [2, 2, 1, 1, 0, 0],
+        [2, 2, 0, 1, 0, 1],
+        [0, 3, 0, 2, 0, 1],
+        [0, 2, 3, 1, 0, 0],
+        [2, 0, 3, 1, 0, 0],
+        [1, 2, 3, 1, 0, 0],
+        [0, 2, 2, 2, 0, 1],
+        [0, 1, 3, 2, 1, 0],
+        [0, 1, 2, 2, 1, 1],
+        [2, 2, 2, 1, 0, 0],
+        [2, 2, 1, 1, 0, 1],
+        [1, 2, 1, 2, 0, 1],
+        [1, 2, 0, 2, 0, 2],
+        [1, 1, 2, 2, 1, 0],
+        [1, 1, 1, 2, 1, 1],
+        [0, 3, 1, 2, 0, 1],
+        [0, 3, 0, 2, 0, 2],
+        [2, 1, 3, 1, 0, 0],
+        [2, 0, 3, 1, 1, 0],
+        [1, 1, 2, 2, 0, 1],
+        [1, 0, 3, 2, 1, 0],
+        [1, 3, 0, 2, 0, 1],
+        [0, 0, 4, 2, 1, 0],
+        [4, 1, 2, 0, 0, 0],
+        [2, 2, 3, 1, 0, 0],
+        [1, 2, 2, 2, 0, 1],
+        [1, 1, 3, 2, 1, 0],
+        [1, 1, 2, 2, 1, 1],
+        [0, 2, 1, 3, 0, 2],
+        [0, 1, 2, 3, 1, 1],
+        [1, 3, 1, 2, 0, 1],
+        [1, 3, 0, 2, 0, 2],
+        [0, 3, 0, 3, 0, 2],
+        [0, 2, 1, 3, 1, 1],
+        [0, 1, 4, 2, 1, 0],
+        [0, 3, 2, 2, 0, 1],
+        [1, 0, 4, 2, 1, 0],
+        [3, 1, 3, 1, 0, 0],
+        [3, 0, 3, 1, 1, 0],
+        [3, 2, 2, 1, 0, 0],
+        [3, 2, 1, 1, 0, 1],
+        [1, 3, 2, 2, 0, 1],
+        [0, 3, 1, 3, 0, 2],
+        [0, 2, 2, 3, 1, 1],
+        [0, 2, 1, 3, 1, 2],
+        [1, 1, 4, 2, 1, 0],
+        [0, 1, 3, 3, 1, 1],
+        [3, 2, 3, 1, 0, 0],
+        [2, 2, 2, 2, 0, 1],
+        [2, 1, 3, 2, 1, 0],
+        [2, 1, 2, 2, 1, 1],
+        [2, 3, 1, 2, 0, 1],
+        [2, 3, 0, 2, 0, 2],
+        [0, 4, 0, 3, 0, 2],
+        [2, 0, 4, 2, 1, 0],
+        [0, 2, 3, 3, 1, 1],
+        [2, 3, 2, 2, 0, 1],
+        [1, 3, 1, 3, 0, 2],
+        [1, 2, 2, 3, 1, 1],
+        [1, 2, 1, 3, 1, 2],
+        [0, 4, 1, 3, 0, 2],
+        [2, 1, 4, 2, 1, 0],
+        [1, 1, 3, 3, 1, 1],
+        [1, 4, 0, 3, 0, 2],
+        [4, 2, 3, 1, 0, 0],
+        [1, 2, 3, 3, 1, 1],
+        [0, 2, 2, 4, 1, 2],
+        [1, 4, 1, 3, 0, 2],
+        [0, 3, 1, 4, 1, 2],
+        [3, 1, 4, 2, 1, 0],
+        [3, 3, 2, 2, 0, 1],
+        [0, 3, 2, 4, 1, 2],
+        [2, 2, 3, 3, 1, 1],
+        [2, 4, 1, 3, 0, 2],
+        [1, 3, 2, 4, 1, 2],
+      ],
+    )
+
+    basis_nz = basis_lie_highest_weight_nz(:A, 3, [2, 1, 1], [2, 3, 1, 2, 3, 1])
+
+    @test issetequal(
+      [only(exponents(m)) for m in monomials(basis_nz)],
+      [
+        [0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 1],
+        [0, 0, 0, 1, 0, 0],
+        [0, 0, 0, 0, 1, 0],
+        [0, 0, 0, 1, 0, 1],
+        [0, 0, 1, 1, 0, 0],
+        [0, 0, 0, 1, 1, 0],
+        [0, 1, 0, 1, 0, 0],
+        [0, 0, 0, 0, 0, 2],
+        [0, 0, 0, 0, 1, 1],
+        [0, 0, 0, 1, 1, 1],
+        [0, 1, 0, 1, 0, 1],
+        [0, 0, 1, 1, 1, 0],
+        [0, 1, 1, 1, 0, 0],
+        [0, 0, 0, 1, 0, 2],
+        [0, 0, 1, 1, 0, 1],
+        [0, 0, 0, 2, 0, 1],
+        [0, 0, 0, 2, 1, 0],
+        [0, 1, 0, 1, 1, 0],
+        [0, 0, 0, 0, 1, 2],
+        [0, 0, 0, 1, 1, 2],
+        [0, 1, 0, 1, 0, 2],
+        [0, 0, 1, 1, 1, 1],
+        [0, 1, 1, 1, 0, 1],
+        [0, 0, 0, 2, 1, 1],
+        [0, 1, 0, 2, 0, 1],
+        [0, 0, 1, 2, 1, 0],
+        [1, 1, 1, 1, 0, 0],
+        [0, 1, 0, 1, 1, 1],
+        [0, 1, 1, 1, 1, 0],
+        [0, 0, 0, 2, 0, 2],
+        [0, 0, 1, 2, 0, 1],
+        [0, 1, 0, 2, 1, 0],
+        [0, 0, 1, 1, 0, 2],
+        [0, 0, 0, 2, 1, 2],
+        [0, 1, 0, 2, 0, 2],
+        [0, 0, 1, 2, 1, 1],
+        [0, 1, 1, 2, 0, 1],
+        [1, 1, 1, 1, 0, 1],
+        [0, 0, 2, 2, 1, 0],
+        [0, 1, 0, 2, 1, 1],
+        [0, 2, 0, 2, 0, 1],
+        [0, 1, 1, 2, 1, 0],
+        [1, 1, 1, 1, 1, 0],
+        [0, 0, 1, 1, 1, 2],
+        [0, 1, 1, 1, 0, 2],
+        [0, 1, 0, 1, 1, 2],
+        [0, 1, 1, 1, 1, 1],
+        [0, 0, 1, 2, 0, 2],
+        [0, 0, 0, 3, 0, 2],
+        [0, 0, 0, 3, 1, 1],
+        [0, 1, 0, 2, 1, 2],
+        [0, 2, 0, 2, 0, 2],
+        [0, 1, 1, 2, 1, 1],
+        [1, 1, 1, 1, 1, 1],
+        [0, 2, 1, 2, 0, 1],
+        [0, 1, 2, 2, 1, 0],
+        [0, 0, 1, 2, 1, 2],
+        [0, 1, 1, 2, 0, 2],
+        [1, 1, 1, 1, 0, 2],
+        [0, 0, 2, 2, 1, 1],
+        [0, 0, 0, 3, 1, 2],
+        [0, 1, 0, 3, 0, 2],
+        [0, 0, 1, 3, 1, 1],
+        [1, 1, 1, 2, 0, 1],
+        [0, 1, 0, 3, 1, 1],
+        [1, 1, 1, 2, 1, 0],
+        [0, 2, 0, 2, 1, 1],
+        [0, 1, 1, 1, 1, 2],
+        [0, 0, 1, 3, 0, 2],
+        [0, 1, 1, 2, 1, 2],
+        [1, 1, 1, 1, 1, 2],
+        [0, 2, 1, 2, 0, 2],
+        [0, 1, 2, 2, 1, 1],
+        [0, 1, 0, 3, 1, 2],
+        [0, 2, 0, 3, 0, 2],
+        [0, 1, 1, 3, 1, 1],
+        [1, 1, 1, 2, 1, 1],
+        [1, 2, 1, 2, 0, 1],
+        [1, 1, 2, 2, 1, 0],
+        [0, 2, 0, 2, 1, 2],
+        [0, 2, 1, 2, 1, 1],
+        [0, 0, 1, 3, 1, 2],
+        [0, 1, 1, 3, 0, 2],
+        [1, 1, 1, 2, 0, 2],
+        [0, 0, 2, 3, 1, 1],
+        [0, 2, 0, 3, 1, 1],
+        [0, 0, 2, 2, 1, 2],
+        [0, 0, 0, 4, 1, 2],
+        [0, 1, 1, 3, 1, 2],
+        [1, 1, 1, 2, 1, 2],
+        [0, 2, 1, 3, 0, 2],
+        [1, 2, 1, 2, 0, 2],
+        [0, 1, 2, 3, 1, 1],
+        [1, 1, 2, 2, 1, 1],
+        [0, 2, 0, 3, 1, 2],
+        [0, 3, 0, 3, 0, 2],
+        [0, 2, 1, 3, 1, 1],
+        [1, 2, 1, 2, 1, 1],
+        [0, 1, 2, 2, 1, 2],
+        [0, 2, 1, 2, 1, 2],
+        [0, 0, 2, 3, 1, 2],
+        [0, 0, 1, 4, 1, 2],
+        [1, 1, 1, 3, 0, 2],
+        [0, 1, 0, 4, 1, 2],
+        [1, 1, 1, 3, 1, 1],
+        [0, 2, 1, 3, 1, 2],
+        [1, 2, 1, 2, 1, 2],
+        [0, 3, 1, 3, 0, 2],
+        [0, 2, 2, 3, 1, 1],
+        [0, 1, 2, 3, 1, 2],
+        [1, 1, 2, 2, 1, 2],
+        [0, 1, 1, 4, 1, 2],
+        [1, 1, 1, 3, 1, 2],
+        [1, 2, 1, 3, 0, 2],
+        [1, 1, 2, 3, 1, 1],
+        [0, 2, 0, 4, 1, 2],
+        [1, 2, 1, 3, 1, 1],
+        [0, 3, 0, 3, 1, 2],
+        [0, 0, 2, 4, 1, 2],
+        [0, 2, 2, 3, 1, 2],
+        [0, 2, 1, 4, 1, 2],
+        [1, 2, 1, 3, 1, 2],
+        [1, 3, 1, 3, 0, 2],
+        [1, 2, 2, 3, 1, 1],
+        [0, 3, 1, 3, 1, 2],
+        [0, 1, 2, 4, 1, 2],
+        [1, 1, 2, 3, 1, 2],
+        [0, 3, 0, 4, 1, 2],
+        [1, 1, 1, 4, 1, 2],
+        [0, 2, 2, 4, 1, 2],
+        [1, 2, 2, 3, 1, 2],
+        [0, 3, 1, 4, 1, 2],
+        [1, 3, 1, 3, 1, 2],
+        [1, 1, 2, 4, 1, 2],
+        [1, 2, 1, 4, 1, 2],
+        [0, 3, 2, 4, 1, 2],
+        [1, 2, 2, 4, 1, 2],
+        [1, 3, 1, 4, 1, 2],
+        [1, 3, 2, 4, 1, 2],
+      ],
+    )
+  end
+
+  @testset "Check dimension" begin
+    @testset "Monomial order $monomial_ordering" for monomial_ordering in
+                                                     (:lex, :revlex, :degrevlex)
+      check_dimension(:A, 3, [1, 1, 1], monomial_ordering)
+      #check_dimension(:B, 3, [2,1,0], monomial_ordering)
+      #check_dimension(:C, 3, [1,1,1], monomial_ordering)
+      #check_dimension(:D, 4, [3,0,1,1], monomial_ordering)
+      #check_dimension(:F, 4, [2,0,1,0], monomial_ordering)
+      #check_dimension(:G, 2, [1,0], monomial_ordering)
+      #check_dimension(:G, 2, [2,2], monomial_ordering)
+    end
+  end
+end
diff --git a/experimental/BasisLieHighestWeight/test/NewMonomial-test.jl b/experimental/BasisLieHighestWeight/test/NewMonomial-test.jl
new file mode 100644
index 000000000000..91dbbfb6c38c
--- /dev/null
+++ b/experimental/BasisLieHighestWeight/test/NewMonomial-test.jl
@@ -0,0 +1,30 @@
+@testset "Test NewMonomial" begin
+  weight = BasisLieHighestWeight.weight
+  calc_vec = BasisLieHighestWeight.calc_vec
+
+  ZZx, _ = PolynomialRing(ZZ, 2)
+  x = gens(ZZx)
+  mon1 = ZZx(1)
+  mon2 = x[1]^2 * x[2]
+  weights = [[ZZ(1), ZZ(1)], [ZZ(2), ZZ(1)]]
+  A = sparse_matrix(ZZ, 2, 2) # [0, 1; 2, 1]
+  setindex!(A, sparse_row(ZZ, [2], [ZZ(1)]), 1)
+  setindex!(A, sparse_row(ZZ, [1, 2], [ZZ(2), ZZ(1)]), 2)
+  B = sparse_matrix(ZZ, 2, 2) # [1, 0; 2, 0]
+  setindex!(B, sparse_row(ZZ, [1], [ZZ(1)]), 1)
+  setindex!(B, sparse_row(ZZ, [2], [ZZ(2)]), 2)
+  matrices_of_operators = [A, B]
+  v0 = sparse_row(ZZ, [1], [1]) # [1, 0]
+
+  mon2_vec = sparse_row(ZZ, [1, 2], [2, 2])
+
+  @testset "weight" begin
+    @test isequal(weight(mon1, weights), [ZZ(0), ZZ(0)])
+    @test isequal(weight(mon2, weights), [ZZ(4), ZZ(3)])
+  end
+
+  @testset "calc_vec" begin
+    @test isequal(calc_vec(v0, mon1, matrices_of_operators), v0)
+    @test isequal(calc_vec(v0, mon2, matrices_of_operators), mon2_vec)
+  end
+end
diff --git a/experimental/BasisLieHighestWeight/test/RootConversion-test.jl b/experimental/BasisLieHighestWeight/test/RootConversion-test.jl
new file mode 100644
index 000000000000..f4f8fb66de35
--- /dev/null
+++ b/experimental/BasisLieHighestWeight/test/RootConversion-test.jl
@@ -0,0 +1,28 @@
+
+@testset "Test RootConversion" begin
+  w_to_alpha = BasisLieHighestWeight.w_to_alpha
+  alpha_to_w = BasisLieHighestWeight.alpha_to_w
+
+  function test_inverse_alpha_w(L, weight)
+    @test w_to_alpha(L, alpha_to_w(L, weight)) == weight # alpha -> w -> alpha
+    @test alpha_to_w(L, w_to_alpha(L, weight)) == weight # w -> alpha -> w
+  end
+
+  @testset "Dynkin type $dynkin" for dynkin in (:A, :B, :C, :D, :E, :F, :G)
+    @testset "n = $n" for n in 1:10
+      if (
+        !(dynkin == :B && n < 2) &&
+        !(dynkin == :C && n < 2) &&
+        !(dynkin == :D && n < 4) &&
+        !(dynkin == :E && !(n == 6 || n == 7 || n == 8)) &&
+        !(dynkin == :F && n != 4) &&
+        !(dynkin == :G && (n != 2))
+      )
+        weight = [rand(QQ, -10:10) for _ in 1:n]
+        print(".")
+        L = BasisLieHighestWeight.lie_algebra(dynkin, n)
+        test_inverse_alpha_w(L, weight)
+      end
+    end
+  end
+end
diff --git a/experimental/BasisLieHighestWeight/test/runtests.jl b/experimental/BasisLieHighestWeight/test/runtests.jl
new file mode 100644
index 000000000000..b07a9d40abbe
--- /dev/null
+++ b/experimental/BasisLieHighestWeight/test/runtests.jl
@@ -0,0 +1,3 @@
+include("NewMonomial-test.jl")
+include("RootConversion-test.jl")
+include("MainAlgorithm-test.jl")