diff --git a/docs/Build.jl b/docs/Build.jl
index 9977684fef..fd45db4e50 100644
--- a/docs/Build.jl
+++ b/docs/Build.jl
@@ -44,6 +44,7 @@ pages = [
                      "misc/conjugacy.md",
                      ],
          "Extra features" => ["features/macros.md",
+                              "features/mset.md",
                              ],
          "Examples" => "examples.md",
          "References" => "references.md",
diff --git a/docs/src/features/mset.md b/docs/src/features/mset.md
new file mode 100644
index 0000000000..dabfa94ba2
--- /dev/null
+++ b/docs/src/features/mset.md
@@ -0,0 +1,99 @@
+# Multi-sets and sub-set iterators
+
+```@meta
+CurrentModule = Hecke
+DocTestSetup = quote
+    using Hecke
+end
+```
+
+## Multi-sets
+
+### Type and constructors
+
+Objects of type `MSet` consists of a dictionary whose keys are the elements in
+the set, and the values are their respective multiplicity.
+
+```@docs
+MSet
+```
+
+We can create multi-sets from any finite iterator, dictionary or pair of lists
+with the appropriate conditions.
+
+```@docs
+multiset
+```
+
+### Functions
+
+One can iterate over an `MSet` as on a regular `Set`. Here is moreover a list
+of functions defined for collections of objects which are currently available
+for `MSet`:
+
+* `==`
+* `all`
+* `any`
+* `copy`
+* `delete!`
+* `eltype`
+* `filter`
+* `filter!`
+* `in`
+* `intersect`
+* `intersect!`
+* `isempty`
+* `issubset`
+* `length`
+* `pop!`
+* `push!`
+* `setdiff`
+* `setdiff!`
+* `similar`
+* `unique`
+* `union`
+* `union!`
+* ...
+
+Note that `pop!` and `delete!` for `MSet` are available but have a different behaviour.
+For an element `x` in an multi-set `M <: MSet`, then `pop!(M, x)` will remove
+*one* instance of `x` in `M` - in particular `multiplicity(M, x)` will drop by
+$1$. Much stronger, `delete!(M, x)` will remove *all* instances of `x` in `M` and
+so `multiplicity(M, x)` will be $0$.
+
+While `unique` will return the keys of the underlying dictionary, one can access
+the values (i.e. the multiplicities of the elements in the multi-set) via the
+following functions:
+
+```@docs
+multiplicities(::MSet)
+multiplicity(::MSet{T}, ::T) where {T}
+```
+
+Finally, the sum and difference for `MSet` are also available. Difference is
+given by the complements of sets and the sum is given by disjoint union of sets.
+
+```@docs
+Base.:(+)(::MSet, ::MSet)
+Base.:(-)(::MSet, ::MSet...)
+```
+
+## Sub-set iterators
+
+### Sub-multi-sets
+
+```@docs
+subsets(::MSet{Int})
+```
+
+### Sub-sets
+
+```@docs
+subsets(::Set{Int})
+```
+
+### Sub-sets of a given size
+
+```@docs
+subsets(::Set, ::Int)
+```
diff --git a/docs/src/function_fields/degree_localization.md b/docs/src/function_fields/degree_localization.md
index b100a09dcd..7e01d9d8a3 100644
--- a/docs/src/function_fields/degree_localization.md
+++ b/docs/src/function_fields/degree_localization.md
@@ -13,7 +13,7 @@ Given $k(x)$ a (univariate) rational function field, there are two rings of inte
 both of which are Euclidean:
 
 * $k[x]$
-* $k_\infty(x) = \{a/b | a, b \in k[x] \;\;\mbox{where}\;\; \deg(a) \leq \deg(b)\}
+* $k_\infty(x) = \{a/b | a, b \in k[x] \;\;\mbox{where}\;\; \deg(a) \leq \deg(b)\}$
 
 The second of these rings is the localization of $k[1/x]$ at $(1/x)$ inside the rational
 function field $k(x)$, i.e. the localization of the function field at the point at
diff --git a/examples/NFDB.jl b/examples/NFDB.jl
index 34b6d9bf1f..16bc3796c5 100644
--- a/examples/NFDB.jl
+++ b/examples/NFDB.jl
@@ -994,10 +994,12 @@ end
 #
 ################################################################################
 
-function Base.merge!(R::NFDB{1}, D1::NFDB{1})
+function Base.merge!(R::NFDB{1}, D1::NFDB{1}; skip_update = false)
   sizehint!(R.fields, length(R) + length(D1))
   append!(R.fields, D1.fields)
-  update_properties!(R)
+  if !skip_update
+    update_properties!(R)
+  end
   return R
 end
 
@@ -1014,10 +1016,13 @@ function Base.merge(D::Vector{NFDB{1}})
   end
 
   R = NFDB{1}()
+  sizehint!(R.fields, sum(length(d) for d in D))
   for i in 1:length(D)
-    merge!(R, D[i])
+    merge!(R, D[i], skip_update = true)
   end
 
+  update_properties!(R)
+
   return R
 end
 
diff --git a/src/Deprecations.jl b/src/Deprecations.jl
index 8f99968897..1559abce94 100644
--- a/src/Deprecations.jl
+++ b/src/Deprecations.jl
@@ -125,6 +125,7 @@
 
 @deprecate mul(A::SRow{T}, B::SMat{T}) where T A*B
 
+@deprecate field_of_fractions(O::GenOrd) function_field(O::GenOrd)
 
 # Things that moved to Nemo
 
diff --git a/src/FunField/DegreeLocalization.jl b/src/FunField/DegreeLocalization.jl
index 1ced56ee26..b5943c83c7 100644
--- a/src/FunField/DegreeLocalization.jl
+++ b/src/FunField/DegreeLocalization.jl
@@ -89,7 +89,7 @@ end
     in(a::Generic.RationalFunctionFieldElem{T}, R::KInftyRing{T}) where T <: FieldElement
 
 Return `true` if the given element of the rational function field is an
-element of `k_\infty(x)`, i.e. if `degree(numerator) <= degree(denominator)`.
+element of $k_\infty(x)$, i.e. if `degree(numerator) <= degree(denominator)`.
 """
 function in(a::Generic.RationalFunctionFieldElem{T}, R::KInftyRing{T}) where T <: FieldElement
   if parent(a) != function_field(R)
diff --git a/src/FunField/Divisor.jl b/src/FunField/Divisor.jl
index ff44c5d7b3..84cdd633db 100644
--- a/src/FunField/Divisor.jl
+++ b/src/FunField/Divisor.jl
@@ -171,15 +171,10 @@ function ideals(D)
 end
 
 @doc raw"""
-    field_of_fractions(O::GenOrd) -> FunctionField
     function_field(O::GenOrd) -> FunctionField
 
 Return the function field of O.
 """
-function field_of_fractions(O::GenOrd)
-  return O.F
-end
-
 function function_field(O::GenOrd)
   return O.F
 end
diff --git a/src/Misc/MSet.jl b/src/Misc/MSet.jl
index 2385328758..f9f9b5dc42 100644
--- a/src/Misc/MSet.jl
+++ b/src/Misc/MSet.jl
@@ -1,14 +1,32 @@
-export MSet, multiplicities, multiplicity, subsets
+export MSet
+export multiplicities
+export multiplicity
+export multiset
+export subsets
+
+###############################################################################
+#
+#  Multi-sets
+#
+###############################################################################
+
+### Type and constructors
 
 @doc raw"""
+    MSet{T} <: AbstractSet{T}
+
 Type for a multi-set, i.e. a set where elements are not unique, they
 (can) have a multiplicity. `MSet`s can be created from any finite iterator.
 
-# Example
+# Examples
 ```jldoctest
 julia> MSet([1,1,2,3,4,4,5])
-MSet(5, 4 : 2, 2, 3, 1 : 2)
-
+MSet{Int64} with 7 elements:
+  5
+  4 : 2
+  2
+  3
+  1 : 2
 ```
 `4 : 2` means the element `4` has multiplicity `2`, i.e. was included twice.
 """
@@ -16,26 +34,159 @@ mutable struct MSet{T} <: AbstractSet{T}
   dict::Dict{T,Int}
 
   MSet{T}() where {T} = new{T}(Dict{T,Int}())
-  MSet{T}(itr) where {T} = union!(new{T}(Dict{T,Int}()), itr)
+
+  function MSet{T}(itr) where {T}
+    s = new{T}(Dict{T, Int}())
+    for x in itr
+      push!(s, x)
+    end
+    return s
+  end
+
+  MSet{T}(d::Dict{T, Int}) where {T} = new{T}(d)
+  MSet{T}(l::Vector{T}, m::Vector{Int}) where {T} = MSet{T}(Dict(zip(l, m)))
 end
 
 MSet() = MSet{Any}()
 MSet(itr) = MSet{eltype(itr)}(itr)
+MSet(d::Dict{T, Int}) where {T} = MSet{T}(d)
+MSet(l::Vector{T}, m::Vector{Int}) where {T} = MSet{T}(l, m)
+
+@doc raw"""
+    multiset(iter) -> MSet{eltype(iter)}
+    multiset(d::Dict{T, Int}) -> MSet{T}
+    multiset{l::Vector{T}, m::Vector{Int}} -> MSet{T}
+
+Given either:
+ - a finite iterator `iter`;
+ - a dictionary `d` whose values are positive integers;
+ - a list `l` and a list of positive integers `m` of same length as `l`;
+return the asscciated multi-set `M`.
+
+# Examples
+```jldoctest
+julia> str = "A nice sentence"
+"A nice sentence"
+
+julia> multiset(str)
+MSet{Char} with 15 elements:
+  'n' : 3
+  'A'
+  'c' : 2
+  'i'
+  'e' : 4
+  's'
+  't'
+  ' ' : 2
+
+julia> multiset(Int[x^3%8 for x = 1:50])
+MSet{Int64} with 50 elements:
+  0 : 25
+  5 : 6
+  7 : 6
+  3 : 6
+  1 : 7
+
+julia> d = Dict("a" => 4, "b" => 1, "c" =>9)
+Dict{String, Int64} with 3 entries:
+  "c" => 9
+  "b" => 1
+  "a" => 4
+
+julia> multiset(d)
+MSet{String} with 14 elements:
+  "c" : 9
+  "b"
+  "a" : 4
+
+julia> multiset(["a", "b", "c"], [4, 1, 9])
+MSet{String} with 14 elements:
+  "c" : 9
+  "b"
+  "a" : 4
+```
+"""
+multiset(itr) = MSet(itr)
+
+function multiset(d::Dict{T, Int}) where {T}
+  @req minimum(values(d)) > 0 "The values of d must be positive integers"
+  return MSet{T}(d)
+end
+
+function multiset(l::Vector{T}, m::Vector{Int}) where {T}
+  @req length(m) == length(l) "Lists must have the same length"
+  @req minimum(m) > 0 "Multiplicity list must consist of positive integers"
+  return MSet{T}(l, m)
+end
+
+@doc raw"""
+    multiset(T::Type) -> MSet{T}
 
+Create an empty multi-set `M` with elements of type `T`.
 
+# Examples
+```jldoctest
+julia> multiset(QQFieldElem)
+MSet{QQFieldElem}()
+
+julia> multiset()
+MSet{Any}()
+```
+"""
+multiset() = MSet()
+
+multiset(T::DataType) = MSet{T}()
+
+@doc raw"""
+    similar(M::MSet{T}) -> MSet{T}
+    similar(M::MSet, T::Type) -> MSet{T}
+
+Create an empty multi-set with elements of type `T`.
+
+# Examples
+```jldoctest
+julia> m = multiset([1,1,2,3,4,4,5])
+MSet{Int64} with 7 elements:
+  5
+  4 : 2
+  2
+  3
+  1 : 2
+
+julia> similar(m)
+MSet{Int64}()
+
+julia> similar(m, QQFieldElem)
+MSet{QQFieldElem}()
+```
+"""
 Base.similar(::MSet{T}) where {T} = MSet{T}()
 Base.similar(::MSet, T::Type) = MSet{T}()
 
-#TODO: compact printing, remove trailing , ... the works...
-function Base.show(io::IO, ::MIME"text/plain", s::MSet)
-    print(io,"MSet")
-    if isempty(s)
-        print(io,"{",eltype(s),"}()")
-        return
+Base.copy(s::MSet) = union!(similar(s), s)
+
+# Only for internal use
+dict(s::MSet) = s.dict
+
+### Show methods
+
+# We try to adopt the same conventions as in Oscar, so one-line printing should
+# stay in one line, and we do not give details about what is in the MSet: the
+# detailled printing will take care of it
+function Base.show(io::IO, s::MSet{T}) where {T}
+  if isempty(s)
+    print(io, "MSet{$(eltype(s))}()")
+  elseif get(io, :supercompact, false)
+    io = pretty(io)
+    print(io, "Multi-set with ", ItemQuantity(length(s), "element"))
+    if !(T == Any)
+      print(io, " of type $T")
     end
-    print(io,"(")
+  else
+    print(io, "MSet(")
     first = true
-    for (k,v) = s.dict
+    d = dict(s)
+    for (k, v) in d
       first || print(io, ", ")
       first = false
       if v > 1
@@ -44,21 +195,103 @@ function Base.show(io::IO, ::MIME"text/plain", s::MSet)
         print(io, "$k")
       end
     end
-    print(io,")")
+    print(io, ")")
+  end
 end
 
+# We use the one-line (hopefully implemented!) printing system for the elements
+# in
+function Base.show(io::IO, ::MIME"text/plain", s::MSet)
+  print(io,"MSet{",eltype(s),"}")
+  if isempty(s)
+    print(io,"()")
+  else
+    io = pretty(io)
+    szh, szw = displaysize()
+    szh -= 5
+    szw -= 10
+    print(io, " with ", ItemQuantity(length(s), "element"), ":")
+    print(io, Indent())
+    d = dict(s)
+    rmax = maximum(ndigits(k) for k in values(d))
+    offmax = szw - (rmax + 3)
+    if length(d) <= szh
+      lmax = min(maximum(length(sprint(show, a)) for a in keys(d)), offmax)
+      for (k, v) in d
+        pk = sprint(show, k)
+        lk = length(pk)
+        println(io)
+        if lk > offmax
+          print(io, pk[1:offmax-length(" \u2026")], " \u2026")
+        else
+          print(io, pk)
+        end
+        lk = min(offmax, lk)
+        if v > 1
+          print(io, " "^(lmax-lk+1), ": $v")
+        end
+      end
+    else
+      un = collect(keys(d))[1:szh]
+      lmax = min(maximum(length(sprint(show, a)) for a in un), offmax)
+      for k in un
+        println(io)
+        pk = sprint(show, k)
+        lk = length(sprint(show, k))
+        v = d[k]
+        if lk > offmax
+          print(io, pk[1:offmax-length(" \u2026")], " \u2026")
+        else
+          print(io, pk)
+        end
+        lk = min(offmax, lk)
+        if v > 1
+          print(io, " "^(lmax-lk+1), ": $v")
+        end
+      end
+      println(io)
+      print(io, "\u22ee")
+    end
+  end
+end
+
+### Iteration
+
 Base.isempty(s::MSet) = isempty(s.dict)
 Base.length(s::MSet) = sum(values(s.dict))
 Base.IteratorSize(::Type{MSet}) = Base.HasLength()
 Base.IteratorEltype(::Type{MSet}) = Base.HasEltype()
 Base.eltype(::Type{MSet{T}}) where {T} = T
-Base.in(x, s::MSet) = haskey(s.dict, x)
+Base.in(x::T, s::MSet{T}) where {T} = haskey(dict(s), x)
+Base.in(x, s::MSet{T}) where {T} = haskey(dict(s), convert(T, x))
+
+function Base.iterate(s::MSet)
+  I = iterate(s.dict)
+  I === nothing && return I
+  return I[1][1], (I[1], I[2], 1)
+end
+
+function Base.iterate(s::MSet, state)
+  if state[3] < state[1][2]
+    return state[1][1], (state[1], state[2], state[3]+1)
+  else
+    I = iterate(s.dict, state[2])
+    I === nothing && return I
+    val, st = I
+    return (val[1], (val, st, 1))
+  end
+end
 
-function Base.push!(s::MSet, x, mult::Int=1)
-  add_to_key!(s.dict, x, mult)
+### MSets operations
+
+function Base.push!(s::MSet{T}, x, mult::Int=1) where {T}
+  @req promote_type(T, typeof(x)) == T "Cannot coerce element"
+  y = x isa T ? x : T(x)
+  add_to_key!(s.dict, y, mult)
 end
 
 function Base.pop!(s::MSet{T}, x) where {T}
+  @req promote_type(T, typeof(x)) == T "Cannot coerce element"
   y = x isa T ? x : T(x)
   y in s || throw(KeyError(y))
   add_to_key!(s.dict, y, -1)
@@ -66,52 +299,226 @@ function Base.pop!(s::MSet{T}, x) where {T}
 end
 
 function Base.pop!(s::MSet{T}, x, default) where {T}
+  @req promote_type(T, typeof(x)) == T "Cannot coerce element"
   y = x isa T ? x : T(x)
-  return y in s ? pop!(s, y) : (default isa T ? default : T(default))
+  return y in s ? pop!(s, y) : (default isa T ? default : convert(T, default))
 end
+
 Base.pop!(s::MSet) = (val = iterate(s.dict)[1][1]; pop!(s, val))
 
 function Base.delete!(s::MSet{T}, x) where {T}
+  @req promote_type(T, typeof(x)) == T "Cannot coerce element"
   y = x isa T ? x : T(x)
   delete!(s.dict, y)
   return s
 end
 
-Base.copy(s::MSet) = union!(similar(s), s)
+Base.setdiff(s::MSet, itrs...) = setdiff!(copy(s), itrs...)
 
-function Base.iterate(s::MSet)
-  I = iterate(s.dict)
-  I === nothing && return I
-  return I[1][1], (I[1], I[2], 1)
+function Base.setdiff!(s::MSet, itrs...)
+  for x in itrs
+    setdiff!(s, x)
+  end
+  return s
 end
 
-function Base.iterate(s::MSet, state)
-  if state[3] < state[1][2]
-    return state[1][1], (state[1], state[2], state[3]+1)
-  else
-    I = iterate(s.dict, state[2])
-    I === nothing && return I
-    val, st = I
-    return (val[1], (val, st, 1))
+function Base.setdiff!(s::MSet, itr)
+  for x in itr
+    pop!(s, x, x)
   end
+  return s
+end
+
+@doc raw"""
+    (-)(s::MSet, itrs...) -> MSet
+
+Return the multi-set associated to the complement in `s` of the collections
+in `itrs`.
+
+Alias for `setdiff(s, itrs...)`.
+
+# Examples
+```jldoctest
+julia> m = multiset("A very nice sentence")
+MSet{Char} with 20 elements:
+  'n' : 3
+  'e' : 5
+  'A'
+  'y'
+  'i'
+  'r'
+  's'
+  't'
+  ' ' : 3
+  'c' : 2
+  'v'
+
+julia> n = multiset("A nice sentence")
+MSet{Char} with 15 elements:
+  'n' : 3
+  'A'
+  'c' : 2
+  'i'
+  'e' : 4
+  's'
+  't'
+  ' ' : 2
+
+julia> n-m
+MSet{Char}()
+
+julia> m-n
+MSet{Char} with 5 elements:
+  'e'
+  'y'
+  'r'
+  ' '
+  'v'
+```
+"""
+Base.:(-)(s::MSet, itrs...) = setdiff(s, itrs...)
+
+function Base.unique(s::MSet)
+  return collect(keys(dict(s)))
+end
+
+function Base.issubset(s1::MSet{T}, s2::MSet{U}) where {T, U}
+  @req promote_type(T, U) == U "Cannot compare multi-sets"
+  for (x, k) in dict(s1)
+    y = convert(U, x)
+    !haskey(dict(s2), y) && return false
+    k > multiplicity(s2, y) && return false
+  end
+  return true
+end
+
+@doc raw"""
+    (+)(s::MSet, itrs...) -> MSet
+
+Return the multi-sets associated to the disjoint union of `s` and the
+collections of objects in `itrs`.
+
+# Examples
+```jldoctest
+julia> m = multiset("A nice sentence")
+MSet{Char} with 15 elements:
+  'n' : 3
+  'A'
+  'c' : 2
+  'i'
+  'e' : 4
+  's'
+  't'
+  ' ' : 2
+
+julia> n = multiset("A very nice sentence")
+MSet{Char} with 20 elements:
+  'n' : 3
+  'e' : 5
+  'A'
+  'y'
+  'i'
+  'r'
+  's'
+  't'
+  ' ' : 3
+  'c' : 2
+  'v'
+
+julia> m + n
+MSet{Char} with 35 elements:
+  'n' : 6
+  'e' : 9
+  'A' : 2
+  's' : 2
+  'i' : 2
+  't' : 2
+  'y'
+  'r'
+  ' ' : 5
+  'c' : 4
+  'v'
+```
+"""
+function Base.:(+)(s1::MSet, s2::MSet)
+  T = Base.promote_eltype(s1, s2)
+  s = similar(s1, T)
+  d = dict(s)
+  for (x, k) in dict(s1)
+    add_to_key!(d, convert(T, x),  k)
+  end
+
+  for (y, k) in dict(s2)
+    add_to_key!(d, convert(T, y), k)
+  end
+  return s
+end
+
+function Base.:(+)(s::MSet, itrs...)
+  s2 = s + multiset(itrs[1])
+  return (+)(s2, itrs[2:end]...)
 end
 
 Base.union(s::MSet) = copy(s)
-function Base.union(s::MSet, sets::Set...)
-    u = MSet{join_eltype(s, sets...)}()
-    union!(u,s)
-    for t in sets
-        union!(u,t)
+
+function Base.union(s::MSet, itrs...)
+  T = Base.promote_eltype(s, itrs...)
+  return union!(similar(s, T), s, itrs...)
+end
+
+function Base.union!(s1::MSet{T}, s2::MSet{U}) where {T, U}
+  @req promote_type(T, U) == T "Cannot coerce elements"
+  d = dict(s1)
+  for (x, k) in d
+    fi1 = filter(isequal(x), keys(dict(s2)))
+    if !isempty(fi1)
+      k = max(k, multiplicity(s2, first(fi1)))
+      d[x] =  k
+    end
+  end
+
+  for (y, k) in dict(s2)
+    fi2 = filter(isequal(y), keys(d))
+    if isempty(fi2)
+      d[convert(T, y)] = k
     end
-    return u
+  end
+  return s1
 end
-Base.union!(s::MSet, xs) = (for x=xs; push!(s,x); end; s)
-Base.union!(s::MSet, xs::AbstractArray) = (for x=xs; push!(s,x); end; s)
 
+function Base.union!(s::MSet, itrs...)
+  union!(s, multiset(itrs[1]))
+  return union!(s, itrs[2:end]...)
+end
+
+function Base.intersect(s::MSet, itrs...)
+  T = Base.promote_eltype(s, itrs...)
+  return intersect!(union!(similar(s, T), s), itrs...)
+end
+
+function Base.intersect!(s1::MSet{T}, s2::MSet) where {T}
+  val = intersect(keys(dict(s1)), keys(dict(s2)))
+  @req promote_type(T, typeof.(val)...) == T "Cannot coerce elements"
+  d = dict(s1)
+  for (x, k) in d
+    if !(x in val)
+      delete!(s1, x)
+    else
+      y = first(filter(y -> x == y[1], dict(s2)))
+      d[x] = min(k, y[2])
+    end
+  end
+  return s1
+end
+
+function Base.intersect!(s::MSet, itrs...)
+  s2 = intersect!(s, multiset(itrs[1]))
+  return intersect!(s2, itrs[2:end]...)
+end
 
 function Base.filter(pred, s::MSet)
   t = similar(s)
-  for (x, m) in s.dict
+  for (x, m) in dict(s)
     if pred(x)
       push!(t, x, m)
     end
@@ -119,48 +526,90 @@ function Base.filter(pred, s::MSet)
   return t
 end
 
+function Base.filter!(pred, s::MSet)
+  for x in keys(dict(s))
+    if !pred(x)
+      delete!(s, x)
+    end
+  end
+  return s
+end
 
 @doc raw"""
-    multiplicities(s::MSet)
+    multiplicities(s::MSet{T}) -> ValueIterator{Dict{T, Int}}
+
+Return an iterator for the multiplicities of all the elements in `s`.
 
-Return an iterator for the multiplicities of all the elements.
+# Examples
+```jldoctest
+julia> m = multiset([1,1,2,3,4,4,5])
+MSet{Int64} with 7 elements:
+  5
+  4 : 2
+  2
+  3
+  1 : 2
+
+julia> mult_m = multiplicities(m)
+ValueIterator for a Dict{Int64, Int64} with 5 entries. Values:
+  1
+  2
+  1
+  1
+  2
+
+julia> collect(mult_m)
+5-element Vector{Int64}:
+ 1
+ 2
+ 1
+ 1
+ 2
+```
 """
 function multiplicities(s::MSet)
-  return values(s.dict)
+  return values(dict(s))
 end
 
 @doc raw"""
-    multiplicity(s::MSet, x)
+    multiplicity(s::MSet{T}, x::T) -> Int
+
+Return the multiplicity of the element `x` in the multi-set `s`. If `x` is not
+in `s`, return 0.
 
-The multiplicity of the element `x`` in the multi-set `s - or zero if
-  `x` is not in `s`,
+# Examples
+```jldoctest
+julia> m = multiset([1,1,2,3,4,4,5])
+MSet{Int64} with 7 elements:
+  5
+  4 : 2
+  2
+  3
+  1 : 2
+
+julia> multiplicity(m, 2)
+1
+
+julia> multiplicity(m, 6)
+0
+```
 """
 function multiplicity(s::MSet{T}, x::T) where {T}
   y = x isa T ? x : T(x)
-  if haskey(s.dict, y)
-    return s.dict[y]
+  if haskey(dict(s), y)
+    return dict(s)[y]
   else
     return 0
   end
 end
 
-function Base.unique(s::MSet)
-  return collect(keys(s.dict))
-end
-
-function Base.setdiff(s::MSet, itrs...)
-  s = copy(s)
-  for i = itrs
-    if i in s
-      pop!(s, i)
-    end
-  end
-  return s
-end
+###############################################################################
+#
+#  Sub-set iterators
+#
+###############################################################################
 
-############################################
-# subsets iterator
-############################################
+### Sub-multi-sets
 
 struct MSubSetItr{T}
   b::Vector{T}
@@ -169,23 +618,17 @@ struct MSubSetItr{T}
 end
 
 @doc raw"""
-    subsets(s::MSet)
+    subsets(s::MSet) -> MSubSetIt{T}
 
-An iterator for all sub-multi-sets of `s`.
+Return an iterator on all sub-multi-sets of `s`.
 """
 function subsets(s::MSet{T}) where T
   # subsets are represented by integers in a weird base
   # the distinct elements are b0...bn with mult mi
   # subset (bi, ni) -> sum ni gi where gi = prod (mj+1)
-  b = collect(unique(s))
-  m = [s.dict[x] for x = b]
-  #= not needed for the iterator
-  g = [1]
-  for i=2:length(b)
-    push!(g, g[end]*(m[i]+1))
-  end
-  =#
-  return MSubSetItr{T}(b, m, length(m) == 0 ? 1 : prod(x+1 for x=m))
+  b = unique(s)
+  m = Int[multiplicity(s, x) for x in b]
+  return MSubSetItr{T}(b, m, length(m) == 0 ? 1 : prod(x+1 for x in m))
 end
 
 function int_to_elt(M::MSubSetItr{T}, i::Int) where T
@@ -218,28 +661,37 @@ Base.IteratorSize(::Type{MSubSetItr}) = Base.HasLength()
 Base.IteratorEltype(::Type{MSubSetItr}) = Base.HasEltype()
 Base.eltype(::Type{MSubSetItr{T}}) where {T} = MSet{T}
 
-function Base.show(io::IO, M::MSubSetItr)
-  println(io, "subset iterator of length $(M.length) for $(M.b) with multiplicities $(M.m)")
+function Base.show(io::IO, M::MSubSetItr{T}) where {T}
+  print(io, "Iterator for subsets of a multi-set")
+  if !get(io, :supercompact, false) && T != Any
+    print(io, " with elements of type $T")
+  end
+end
+
+function Base.show(io::IO, ::MIME"text/plain", M::MSubSetItr)
+  io = pretty(io)
+  println(io, "Iterator for subsets")
+  println(io, Indent(), "of ", MSet(M.b, M.m))
+  print(io, Dedent(), "of length ", M.length)
 end
 
 #... to be completed from base/Set.jl ...
 
-#subsets for Set
+### Arbitrary sub-sets
+
 struct SubSetItr{T}
   b::Vector{T}
   length::Int
 end
 
 @doc raw"""
-    subsets(s::Set)
-    subsets(s::Set, k::Int)
+    subsets(s::Set) -> SubSetItr{T}
 
-An iterator for all sub-sets of `s`. In the second for only
-subsets of size `k` are listed.
+Return an iterator for all sub-sets of `s`.
 """
 function subsets(s::Set{T}) where T
   # subsets are represented by integers in base 2
-  b = collect(unique(s))
+  b = unique(s)
   return SubSetItr{T}(b, 2^length(b))
 end
 
@@ -274,11 +726,22 @@ Base.IteratorSize(::Type{SubSetItr}) = Base.HasLength()
 Base.IteratorEltype(::Type{SubSetItr}) = Base.HasEltype()
 Base.eltype(::Type{SubSetItr{T}}) where {T} = Set{T}
 
-function Base.show(io::IO, M::SubSetItr)
-  println(io, "subset iterator of length $(M.length) for $(M.b)")
+function Base.show(io::IO, M::SubSetItr{T}) where {T}
+  print(io, "Iterator for subsets of a set")
+  if !get(io, :supercompact, false) && T != Any
+    print(io, " with elements of type $T")
+  end
+end
+
+function Base.show(io::IO, ::MIME"text/plain", M::SubSetItr)
+  io = pretty(io)
+  println(io, "Iterator for subsets")
+  println(io, Indent(), "of ", Set(M.b))
+  print(io, Dedent(), "of length ", M.length)
 end
 
-#only subsets of a given size
+### Sub-sets of a given size
+
 struct SubSetSizeItr{T}
   b::Vector{T}
   k::Int #subsets of size k only
@@ -286,12 +749,17 @@ struct SubSetSizeItr{T}
   length::Int
 end
 
+@doc raw"""
+    subsets(s::Set, k::Int) -> SubSetSizeItr{T}
+
+Return an iterator on all sub-sets of size `k` of `s`.
+"""
 function subsets(s::Set{T}, k::Int) where T
   # subsets are represented by integers in the Combinatorial_number_system
   # https://en.wikipedia.org/wiki/Combinatorial_number_system
   # this iterator could indexed, the other one below not
   # maybe use "The coolest way to generate combinations" instead
-  b = collect(unique(s))
+  b = unique(s)
   m = Int(binomial(length(b), k))
   C = Vector{Vector{Int}}()
   while k >= 1
@@ -309,11 +777,9 @@ function subsets(s::Set{T}, k::Int) where T
     push!(C, B)
     k -=1
   end
-
   return SubSetSizeItr{T}(b, length(C), C, m)
 end
 
-
 function int_to_elt(M::SubSetSizeItr{T}, i::Int) where T
   s = Set{T}()
   j = 1
@@ -326,7 +792,6 @@ function int_to_elt(M::SubSetSizeItr{T}, i::Int) where T
   while length(s) < M.k
     push!(s, M.b[length(s)])
   end
-
   return s
 end
 
@@ -353,6 +818,16 @@ function Base.getindex(S::SubSetSizeItr, i::Int)
   return Hecke.int_to_elt(S, i)
 end
 
-function Base.show(io::IO, M::SubSetSizeItr)
-  println(io, "subset iterator of length $(M.length) for $(M.b) and subsets of size $(M.k)")
+function Base.show(io::IO, M::SubSetSizeItr{T}) where {T}
+  print(io, "Iterator for subsets of size $(M.k) of a set")
+  if !get(io, :supercompact, false) && T != Any
+    print(io, " with elements of type $T")
+  end
+end
+
+function Base.show(io::IO, ::MIME"text/plain", M::SubSetSizeItr)
+  io = pretty(io)
+  println(io, "Iterator for subsets of size $(M.k)")
+  println(io, Indent(), "of ", Set(M.b))
+  print(io, Dedent(), "of length ", M.length)
 end
diff --git a/src/Misc/Matrix.jl b/src/Misc/Matrix.jl
index 6a1c40e023..70e00be7c4 100644
--- a/src/Misc/Matrix.jl
+++ b/src/Misc/Matrix.jl
@@ -1468,3 +1468,20 @@ function solve(L::LinearSolveCtx, b::Vector)
   #end
   return fl, deepcopy(w)
 end
+
+################################################################################
+#
+#  Determinant of triangular matrix
+#
+################################################################################
+
+# Compute the determinant of a (lower-left or upper-right) triangular matrix by
+# multiplying the diagonal entries. Nothing is checked.
+function _det_triangular(M::MatElem)
+  R = base_ring(M)
+  d = one(R)
+  for i in 1:nrows(M)
+    mul!(d, d, M[i, i])
+  end
+  return d
+end
diff --git a/src/NumFieldOrd/NfOrd/NfOrd.jl b/src/NumFieldOrd/NfOrd/NfOrd.jl
index 6b08093e6c..007436e0f3 100644
--- a/src/NumFieldOrd/NfOrd/NfOrd.jl
+++ b/src/NumFieldOrd/NfOrd/NfOrd.jl
@@ -747,9 +747,8 @@ checked by computing minimal polynomials. If `isbasis` is set, then elements are
 assumed to form a $\mathbf{Z}$-basis. If `cached` is set, then the constructed
 order is cached for future use.
 """
-function Order(::S, a::Vector{T}; check::Bool = true, isbasis::Bool = false,
+function Order(K::S, a::Vector{T}; check::Bool = true, isbasis::Bool = false,
                cached::Bool = false) where {S <: NumField{QQFieldElem}, T <: NumFieldElem{QQFieldElem}}
-  K = parent(a[1])
   @assert all(x->K == parent(x), a)
   if isbasis
     if check
@@ -769,9 +768,10 @@ end
 
 function Order(K, a::Vector; check::Bool = true, isbasis::Bool = false,
                cached::Bool = true)
-  local b
+  local b::Vector{elem_type(K)}
   try
     b = map(K, a)
+    b = convert(Vector{elem_type(K)}, b)
   catch
     error("Cannot coerce elements from array into the number field")
   end
@@ -886,11 +886,6 @@ function any_order(K::NfAbsNS)
   g = gens(K)
   for i in 1:ngens(K)
     f = denominator(K.pol[i]) * K.pol[i]
-    @show f
-    @show isone(coeff(f, 1))
-    @show coeff(f, 1)
-    @show typeof(f)
-    @show g[i]
     if isone(coeff(f, 1))
       normalized_gens[i] = g[i]
     else
@@ -898,8 +893,6 @@ function any_order(K::NfAbsNS)
     end
   end
 
-  @show normalized_gens
-
   b = Vector{NfAbsNSElem}(undef, degree(K))
   ind = 1
   it = cartesian_product_iterator([1:degrees(K)[i] for i in 1:ngens(K)], inplace = true)
@@ -999,114 +992,116 @@ The equation order of the number field.
 """
 equation_order(M::NfAbsOrd) = equation_order(nf(M))
 
+# Construct the smallest order of K containing the elements in elt.
+# If check == true, it is checked whether the given elements in elt are integral
+# and whether the constructed order is actually an order.
+# Via extends one may supply an order which will then be extended by the elements
+# in elt.
 function _order(K::S, elt::Vector{T}; cached::Bool = true, check::Bool = true, extends = nothing) where {S <: NumField{QQFieldElem}, T}
-  #=
-    check == true: the elements are known to be integral
-    extends !== nothing: then extends is an order, which we are extending
-  =#
+  elt = unique(elt)
   n = degree(K)
 
-  extending = false
-
-  local B::FakeFmpqMat = FakeFmpqMat()
-
   if extends !== nothing
     extended_order::order_type(K) = extends
     @assert K === nf(extended_order)
-    extend = true
 
     if is_maximal_known_and_maximal(extended_order) || length(elt) == 0
       return extended_order
     end
-    #in this case we can start with phase 2 directly as we have mult. closed
-    #module to start with, so set everything up for it...
     B = basis_matrix(extended_order)
     bas = basis(extended_order, K)
-    phase = 2
+    full_rank = true
+    m = _det_triangular(numerator(B, copy = false))//denominator(B, copy = false)
   else
+    if isempty(elt)
+      elt = elem_type(K)[one(K)]
+    end
     bas = elem_type(K)[one(K)]
-    phase = 1
+    B = basis_matrix(bas, FakeFmpqMat) # trivially in lower-left HNF
+    full_rank = false
   end
 
+  dummy_vector = elem_type(K)[K()]
+  function in_span_of_B(x::T)
+    if mod(denominator(B, copy = false), denominator(x)) == 0
+      dummy_vector[1] = x
+      C = basis_matrix(dummy_vector, FakeFmpqMat)
+      return is_zero_mod_hnf!(div(denominator(B, copy = false), denominator(x))*numerator(C, copy = false), numerator(B, copy = false))
+    end
+    return false
+  end
 
   for e in elt
-#    @show findall(isequal(e), elt)
-    if phase == 2
-      if denominator(B) % denominator(e) == 0
-        C = basis_matrix([e], FakeFmpqMat)
-        fl, _ = can_solve_with_solution(B.num, div(B.den, denominator(e))*C.num, side = :left)
-#        fl && println("elt known:", :e)
-        fl && continue
-      end
-    end
+    # Check if e is already in the multiplicatively closed module generated by
+    # the previous elements of elt
+    in_span_of_B(e) && continue
+
+    # Multiply powers of e to the existing basis elements
     if check
       f = minpoly(e)
-      isone(denominator(f)) || error("data does not define an order, $e is non-integral")
-      df = degree(f)-1
+      isone(denominator(f)) || error("The elements do not define an order: $e is non-integral")
+      df = degree(f) - 1
     else
-      df = n-1
+      df = n - 1
     end
-    f = one(K)
-    for i=1:df
-      mul!(f, f, e)
-      if phase == 2 # don't understand this part
-        if denominator(B) % denominator(f) == 0
-          C = basis_matrix(elem_type(K)[f], FakeFmpqMat)
-          fl = is_zero_mod_hnf!(div(B.den, denominator(f))*C.num, B.num)
-#          fl && println("inner abort: ", :e, " ^ ", i)
-          fl && break
-        end
-      end
-      if phase == 1
-        # [1] -> [1, e] -> [1, e, e, e^2] -> ... otherwise
-        push!(bas, deepcopy(f))
-      else
-        b = elem_type(K)[e*x for x in bas]
-        append!(bas, b)
+
+    start = 1
+    # We only multiply the elements of index start:length(bas) by e .
+    # Example: bas = [a_1, ..., a_k] with a_1 = 1. Then
+    # new_bas := [e, e*a_2, ..., e*a_k] and we append this to bas and set
+    # start := k + 1. In the next iteration, we then have
+    # new_bas := [e^2, e^2*a_2, ..., e^2*a_k] (assuming that there was no
+    # reduction of the basis in between).
+    for i in 1:df
+      new_bas = elem_type(K)[]
+      for j in start:length(bas)
+        t = e*bas[j]
+        in_span_of_B(t) && continue
+        push!(new_bas, t)
       end
+      isempty(new_bas) && break
+      start = length(bas) + 1
+      append!(bas, new_bas)
+
       if length(bas) >= n
+        # HNF reduce the basis we have so far, if B is already of full rank,
+        # we can do this with the modular algorithm
         B = basis_matrix(bas, FakeFmpqMat)
-        if extending
-          # We are extending extended_order, which has basis matrix M/d
-          # Thus we know that B.den/d * M \subseteq <B.num>
-          # So we can take B.den/d * largest_elementary_divisor(M) as the modulus
-          B = hnf_modular_eldiv(B, B.den, shape = :lowerleft)
+        if full_rank
+          # We have M \subseteq B, where M is a former incarnation of B.
+          # So we have B.den * M.num/M.den \subseteq B.num \subseteq Z^n, so
+          # M.d divides B.den and we can choose (B.den/M.den)*det(M.num) as
+          # modulus for the HNF of B.num.
+          mm = ZZ(m*denominator(B, copy = false))
+          hnf_modular_eldiv!(B, mm, shape = :lowerleft)
+          B = sub(B, nrows(B) - n + 1:nrows(B), 1:n)
+
+          # Check if we have a better modulus
+          new_m = _det_triangular(numerator(B, copy = false))//denominator(B, copy = false)
+          if new_m < m
+            m = new_m
+          end
         else
           hnf!(B)
+          k = findfirst(k -> !is_zero_row(B, k), nrows(B) - n + 1:nrows(B))
+          B = sub(B, nrows(B) - n + k:nrows(B), 1:n)
+          if nrows(B) == n
+            full_rank = true
+            m = _det_triangular(numerator(B, copy = false))//denominator(B, copy = false)
+          end
         end
-        rk = nrows(B) - n + 1
-        while is_zero_row(B, rk)
-          rk += 1
-        end
-        B = sub(B, rk:nrows(B), 1:n)
-        phase = 2
-        bas = elem_type(K)[ elem_from_mat_row(K, B.num, i, B.den) for i = 1:nrows(B) ]
+        bas = elem_type(K)[ elem_from_mat_row(K, numerator(B, copy = false), i, denominator(B, copy = false)) for i = 1:nrows(B) ]
+        start = 1
         if check
           @assert isone(bas[1])
         end
       end
     end
   end
-
-  if length(bas) > n # == n can only happen here after an hnf was computed
-                     # above. Don't quite see how > n can happen here either
-    B = basis_matrix(bas, FakeFmpqMat)
-    hnf!(B)
-    rk = nrows(B) - n + 1
-    if is_zero_row(B.num, rk)
-      error("data does not define an order: dimension to small")
-    end
-    B = sub(B, rk:nrows(B), 1:n)
-    bas = elem_type(K)[ elem_from_mat_row(K, B.num, i, B.den) for i = 1:nrows(B) ]
-  end
-
-  if !isdefined(B, :num)
-    error("data does not define an order: dimension to small")
+  if length(bas) < n
+    error("The elements do not define an order: rank too small")
   end
-
-  # Make an explicit check
-  @hassert :NfOrd 1 defines_order(K, B)[1]
-  return Order(K, B, cached = cached, check = check)
+  return Order(K, B, cached = cached, check = check)::order_type(K)
 end
 
 ################################################################################
diff --git a/test/Misc.jl b/test/Misc.jl
index 59ee89b5b3..c8795f52d6 100644
--- a/test/Misc.jl
+++ b/test/Misc.jl
@@ -12,4 +12,5 @@
   include("Misc/UnitsModM.jl")
   include("Misc/Matrix.jl")
   include("Misc/PIDIdeal.jl")
+  include("Misc/MSet.jl")
 end
diff --git a/test/Misc/MSet.jl b/test/Misc/MSet.jl
new file mode 100644
index 0000000000..b37b63bcaf
--- /dev/null
+++ b/test/Misc/MSet.jl
@@ -0,0 +1,131 @@
+@testset "Multi-sets" begin
+  str = "A nice sentence"
+  m = @inferred multiset(str)
+  io = IOBuffer()
+  show(io, MIME"text/plain"(), m)
+  @test length(String(take!(io))) == 92
+  show(io, m)
+  @test length(String(take!(io))) == 44
+  print(IOContext(io, :supercompact => true), m)
+  @test length(String(take!(io))) == 39
+
+  M = MSet(root_lattice(:A, i) for j in 1:10 for i in 1:100)
+  show(io, MIME"text/plain"(), M)
+  s = String(take!(io))
+
+  M = MSet{String}("$i"^100 for j in 1:4 for i in 1:130)
+  show(io, MIME"text/plain"(), M)
+  s = String(take!(io))
+
+  m = MSet{Int}()
+  show(io, MIME"text/plain"(), m)
+  @test length(String(take!(io))) == 13
+
+  m = @inferred multiset(Int[x^3%8 for x = 1:50])
+  @test !isempty(m)
+  @test length(m) == 50
+  @test eltype(m) == Int
+  @test !(2 in m)
+
+  m2 = copy(m)
+  @test m2 == m
+  @test unique(m) == unique(collect(m))
+
+  push!(m, 2, 4)
+  pop!(m, 2)
+  @test (2 in m)
+  @test multiplicity(m, 2) == 3
+
+  delete!(m, 2)
+  @test multiplicity(m, 2) == 0
+
+  @test pop!(m, 2, 0) == 0
+  k = pop!(m2)
+  @test multiplicity(m2, k) == multiplicity(m, k) - 1
+
+  m = @inferred multiset(Dict("a" => 4, "b" => 1, "c" => 9))
+  lis = @inferred collect(m)
+  @test length(m) == length(lis)
+
+  m2 = @inferred m + m
+  for i in m
+    @test multiplicity(m2, i) == 2*multiplicity(m, i)
+  end
+
+  m3 = @inferred m-m
+  @test length(m3) == 0
+
+  @test union(m) == m
+  @test length(filter(x -> multiplicity(m, x) != 1, m)) == length(m) - 1
+
+  m = multiset(Int[x^3%8 for x = 1:50])
+  filter!(isodd, m)
+  @test length(m) == 25
+
+  val = @inferred multiplicities(m)
+  @test all(x -> x > 1, val)
+
+  @test isempty(setdiff(m, m))
+
+  m = MSet(Dict("a" => 4, "b" => 1, "c" => 9))
+  @test length(setdiff!(m, unique(m), unique(m))) == 9
+  m3 = m + m + m
+  @test length(m3) == 27
+  m = multiset(Int[x^3%8 for x = 1:50])
+  @test length(union(m, m3)) == 77
+  @test union(m3, m3, m3) == m3
+
+  @test_throws ArgumentError union!(m, m3)
+
+  m1 = multiset(fill(3,4))
+  m2 = multiset(fill(2,6))
+  m3 = multiset(Int[2,2,3,3,4,4])
+  m4 = multiset(Int[3,4,4])
+  @test isempty(intersect(m1, m2, m3))
+  intersect!(m3, m1, m4)
+  @test length(m3) == 1
+end
+
+@testset "Sub-multi-set iterator" begin
+  m = @inferred multiset(["a", "b", "c"], [4, 1, 9])
+  M = @inferred subsets(m)
+  io = IOBuffer()
+  show(io, MIME"text/plain"(), M)
+  @test length(String(take!(io))) == 61
+  show(io, M)
+  @test length(String(take!(io))) == 64
+  print(IOContext(io, :supercompact => true), M)
+  @test length(String(take!(io))) == 35
+  @test eltype(M) == typeof(m)
+  @test length(collect(M)) == length(M)
+
+  n = collect(M)[end]
+  @test union(m, n) == m
+  @test intersect(m, n) == n
+end
+
+@testset "Sub-set iterators" begin
+  s = Set(String["a", "b", "c"])
+  S = @inferred subsets(s)
+  io = IOBuffer()
+  show(io, MIME"text/plain"(), S)
+  @test length(String(take!(io))) == 58
+  show(io, S)
+  @test length(String(take!(io))) == 58
+  print(IOContext(io, :supercompact => true), S)
+  @test length(String(take!(io))) == 29
+  @test eltype(S) == typeof(s)
+  @test length(collect(S)) == length(S)
+
+  S = @inferred subsets(s, 2)
+  io = IOBuffer()
+  show(io, MIME"text/plain"(), S)
+  @test length(String(take!(io))) == 68
+  show(io, S)
+  @test length(String(take!(io))) == 68
+  print(IOContext(io, :supercompact => true), S)
+  @test length(String(take!(io))) == 39
+  @test eltype(S) == typeof(s)
+  @test length(collect(S)) == length(S)
+  @test length(S[1]) == 2
+end
diff --git a/test/NfOrd/NfOrd.jl b/test/NfOrd/NfOrd.jl
index f0a2d2fedc..c3e7da1138 100644
--- a/test/NfOrd/NfOrd.jl
+++ b/test/NfOrd/NfOrd.jl
@@ -64,7 +64,7 @@
     #@test O7 == O77
     #@test !(O7 === O77)
 
-    O8 = Order(K6, [a1])
+    O8 = Order(K1, [a1])
     @test O8 == EquationOrder(K1)
 
     @test_throws ErrorException Order(K1, [a1, a1, a1], isbasis = true)
diff --git a/test/NumFieldOrd/NumFieldOrd.jl b/test/NumFieldOrd/NumFieldOrd.jl
index d7ddc37d7f..b93a568b9e 100644
--- a/test/NumFieldOrd/NumFieldOrd.jl
+++ b/test/NumFieldOrd/NumFieldOrd.jl
@@ -124,5 +124,20 @@ end
   @test extend(R, []) == R
   @test extend(R, [1//2 + a//2]) == maximal_order(K)
   @test extend(maximal_order(R), [a]) == maximal_order(R)
+
+  K, a = number_field(x, "a")
+  @test Order(K, [1]) == equation_order(K)
+  @test Order(K, []) == equation_order(K)
+
+  K, a = NumberField(x^4 - 10*x^2 + 1, "a")
+  x = 1//2*a^3 - 9//2*a # sqrt(2)
+  y = 1//2*a^3 - 11//2*a # sqrt(3)
+  O = Order(K, [x, y, x*y])
+  @test O == Order(K, [x, y])
+  @test O == Order(K, [x, y], check = false)
+  z = 1//4*a^3 + 1//4*a^2 + 3//4*a + 3//4
+  OO = Hecke._order(K, [z], extends = O)
+  @test is_maximal(OO)
+  @test_throws ErrorException Order(K, [x])
 end