# FFreduce
#
# Author: Howard Cheng and George Labahn
# 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.
#
# ---> FFreduce(F, Z, sigma, delta, n)
#
# Implements the FFreduce algorithm in
#
# Beckermann, B., Cheng, H. and Labahn, G. "Fraction-free Row Reduction 
#   of Matrices of Ore Polynomials." Journal of Symbolic Computation, 
#   41(5), pages 513-543, 2006. 
#
# Input:
#
#    F    : a matrix of skew polynomials over D[Z;sigma,delta], where the 
#           powers of Z are assumed to be on the right, where D is an 
#           appropriate integral domain
#    Z    : the indeterminate
#    sigma: the automorphism, implemented as a function: (D,n) -> D.
#    delta: the derivation, implemented as a function: (D,n) -> D.
#    n    : the variable in the domain D (e.g. Q[n,q^n])
#
# Output:
#
#    M and R: where
#
#          M . F = R . Z^omega;
#
#          the trailing coefficient of R has the same rank as F;
#          the number of rows in M and R = rank of F;
#          there exists an M^* such that  M^* . M = Z^p for some p
#
#   omega: 
#          the final order achieved
#   mu: 
#          the final degree constraints achieved
#

FFreduce := proc(F :: Matrix, Z :: name, sigma, delta, n :: name)
  local M, R, j, k, l, v, m, s, lambda, Lambda, omega, r, pi, p, c,
        mu, t, d, rho, i, N, pivotrows;

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

  # Initialization
  m,s := op(1,F);
  M := Matrix(m,m,(i,j) -> if i=j then 1 else 0 fi);
  R := Matrix(F);
  d := 1;
  omega := Vector(s);
  mu := Vector(m);
  rho := 0;
  N := max(seq(seq(degree(R[i,j],Z),i=1..m),j=1..s));

  while omega[1] < m*N+1 do
    rho := 0;
    pivotrows := NULL;
    for lambda from 1 to s do
      r := Vector(m, [seq(coeff(R[l,lambda], Z, 0),l=1..m)]);
      Lambda := select(l -> evalb(r[l] <> 0), [seq(l,l=1..m)]);
      if nops(Lambda) > 0 then
        pi := min(op(select(l -> evalb(mu[l] = min(seq(mu[v], v=Lambda))), 
                            Lambda)));

        p := Vector(m, l -> if l <> pi then
                              coeff(M[pi,l],Z,mu[l]-1);
                            else
                              coeff(M[pi,l],Z,mu[l]);
                            fi);

        # increase order for l = 1, ..., m, l <> pi
        for l to m do
          if l <> pi then
            for j to m do
              t := collect(expand(r[pi]*M[l,j] - r[l]*M[pi,j]),Z);
              M[l,j] := add(normal(coeff(t,Z,k)/d)*Z^k,k=0..degree(t,Z));
            od;
            for j to s do
              t := collect(expand(r[pi]*R[l,j] - r[l]*R[pi,j]),Z);
              R[l,j] := add(normal(coeff(t,Z,k)/d)*Z^k,k=0..degree(t,Z));
            od;
          fi;
        od;

        # increase order and adjust degree constraints for row pi
        d := sigma(d,n);
        for j to m do
          t := collect(expand(r[pi]*(ApplyCoeff(M[pi,j],Z,sigma,n)*Z
                              +ApplyCoeff(M[pi,j],Z,delta,n))
                      - delta(r[pi],n)*M[pi,j]
                      - add(sigma(p[l],n)*M[l,j],l=1..pi-1) 
                      - add(sigma(p[l],n)*M[l,j],l=pi+1..m)),Z);
          M[pi,j] := add(normal(coeff(t,Z,k)/d)*Z^k,k=0..degree(t,Z));

        od;
        for j to s do
          t := collect(expand(r[pi]*(ApplyCoeff(R[pi,j],Z,sigma,n)*Z 
                              +ApplyCoeff(R[pi,j],Z,delta,n))
                      - delta(r[pi],n)*R[pi,j]
                      - add(sigma(p[l],n)*R[l,j],l=1..pi-1) 
                      - add(sigma(p[l],n)*R[l,j],l=pi+1..m)),Z);
          R[pi,j] := add(normal(coeff(t,Z,k)/d)*Z^k,k=0..degree(t,Z));
        od;
	    
        d := r[pi];
        mu[pi] := mu[pi]+1;
        rho := rho+1;
        pivotrows := pivotrows, pi;
      fi;

      # adjust residual in column lambda
      for l to m do
        R[l,lambda] := collect(R[l,lambda]/Z,Z);
      od;

      omega[lambda] := omega[lambda] + 1;
    od;
  od;


  return M, R, omega, mu;
end:

#
# ---> ApplyCoeff(p, sigma, n)
#
# Apply sigma to each coefficient of p
#
ApplyCoeff := proc(p, Z, sigma, n)
  local k;
  if p = 0 then
    return p;
  else
    return add(sigma(coeff(p,Z,k),n)*Z^k, k=ldegree(p,Z)..degree(p,Z));
  fi;
end: