// ****************************************************************************
// This file contains several functions regarding the solution of univariate polynomials of high degree
// issquarefree(pol) -- checks, if the univariate polynomial pol is square free
// makesquarefree(pol) -- returns a set B of the square free components of pol; uses makesquarefreeround as a subroutine
// DDF(pol) -- Distict Degree Factoring Algorithm for monic square free polynomials
// FDF(pol) -- Fixed Degree Factoring Algorithm
// Berlekamp(pol) and CantorZassenhaus(pol) are there to illustrate the algorithms
// they may fail in extreme cases where the Magma function Factorization(pol) might succeed
//
/* For usage see the example below with polynomials given in the book
   q:= 2^7;
   GF:=GaloisField(q);
   Pol<X>:=PolynomialRing(GF);
   pol:=X^6+X^5+X^4+X^3+1;
   load "Univariate";
   Berlekamp(pol);         // Berlekamp's Algorithm
   CantorZassenhaus(pol);  // Cantor-Zassenhaus Algorithm  
   Factorization(pol);     // Magma function to check the two algorithms

   clear ;
   q:= 3;
   GF:=GaloisField(q);
   Pol<X>:=PolynomialRing(GF);
   pol:= X^7+2*X^5+X^3+2*X ;
   load "Univariate";
   Berlekamp(pol);         // Berlekamp's Algorithm
   CantorZassenhaus(pol);  // Cantor-Zassenhaus Algorithm  
   Factorization(pol);     // Magma function to check the two algorithms



*********************************************************************/
printf "Solving Univariate Polynomials of High Degree ----------- \n \n";

function issquarefree(pol)
  fp:=Derivative(pol);
  tr := Gcd(pol,fp) eq 1 ;
  return tr;
end function;

function makesquarefreeround (pol)
  B1:=[];
    fp:=Derivative(pol);
    d:=Gcd(fp,pol);
    if d eq pol then
      h:=pol mod (X^q-X); // set X^q=X
      for i:=1 to q do
        B1:=B1 cat [h];
      end for;
    elif d ne 1 then
      B1:=B1 cat [d,Pol!(pol/d)];
    end if;
  return B1;
end function;

function makesquarefree (pol)
  B:=[pol];
  repeat
    sqf:=true;
    for i:=1 to #B do
      if issquarefree(B[i]) eq false then
        sqf:=false;
        B1:=makesquarefreeround(B[i]);
        B:=Remove(B,i);
        B:=B cat B1; 
      end if;
    end for;
  until sqf ;
  return B;
end function;
 
function DDF(pol) //let pol be a monic squarefree polynomial
  B:=[];
  i:=1;
  repeat
    R:=Gcd(X^(q^i)-X,pol);
    B:=B cat [R];
    pol:= Pol!(pol/R);
    i:=i+1;
  until pol eq 1 or i gt Degree(pol);
  return B;
end function;

function FDF(pol,d); 
  B:=[];
  repeat
    a:=X^d;
    for i:=0 to d-1 do
      a:=a+Random(GF)*X^i;
    end for;
    g:=Gcd(pol,a);
    if g ne 1 then
      B:= B cat [g];
     pol:=Pol!(pol/g);
     else
      t:=a^Integers()!((q^d-1)/2);
      g:=Gcd(t-1,pol);  
      if Degree(g) eq d then
        B:=B cat [g];
        pol:=Pol!(pol/g);
      end if;
    end if;
  until Degree(pol) eq d;
  B:=B cat [pol];
  return B;
end function;


function Bround(pol);
  B:=[];
  d:=Degree(pol);
  Qf:=ZeroMatrix(GF,d,d);
  a := 1;
  for i:=1 to d do
    for j:=1 to d do
      Qf[i][j]:=MonomialCoefficient(a,X^(j-1));
    end for;
    a:=a*X^q mod pol;
  end for;
  //print "BMatrix \n%o \n\n", Qf;
  M:=Qf-IdentityMatrix(GF,d);
  //print "Qf-I \n%o \n\n", M;
  ker:=Kernel(M);
  printf"Ker(Qf-I)= %o", ker;
  for k:=2 to #Basis(ker) do
    fk:=Basis(ker)[k];
    fp:=Pol!0;
    for j:=1 to d do
      fp:=fp+fk[j]*X^(j-1);
    end for;
    fp;
    for s in GF do
      g:=Gcd(pol,fp-s);
      if g ne 1 then
        //printf "Found Factor: %o \n", g;
        B:=B cat [g];
        pol:=Pol!(pol/g);
      end if;
    end for;
  end for;
  if pol ne 1 then
    B:=B cat [pol];
    //printf "Found Factor: %o \n", pol;
  end if;
  return B;
end function;

function CantorZassenhaus(pol)
  A:=makesquarefree(pol); // A is a set containing the square free factors of pol
  B:=[];
  C:=[];
  for i:=1 to #A do
    B[i]:=DDF(A[i]); // B[i] is a set containing the distinct degree factors of the polynomial A[i] 
    for j:=1 to #B[i] do
      if Degree(B[i][j]) gt j then
	C:=C cat FDF(B[i][j],j);
      elif B[i][j] ne 1 then
        C:=C cat [B[i][j]];
      end if;
    end for;
  end for;
  return C;
end function;

function Berlekamp(pol) // monic squarefree polynomial
  B:=[pol];
  for i:=1 to #B do
    poly:=B[i];
    if not IsIrreducible(poly) then
      B1:=Bround(poly);
    else
      B1:=poly;
    end if;   
  end for;
  return B1;
end function;