# Modreduce
#
# Author: Howard Cheng
#
# License: CC BY-SA (http://creativecommons.org/licenses/by-sa/4.0/)
#
# The software is provided "AS IS" and there is no warranty of any kind.
#
# Implements the modular algorithm in
#
# Cheng, H. and Labahn, G. "Modular Computation for Matrices of Ore 
#   Polynomials." Computer Algebra 2006: Latest Advances in Symbolic 
#   Algorithms, pages 43-66, 2007. 
#
# Input:
#
#    F    : a matrix of skew polynomials over Z[t][Z;sigma,delta], where the
#           powers of Z are assumed to be on the right
#    Z    : the indeterminate
#    sigma: the automorphism, implemented as a function: (Z[t], t) -> Z[t]
#    delta: the derivation, implemented as a function: (Z[t], t) -> Z[t]
#    t    : the variable in the domain Z[t]
#
# Output:
#
#    M and R: where
#
#          M . F = R . Z^omega;
#
#    - the number of nonzero rows in R = rank of F;
#    - the rows in M corresponding to the zero rows in R form a basis for
#      the left nullspace of F
#
#   omega:
#          the final order achieved (always (mN+1) \cdot \vec e)
#   mu:
#          the final degree constraints achieved
#

Modreduce := proc(F :: Matrix, Z :: name, sigma, delta, t :: name)
  local M, R, m, n, N, omega, mu, J, 
        M_p, R_p, omega_p, mu_p, J_p, bound,
        smu, smu_p, tau, p, prod, i, j, unchanged, q, newM, newR, T, kappa,
        row;

  option `Copyright (c) 2005 by Howard Cheng`;

  # replace mod with symmetric mod
  `mod` := mods;

  # Initialization
  m, n := LinearAlgebra:-Dimensions(F);
  N := max(seq(seq(degree(F[i,j],Z),i=1..m),j=1..n));
  T, kappa := ComputeTKappa(F, Z, sigma, delta, t);

  userinfo(3, Modreduce, `T = `, T);
  userinfo(3, Modreduce, `kappa = `, kappa);
  
  # compute bound for early termination

  bound := ((m*N+1)*n)*kappa*(T+1);

  tau := 0;
  prod := 1;
  p := 30000;   # We start with a relatively large number
  while prod < bound do
    tau := tau+1;
    p := nextprime(p);
    prod := prod * p;
  od;

  userinfo(3, Modreduce, `bound = `, bound);
  userinfo(3, Modreduce, `tau = `, tau);

  # initialize before entering loop
  p := 30000;
  unchanged := 0;
  omega := Vector[row](n);
  mu := Vector[row](m);
  J := [];
  q := 1;                 # modulus for what has been reconstructed so far  

  # now do the loop
  while true do
    p := nextprime(p);

    userinfo(3, Modreduce, `Processing p = `, p);

    # check for unlucky primes
    if irem(lcoeff(sigma(t, t), t), p) = 0 then
      userinfo(3, Modreduce, `  unlucky: hc(sigma(t)) = 0 mod p`);
      next;
    fi;

    M_p, R_p, omega_p, mu_p, J_p := ReduceZpt(F, m, n, p, N, T, Z, 
                                              sigma, delta, t);

    if q > 1 and not LinearAlgebra:-Equal(omega, omega_p) then
      error `omega is not the same for different primes!  omega = `, omega,
            `omega_p = `, omega_p;
    fi;

    # now check for unlucky homomorphisms
    if (q = 1) or
       (LinearAlgebra:-Equal(mu, mu_p) and LinearAlgebra:-Equal(J, J_p)) then
      # seems to be lucky
      userinfo(3, Modreduce, `  seems to be lucky`);
      
      if q = 1 then
        mu, J, M, R, omega := mu_p, J_p, M_p, R_p, omega_p;
      else
        newM := MatrixCRA(M, M_p, q, p, Z);
        newR := MatrixCRA(R, R_p, q, p, Z);
        M_p, R_p := 0, 0;

        # check for early termination
        if LinearAlgebra:-Equal(M, newM) and LinearAlgebra:-Equal(R, newR) then
          unchanged := unchanged + 1;
          newM, newR := 0, 0;
        else
          unchanged := 0;
          M, R := newM, newR;
        fi;
        if unchanged >= tau then
          userinfo(3, Modreduce, `  early termination`);
          break;
        fi;
      fi;
      q := q*p;
    else
      # one of the is unlucky, figure out which one
      userinfo(3, Modreduce, `  unlucky prime(s) detected`);
      smu := add(mu[i], i=1..m);
      smu_p := add(mu_p[i], i=1..m);
      if smu_p < smu or 
         (smu_p = smu and (LexSmaller(J, J_p) or Closer(mu, mu_p, m))) then
        userinfo(3, Modreduce, `  p = `, p, ` is unlucky`);
      else
        userinfo(3, Modreduce, `  previous primes are unlucky`);
        mu, J, M, R, q := mu_p, J_p, M_p, R_p, p;
        unchanged := 0;
      fi;
    fi;
  od;

  for row from 1 to m do
    LinearAlgebra:-RowOperation(R, row, 1/Z^omega[row], inplace=true)
  od;
  LinearAlgebra:-Map(expand, R);
  return M, R, omega, mu;
