Use 'FFLU' in SolveCtx for fraction fields

src/Matrix.jl

@@ -1578,6 +1578,51 @@ function *(P::Perm, x::MatrixElem{T}) where T <: NCRingElement
    return z
 end
 
+@doc raw"""
+    *(x::MatrixElem{T}, P::Perm) where T <: NCRingElement
+
+Apply the pemutation $P$ to the columns of the matrix $x$ and return the result.
+
+# Examples
+
+```jldoctest; setup = :(using AbstractAlgebra)
+julia> R, t = polynomial_ring(QQ, "t")
+(Univariate polynomial ring in t over rationals, t)
+
+julia> S = matrix_space(R, 3, 3)
+Matrix space of 3 rows and 3 columns
+  over univariate polynomial ring in t over rationals
+
+julia> G = SymmetricGroup(3)
+Full symmetric group over 3 elements
+
+julia> A = S([t + 1 t R(1); t^2 t t; R(-2) t + 2 t^2 + t + 1])
+[t + 1       t             1]
+[  t^2       t             t]
+[   -2   t + 2   t^2 + t + 1]
+
+julia> P = G([1, 3, 2])
+(2,3)
+
+julia> B = A*P
+[t + 1             1       t]
+[  t^2             t       t]
+[   -2   t^2 + t + 1   t + 2]
+
+```
+"""
+function *(x::MatrixElem{T}, P::Perm) where T <: NCRingElement
+   z = similar(x)
+   m = nrows(x)
+   n = ncols(x)
+   for i = 1:m
+      for j = 1:n
+         z[i, P[j]] = x[i, j]
+      end
+   end
+   return z
+end
+
 ###############################################################################
 #
 #   LU factorisation

src/Solve.jl

@@ -3,16 +3,22 @@ module Solve
 using AbstractAlgebra
 
 import AbstractAlgebra:
-  base_ring,
+  @attributes,
   _can_solve_with_solution_fflu,
   _can_solve_with_solution_interpolation,
+  _solve_fflu_precomp,
+  base_ring,
+  fflu!,
   Generic,
+  get_attribute,
+  has_attribute,
   kernel,
   matrix,
   nrows,
   ncols,
   PrettyPrinting,
-  rank
+  rank,
+  set_attribute!
 
 import Base: show
 
