Some corrections for algo remplacing LLL

framework/geometry.sage

@@ -75,6 +75,17 @@ def cannonical_param(v):
     i = [x != 0 for x in v[0]].index(True)
     return i, v[0, i]
 
+def xgcd_of_list(a):
+    (g, s, t) = xgcd(a[0], a[1])
+    bezouts = [s, t]
+
+    for v in a[2:]:
+        (g, s_, t_) = xgcd(g, v)
+        bezouts = [b*s_ for b in bezouts]
+        bezouts.append(t_)
+
+    return (g, bezouts)
+
 def remove_linear_dependencies(B, dim=None):
     nrows = B.nrows()
     if dim is None or nrows > dim:
@@ -85,22 +96,46 @@ def remove_linear_dependencies(B, dim=None):
             r = K.dimensions()[0]
         else:
             r = nrows-dim
-        
-        if r == 1:
+
+        if r == 1 and False:
+            print("Use Better Algo")
             # Find a linear dependency
             if K is None:
                 K = B.left_kernel().basis_matrix()
             assert K.dimensions()[0] == 1
             combinaison = K[0]
-            
-            # Detect the redundant vector
+            combinaison *= lcm([v.denominator() for v in combinaison])
+            print(combinaison)
+
             pivot, pivot_value = None, None
             for ind, value in enumerate(combinaison):
-                if abs(value) > 0 and (pivot is None or abs(value)<abs(pivot_value)):
+                if abs(value) == 1:
                     pivot, pivot_value = ind, value
-            
+                    break
+
+            if pivot_value is None:
+                print('Complex case')
+                for ind, value in enumerate(combinaison):
+                    if abs(value) > 0 and gcd([v for i,v in enumerate(combinaison) if i!=ind]) == 1:
+                        if pivot is None or abs(value)<abs(pivot_value):
+                            pivot, pivot_value = ind, value
+
+                if pivot_value < 0:
+                    combinaison = vector(ZZ, [-v for v in combinaison])
+
+                _, bezouts = xgcd_of_list([v for i,v in enumerate(combinaison) if i!=pivot])
+
+                factor = combinaison[pivot]-1
+                for i in range(len(combinaison)):
+                    ind = i if i < pivot else i-1
+                    if i != pivot:
+                        B[i] += factor*bezouts[ind]*B[pivot]
+                combinaison[pivot] += -factor
+                assert (combinaison*B).is_zero(), 'It is not a linear dependency anymore !'
+
+            assert abs(combinaison[pivot]) == 1, f'Error abs(combinaison[pivot]) == {abs(combinaison[pivot])} != 1'
             B = B[[i for i in range(B.dimensions()[0]) if i != pivot]]
-            
+
         else:
             B = B.LLL()
             B = B[r:]
@@ -129,7 +164,7 @@ def lattice_orthogonal_section(D, V, maintains_basis=True):
     # Eliminate linear dependencies
     if maintains_basis:
         D = remove_linear_dependencies(D)
-    
+
     # Go back to the primal
     return D
 
@@ -137,7 +172,7 @@ def lattice_orthogonal_section(D, V, maintains_basis=True):
 def lattice_project_against(B, V, maintains_basis=True):
     """
     Compute the projection of the lattice L(B) orthogonally to Span(V). All vectors if V
-    (or at least their projection on Span(B)) must belong to L(B). 
+    (or at least their projection on Span(B)) must belong to L(B).
     Algorithm:
     - project V onto Span(B)
     - project the basis onto orth(V)
@@ -200,7 +235,7 @@ def is_diagonal(M):
 
 
 def logdet(M, exact=False):
-    """ 
+    """
     Compute the log of the determinant of a large rational matrix,
     tryping to avoid overflows.
     """
@@ -221,7 +256,7 @@ def logdet(M, exact=False):
 
 
 def degen_inverse(S, B=None):
-    """ Compute the inverse of a symmetric matrix restricted 
+    """ Compute the inverse of a symmetric matrix restricted
     to its span
     """
     # Get an orthogonal basis for the Span of B
@@ -294,13 +329,13 @@ def square_root_inverse_degen(S, B=None):
 
 
 def build_standard_substitution_matrix(V, pivot=None, data=None, output_data=False):
-    _, pivot = V.nonzero_positions()[0] if (pivot is None) else (None, pivot) 
+    _, pivot = V.nonzero_positions()[0] if (pivot is None) else (None, pivot)
     assert V[0,pivot] != 0, 'The value of the pivot must be non-zero.'
     dim = V.ncols()
-    
+
     V1 = - V[0,:pivot] / V[0,pivot]
     V2 = - V[0,pivot+1:] / V[0,pivot]
-    
+
     Gamma = zero_matrix(QQ, dim,dim-1)
     Gamma[:pivot,:pivot] = identity_matrix(pivot)
     Gamma[pivot,:pivot] = V1
@@ -315,11 +350,11 @@ def build_standard_substitution_matrix(V, pivot=None, data=None, output_data=Fal
         normalization_matrix[pivot:,pivot:] = identity_matrix(dim-pivot-1) - V2.T*V2 / norm2
         normalization_matrix[:pivot,pivot:] = - V1.T*V2 / norm2
         normalization_matrix[pivot:,:pivot] = - V2.T*V1 / norm2
-        
+
         data['det'] = norm2
         data['normalization_matrix'] = normalization_matrix
 
     if output_data:
         return Gamma, data
     else:
-        return Gamma
+        return Gamma