end:

#
# ReduceZpt
#
# Computes the order basis and residual over Zp by modular computation.
# This function only takes care of the parts for detecting unlucky
# homomorphisms and for reconstructing the results.
#
# Input:
#
#    F    : a matrix of skew polynomials over Z[t][Z;sigma,delta], where the
#           powers of Z are assumed to be on the right
#    m, n : dimensions of matrix F
#    p    : the prime p
#    N    : degree(F, Z)
#    T    : a bound on the degree(Z^l * F, t) for all 0 <= l <= mN+1
#    Z    : the indeterminate
#    sigma: the automorphism, implemented as a function: (Z[t], t) -> Z[t]
#    delta: the derivation, implemented as a function: (Z[t], t) -> Z[t]
#    t    : the variable in the domain Z[t]
#
# Output:
#
#    M and R: where
#
#          M . F = R . Z^omega;
#
#    - the number of nonzero rows in R = rank of F;
#    - the rows in M corresponding to the zero rows in R form a basis for
#      the left nullspace of F
#
#   omega:
#          the final order achieved (always (mN+1) \cdot \vec e)
#   mu:
#          the final degree constraints achieved
#   J:
#          the set of columns used in elimination
#

ReduceZpt := proc(F :: Matrix, m, n, p, N, T, Z :: name, sigma, delta, 
                  t :: name)
  local M, R, newM, newR, KM, unchanged, alpha, omega, mu, J, q, deg, 
        M_a, R_a, omega_a, mu_a, J_a, i, j, smu, smu_a, Fc, Fv, zv;

  option `Copyright (c) 2005 by Howard Cheng`;

  # first we allocate a matrix to minimize memory management:
  #   (mN+1)n + m rows:
  #     - start with m rows, and each of the (mN+1)n steps adds at most one
  #       row
  #   ((m+1)N+2)n columns for K:
  #     - need to go to order (mN+1), but we had N+1 coefficients to start
  #   (mN+1)n + m columns for M
  KM := LinearAlgebra:-Modular:-Create(p, (m*N+1)*n+m, 
                                       ((m+1)*N+2)*n + (m*N+1)*n+m,
                                       integer[kernelopts(wordsize)/8]);

  # initialize before entering loop
  unchanged := 0;
  omega := Vector[row](n);
  mu := Vector[row](m);
  J := [];
  deg := 0;                 # degree of what has been reconstructed so far  
  q := 1;
  zv := Vector[row](n);

  # Fc is a matrix of the shifts of the rows of F so we don't have to
  # recompute this over different evaluation points.  Fv stores the number
  Fc := [seq(Matrix((m*N+1)*n+m, ((m+1)*N+2)*n), i=1..m)];
  for i to m do
    for j to ((m+1)*N+2)*n do
      Fc[i][1,j] := coeff(F[i, modp(j-1, n)+1], Z, iquo(j-1, n));
    od;
  od;

  Fv := Vector[row](m, 1);

  # now do the loop: try all evaluation points until we are done
  for alpha from 0 to p-1 do
    userinfo(4, Modreduce, `Processing alpha = `, alpha);

    M_a, R_a, omega_a, mu_a, J_a := ReduceZp('Fc', Fv, KM, m, n, p, alpha, N, Z, 
                                             sigma, delta, t);
    if LinearAlgebra:-Equal(omega_a, zv) then
      # unlucky, just do the next one
      next;
    fi;

    if deg > 0 and not LinearAlgebra:-Equal(omega, omega_a) then
      print(`omega = `, omega);
      print(`omega_a = `, omega_a);
      error `omega is not the same for different points!  omega = `, omega,
            `omega_a = `, omega_a;
    fi;

    # now check for unlucky homomorphisms
    if (deg = 0) or
       (LinearAlgebra:-Equal(mu, mu_a) and LinearAlgebra:-Equal(J, J_a)) then
      # seems to be lucky
      userinfo(4, Modreduce, `  seems to be lucky`);
      
      if deg = 0 then
        mu, J, M, R, omega := mu_a, J_a, M_a, R_a, omega_a;
      else
        newM := MatrixInterp(M, M_a, q, alpha, Z, t, p);
        newR := MatrixInterp(R, R_a, q, alpha, Z, t, p);
        M_a, R_a := 0, 0;

        # check for early termination
        if LinearAlgebra:-Equal(M, newM) and LinearAlgebra:-Equal(R, newR) then
          unchanged := unchanged + 1;
          newM, newR := 0, 0; 
        else
          unchanged := 0;
          M, R := newM, newR;
        fi;
        if unchanged >= T then
          userinfo(4, Modreduce, `  early termination`);
          return M, R, omega, mu, J;
        fi;
      fi;
      deg := deg+1;
      q := q * (t-alpha) mod p;
    else
      # one of the is unlucky, figure out which one
      userinfo(4, Modreduce, `  unlucky evaluation(s) detected`);
      smu := add(mu[i], i=1..m);
      smu_a := add(mu_a[i], i=1..m);

      if smu_a < smu or 
         (smu_a = smu and (LexSmaller(J, J_a) or Closer(mu, mu_a, m))) then
        userinfo(4, Modreduce, `  alpha = `, alpha, ` is unlucky`);
      else
        userinfo(4, Modreduce, `  previous points are unlucky`);
        mu, J, M, R, deg, q := mu_a, J_a, M_a, R_a, 1, t-alpha;
        unchanged := 0;
      fi;
    fi;
  od;

  # we have tried all evaluation points
  error `not enough lucky evaluation points mod p = `, p;
end:

#
# ReduceZp
#
# Computes the order basis and residual over Zp by performing Gaussian
# elimination on the striped Krylov matrix that is built incrementally.
#
# Input:
#
#    Fc   : a list of matrices, each containing the shifts of the coefficients
#           of a stripe of F(Z)
#    Fv   : a vector giving the number of rows of Fc computed in each stripe
#    KM   : an ((mN+1)n+m) x (((m+1)N+2)n + (mN+1)n+m) matrix for storing the 
#           striped Krylov matrix and augmented transformation matrix
#    m, n : dimensions of matrix F
#    p    : the prime p
#    alpha: the evaluation point
#    N    : degree(F, Z)
#    Z    : the indeterminate
#    sigma: the automorphism, implemented as a function: (Z[t], t) -> Z[t]
#    delta: the derivation, implemented as a function: (Z[t], t) -> Z[t]
#    t    : the variable in the domain Z[t]
#
# Output:
#
#    M and R: where
#
#          M . F = R . Z^omega;
#
#    - the number of nonzero rows in R = rank of F;
#    - the rows in M corresponding to the zero rows in R form a basis for
#      the left nullspace of F
#
#   omega:
#          the final order achieved (always (mN+1) \cdot \vec e)
#   mu:
#          the final degree constraints achieved
#   J:
#          the set of columns used in elimination
#    

ReduceZp := proc(Fc :: name, Fv :: Vector, KM :: Matrix, m, n, p, alpha, N, 
                 Z :: name, sigma, delta, t :: name)
  local i, j, k, l, pi, omega, mu, J, rsize, csize, AddRow, 
        M, R, rank, found, mult, inv, det, rowindex;

  # add the next row of stripe i to K, and advance the row
  AddRow := proc(i)
    local j;

    # see if the row we need is already computed, it should only be
    # off by 1.
    if Fv[i] < mu[i]+1 then
      for j from csize to n+1 by -1 do
        Fc[i][mu[i]+1, j] := modp(sigma(Fc[i][mu[i],j-n], t) + 
                                  delta(Fc[i][mu[i],j], t), p);
      od;
      for j to n do
        Fc[i][mu[i]+1, j] := modp(delta(Fc[i][mu[i],j], t), p);
      od;
      Fv[i] := Fv[i]+1;
    fi;

    rsize := rsize+1;
    rowindex[i][mu[i]] := rsize;
    for j to csize do
      KM[rsize, j] := modp(subs(t=alpha, Fc[i][mu[i]+1,j]), p);
    od;
    for j to (m*N+1)*n+m do
      KM[rsize, j+csize] := 0;
    od;
    KM[rsize, rsize+csize] := 1;
  end:

  J := { NULL };
  mu := Vector[row](m);
  omega := Vector[row](n);
  rsize := 0;                    # number of rows currently in use in K
  csize := ((m+1)*N+2)*n;        # number of cols currently in use in K
  det := 1;
  rank := 0;

  # indexed by i and l, gives the index of the row in K containing the row
  # for the i-th stripe of Z^l * F.  Allocated for efficiency.
  rowindex := [seq(array(0..(m*N+1)*n+m), i=1..m)];

  # Initialize KM
  for i to m do
    AddRow(i);
  od;
  # now start eliminating each column
  for j to (m*N+1)*n do
    # look for pivot in the current rows
    found := false;
    for l to m do
      if KM[rowindex[l][mu[l]], j] <> 0 and 
         (not found or mu[l] < mu[pi]) then
        pi := l;
        found := true;
      fi;
    od;

    # if column is linearly dependent, do nothing
    if not found then
      userinfo(4, Modreduce, `  column `, j, ` dependent`);
      next;
    fi;

