
#this file contains methods for computing with polynomials with n binary variables, they are needed for the computation of isogenies for curves over fields of characteristic 2


   _Coefficients:=function(f) 
      return Copy(Base(f).coeffs);
   end;

   _Variables:=function(f) 
      return Copy(Base(f).vars); 
   end;

   _CoefficientRing:= function(PA)
      return Base(PA);
   end;

   # this function is only for forwarding
   MultPol:=function(a)
      return a;
   end;

   _alg_mult_ops:=rec();

   _alg_mult_ops.Print:=function(al)
      Print("Polynomial Algebra over ",Base(al)," in ", al.vars," variables ");
   end;

   _elt_alg_mult_pol_ops:=rec(); 

   Merge:=function(L,c,left,middle,right)
     local i,j,m,n,L1,L2,c1,c2,L,c;
     L1:=[];
     L2:=[];
     c1:=[];
     c2:=[];
     m:=middle-left+1;
     n:=right-middle;
     for i in [1..m] do
        L1[i]:=L[left+i-1];
        c1[i]:=c[left+i-1];
     od;
     for i in [1..n] do
        L2[i]:=L[middle+i];
        c2[i]:=c[middle+i];
     od;
     i:=1;
     j:=1;
     while i<=m and j<=n do
        if L1[i]<L2[j] then
           L[left+i+j-2]:=L1[i];
           c[left+i+j-2]:=c1[i];
           i:=i+1;
        else 
           L[left+i+j-2]:=L2[j];
           c[left+i+j-2]:=c2[j];
           j:=j+1;
        fi;
     od;
     if j = n+1 then
        while i <=m do
           L[left+i+j-2]:= L1[i];
           c[left+i+j-2]:= c1[i];
           i:=i+1;
        od;
     fi;
   end;

   MergeSortRec:=function(L,c,first,ende)
     local middle;
     if first = ende then
        return;
     fi;
     middle:= Floor((first+ende)/2);
     MergeSortRec(L,c,first,middle);
     MergeSortRec(L,c,middle+1,ende);
     Merge(L,c,first,middle,ende);
   end;

  
   MergeSort:=function(L,c)
     return MergeSortRec(L,c,1,Length(L));
   end;

   _Reduction_:=function(f)
     local i,j,L,c;
     L:=Base(f).vars;
     c:=Base(f).coeffs;
     MergeSort(L,c);
     i:=0;
     while i <Length(L)-1 do
        i:=i+1;
        j:=i;
        while j<Length(L) do
           j:=j+1;
           if L[i]=L[j] then
              Remove_(L,j);
              c[i]:=c[i]+c[j];
              Remove_(c,j);
              j:=j-1;
           fi;
        od;
     od;
     i:=0;
     while i <Length(c) do
        i:=i+1;
        if c[i]=0 then
           Remove_(c,i);
           Remove_(L,i);
           i:=i-1;
        fi;
     od;
     if L=[] and c = [] then
        Base(f).vars:=[[0]];
        Base(f).coeffs:=[0];
     else
        Base(f).vars:=Set(L);
        Base(f).coeffs:=c;
     fi;
   end;


   _Coerce:=function(F,f)
     local L,c;
     L:=_Variables(f);
     c:=Coefficients(f);
     if Length(L)=1 and L[1][1]=0 then
        return Coerce(F,c[1]);
     else 
        return Error("cannot coerce ",f," into ",F);
     fi;
   end;

 
   MultPolAlg:=function(K,n)
     local al;
     al:=rec(base:=K,type:=alg^mult,operations:=_alg_mult_ops,vars:=n);
     return al;
   end;


   _MultPol:=function(PA,L,c)
     local f,i;
     if L=[] then
        if not IsList(c) then
           return rec(base:=rec(vars:=[[0]],coeffs:=[c]),type:=elt-alg^mult,parent:=PA,operations:=_elt_alg_mult_pol_ops);
        else 
           return Error("not enough information for creating a polynomial");
        fi;
     fi;
     if Is(Type(L[1]),list) then
        for i in [1..Length(L)] do
           L[i]:=Set(L[i]);
           if L[i][Length(L[i])]>PA.vars then
              return Error("too many variables for ",PA);
           fi;
        od;
     else 
        L:=[Set(L)];
        if L[1][Length(L[1])] > PA.vars then
           return Error("too many variables for ",PA);
        fi;
        c:=[c];
     fi;
     for i in [1..Length(c)] do
        c[i]:=Coerce(Base(PA),c[i]);
     od;
     if not Length(L) = Length(c) then
        if Length(L) > Length(c) then
           return Error("not enough coefficients");
        else return Error("too many coefficients");
        fi;
     fi;
     f:=rec(base:=rec(vars:=L,coeffs:=c),type:=elt-alg^mult,parent:=PA,operations:=_elt_alg_mult_pol_ops);
     _Reduction_(f);
     return f;
   end;


   _MultPol2:=function(PA,L,c,d)
     local f,i;
     if Is(Type(L[1]),list) then
        for i in [1..Length(L)] do
           L[i]:=Set(L[i]);
           if L[i][Length(L[i])]>PA.vars then
              return Error("too many variables for ",PA);
           fi;
        od;
     else L:=[Set(L)];
        if L[1][Length(L[1])] > PA.vars then
           return Error("too many variables for ",PA);
        fi; 
        c:=[c];
     fi;
     for i in [1..Length(c)] do
        c[i]:=Coerce(Base(PA),c[i]);
     od;
     if not Length(L) = Length(c) then
        if Length(L) > Length(c) then
           return Error("not enough coefficients");
        else return Error("too many coefficients");
        fi;
     fi;
     f:=rec(base:=rec(vars:=L,coeffs:=c),type:=elt-alg^mult,parent:=PA,operations:=_elt_alg_mult_pol_ops,denom:=d);
     _Reduction_(f);
     return f;
   end;



   _elt_alg_mult_pol_ops.Print:= function(f)
     local L,c,i,j,den;
     L:=_Variables(f);
     c:=Coefficients(f);
     if "denom" in RecFields(f) then
        if f.denom<>1 then
           den:=true;
           Print("(");
        else den:=false;
        fi;
     else den := false;
     fi;
     for i in [1..Length(L)] do
        Print(c[i]);
        if L[i][1] <> 0 then
           Print("*");
           for j in [1..Length(L[i])] do
              Print("x",L[i][j]);
              if not j = Length(L[i]) then
                 Print("*");
              fi;
           od;
        fi;
        if not i = Length(L) then
           Print(" + ");
        fi;
     od;
     if den then
        Print(")/(");
        Print(f.denom);
        Print(")");
     fi;
   end;

   _elt_alg_mult_pol_ops.\+:= function(f1,f2)
     local L1,L2,c1,c2,g,denom,f11,f21;
     if Is(Type(f1),elt-alg^mult) and Is(Type(f2),elt-alg^mult) then
        if "denom" in RecFields(f1) and "denom" in RecFields(f2) then
           f11:=f1*f2.denom;
           f21:=f2*f1.denom;
           denom:=f1.denom*f2.denom;
        elif "denom" in RecFields(f1) then
           f11:=f1;
           f21:=f2*f1.denom;
           denom:=f1.denom;
        elif "denom" in RecFields(f2) then
           f11:=f1*f2.denom;
           f21:=f2;
           denom:=f2.denom;
        else 
           f11:=f1;
           f21:= f2;
           denom:=1;
        fi;
        L1:=_Variables(f11);
        L2:=_Variables(f21);
        c1:=Coefficients(f11);
        c2:=Coefficients(f21);
        if denom <>1 then
           g:=MultPol(Parent(f1),Append(L1,L2),Append(c1,c2),denom);
        else 
           g:=MultPol(Parent(f1),Append(L1,L2),Append(c1,c2));
        fi;
        return g;
     elif Is(Type(f1),elt-alg^mult) then
        if not Is(Type(f2),elt-fld^fin) and not  Is(Type(f2),elt-ord^rat) then
           return Error("wrong type in addition");
        fi;
        if "denom" in RecFields(f1) then
           f21:=f1.denom*f2;
           return f21+f1;
        else return MultPol(Parent(f1),Append(_Variables(f1),[[0]]),Append(Coefficients(f1),[f2]));
        fi;
     elif Is(Type(f2),elt-alg^mult) then
        if not Is(Type(f1),elt-fld^fin) and not Is(Type(f1),elt-ord^rat) then
           return Error("wrong type in addition");
        fi;
        if "denom" in RecFields(f2) then
           f11:=f2.denom*f1;   
           return f11+f2;
        else
           return MultPol(Parent(f2),Append(_Variables(f2),[[0]]),Append(Coefficients(f2),[f1]));
        fi;
     fi;
   end;

   _elt_alg_mult_pol_ops.\^:= function(f1,n)
     local L1,L2,c1,c2,g,i;
     if Is(Type(f1),elt-alg^mult) and Is(Type(n),elt-ord^rat) then
        g:=f1;
        for i in [2..n] do
           g:=g*f1;
        od;
        return g;
     else return Error("wrong type in ^");
     fi;
   end;

   _elt_alg_mult_pol_ops.\-:= function(f1,f2)
     return f1+f2; #nur fuer den Fall K=GF(2,n)
   end;

   _elt_alg_mult_pol_ops.\/:=function(f1,f2)
     local L,c,PA,i;
     PA:=Parent(f1);
     if Is(Type(f2),elt-ord^rat) or Is(Type(f2),elt-fld^fin) then
        f2:= Coerce(PA,f2);
        L:= _Variables(f1);
        c:= Coefficients(f1);
        for i in [1..Length(c)] do
           c[i]:= c[i]/f2;
        od;
        return MultPol(PA,L,c);
     else return Error("wrong type in / ");
     fi;
   end;

   _elt_alg_mult_pol_ops.\=:= function(f1,f2)
     local L,c;
     if Is(Type(f1),elt-alg^mult) and Is(Type(f2),elt-alg^mult) then
        return Coefficients(f1)=Coefficients(f2) and _Variables(f1)=_Variables(f2);
     elif Is(Type(f1),elt-fld^fin) or Is(Type(f1),elt-ord^rat) then
        L:=_Variables(f2);
        if Length(L) >1 then
           return false;
        elif Length(L[1])>1 then
           return false;
        elif L[1][1] <> 0 then
           return false;
        else
           c:=Coefficients(f2);
           f1:=Coerce(Parent(f2),f1);
           return c[1] = f1;
        fi;
     elif Is(Type(f2),elt-fld^fin) or Is(Type(f2),elt-ord^rat) then
        L:=_Variables(f1);
        if Length(L) >1 then
           return false;
        elif Length(L[1])>1 then
           return false;
        elif L[1][1] <> 0 then
           return false;
        else
           c:=Coefficients(f1);
           f1:=Coerce(Parent(f1),f2);
           return c[1] = f2;
        fi;
     else return Error("wrong type in = ");
     fi;
   end;

   _elt_alg_mult_pol_ops.\*:= function(f1,f2)
     local PA,f12,g,L,L1,L2,c,c1,c2,i,j,pos,dL1,dL2,dc1,dc2,dL,dc,denom;
     if Is(Type(f1),elt-ord^rat) or Is(Type(f1),elt-fld^fin) then
        PA:=Parent(f2);
        f1:=Coerce(PA,f1);
        c:=Coefficients(f2);
        L:=_Variables(f2);
        for i in [1..Length(c)] do
           c[i]:= f1*c[i];
        od;
        if "denom" in RecFields(f2) then
           return MultPol(PA,L,c,f2.denom);
        else return MultPol(PA,L,c);
        fi;
     elif Is(Type(f2),elt-ord^rat) or Is(Type(f2),elt-fld^fin) then
        PA:=Parent(f1);
        f2:=Coerce(PA,f2);
        c:=Coefficients(f1);
        L:=_Variables(f1);
        for i in [1..Length(c)] do
           c[i]:= f2*c[i];
        od; 
        if "denom" in RecFields(f1) then
           return MultPol(PA,L,c,f1.denom);
        else return MultPol(PA,L,c);
        fi;
     elif Is(Type(f1),elt-alg^mult) and Is(Type(f2),elt-alg^mult) then
        PA:= Parent(f1);
        L1:= _Variables(f1);
        L2:= _Variables(f2);
        c1:= Coefficients(f1);
        c2:= Coefficients(f2);
        L:=[];
        c:=[];
        for i in [1..Length(L1)] do
           for j in [1..Length(L2)] do
              Append_(L,[Set(Append(L1[i],L2[j]))]);
              Append_(c,[c1[i]*c2[j]]);
              if Length(Last(L))>1 then
                 pos:=Position(Last(L),0);
                 if pos <> FAILURE then
                    Remove_(Last(L),pos);
                 fi;
              fi;
           od;
        od;
        if "denom" in RecFields(f1) and "denom" in RecFields(f2) then
           if f1.denom = 1 and f2.denom =1 then
              return MultPol(PA,L,c);
           elif f1.denom <>1 then
              if f2.denom = 1 then
                 return MultPol(PA,L,c,f1.denom);
              else
                 dL1:=_Variables(f1.denom);
                 dL2:=_Variables(f2.denom);
                 dc1:=Coefficients(f1.denom);
                 dc2:=Coefficients(f2.denom);
                 dL:=[];
                 dc:=[];
                 for i in [1..Length(dL1)] do
                    for j in [1..Length(dL2)] do
                       Append_(dL,[Set(Append(dL1[i],dL2[j]))]);
                       Append_(dc,[dc1[i]*dc2[j]]);
                       if Length(Last(dL))>1 then
                          pos:=Position(Last(dL),0);
                          if pos <> FAILURE then
                             Remove_(Last(dL),pos);
                          fi;
                       fi;
                    od;
                 od;
                 denom:=MultPol(PA,dL,dc);
                 return MultPol(PA,L,c,denom);
              fi;
           else
              return MultPol(PA,L,c,f2.denom);
           fi;
        elif "denom" in RecFields(f1) then
           return MultPol(PA,L,c,f1.denom);
        elif "denom" in RecFields(f2) then
           return MultPol(PA,L,c,f2.denom);
        else 
           return  MultPol(PA,L,c);
        fi;
     else return Error("wrong types in multiplication");
     fi;
   end;

   IsConstant:=function(g)
     if Is(Type(g),elt-ord^rat) or Is(Type(g),elt-fld^fin) then 
        return true;
     else 
        return (Length(_Variables(g)) =1 and Last(_Variables(g)[1]) =0);
     fi;
   end;

   Elimination_:=function(g,k,f)
     local i,L,c,j,pos,g2,PA,denom,den,L2,gd;
     if IsConstant(g) then 
        return g;
     fi;
     PA:=Parent(g);
     L:=_Variables(g);
     c:=Coefficients(g);
     g2:=MultPol(PA,[0],0);
     j:=0;
     denom:=("denom" in RecFields(g));      
     for j in [1..Length(L)] do
        pos:=Position(L[j],k);
        if pos = FAILURE then
           g2:=g2+MultPol(PA,L[j],c[j]);
        else
           L2:=Remove(L[j],pos);
           if L2=[] then
              g2:=g2+f*c[j];
           else 
              g2:=g2+MultPol(PA,L2,c[j])*f;
           fi;
        fi; 
     od;
     if denom then 
        den:=g.denom;
        L:=_Variables(den);
        c:=Coefficients(den);
        gd:=MultPol(PA,[0],0);
        j:=0;
        for j in [1..Length(L)] do
           pos:=Position(L[j],k);
           if pos = FAILURE then
              gd:=gd+MultPol(PA,L[j],c[j]);
           else
              L2:=Remove(L[j],pos);
              if L2=[] then
                 gd:=gd+f*c[j];
              else 
                 gd:=gd+MultPol(PA,L2,c[j])*f;
              fi;
           fi; 
        od;
        if IsConstant(gd) then 
           return g2/Coerce(Base(PA),gd);
        else 
           return MultPol(PA,_Variables(g2),Coefficients(g2),gd);
        fi;
     fi;
     return g2;
   end;


InstallMethod(rec(
    kind:="FUNCTION",
    name:="Coerce",
    sin:=[[fld^fin,"F"],[elt-alg^mult,"f"]],
    sou:=[[elt-fld^fin,"f2"]],
    short:=
      "Coerces f into F, if possible ",
    ex:=[],
    see:=[]), _Coerce);

InstallMethod(
  rec(
  kind := "FUNCTION",
  name := "Coefficients",
  sin := [[elt-alg^mult]],
  sou := [[list]]
  ),_Coefficients);

InstallMethod(
  rec(
  kind := "FUNCTION",
  name := "CoefficientRing",
  sin := [[alg^mult]],
  sou := [[fld^fin]]
  ),_CoefficientRing);


InstallMethod(
  rec(
  kind := "FUNCTION",
  name := "Variables",
  sin := [[elt-alg^mult]],
  sou := [[list]]
  ),_Variables);

InstallMethod(rec(
    kind:="FUNCTION",
    name:="MultPol",
    sin:=[[alg^mult,"PAMult"],[list,"vars"],[list,"coeffs"],[elt-alg^mult,"denom"]],
    sou:=[[elt-alg^mult,"f"]],
    short:=
      "Returns an element of PAMult with variables vars, coefficients coeffs and denominator denom ",
    ex:=[],
    see:=[]), _MultPol2);

