Bugs in Order(...)

[fagu]

Here are some bugs in the current Order(...) function, a discussion of the algorithm and some half-baked ideas. I could try to implement a fix, but wanted to first ask for opinions about what the algorithm should actually be...

Issues

There are several (more or less separate) bugs in the Order(...) function that in a given number field K should compute the order generated by a given list of elements.

using Oscar
R,x=QQ["x"]

Issue 1: Running K,w=NumberField(x,"w"); Order(K,[1]) should produce the order Z in the number field Q, but gives:

ERROR: data does not define an order: dimension to small
Stacktrace:
 [1] _order(K::AnticNumberField, elt::Vector{nf_elem}; cached::Bool, check::Bool, extends::Nothing)
   @ Hecke ~/.julia/packages/Hecke/D2kay/src/NumFieldOrd/NfOrd/NfOrd.jl:1097
 [2] _order
   @ ~/.julia/packages/Hecke/D2kay/src/NumFieldOrd/NfOrd/NfOrd.jl:995 [inlined]
 [3] Order(::AnticNumberField, a::Vector{nf_elem}; check::Bool, isbasis::Bool, cached::Bool)
   @ Hecke ~/.julia/packages/Hecke/D2kay/src/NumFieldOrd/NfOrd/NfOrd.jl:766
 [4] Order
   @ ~/.julia/packages/Hecke/D2kay/src/NumFieldOrd/NfOrd/NfOrd.jl:750 [inlined]
 [5] Order(K::AnticNumberField, a::Vector{Int64}; check::Bool, isbasis::Bool, cached::Bool)
   @ Hecke ~/.julia/packages/Hecke/D2kay/src/NumFieldOrd/NfOrd/NfOrd.jl:778
 [6] Order(K::AnticNumberField, a::Vector{Int64})
   @ Hecke ~/.julia/packages/Hecke/D2kay/src/NumFieldOrd/NfOrd/NfOrd.jl:770
 [7] top-level scope
   @ REPL[7]:1

The problem is that although bas is initialized to the basis 1 of Q, the matrix B still has zero rows:

Hecke.jl/src/NumFieldOrd/NfOrd/NfOrd.jl

Lines 1027 to 1028 in d23dfff

	bas = elem_type(K)[one(K)]
	phase = 1

Issue 2: Running K,w=NumberField(x,"w"); Order(K,[]) should also produce the order Z in the number field Q, but gives:

ERROR: StackOverflowError:
Stacktrace:
 [1] Order(K::AnticNumberField, a::Vector{Any}; check::Bool, isbasis::Bool, cached::Bool)
   @ Hecke ~/.julia/packages/Hecke/D2kay/src/NumFieldOrd/NfOrd/NfOrd.jl:773
 [2] Order(K::AnticNumberField, a::Vector{Any}; check::Bool, isbasis::Bool, cached::Bool) (repeats 16360 times)
   @ Hecke ~/.julia/packages/Hecke/D2kay/src/NumFieldOrd/NfOrd/NfOrd.jl:778
 [3] Order(K::AnticNumberField, a::Vector{Any})
   @ Hecke ~/.julia/packages/Hecke/D2kay/src/NumFieldOrd/NfOrd/NfOrd.jl:770
 [4] top-level scope
   @ REPL[8]:1

I think this is a type confusion. Somehow, map(K, []) returns an object of type Vector{Any} instead of Vector{nf_elem}. (Why?) That's why this function

Hecke.jl/src/NumFieldOrd/NfOrd/NfOrd.jl

Lines 770 to 779 in d23dfff

	function Order(K, a::Vector; check::Bool = true, isbasis::Bool = false,
	cached::Bool = true)
	local b
	try
	b = map(K, a)
	catch
	error("Cannot coerce elements from array into the number field")
	end
	return Order(K, b, check = check, cached = cached, isbasis = isbasis)
	end

keeps calling itself instead of calling

Hecke.jl/src/NumFieldOrd/NfOrd/NfOrd.jl

Lines 750 to 751 in d23dfff

	function Order(::S, a::Vector{T}; check::Bool = true, isbasis::Bool = false,
	cached::Bool = false) where {S <: NumField{QQFieldElem}, T <: NumFieldElem{QQFieldElem}}

Issue 3: Running K,w=NumberField(x^4-10*x^2+1,"w"); a=1//2*w^3-9//2*w; b=1//2*w^3-11//2*w; Order(K,[a,b]) should compute the field $K=\mathbb Q(\sqrt2,\sqrt3)$, the elements $a=\sqrt2$, $b=\sqrt3$, and return the order $\mathbb Z[a,b]$ of $K$. Instead, it fails:

ERROR: data does not define an order: dimension to small
Stacktrace:
 [1] _order(K::AnticNumberField, elt::Vector{nf_elem}; cached::Bool, check::Bool, extends::Nothing)
   @ Hecke ~/.julia/packages/Hecke/D2kay/src/NumFieldOrd/NfOrd/NfOrd.jl:1097
 [2] _order
   @ ~/.julia/packages/Hecke/D2kay/src/NumFieldOrd/NfOrd/NfOrd.jl:995 [inlined]
 [3] Order(::AnticNumberField, a::Vector{nf_elem}; check::Bool, isbasis::Bool, cached::Bool)
   @ Hecke ~/.julia/packages/Hecke/D2kay/src/NumFieldOrd/NfOrd/NfOrd.jl:766
 [4] Order(::AnticNumberField, a::Vector{nf_elem})
   @ Hecke ~/.julia/packages/Hecke/D2kay/src/NumFieldOrd/NfOrd/NfOrd.jl:750
 [5] top-level scope
   @ REPL[34]:1

(Both Order(K,[a,b,a*b]) and Order(K,[a,b],check=false) actually succeed and produce the correct output.)
In this function, when check=true, the function _order(...) actually never reaches "phase 2". It never computes the product a*b, only produces the Z-module generated by 1,a,b, which is indeed not an order. I'm confused about the purpose of the whole phase management. (See below.)

Issue 4: This coincidentally doesn't produce any serious problems, but I think there's a mistyped variable name. Compare these lines:

Hecke.jl/src/NumFieldOrd/NfOrd/NfOrd.jl

Line 1009 in d23dfff

extending = false

Hecke.jl/src/NumFieldOrd/NfOrd/NfOrd.jl

Line 1016 in d23dfff

extend = true

This branch, which would use the modular Hermite normal form optimization, is never taken:

Hecke.jl/src/NumFieldOrd/NfOrd/NfOrd.jl

Line 1069 in d23dfff

if extending

Regardless of this typo, I don't understand why this optimization should be used only when we're extending an order. (See below.)

There are also some unnecessary checks in the code, as was already pointed out by @fieker in inline comments:

Hecke.jl/src/NumFieldOrd/NfOrd/NfOrd.jl

Line 1052 in d23dfff

if phase == 2 # don't understand this part

Hecke.jl/src/NumFieldOrd/NfOrd/NfOrd.jl

Lines 1091 to 1092 in d23dfff

	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

Fixed algorithm

I don't know if the current algorithm is documented anywhere, so here's an outline in pseudocode for what I think the function was trying to do to compute the order generated by the elements in a given vector gens, where the number field K has degree n:

A and B will be Z-submodules of K.
A := 1*Z (or whatever order we're extending)
for g in gens:
    B := A
    for i in 1,2,...:
        if i == n and not check:
            break
        if g^i in B:
            break
        if i == n:
            raise error(g is not an algebraic integer)
        B := B + g^i * A
        # At this point, B = A + g*A + ... + g^i*A
    A := B
    # At this point, A is the subring of K generated by the generators up to g.
return A

The Z-submodules of K are represented by Hermite-reduced bases. One can use a modular HNF algorithm (the optimization mentioned above) as soon as B has full rank. (Using "the determinant of B" as the modulus.)

If you append k rows to an $n\times n$ Hermite reduced matrix and apply Hermite reduction to the result, I think the algorithm might do about $n^3 + kn^2$ operations in Z. (The summand $n^3$ arises because all diagonal entries might decrease, and then all the off-diagonal entries have to be reduced.)

Hence, it is faster by a factor of ~n to append n rows and Hermite reduce once than to n times append one row and reduce. The attempt to implement such an optimization is the root of all evil in the current implementation. :)

One easy way to solve this issue (if one is willing to sacrifice a constant factor in the running time) should be to represent a Z-module by a list of generators which is not necessarily Hermite reduced. Only apply Hermite reduction after line 12 if we have accumulated $\geq2n$ generators, and once at the end of the algorithm. It's okay to just skip the checks in lines 8-11 if the current list of generators of B hasn't been reduced, yet.

More efficient methods?

I feel like there must be more efficient algorithms for this problem, at least when the number r of generators is large. Consider the total number N of rows in all calls to the Hermine normalization function hnf.

If you run the above algorithm (or the currently implemented one) with r >= 2 generators, I suppose N could be on the order of $rn^2$. (After the first generator, A can have rank n, so we then have to append ~n rows, for i=1,...,n, for each generator.)

Here's an algorithm with $N=\mathcal O(r+n^2\log(rn))$:

A := Z-module generated by 1 and the given generators
for i in 1,2,...:
    B := A*A
    if B == A:
        break
    if 2^(i-1) >= r*n:
        raise error(not all algebraic integers)
    A := B
    # At this point, A is the Z-module generated by all products g_1 ... g_k of k <= 2^i of the given generators.
return A

(Computing A*A involves a matrix with $\leq n^2$ rows.)

[fagu]

Another idea: First, run one of the algorithms with vector spaces over Q to find a set of n linearly independent elements of the order. That gives you a finite index submodule of the order, allowing you to use modular Hermite reduction from the beginning.
It could then make sense to relax the notion of Hermite reduced matrix as follows: entries don't need to be reduced modulo the "diagonal" entry in the same column, only modulo the determinant bound. This should get the number of operations in Hermite reduction down to $\mathcal O(kn^2)$, so it might be okay to reduce each time you add a vector. (?)

[fagu]

Maybe there's a probabilistic (Las Vegas) algorithm using the following lemma that's faster than the one described above under "More efficient methods?" by a factor of $\log(rn)$, and an even faster Monte Carlo algorithm.

Lemma. Let G be a finite abelian group of order s. Then, k random elements of G generate the group G with probability $\geq 1 - s/2^k$.
Proof See Hafner, McCurley: A rigorous subexponential algorithm for computation of class groups.

So if you have a rank n module $A$ containing $1$ and you want to compute $A\cdot A$, knowing that $[A\cdot A:A]\leq s$, instead of adding all $\sim n^2$ products of pairs of basis elements of $A$, adding $k$ random elements of $A\cdot A\mod s$ to $A$ will produce $A\cdot A$ with probability $\geq 1 - s/2^k$. So adding $\log_2s + \mathcal O(1)$ elements is enough to give $A\cdot A$ with large probability.

To obtain a Las Vegas algorithm with guaranteed correctness, at the end of the algorithm, check that $A\cdot A=A$ by adding all products of pairs of basis elements. If that fails, you could restart the algorithm.

I haven't tried to rigorously analyze the running time. Of course, one would need to also keep track of the number of digits of all numbers involved in the computation.

[thofma]

Thanks for the report and the analysis. I have a feeling phase = 1 refers to "not full rank yet" and phase = 2 to "full rank". Why some of the things should be done only in the first case, I do not know.

I will try to cook up a simpler implementation to fix the bugs that you reported.

I am also currently collecting all Order(K, [...]) calls that occur in the test suite. This should be helpful when trying to improve the algorithm.

[fieker]

On Mon, Sep 25, 2023 at 01:00:41AM -0700, Tommy Hofmann wrote: Thanks for the report and the analysis. I have a feeling `phase = 1` refers to "not full rank yet" and `phase = 2` to "full rank". Why some of the things should be done only in the first case, I do not know.

In phase 1 we're adding also powers of elements (possibly all n in one go), in phase only products are added. Rationale: if you're creating an equation order, adding (and HNFing) 1, a, then a^2, ... is painful. Keeping track of redundant products is also hard

…

I will try to cook up a simpler implementation to fix the bugs that you reported. I am also currently collecting all `Order(K, [...])` calls that occur in the test suite. This should be helpful when trying to improve the algorithm. -- Reply to this email directly or view it on GitHub: #1223 (comment) You are receiving this because you were mentioned. Message ID: ***@***.***>