#    userinfo(4, Modreduce, `  column `, j, ` pivot `, pi);	
    # swap the pivot row
    if rowindex[pi][mu[pi]] <> rank+1 then
      # look for a row which is in rank+1 position
      found := false;
      for i to m while not found do
        for k from 0 to mu[i] while not found do
          if rowindex[i][k] = rank+1 then
            found := true;
            det := modp(-det, p);
            LinearAlgebra:-Modular:-Swap(p, KM, rowindex[i][k], 
                                            KM, rowindex[pi][mu[pi]]);
            LinearAlgebra:-Modular:-Swap(p, KM, 1..-1, 
                                            rowindex[i][k]+csize, 
                                            KM, 1..-1, 
                                            rowindex[pi][mu[pi]]+csize);
            rowindex[i][k] := rowindex[pi][mu[pi]];
            rowindex[pi][mu[pi]] := rank+1;
          fi;
        od;
      od;
      if not found then
        error `cannot find a row to swap`;
      fi;
    fi;
    
    # now add the next row, and eliminate it using previous pivots
    mu[pi] := mu[pi]+1;
    AddRow(pi);
    for i from pi+1 to m do
      if type(mu[i]+1, odd) then
        det := modp(-det, p);
      fi;
    od;
    for i to rank do
      found := false;
      for k from i to j-1 while not found do
        if KM[i,k] <> 0 then
          found := true;
          if KM[rsize,k] = 0 then
            next;
          fi;
          mult := modp(-KM[rsize,k]/KM[i,k], p);
          LinearAlgebra:-Modular:-AddMultiple(p, mult, 
                                              KM, rsize, k..rsize+csize,
                                              KM, i, k..rsize+csize, 
                                              KM, rsize, k..rsize+csize,
                                              'sparse');
        elif KM[rsize,k] <> 0 then
          # pivot rows are used out of sequence.  This doesn't necessarily
          # mean that we have an unlucky evaluation point and we may be
          # able to fix it.  But since this doesn't happen often we will
          # just return from here with some values that's guarantee to be
          # unlucky (because mu = 0)
          userinfo(4, Modreduce, `   Pivot rows out of sequence.`);
          return Matrix(m, m), Matrix(m, n), Vector[row](n), 
                 Vector[row](m), Vector[row](1);
        fi;
      od;
    od;

    # eliminate the remaining rows using pivot row
    det := modp(det * KM[rowindex[pi][mu[pi]-1], j], p);
    inv := modp(-KM[rowindex[pi][mu[pi]-1], j]^(-1), p);
    for i from rank+2 to rsize do
      mult := modp(inv * KM[i,j], p);
      LinearAlgebra:-Modular:-AddMultiple(p, mult, KM, i, j..rsize+csize,
                                          KM, rank+1, j..rsize+csize, 
                                          KM, i, j..rsize+csize,
                                          'sparse');
    od;

    rank := rank+1;
    J := J union {j};
    omega[modp(j-1, n)+1] := omega[modp(j-1, n)+1] + 1;
  od;

  # normalize answer to determinant
  for i to m do
    k := rowindex[i][mu[i]];
    mult := modp(det / KM[k, k+csize], p);
    LinearAlgebra:-Modular:-Multiply(p, mult, KM, k, KM, k);
  od;

  # now we extract M and R
  M := Matrix(m, m,
              (i,j) -> add(KM[rowindex[i][mu[i]], 
                              rowindex[j][k]+csize]*Z^k, 
                           k=0..mu[j]));
  R := Matrix(m, n, 
              (i,j) -> add(KM[rowindex[i][mu[i]], j+k*n]*Z^k, 
                           k=0..csize/n-1));
  return M, R, omega, mu, Vector[row](sort([op(J)]));
end:


#
# ComputeTKappa
#
# Compute the bounds T and kappa as defined in the paper.
#