InstallMethod(rec(
    kind:="FUNCTION",
    name:="MultPol",
    sin:=[[alg^mult,"PAMult"],[list,"vars"],[elt-fld^fin,"coeffs"],[elt-alg^mult,"denom"]],
    sou:=[[elt-alg^mult,"f"]],
    short:=
      "Returns an element of PAMult with variables vars and coefficients coeffs ",
    ex:=[],
    see:=[]), _MultPol2);

InstallMethod(rec(
    kind:="FUNCTION",
    name:="MultPol",
    sin:=[[alg^mult,"PAMult"],[list,"vars"],[elt-ord^rat,"coeffs"],[elt-alg^mult,"denom"]],
    sou:=[[elt-alg^mult,"f"]],
    short:=
      "Returns an element of PAMult with variables vars and coefficients coeffs ",
    ex:=[],
    see:=[]), _MultPol2);

InstallMethod(rec(
    kind:="FUNCTION",
    name:="MultPol",
    sin:=[[alg^mult,"PAMult"],[list,"vars"],[list,"coeffs"]],
    sou:=[[elt-alg^mult,"f"]],
    short:=
      "Returns an element of PAMult with variables vars and coefficients coeffs ",
    ex:=[],
    see:=[]), _MultPol);

InstallMethod(rec(
    kind:="FUNCTION",
    name:="MultPol",
    sin:=[[alg^mult,"PAMult"],[list,"vars"],[elt-fld^fin,"coeffs"]],
    sou:=[[elt-alg^mult,"f"]],
    short:=
      "Returns an element of PAMult with variables vars and coefficients coeffs ",
    ex:=[],
    see:=[]), _MultPol);

InstallMethod(rec(
    kind:="FUNCTION",
    name:="MultPol",
    sin:=[[alg^mult,"PAMult"],[list,"vars"],[elt-ord^rat,"coeffs"]],
    sou:=[[elt-alg^mult,"f"]],
    short:=
      "Returns an element of PAMult with variables vars and coefficients coeffs ",
    ex:=[],
    see:=[]), _MultPol);