# Dixon/Bezout Resultant Package
# Date: Sep, 2003 Version 1.5
# Author: Arthur D. Chtcherba
# Email: cherba@cs.panam.edu

# Implements:  (See accompaning manual for details)
#  Dixon/Multivariate Bezoutian method to compute
#  the resultant of a polynomial system, introduced by Kapur/Saxena and Yang 1994
#
#  Dixon Dialytic method of A.D. Chtcherba and D. Kapur
#---------------------------------------------------------------------------------------------------------
#  Exported Methods:
#    mtx := DixonMatrix(pols, vars, order, 'rowH', 'colH', [roworder], [colorder])
#     input:  
#        pols - (d+1) polynomials
#        vars - d variables
#        order - specify order function of monomials which determines order arrangement of rows and columns
#                Note: maple has build in orders like: plex, tdeg, ...
#        rowH - name to assign
#        colH - name to assign
#        roworder - optional, specifies monomial order, to arrange matrix rows, 
#                   possible orders include plex, tdeg, ...
#        colorder - optional, specifies monomial order, to arrange matrix columns
#                   possible orders include plex, tdeg, ...
#     returns
#         mtx  - polynomial matrix
#         rowH - "row header", monomials (in terms of new variables "vars[i]_") labeling rows of the matrix
#         colH - "column header", monomials (in terms of vars[i]) labeling rows of the matrix
#     Implementation:
#       Calls function "DixonPolynomial", and then arranges terms in a polynomial matrix.
#----------------------------------------------------------------------------------------------------------
#  dixPol := DixonPolynomial(pols, vars, [newVars])
#     input:
#        pols - list of (d+1) polynomials (different order just changes sign)
#        vars - list d variables (different order produce completely different polynomials !)
#        newVars - a name for output or a list of new variables
#     output:
#        dixPol - expanded Dixon polynomial
#        newVars - a list of new variables if was given as name in input, otherwise unchanged
#     Implementation:
#        Expands cancelation matrix, by brute force, i.e. using maple LinearAlgebra, det function.
#        Slow for very big with many variables polynomials.
#     Example:
#       dixPol := DixonPolynomial(myPolSystem, [x,y]);
#      or
#       dixPol := DixonPolynomial(myPolSystem, [x,y],[alpha,beta]);
#      or
#      > dixPol := DixonPolynomial(myPolSystem, [x,y], newVar):
#      > newVars;
#                       [x_, y_]
#----------------------------------------------------------------------------------------------------------
#  mults:=DixonMults(g:polynom, pols::list(polynom), vars::list(symbol) T::name, rowH::name, colH::name, [order], [groupby])
#    Description:
#        Based on formula "g Theta = T x M", where Theta is Dixon matrix, T is some transformation matrix
#        and M is Sylvester-type matrix, compute T, and multipliers used for constructing M.
#    Inputs: 
#        g - arbitrary polynomial, often single monomial
#        pols - system of d+1 polynomials
#        vars - d - variables in the system
#        T - name to assign
#        rowH - name to assign
#        colH -  name to assign
#        order - optional, can be any monomial order default plex.
#                  Specifies how to arrange rows of matrix T.
#        groupby - optional, possible values `polynomial` or `multiplier`. Default `multiplier`.
#                  Specifies how to arrange columns of matrix T.
#    Outputs:
#        mults - a list of list of monomial multipliers for each of the polynomials to construct Sylvester-type matrix.
#        T -  transformation matrix.
#        rowH - (list of manomials) row header for the transformation matrix.
#        colH - (list of pol names) column header for the transformation matrix.
#
##################################################################################
Dixon := module()

export DixonPolynomial, 
       DixonMatrix, 
       DixonDialytic,
       DixonMults,
       DixonSupport,       
# Auxilary useful functions
       Polynomial2Matrix,  # Translate polynomial into matrix, by decomposing it based on two sets of variables.
       SizeDialyticMatrix,
       RSCCriteria, RSC,
       ConstructDialyticMatrix;
#*********************************************************************************
DixonMatrix := proc(pols::list(polynom), vars::list(name), rowH::symbol, colH::symbol)
local dixPol, newVars, roworder, colorder;
  if (nargs > 4) then
    roworder := args[5];
    if (nargs > 5) then
      colorder := args[6];
    else
      colorder := tdeg;
    fi;
  else
    roworder := tdeg;
    colorder := tdeg;
  fi;
  dixPol := DixonPolynomial(pols, vars, 'newVars');
  Polynomial2Matrix(dixPol, newVars, vars, rowH, colH, roworder, colorder);
end proc:
#*********************************************************************************
DixonPolynomial := proc(pols::list(polynom), vars::list(symbol))
local new_vars, nVars, nPols, CM, i, j, dv, dixp;

  # Do parameter checking
  if (nargs > 3) then
     ERROR("Too many arguments");
  else 
    if (nargs < 2) then
      ERROR("Too few arguments");
    fi;
  fi;
 
  nPols := nops(pols);
  nVars := nops(vars);

  if not(nPols = nVars + 1) then
     ERROR(`Number of equations should be number of variables plus one, but number of equations is %1 and number of vairables is %2`, nPols, nVars);
  fi;
  
  if (nargs = 3) then
    if type(args[3],symbol) then
      new_vars := [seq(cat(vars[i],`_`),i=1..nVars)];
      assign(args[3], new_vars);
    else
      if type(args[3],list(symbol)) then
        new_vars := args[3];
      else
        ERROR("Optional argument is expected to be of type symbol or list of symbols")
      fi;
    fi;
  fi;

  # Construct Cancelation Matrix
  CM := Matrix(nPols,nPols);

  for i from 1 to nPols do
    for j from 1 to nPols do
      if (i = 1) then
        CM[i,j] := pols[j];
      else
        CM[i,j] := subs(vars[i-1]=new_vars[i-1], CM[i-1,j]);
      fi;
    od;
  od;

  # Now, divide the determinant by the product of (vars[i-1] - new_vars[i-1]).
  # This is accomplished by subtracting ith row from the (i+1)^th, and
  # dividing the result by (vars[i-1] - new_vars[i-1]).

  for i from nPols by -1 to 2 do
    dv := (vars[i-1] - new_vars[i-1]);
    for j from 1 to nPols do
      CM[i, j] := simplify((CM[i, j] - CM[i-1, j])/dv);
    od;
  od;

  dixp := expand(Determinant(CM), method=fracfree):  # Dixon Polynomial
end proc:  
#**********************************************************************************************  
Polynomial2Matrix := proc(pol::polynom, rowVars::list(name), colVars::list(name), rowH::name, colH::name, roworder, colorder)
local rowPols, nRows, rowMon, rowMonS, row, col, colMon, nCols, M, i,j,polCoeffs,polColMon;
  rowPols := [coeffs(pol, rowVars, 'rowMon')];
  nRows := nops(rowPols);
  rowMon := [rowMon];

  coeffs(pol, colVars,'colMon');
  colMon := [colMon];
  nCols := nops(colMon);

  if (nCols > 1) then
    colMon := convert(sort(convert(colMon,`+`),colVars, colorder),list);
  fi;
  if (nRows > 1) then
    rowMonS := convert(sort(convert(rowMon,`+`), rowVars, roworder),list);
  else
    rowMonS := rowMon;
  fi;

  M := Matrix(nRows, nCols);
  for i from 1 to nRows do
    member(rowMonS[i], rowMon, 'row');
    polCoeffs := [coeffs(rowPols[row], colVars,'polColMon')];
    polColMon := [polColMon];

    for j from 1 to nops(polColMon) do 
       member(polColMon[j], colMon, 'col');
       M[row,col]:=polCoeffs[j];
    od:
  od:

  rowH := rowMon;
  colH := colMon;
  M;
end:
#**********************************************************************************************
DixonMults := proc(g::polynom, pols::list(polynom), vars::list(symbol), T::symbol, rowH::symbol, colH::symbol)
local groupby, order, nPols, mon, mult, rowNames, TMtable, nCols, 
  i, p, tm, j, k, TM, rc, colNames, cheaderNames, m, newPols;


  groupby := `polynomial`;  # Group matrix columns by polynomial, not by multiplier
  order := plex;            # Order matrix row by monomials using lex order
  
  # See if optional arguments for groupby and order have been specified  
  if (nargs > 5) then
    order := args[6];
    if (nargs > 6) then
      groupby := args[7];
    fi;
  fi;

  nPols := nops(pols);  # Number of polynomials in the system
  
  mult := [seq([], i=1..nPols)];  # Dialytic multiplies
  rowNames := {};                 # Row names/labels for matrix T
  TMtable := table();             # Table used to construct matrix T

  nCols := 0;     # Count number of columns
  
  for i from 1 to nPols do 
     newPols := subsop(i=g, pols);   # Substitute ith polynomial by g
     tm := DixonMatrix(newPols, vars, 'rowHi', 'colHi');  # Construct corresponding Dixon matrix
     # Add matrix to the table for for future construction of matrix T
     for j from 1 to nops(rowHi) do
       for k from 1 to nops(colHi) do
          if (tm[j,k]<>0) then 
             TMtable[i,rowHi[j],k] := tm[j,k];
          fi;
       od;
     od;
     mult[i] := colHi;              # Dixon matrix column labels are multipliers
     nCols := nCols + nops(colHi);  # Count number of new columns
     rowNames := rowNames union convert(rowHi,set); # Keep track of all rows
  od;
  #------------------------------------------------
  # Reconstruct matrix T from the table
  rowNames := convert(rowNames, list);
  TM := Matrix(nops(rowNames),nCols,0);
  rc := 0;
  cheaderNames := [];

  if (groupby = `multiplier`) then
     colNames := {}; 
     # Compute Union of All multiplies
     for i from 1 to nPols do 
       colNames := colNames union convert(mult[i],set); 
     od;
     # Sort them
     if (nops(colNames)>1) then 
       colNames := convert(sort(convert(colNames, `+`), vars, order), list); 
     else 
       colNames := convert(colNames, list); 
     fi;
     # Fill in the matrix
     for m in colNames do
       for i from 1 to nPols do
         if (member(m, mult[i], 'pos')) then
           rc := rc+1;
           cheaderNames := [op(cheaderNames), m*cat(`f`,(i-1))];
           for k from 1 to nops(rowNames) do
             if (assigned(TMtable[i, rowNames[k], pos])) then
               TM[k,rc] := TMtable[i, rowNames[k], pos];
             fi;
           od;
         fi;
       od;
     od;
   else
     # Fill in the matrix polynomial by polynomial
     rc:=1;
     for i from 1 to nPols do
        for j from 1 to nops(mult[i]) do
           cheaderNames := [op(cheaderNames), (mult[i][j])*cat(`f`,(i-1))];
           for k from 1 to nops(rowNames) do
             # print("Pol: ",i,"Row: ", rowNames[k],"Mult: ",j);
             if (assigned(TMtable[i, rowNames[k], j])) then
               TM[k, rc] := TMtable[i, rowNames[k], j];
             fi;
          od;
          rc := rc + 1;
       od;        
     od;
  fi;
  # Return Computed Values
  rowH := rowNames;
  colH := cheaderNames;
  T := TM;
  mult;
end:
#-------------------------------------------------------------------------------------------
# Given a set of vectors of cardinality (d+1) and list of (d+1) polynomials construct
# a matrix where each polynomial is translated by its own set of vectors (monomials)
# i'th polynomial is named f(i-1), hence one gets f0,...,fd (d -dimension, number of variables)
# Returns:
# Matrix
# bt - row names
# xt - column names
# order - defines monomial order
# groupby - is either `polynomial` or `multiplier`
ConstructDialyticMatrix :=proc(B::list, pols::list(polynom), vars::list(name), bt::name, xt::name)
local nRows,  
      tb, 
      cf, p,
      colH,    # Column Header
      rowH,    # Row header
      nPols,
      i,j,k,
      t,
      mon,
      nTerms,
      M, groupby, order, mult;

  if (nops(B) <> nops(pols)) then
     ERROR("Number of multplier sets does not match number of polynomials");
  fi;

  groupby := `polynomial`;
  order := plex;
  if (nargs > 5) then
    order := args[6];
    if (nargs > 6) then
      groupby := args[7];
    fi;
  fi;

  tb := table();
  colH := {};
  nPols := nops(pols);
  nRows := 0;
  rowH := [];
  if (groupby = `multiplier`) then
    mult := {};
    for i from 1 to nPols do
      mult := mult union {op(B[i])};
    od;
    if (nops(mult) > 1) then
        mult := convert(sort(convert(mult, `+`), vars, plex), list);
    fi;
    for i from 1 to nops(mult) do
      mon := mult[i];
      for j from 1 to nPols do
        if (member(mon, B[j])) then
          nRows := nRows + 1;
          rowH := [op(rowH), mon*cat(`f`,j-1)];
          p := expand(mon*pols[j]);
          cf := [coeffs(p, vars, 't')];
          colH := colH union {t};
          t := [t];
          for k from 1 to nops(t) do
            tb[nRows,t[k]] := cf[k];
          od;
        fi;
      od;

    od;

  else
    for i from 1 to nPols do
      for j from 1 to nops(B[i]) do 
        mon := B[i][j];
        nRows := nRows + 1;
        rowH := [op(rowH), mon*cat(`f`,i-1)];
        p := expand(mon*pols[i]);
        cf := [coeffs(p, vars, 't')];
        colH := colH union {t};
        t := [t];
        for k from 1 to nops(t) do
          tb[nRows,t[k]] := cf[k];
        od;
      od;
    od;
  fi;

  # Sort column header in ascending total degree order
  nTerms := nops(colH);

  if (nTerms > 1) then
    colH := convert(sort(convert(colH, `+`), vars, order), list);
  else
    colH := convert(colH, list);
  fi;

  M := Matrix(nRows, nTerms);
  for i from 1 to nRows do
   for j from 1 to nTerms do
     if (assigned(tb[i, colH[j]])) then
       M[i,j] := tb[i, colH[j]];
     else
       M[i,j] := 0;
     fi;
   od;
  od;
  xt := colH;
  bt := rowH;
  M;
end:
#**********************************************************************************************
#-------------------------------------------------------------------------------------------
# Given a list of list of monomials of cardinality (d+1) and list of (d+1) polynomials compute
# the size of a dialytic matrix where each polynomial is "translated" by its set of monomials.
# i'th polynomial is named f(i-1), hence one gets f0,...,fd (d -dimension, number of variables)
# Returns:
# tuple [#rows, #cols]
SizeDialyticMatrix :=proc(B::list, pols::list(polynom), vars::list(name), bt::name, xt::name)
local nRows,  
      tb, 
      cf, p,
      colH,    # Column Header
      rowH,    # Row header
      nPols,
      i,j,k,
      t,
      mon,
      nTerms,
      M, groupby, order, mult;

  groupby := `polynomial`;
  order := plex;
  if (nargs > 5) then
    order := args[6];
    if (nargs > 6) then
      groupby := args[7];
    fi;
  fi;

  tb := table();
  colH := {};
  nPols := nops(pols);
  nRows := 0;
  rowH := [];
  if (groupby = `multiplier`) then
    mult := {};
    for i from 1 to nPols do
      mult := mult union {op(B[i])};
    od;
    mult := convert(sort(convert(mult, `+`), vars, plex), list);
    for i from 1 to nops(mult) do
      mon := mult[i];
      for j from 1 to nPols do
        if (member(mon, B[j])) then
          nRows := nRows + 1;
          rowH := [op(rowH), mon*(`f`.(j-1))];
          p := expand(mon*pols[j]);
          cf := [coeffs(p, vars, 't')];
          colH := colH union {t};
        fi;
      od;

    od;

  else
    for i from 1 to nPols do
      for j from 1 to nops(B[i]) do 
        mon := B[i][j];
        nRows := nRows + 1;
        rowH := [op(rowH), mon*(`f`.(i-1))];
        p := expand(mon*pols[i]);
        cf := [coeffs(p, vars, 't')];
        colH := colH union {t};
      od;
    od;
  fi;

  # Sort column header in ascending total degree order
  nTerms := nops(colH);

  if (nTerms > 1) then
    colH := convert(sort(convert(colH, `+`), vars, order), list);
  else
    colH := convert(colH, list);
  fi;

  xt := colH;
  bt := rowH;
  [nops(rowH),nops(colH)];
end:
#===========================================================================
RSCCriteria := proc(M::Matrix, varLst::list(name))
  local rows, cols, i, j, rM, rN, N, Ms, lst, iCols, p, v, substLst;

  (rows, cols) := Dimension(M);

  # Collect the names of all variables apearing in the matrix

  #printf(`Matrix is %d x %d\n`, rows, cols);

  #varLst := {};
  #for i from 1 to rows do
  #  for j from 1 to cols do
  #    varLst := varLst union indets(M[i,j]);
  #  od;
  #od;

  # Create Substitution List
  substLst := []:
  p := 137;
  for v in varLst do 
    p := nextprime(p);
    substLst := [op(substLst),v=p]: 
  od:

  # Create Numberic Matrix
  Ms := Matrix(rows,cols);
  for i from 1 to rows do
    for j from 1 to cols do
      Ms[i,j] := subs(op(substLst),M[i,j]);
    od:
  od:

  rM := Rank(Ms);
  lst := [seq(i, i=1..cols)];
  printf("The matrix has rank %d\n", rM);
  iCols := [];
  for i from 1 to cols do
     N := SubMatrix(Ms, 1..rows, subsop(i=NULL, lst));
     rN := Rank(N);
     printf("M\\C_%d rank is %d\n", i, rN);
     if (rN < rM) then
       iCols := [op(iCols), i];
     fi;
  od;
  if (nops(iCols) = 0) then
    printf("Failed\n");
  else
    printf("# of indep cols %d out of %d\n", nops(iCols), cols);
  fi;
  iCols;
end:
#====================================================================
# RankSubMinorConstruction

RSC := proc(M::Matrix, varLst::list(name))
  local rows, cols, i, j, rM, rN, N, Um, Ms, lst, iCols, p, v, substLst, newLst;

  (rows, cols) := Dimension(M);

  # Collect the names of all variables apearing in the matrix

  #printf(`Matrix is %d x %d\n`, rows, cols);

  #varLst := {};
  #for i from 1 to rows do
  #  for j from 1 to cols do
  #    varLst := varLst union indets(M[i,j]);
  #  od;
  #od;

  # Create Substitution List
  substLst := []:
  p := 137;
  for v in varLst do 
    p := nextprime(p);
    substLst := [op(substLst),v=p]: 
  od:

  # Create Numberic Matrix
  Ms := Matrix(rows,cols);
  for i from 1 to rows do
    for j from 1 to cols do
      Ms[i,j] := subs(op(substLst),M[i,j]);
    od:
  od:

  Um := LUDecomposition(Ms, output='U');

  lst := [seq([], i=1..cols)];
  # printf("The matrix has rank %d\n", rM);
  for i from 1 to cols do
    for j from min(i,rows) to 1 by -1 do
      if (Um[j,i] <> 0) then
        lst[j] := [op(lst[j]), i];
        break;
      fi;
    od;
  od;
  newLst := [];
  for i from 1 to cols do
    if (nops(lst[i])>0) then newLst := [op(newLst), lst[i]]; fi;
  od;
  newLst;
end:
end module: