solve for general zzModMatrix

[martinra]

Currently, the code in Flint is only guaranteed to work for prime modulus. I needed more and have implemented it. Currently, I don't want to take the time to polish this up (but if nobody comes before me, I'll do this once my current priority is completed in N months for some large N). In the meantime, I'll leave the code for reference and for discussion. Note that there is a similar issue about matrix inversion #956, which can be resolved with this code also.

Questions that I have:

• Flint sometimes succeeds solving general systems and then is much faster. Should we first try to call flint and only run an alternative if this does not work?
• Is is desirable at all to have this in Nemo or should this go to flint. In N+N' months I could also attend to this.

function AbstractAlgebra.Solve._can_solve_internal_no_check(
    A::zzModMatrix,
    b::zzModMatrix,
    task::Symbol;
    side::Symbol = :left
    )
  if side === :left
    fl, sol, K = AbstractAlgebra.Solve._can_solve_internal_no_check(
                     transpose(A), transpose(b), task, side = :right)
    return fl, transpose(sol), transpose(K)
  end

  R = base_ring(A)
  ps = [(UInt64(p),m) for (p,m) in factor(UInt64(characteristic(R)))]

  if length(ps) == 1 && ps[1][2] == 1
    x = similar(A, ncols(A), ncols(b))
    # This is only guaranteed to terminate sucessfully if the characteristic is
    # prime
    fl = ccall((:nmod_mat_can_solve, libflint), Cint,
               (Ref{zzModMatrix}, Ref{zzModMatrix}, Ref{zzModMatrix}),
               x, A, b)
    if task === :only_check || task === :with_solution
      return Bool(fl), x, zero(A, 0, 0)
    end
    return Bool(fl), x, kernel(A, side = :right)
  else
    @assert task === :only_check || task === :with_solution

    r = nrows(A)
    ca = ncols(A)
    cb = ncols(b)
    caprev = 0
    Ks = []
    xorig = zero(b, ca, cb)

    modacc = UInt64(1)
    while true
      p = first(ps)[1]
      if ps[1][2] == 1
        popfirst!(ps)
      else
        ps[1] = (p,ps[1][2]-1)
      end
      F,_ = quo(ZZ,p)

      Ared = matrix(F, [lift(A[ix,jx]) for ix in 1:r, jx in 1:ca])
      bred = matrix(F, [lift(b[ix,jx]) for ix in 1:r, jx in 1:cb])
      flag, xred, Kred = can_solve_with_solution_and_kernel(Ared, bred; side)

      if !flag
        return false, xorig, zero(A, 0, 0)
      end

      x = matrix(R, [lift(xred[ix,jx]) for ix in 1:ca, jx in 1:cb])
      K = matrix(R, [lift(Kred[ix,jx]) for ix in 1:ca, jx in 1:ncols(Kred)])
      if modacc == 1
        xorig = x
      else
        xcontracted = x
        for Kprev in Ks
          xcontracted = 
            ( xcontracted[1:end-ncols(Kprev),:]
            + Kprev * x[end-ncols(Kprev)+1:end,:]
            )
        end
        xorig = xorig + modacc*xcontracted
      end

      if flag && isempty(ps)
        return true, xorig, zero(A, 0, 0)
      end
      
      # A heuristic to keep the matrix A from growing excessively.
      Knz = [jx for jx in 1:ncols(K) if any(!iszero, K[caprev+1:ca,jx])]
      K = K[:,Knz]
      AK = A*K
      AKnz = [jx for jx in 1:ncols(AK) if any(!iszero, AK[:,jx])]
      K = K[:,AKnz]
      AK = AK[:,AKnz]

      b = map(a -> divexact(a,p), b - A * x)
      A = hcat(A, map(a -> divexact(a,p), AK))
      caprev = ca
      ca = ncols(A)
      pushfirst!(Ks, K)

      modacc *= p
    end

    error("unreachable")
  end
end

A little test

mod = 13^4 * 23^2
R,_ = quo(ZZ, mod)
a = matrix(R, [rand(1:mod) for ix in 1:4, jx in 1:7])
x = matrix(R, [rand(1:mod) for ix in 1:3, jx in 1:4])
b = x*a
fl, xsol = can_solve_with_solution(a, b; side = :left)
b == xsol*a

[fredrik-johansson]

Is is desirable at all to have this in Nemo or should this go to flint. In N+N' months I could also attend to this.

We definitely want to have this in Flint!

[thofma]

@joschmitt has been implemented solving stuff for arbitrary PIRs that impelement the Howell form. Since this exists, we should just make use of it.

[thofma]

It is basically this function here: https://github.com/thofma/Hecke.jl/blob/0cf400dfaa770a309c2f0a2c57cb881a8e53fce5/src/NumFieldOrd/NfOrd/StrongEchelonForm.jl#L276 This works at the moment for quotient rings of rings of integers or valuation rings of local fields. But we discussed moving this to up to AbstractAlgebra so that any ring implementing howell_form can hook into this.

[thofma]

@martinra we will cook something up, which will make it work for zzModMatrix and ZZModMatrix (and maybe many more rings!).

[martinra]

Thanks! That's fantastic.

[thofma]

Should be working now. Let us know if something does not work.