ComputeTKappa := proc(F :: Matrix, Z :: name, sigma, delta, t :: name)
  local m, n, N, T, kappa, F2, l, i, j;

  option `Copyright (c) 2005 by Howard Cheng`;

  m, n := LinearAlgebra:-Dimensions(F);
  N := max(seq(seq(degree(F[i,j],Z),i=1..m),j=1..n));
  F2 := copy(F);

  T, kappa := 0, 0;

  T := max(T, seq(seq(degree(collect(F2[i,j], t), t), i=1..m), j=1..n));
  kappa := max(kappa, seq(seq(norm(expand(F2[i,j]), infinity), i=1..m), 
                          j=1..n));
  for l from 1 to m*N+1 do
    LinearAlgebra:-Map(f -> MultOreZ(f, Z, t, sigma, delta, m*N+1), F2);
    T := max(T, seq(seq(degree(collect(F2[i,j], t), t), i=1..m), j=1..n));
    kappa := max(kappa, seq(seq(norm(expand(F2[i,j]), infinity), i=1..m), 
                            j=1..n));
  od;
  return T, kappa;
end:

#
# MultOreZ
#
# Multiply an Ore polynomial f by Z (as defined by sigma and delta), and 
# truncate results to Z^(order-1).
# It is assumed that f is collected in Z already.
#

MultOreZ := proc(f, Z, t, sigma, delta, order)
  local r, i;

  option `Copyright (c) 2005 by Howard Cheng`;

  r := collect(add(sigma(coeff(f, Z, i),t)*Z^(i+1) + 
                   delta(coeff(f, Z, i),t)*Z^i, 
                   i=0..min(degree(f, Z), order-1)), Z);

  if degree(r, Z) >= order then
    return r - lcoeff(r, Z)*Z^order;
  else
    return r;
  fi;
end:

#
# LexSmaller
#
# Determine if one vector is lexicographically smaller than another.
# It is assumed that the two vectors have the same dimensions.
#

LexSmaller := proc(v1, v2)
  local i;

  for i to LinearAlgebra:-Dimension(v1) do
    if v1[i] < v2[i] then
      return true;
    elif v1[i] > v2[i] then
      return false;
    fi;
  od;
  return false;
end:

#
# Closer
#
# Determine if the first vector is lexicographically closer to a "normal"
# path then the second one.  Assumes that the two vectors have the same
# sum of components but different.
#

Closer := proc(v1, v2, m)
  local i, j, d1, d2, w;

  option `Copyright (c) 2005 by Howard Cheng`;

  w := Vector[row](m);
  j := 1;
  while true do
    # this loop should return eventually if the two vectors are
    # different
    d1 := add(max(0,w[i]-v1[i]),i=1..m);
    d2 := add(max(0,w[i]-v2[i]),i=1..m);
    if d1 < d2 then
      return true;
    elif d1 > d2 then
      return false;
    fi;
    w[j] := w[j]+1;
    j := modp(j, m) + 1;
  od;
end:

#
# MatrixCRA
#
# Given a matrix M1 correct mod m1 (product of primes) and a new image Mp
# mod p, return M such that M = M1 mod m1 and M = Mp mod p
#

MatrixCRA := proc(M1 :: Matrix, Mp :: Matrix, m1, p, Z)
  option `Copyright (c) 2005 by Howard Cheng`;

  return Matrix(LinearAlgebra:-Dimensions(M1), 
                (i,j) -> collect(chrem([expand(M1[i,j]), expand(Mp[i,j])], 
                                       [m1, p]), Z));
end:

#
# MatrixInterp
#
# Given a matrix M1 correct mod m1 (product of (t-alpha)'s) and a new image Ma
# mod (x-alpha), return M such that M = M1 mod m1 and M = Ma mod (t-alpha).
# Everything is computed mod p.
#

MatrixInterp := proc(M1 :: Matrix, Ma :: Matrix, m1, alpha, Z :: name, 
                     t :: name, p)
  local inv;
  option `Copyright (c) 2005 by Howard Cheng`;

  inv := modp(expand(subs(t=alpha, m1)^(-1)*m1), p);
  return Matrix(LinearAlgebra:-Dimensions(M1),
                (i,j) -> PolyCRA(M1[i,j], Ma[i,j], inv, alpha, Z, t, p));
end:

#
# PolyCRA
#
# Chinese remaindering/interpolation on each coefficient of 
# the polynomials p1 and p2 (in Z)
#

PolyCRA := proc(p1, p2, inv, alpha, Z :: name, t :: name, p)
  return modp(collect(modp(p1 + (p2 - subs(t=alpha,p1))* inv, p), Z, expand), p);
end: