Squashed commit of the following:

commit 401a711
Author: Johannes Schmitt <schmitt@mathematik.uni-kl.de>
Date:   Fri Oct 6 15:47:13 2023 +0200

    Fix `_order(::NumField, ::Vector{NumFieldElem})`

commit 48c9923
Author: Tommy Hofmann <thofma@gmail.com>
Date:   Mon Sep 25 09:19:26 2023 +0200

    Fixing Order([...])

Author: joschmitt

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

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,20 +886,13 @@ 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
       normalized_gens[i] = coeff(f, 1) * g[i]
     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
 
 ################################################################################

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)

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