Define SqFree(F)
  If Type(F) <> POLY Then Return Error("Input must be a Polynomial"); EndIf;
  Univar:=UnivariateIndetIndex(F);
  If Univar = 0 Then Return Error("NOT a Univariate Polynomial"); EndIf;
  PolyInd:=Indet(Univar);
  P:=Characteristic();
 
  If P = 0 Then Return F/(GCD(F,Der(F,PolyInd))); EndIf;

    S1:=GCD(F,Der(F,PolyInd));

    If S1=1 Then Return F;
    Else
      S11:=S1;

      While S1<>1 Do
	S11:=S1;
	DS1:=Der(S1,PolyInd);

	While DS1<>0 Do
	   S2:= GCD(S1,Der(S1,PolyInd));
	   DS1:=Der(S2,PolyInd);
	   S1:=S2;
	EndWhile;

       Sup:=Support(S1);
       SupNew:=[];
       
       Foreach T In Sup Do
         Li:=Log(T);
         L:=Li*(1/P);
         TNew:=LogToTerm(L);
         Append(SupNew,TNew);
       EndForeach;

       G:=ScalarProduct(Coefficients(S1),SupNew);
       F:=F*G/S11;
       S1:=GCD(F,Der(F,PolyInd));
      EndWhile;

    Return F;
   EndIf;

EndDefine;


Z3::=Z/(3)[x];Z11::=Z/(11)[x];Z5::=Z/(5)[x];

Ex#12

F1:=x^9+2x^8-x^7+2x^5-2x^4+2x^3+x^2+2x-2;
F2:=x^12+2x^11+5x^10-x^9+x^8-4x^7-x^6-4x^5-x^4+5x^3-3x^2-3x+4;
F3:=x^13+x^12-x-1;

OUTPUT

Use Z5;
F1:=x^9+2x^8-x^7+2x^5-2x^4+2x^3+x^2+2x-2;
SqFree(F1);
Squarefree Polynomial is x^4 + 2x^3 - x^2 - x - 1
-------------------------------
Use Z11;
F2:=x^12+2x^11+5x^10-x^9+x^8-4x^7-x^6-4x^5-x^4+5x^3-3x^2-3x+4;
SqFree(F2);
Squarefree Polynomial is x^5 + 2x^3 + x^2 + 2
-------------------------------
Use Z3;
F3:=x^13+x^12-x-1;
SqFree(F3);
Squarefree Polynomial is x^4 - 1
-------------------------------