@@ -58,15 +64,18 @@ Base.similar(M::LazyTransposeMatElem, i::Int, j::Int) = lazy_transpose(similar(d
 #
 ################################################################################
 
-mutable struct SolveCtx{T, MatT, TranspMatT}
+# attribute storage is enable so that the fraction field types can store integral
+# matrices. These don't go in `red` or `red_transp` because they don't have the
+# correct type. Changing the type parameters would be a breaking change.
+@attributes mutable struct SolveCtx{T, MatT, TranspMatT}
   A::MatT # matrix giving the linear system
   red::MatT # reduced/canonical form of A (rref, hnf, lu)
   red_transp::TranspMatT # reduced/canonical form of transpose(A)
   trafo::MatT # transformation: trafo*A == red (not used for lu)
   trafo_transp::TranspMatT # transformation: trafo_transp*transpose(A) == red_transp
                            # (not used for lu)
-  lu_perm::Generic.Perm # permutation used for the lu factorization of A
-  lu_perm_transp::Generic.Perm # permutation used for the lu factorization of transpose(A)
+  lu_perm::Generic.Perm{Int} # permutation used for the lu factorization of A
+  lu_perm_transp::Generic.Perm{Int} # permutation used for the lu factorization of transpose(A)
 
   rank::Int # rank of A
   pivots::Vector{Int} # pivot and non-pivot columns of red
@@ -150,6 +159,31 @@ function _init_reduce_lu(C::SolveCtx{<:FieldElement})
   return nothing
 end
 
+function _init_reduce_fflu(C::SolveCtx{<:Union{FracElem, Rational{BigInt}}})
+  if has_attribute(C, :reduced_matrix_lu)
+    return nothing
+  end
+
+  A = matrix(C)
+  dA = _common_denominator(A)
+  Aint = matrix(parent(dA), nrows(A), ncols(A), [numerator(A[i, j]*dA) for i in 1:nrows(A) for j in 1:ncols(A)])
+  p = one(SymmetricGroup(nrows(A)))
+  r, dLU = fflu!(p, Aint)
+
+  set_rank!(C, r)
+  C.lu_perm = p
+  d = divexact(dA, base_ring(C)(dLU))
+  set_attribute!(C, :reduced_matrix_lu => Aint, :scaling_factor_fflu => d)
+  if r < nrows(A)
+    A2 = p*A
+    A3 = A2[r + 1:nrows(A), :]
+    set_attribute!(C, :permuted_matrix_fflu => A3)
+  else
+    set_attribute!(C, :permuted_matrix_fflu => zero(A, 0, ncols(A)))
+  end
+  return nothing
+end
+
 function _init_reduce(C::SolveCtx{<:RingElement})
   if isdefined(C, :red) && isdefined(C, :trafo)
     return nothing
@@ -187,6 +221,32 @@ function _init_reduce_transpose_lu(C::SolveCtx{<:FieldElement})
   return nothing
 end
 
+function _init_reduce_transpose_fflu(C::SolveCtx{<:Union{FracElem, Rational{BigInt}}})
+  if has_attribute(C, :reduced_matrix_of_transpose_lu)
+    return nothing
+  end
+
+  A = matrix(C)
+  dA = _common_denominator(A)
+  # We transpose A at this point!
+  Aint = matrix(parent(dA), ncols(A), nrows(A), [numerator(A[i, j]*dA) for j in 1:ncols(A) for i in 1:nrows(A)])
+  p = one(SymmetricGroup(nrows(Aint)))
+  r, dLU = fflu!(p, Aint)
+
+  set_rank!(C, r)
+  C.lu_perm_transp = p
+  d = divexact(dA, base_ring(C)(dLU))
+  set_attribute!(C, :reduced_matrix_of_transpose_lu => Aint, :scaling_factor_of_transpose_fflu => d)
+  if r < ncols(A)
+    A2 = A*p
+    A3 = A2[:, r + 1:ncols(A)]
+    set_attribute!(C, :permuted_matrix_of_transpose_fflu => A3)
+  else
+    set_attribute!(C, :permuted_matrix_of_transpose_fflu => zero(A, nrows(A), 0))
+  end
+  return nothing
+end
+
 function _init_reduce_transpose(C::SolveCtx{<:RingElement})
   if isdefined(C, :red_transp) && isdefined(C, :trafo_transp)
     return nothing
@@ -213,21 +273,78 @@ function reduced_matrix_lu(C::SolveCtx)
   return C.red
 end
 
+function reduced_matrix_lu(C::SolveCtx{<:FracElem{T}}) where T
+  _init_reduce_fflu(C)
+  return get_attribute(C, :reduced_matrix_lu)::dense_matrix_type(T)
+end
+
+function reduced_matrix_lu(C::SolveCtx{<:Rational{BigInt}})
+  _init_reduce_fflu(C)
+  return get_attribute(C, :reduced_matrix_lu)::dense_matrix_type(BigInt)
+end
+
 function reduced_matrix_of_transpose_lu(C::SolveCtx)
   _init_reduce_transpose_lu(C)
   return C.red_transp
 end
 
+function reduced_matrix_of_transpose_lu(C::SolveCtx{<:FracElem{T}}) where T
+  _init_reduce_transpose_fflu(C)
+  return get_attribute(C, :reduced_matrix_of_transpose_lu)::dense_matrix_type(T)
+end
+
+function reduced_matrix_of_transpose_lu(C::SolveCtx{<:Rational{BigInt}})
+  _init_reduce_transpose_fflu(C)
+  return get_attribute(C, :reduced_matrix_of_transpose_lu)::dense_matrix_type(BigInt)
+end
+
+# Factor by which any solution needs to be multiplied.
+# This is the chosen denominator of matrix(C) divided by the denominator returned
+# by fflu!(matrix(C)).
+function scaling_factor_fflu(C::SolveCtx{T}) where {T <: Union{FracElem, Rational{BigInt}}}
+  _init_reduce_fflu(C)
+  return get_attribute(C, :scaling_factor_fflu)::T
+end
+
+function scaling_factor_of_transpose_fflu(C::SolveCtx{T}) where {T <: Union{FracElem, Rational{BigInt}}}
+  _init_reduce_transpose_fflu(C)
+  return get_attribute(C, :scaling_factor_of_transpose_fflu)::T
+end
+
+# Let A = matrix(C).
+# Return the matrix (p*A)[rank(A) + 1:nrows(A), :] where p is lu_permutation(C).
+function permuted_matrix_fflu(C::SolveCtx{FldT, MatT}) where {FldT <: Union{FracElem, Rational{BigInt}}, MatT}
+  _init_reduce_fflu(C)
+  return get_attribute(C, :permuted_matrix_fflu)::MatT
+end
+
+# Let A = matrix(C).
+# Return the matrix (A*p)[:, rank(A) + 1:ncols(A)] where p is lu_permutation_of_transpose(C).
+function permuted_matrix_of_transpose_fflu(C::SolveCtx{FldT, MatT}) where {FldT <: Union{FracElem, Rational{BigInt}}, MatT}
+  _init_reduce_transpose_fflu(C)
+  return get_attribute(C, :permuted_matrix_of_transpose_fflu)::MatT
+end
+
 function lu_permutation(C::SolveCtx)
   _init_reduce_lu(C)
   return C.lu_perm
 end
 
+function lu_permutation(C::SolveCtx{<:Union{FracElem, Rational{BigInt}}})
+  _init_reduce_fflu(C)
+  return C.lu_perm
+end
+
 function lu_permutation_of_transpose(C::SolveCtx)
   _init_reduce_transpose_lu(C)
   return C.lu_perm_transp
 end
 
+function lu_permutation_of_transpose(C::SolveCtx{<:Union{FracElem, Rational{BigInt}}})
+  _init_reduce_transpose_fflu(C)
+  return C.lu_perm_transp
+end
+
 function transformation_matrix(C::SolveCtx)
   _init_reduce(C)
   return C.trafo
@@ -248,11 +365,20 @@ end
 
 function AbstractAlgebra.rank(C::SolveCtx{<:FieldElement})
   if C.rank < 0
+    # We should be calling _init_reduce_lu here, but can't because it has to stay
+    # compatible with Nemo.
     _init_reduce(C)
   end
   return C.rank
 end
 
+function AbstractAlgebra.rank(C::SolveCtx{<:Union{FracElem, Rational{BigInt}}})
+  if C.rank < 0
+    _init_reduce_fflu(C)
+  end
+  return C.rank
+end
+
 AbstractAlgebra.nrows(C::SolveCtx) = nrows(matrix(C))
 AbstractAlgebra.ncols(C::SolveCtx) = ncols(matrix(C))
 AbstractAlgebra.base_ring(C::SolveCtx) = base_ring(matrix(C))
@@ -668,6 +794,71 @@ function _can_solve_internal_no_check(A::MatElem{T}, b::MatElem{T}, task::Symbol
   end
 end
 
+function _can_solve_internal_no_check(C::SolveCtx{T}, b::MatElem{T}, task::Symbol; side::Symbol = :left) where T <: Union{FracElem, Rational{BigInt}}
+  # Split up in separate functions to make the compiler happy
+  if side === :right
+    return _can_solve_internal_no_check_right(C, b, task)
+  else
+    return _can_solve_internal_no_check_left(C, b, task)
+  end
+end
+
+function _can_solve_internal_no_check_right(C::SolveCtx{T}, b::MatElem{T}, task::Symbol) where T <: Union{FracElem, Rational{BigInt}}
+  K = base_ring(C)
+  db = _common_denominator(b)
+  bint = matrix(parent(db), nrows(b), ncols(b), [numerator(b[i, j]*db) for i in 1:nrows(b) for j in 1:ncols(b)])
+  fl, y_int = _solve_fflu_precomp(lu_permutation(C), reduced_matrix_lu(C), bint)
+  if !fl
+    return fl, zero(b, 0, 0), zero(b, 0, 0)
+  end
+  # We have fl == true, but we still have to check whether this really is a solution
+  d = scaling_factor_fflu(C)
+  y = change_base_ring(K, y_int)
+  y = y*divexact(d, K(db))
+  if rank(C) < nrows(C)
+    # We have to check whether y is also a solution for the "lower part"
+    # of the system
+    pb = lu_permutation(C)*b
+    pA = permuted_matrix_fflu(C)
+    fl = pA*y == pb[rank(C) + 1:nrows(C), :]
+  end
+  if task === :with_kernel
+    # I don't know how to compute the kernel using an (ff)lu factoring
+    return fl, y, kernel(C, side = :right)
+  else
+    return fl, y, zero(b, 0, 0)
+  end
+end
+
+function _can_solve_internal_no_check_left(C::SolveCtx{T}, b::MatElem{T}, task::Symbol) where T <: Union{FracElem, Rational{BigInt}}
+  K = base_ring(C)
+  db = _common_denominator(b)
+  # bint == transpose(b)*db
+  bint = matrix(parent(db), ncols(b), nrows(b), [numerator(b[i, j]*db) for j in 1:ncols(b) for i in 1:nrows(b)])
+  fl, y_int = _solve_fflu_precomp(lu_permutation_of_transpose(C), reduced_matrix_of_transpose_lu(C), bint)
+  if !fl
+    return fl, zero(b, 0, 0), zero(b, 0, 0)
+  end
+  # We have fl == true, but we still have to check whether this really is a solution
+  d = scaling_factor_of_transpose_fflu(C)
+  # transpose y_int
+  y = matrix(K, ncols(y_int), nrows(y_int), [ K(y_int[i, j]) for j in 1:ncols(y_int) for i in 1:nrows(y_int)])
+  y = y*divexact(d, K(db))
+  if rank(C) < ncols(C)
+    # We have to check whether y is also a solution for the "right hand part"
+    # of the system
+    pb = b*lu_permutation_of_transpose(C)
+    pA = permuted_matrix_of_transpose_fflu(C)
+    fl = y*pA == pb[:, rank(C) + 1:ncols(C)]
+  end
+  if task === :with_kernel
+    # I don't know how to compute the kernel using an (ff)lu factoring
+    return fl, y, kernel(C, side = :left)
+  else
+    return fl, y, zero(b, 0, 0)
+  end
+end
+
 function _can_solve_internal_no_check(A::MatElem{T}, b::MatElem{T}, task::Symbol; side::Symbol = :left) where T <: FracElem{TT} where TT <: PolyRingElem
 
   if side === :left