"Fossies" - the Fresh Open Source Software Archive

Member "KASH3-lib-archindep-2008-07-31/lib/MultPol.k" (3 Sep 2008, 16497 Bytes) of package /linux/misc/old/KASH3-lib-archindep-2008-07-31.tar.gz:


As a special service "Fossies" has tried to format the requested text file into HTML format (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file.

    1 
    2 #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
    3 
    4 
    5    _Coefficients:=function(f) 
    6       return Copy(Base(f).coeffs);
    7    end;
    8 
    9    _Variables:=function(f) 
   10       return Copy(Base(f).vars); 
   11    end;
   12 
   13    _CoefficientRing:= function(PA)
   14       return Base(PA);
   15    end;
   16 
   17    # this function is only for forwarding
   18    MultPol:=function(a)
   19       return a;
   20    end;
   21 
   22    _alg_mult_ops:=rec();
   23 
   24    _alg_mult_ops.Print:=function(al)
   25       Print("Polynomial Algebra over ",Base(al)," in ", al.vars," variables ");
   26    end;
   27 
   28    _elt_alg_mult_pol_ops:=rec(); 
   29 
   30    Merge:=function(L,c,left,middle,right)
   31      local i,j,m,n,L1,L2,c1,c2,L,c;
   32      L1:=[];
   33      L2:=[];
   34      c1:=[];
   35      c2:=[];
   36      m:=middle-left+1;
   37      n:=right-middle;
   38      for i in [1..m] do
   39         L1[i]:=L[left+i-1];
   40         c1[i]:=c[left+i-1];
   41      od;
   42      for i in [1..n] do
   43         L2[i]:=L[middle+i];
   44         c2[i]:=c[middle+i];
   45      od;
   46      i:=1;
   47      j:=1;
   48      while i<=m and j<=n do
   49         if L1[i]<L2[j] then
   50            L[left+i+j-2]:=L1[i];
   51            c[left+i+j-2]:=c1[i];
   52            i:=i+1;
   53         else 
   54            L[left+i+j-2]:=L2[j];
   55            c[left+i+j-2]:=c2[j];
   56            j:=j+1;
   57         fi;
   58      od;
   59      if j = n+1 then
   60         while i <=m do
   61            L[left+i+j-2]:= L1[i];
   62            c[left+i+j-2]:= c1[i];
   63            i:=i+1;
   64         od;
   65      fi;
   66    end;
   67 
   68    MergeSortRec:=function(L,c,first,ende)
   69      local middle;
   70      if first = ende then
   71         return;
   72      fi;
   73      middle:= Floor((first+ende)/2);
   74      MergeSortRec(L,c,first,middle);
   75      MergeSortRec(L,c,middle+1,ende);
   76      Merge(L,c,first,middle,ende);
   77    end;
   78 
   79   
   80    MergeSort:=function(L,c)
   81      return MergeSortRec(L,c,1,Length(L));
   82    end;
   83 
   84    _Reduction_:=function(f)
   85      local i,j,L,c;
   86      L:=Base(f).vars;
   87      c:=Base(f).coeffs;
   88      MergeSort(L,c);
   89      i:=0;
   90      while i <Length(L)-1 do
   91         i:=i+1;
   92         j:=i;
   93         while j<Length(L) do
   94            j:=j+1;
   95            if L[i]=L[j] then
   96               Remove_(L,j);
   97               c[i]:=c[i]+c[j];
   98               Remove_(c,j);
   99               j:=j-1;
  100            fi;
  101         od;
  102      od;
  103      i:=0;
  104      while i <Length(c) do
  105         i:=i+1;
  106         if c[i]=0 then
  107            Remove_(c,i);
  108            Remove_(L,i);
  109            i:=i-1;
  110         fi;
  111      od;
  112      if L=[] and c = [] then
  113         Base(f).vars:=[[0]];
  114         Base(f).coeffs:=[0];
  115      else
  116         Base(f).vars:=Set(L);
  117         Base(f).coeffs:=c;
  118      fi;
  119    end;
  120 
  121 
  122    _Coerce:=function(F,f)
  123      local L,c;
  124      L:=_Variables(f);
  125      c:=Coefficients(f);
  126      if Length(L)=1 and L[1][1]=0 then
  127         return Coerce(F,c[1]);
  128      else 
  129         return Error("cannot coerce ",f," into ",F);
  130      fi;
  131    end;
  132 
  133  
  134    MultPolAlg:=function(K,n)
  135      local al;
  136      al:=rec(base:=K,type:=alg^mult,operations:=_alg_mult_ops,vars:=n);
  137      return al;
  138    end;
  139 
  140 
  141    _MultPol:=function(PA,L,c)
  142      local f,i;
  143      if L=[] then
  144         if not IsList(c) then
  145            return rec(base:=rec(vars:=[[0]],coeffs:=[c]),type:=elt-alg^mult,parent:=PA,operations:=_elt_alg_mult_pol_ops);
  146         else 
  147            return Error("not enough information for creating a polynomial");
  148         fi;
  149      fi;
  150      if Is(Type(L[1]),list) then
  151         for i in [1..Length(L)] do
  152            L[i]:=Set(L[i]);
  153            if L[i][Length(L[i])]>PA.vars then
  154               return Error("too many variables for ",PA);
  155            fi;
  156         od;
  157      else 
  158         L:=[Set(L)];
  159         if L[1][Length(L[1])] > PA.vars then
  160            return Error("too many variables for ",PA);
  161         fi;
  162         c:=[c];
  163      fi;
  164      for i in [1..Length(c)] do
  165         c[i]:=Coerce(Base(PA),c[i]);
  166      od;
  167      if not Length(L) = Length(c) then
  168         if Length(L) > Length(c) then
  169            return Error("not enough coefficients");
  170         else return Error("too many coefficients");
  171         fi;
  172      fi;
  173      f:=rec(base:=rec(vars:=L,coeffs:=c),type:=elt-alg^mult,parent:=PA,operations:=_elt_alg_mult_pol_ops);
  174      _Reduction_(f);
  175      return f;
  176    end;
  177 
  178 
  179    _MultPol2:=function(PA,L,c,d)
  180      local f,i;
  181      if Is(Type(L[1]),list) then
  182         for i in [1..Length(L)] do
  183            L[i]:=Set(L[i]);
  184            if L[i][Length(L[i])]>PA.vars then
  185               return Error("too many variables for ",PA);
  186            fi;
  187         od;
  188      else L:=[Set(L)];
  189         if L[1][Length(L[1])] > PA.vars then
  190            return Error("too many variables for ",PA);
  191         fi; 
  192         c:=[c];
  193      fi;
  194      for i in [1..Length(c)] do
  195         c[i]:=Coerce(Base(PA),c[i]);
  196      od;
  197      if not Length(L) = Length(c) then
  198         if Length(L) > Length(c) then
  199            return Error("not enough coefficients");
  200         else return Error("too many coefficients");
  201         fi;
  202      fi;
  203      f:=rec(base:=rec(vars:=L,coeffs:=c),type:=elt-alg^mult,parent:=PA,operations:=_elt_alg_mult_pol_ops,denom:=d);
  204      _Reduction_(f);
  205      return f;
  206    end;
  207 
  208 
  209 
  210    _elt_alg_mult_pol_ops.Print:= function(f)
  211      local L,c,i,j,den;
  212      L:=_Variables(f);
  213      c:=Coefficients(f);
  214      if "denom" in RecFields(f) then
  215         if f.denom<>1 then
  216            den:=true;
  217            Print("(");
  218         else den:=false;
  219         fi;
  220      else den := false;
  221      fi;
  222      for i in [1..Length(L)] do
  223         Print(c[i]);
  224         if L[i][1] <> 0 then
  225            Print("*");
  226            for j in [1..Length(L[i])] do
  227               Print("x",L[i][j]);
  228               if not j = Length(L[i]) then
  229                  Print("*");
  230               fi;
  231            od;
  232         fi;
  233         if not i = Length(L) then
  234            Print(" + ");
  235         fi;
  236      od;
  237      if den then
  238         Print(")/(");
  239         Print(f.denom);
  240         Print(")");
  241      fi;
  242    end;
  243 
  244    _elt_alg_mult_pol_ops.\+:= function(f1,f2)
  245      local L1,L2,c1,c2,g,denom,f11,f21;
  246      if Is(Type(f1),elt-alg^mult) and Is(Type(f2),elt-alg^mult) then
  247         if "denom" in RecFields(f1) and "denom" in RecFields(f2) then
  248            f11:=f1*f2.denom;
  249            f21:=f2*f1.denom;
  250            denom:=f1.denom*f2.denom;
  251         elif "denom" in RecFields(f1) then
  252            f11:=f1;
  253            f21:=f2*f1.denom;
  254            denom:=f1.denom;
  255         elif "denom" in RecFields(f2) then
  256            f11:=f1*f2.denom;
  257            f21:=f2;
  258            denom:=f2.denom;
  259         else 
  260            f11:=f1;
  261            f21:= f2;
  262            denom:=1;
  263         fi;
  264         L1:=_Variables(f11);
  265         L2:=_Variables(f21);
  266         c1:=Coefficients(f11);
  267         c2:=Coefficients(f21);
  268         if denom <>1 then
  269            g:=MultPol(Parent(f1),Append(L1,L2),Append(c1,c2),denom);
  270         else 
  271            g:=MultPol(Parent(f1),Append(L1,L2),Append(c1,c2));
  272         fi;
  273         return g;
  274      elif Is(Type(f1),elt-alg^mult) then
  275         if not Is(Type(f2),elt-fld^fin) and not  Is(Type(f2),elt-ord^rat) then
  276            return Error("wrong type in addition");
  277         fi;
  278         if "denom" in RecFields(f1) then
  279            f21:=f1.denom*f2;
  280            return f21+f1;
  281         else return MultPol(Parent(f1),Append(_Variables(f1),[[0]]),Append(Coefficients(f1),[f2]));
  282         fi;
  283      elif Is(Type(f2),elt-alg^mult) then
  284         if not Is(Type(f1),elt-fld^fin) and not Is(Type(f1),elt-ord^rat) then
  285            return Error("wrong type in addition");
  286         fi;
  287         if "denom" in RecFields(f2) then
  288            f11:=f2.denom*f1;   
  289            return f11+f2;
  290         else
  291            return MultPol(Parent(f2),Append(_Variables(f2),[[0]]),Append(Coefficients(f2),[f1]));
  292         fi;
  293      fi;
  294    end;
  295 
  296    _elt_alg_mult_pol_ops.\^:= function(f1,n)
  297      local L1,L2,c1,c2,g,i;
  298      if Is(Type(f1),elt-alg^mult) and Is(Type(n),elt-ord^rat) then
  299         g:=f1;
  300         for i in [2..n] do
  301            g:=g*f1;
  302         od;
  303         return g;
  304      else return Error("wrong type in ^");
  305      fi;
  306    end;
  307 
  308    _elt_alg_mult_pol_ops.\-:= function(f1,f2)
  309      return f1+f2; #nur fuer den Fall K=GF(2,n)
  310    end;
  311 
  312    _elt_alg_mult_pol_ops.\/:=function(f1,f2)
  313      local L,c,PA,i;
  314      PA:=Parent(f1);
  315      if Is(Type(f2),elt-ord^rat) or Is(Type(f2),elt-fld^fin) then
  316         f2:= Coerce(PA,f2);
  317         L:= _Variables(f1);
  318         c:= Coefficients(f1);
  319         for i in [1..Length(c)] do
  320            c[i]:= c[i]/f2;
  321         od;
  322         return MultPol(PA,L,c);
  323      else return Error("wrong type in / ");
  324      fi;
  325    end;
  326 
  327    _elt_alg_mult_pol_ops.\=:= function(f1,f2)
  328      local L,c;
  329      if Is(Type(f1),elt-alg^mult) and Is(Type(f2),elt-alg^mult) then
  330         return Coefficients(f1)=Coefficients(f2) and _Variables(f1)=_Variables(f2);
  331      elif Is(Type(f1),elt-fld^fin) or Is(Type(f1),elt-ord^rat) then
  332         L:=_Variables(f2);
  333         if Length(L) >1 then
  334            return false;
  335         elif Length(L[1])>1 then
  336            return false;
  337         elif L[1][1] <> 0 then
  338            return false;
  339         else
  340            c:=Coefficients(f2);
  341            f1:=Coerce(Parent(f2),f1);
  342            return c[1] = f1;
  343         fi;
  344      elif Is(Type(f2),elt-fld^fin) or Is(Type(f2),elt-ord^rat) then
  345         L:=_Variables(f1);
  346         if Length(L) >1 then
  347            return false;
  348         elif Length(L[1])>1 then
  349            return false;
  350         elif L[1][1] <> 0 then
  351            return false;
  352         else
  353            c:=Coefficients(f1);
  354            f1:=Coerce(Parent(f1),f2);
  355            return c[1] = f2;
  356         fi;
  357      else return Error("wrong type in = ");
  358      fi;
  359    end;
  360 
  361    _elt_alg_mult_pol_ops.\*:= function(f1,f2)
  362      local PA,f12,g,L,L1,L2,c,c1,c2,i,j,pos,dL1,dL2,dc1,dc2,dL,dc,denom;
  363      if Is(Type(f1),elt-ord^rat) or Is(Type(f1),elt-fld^fin) then
  364         PA:=Parent(f2);
  365         f1:=Coerce(PA,f1);
  366         c:=Coefficients(f2);
  367         L:=_Variables(f2);
  368         for i in [1..Length(c)] do
  369            c[i]:= f1*c[i];
  370         od;
  371         if "denom" in RecFields(f2) then
  372            return MultPol(PA,L,c,f2.denom);
  373         else return MultPol(PA,L,c);
  374         fi;
  375      elif Is(Type(f2),elt-ord^rat) or Is(Type(f2),elt-fld^fin) then
  376         PA:=Parent(f1);
  377         f2:=Coerce(PA,f2);
  378         c:=Coefficients(f1);
  379         L:=_Variables(f1);
  380         for i in [1..Length(c)] do
  381            c[i]:= f2*c[i];
  382         od; 
  383         if "denom" in RecFields(f1) then
  384            return MultPol(PA,L,c,f1.denom);
  385         else return MultPol(PA,L,c);
  386         fi;
  387      elif Is(Type(f1),elt-alg^mult) and Is(Type(f2),elt-alg^mult) then
  388         PA:= Parent(f1);
  389         L1:= _Variables(f1);
  390         L2:= _Variables(f2);
  391         c1:= Coefficients(f1);
  392         c2:= Coefficients(f2);
  393         L:=[];
  394         c:=[];
  395         for i in [1..Length(L1)] do
  396            for j in [1..Length(L2)] do
  397               Append_(L,[Set(Append(L1[i],L2[j]))]);
  398               Append_(c,[c1[i]*c2[j]]);
  399               if Length(Last(L))>1 then
  400                  pos:=Position(Last(L),0);
  401                  if pos <> FAILURE then
  402                     Remove_(Last(L),pos);
  403                  fi;
  404               fi;
  405            od;
  406         od;
  407         if "denom" in RecFields(f1) and "denom" in RecFields(f2) then
  408            if f1.denom = 1 and f2.denom =1 then
  409               return MultPol(PA,L,c);
  410            elif f1.denom <>1 then
  411               if f2.denom = 1 then
  412                  return MultPol(PA,L,c,f1.denom);
  413               else
  414                  dL1:=_Variables(f1.denom);
  415                  dL2:=_Variables(f2.denom);
  416                  dc1:=Coefficients(f1.denom);
  417                  dc2:=Coefficients(f2.denom);
  418                  dL:=[];
  419                  dc:=[];
  420                  for i in [1..Length(dL1)] do
  421                     for j in [1..Length(dL2)] do
  422                        Append_(dL,[Set(Append(dL1[i],dL2[j]))]);
  423                        Append_(dc,[dc1[i]*dc2[j]]);
  424                        if Length(Last(dL))>1 then
  425                           pos:=Position(Last(dL),0);
  426                           if pos <> FAILURE then
  427                              Remove_(Last(dL),pos);
  428                           fi;
  429                        fi;
  430                     od;
  431                  od;
  432                  denom:=MultPol(PA,dL,dc);
  433                  return MultPol(PA,L,c,denom);
  434               fi;
  435            else
  436               return MultPol(PA,L,c,f2.denom);
  437            fi;
  438         elif "denom" in RecFields(f1) then
  439            return MultPol(PA,L,c,f1.denom);
  440         elif "denom" in RecFields(f2) then
  441            return MultPol(PA,L,c,f2.denom);
  442         else 
  443            return  MultPol(PA,L,c);
  444         fi;
  445      else return Error("wrong types in multiplication");
  446      fi;
  447    end;
  448 
  449    IsConstant:=function(g)
  450      if Is(Type(g),elt-ord^rat) or Is(Type(g),elt-fld^fin) then 
  451         return true;
  452      else 
  453         return (Length(_Variables(g)) =1 and Last(_Variables(g)[1]) =0);
  454      fi;
  455    end;
  456 
  457    Elimination_:=function(g,k,f)
  458      local i,L,c,j,pos,g2,PA,denom,den,L2,gd;
  459      if IsConstant(g) then 
  460         return g;
  461      fi;
  462      PA:=Parent(g);
  463      L:=_Variables(g);
  464      c:=Coefficients(g);
  465      g2:=MultPol(PA,[0],0);
  466      j:=0;
  467      denom:=("denom" in RecFields(g));      
  468      for j in [1..Length(L)] do
  469         pos:=Position(L[j],k);
  470         if pos = FAILURE then
  471            g2:=g2+MultPol(PA,L[j],c[j]);
  472         else
  473            L2:=Remove(L[j],pos);
  474            if L2=[] then
  475               g2:=g2+f*c[j];
  476            else 
  477               g2:=g2+MultPol(PA,L2,c[j])*f;
  478            fi;
  479         fi; 
  480      od;
  481      if denom then 
  482         den:=g.denom;
  483         L:=_Variables(den);
  484         c:=Coefficients(den);
  485         gd:=MultPol(PA,[0],0);
  486         j:=0;
  487         for j in [1..Length(L)] do
  488            pos:=Position(L[j],k);
  489            if pos = FAILURE then
  490               gd:=gd+MultPol(PA,L[j],c[j]);
  491            else
  492               L2:=Remove(L[j],pos);
  493               if L2=[] then
  494                  gd:=gd+f*c[j];
  495               else 
  496                  gd:=gd+MultPol(PA,L2,c[j])*f;
  497               fi;
  498            fi; 
  499         od;
  500         if IsConstant(gd) then 
  501            return g2/Coerce(Base(PA),gd);
  502         else 
  503            return MultPol(PA,_Variables(g2),Coefficients(g2),gd);
  504         fi;
  505      fi;
  506      return g2;
  507    end;
  508 
  509 
  510 InstallMethod(rec(
  511     kind:="FUNCTION",
  512     name:="Coerce",
  513     sin:=[[fld^fin,"F"],[elt-alg^mult,"f"]],
  514     sou:=[[elt-fld^fin,"f2"]],
  515     short:=
  516       "Coerces f into F, if possible ",
  517     ex:=[],
  518     see:=[]), _Coerce);
  519 
  520 InstallMethod(
  521   rec(
  522   kind := "FUNCTION",
  523   name := "Coefficients",
  524   sin := [[elt-alg^mult]],
  525   sou := [[list]]
  526   ),_Coefficients);
  527 
  528 InstallMethod(
  529   rec(
  530   kind := "FUNCTION",
  531   name := "CoefficientRing",
  532   sin := [[alg^mult]],
  533   sou := [[fld^fin]]
  534   ),_CoefficientRing);
  535 
  536 
  537 InstallMethod(
  538   rec(
  539   kind := "FUNCTION",
  540   name := "Variables",
  541   sin := [[elt-alg^mult]],
  542   sou := [[list]]
  543   ),_Variables);
  544 
  545 InstallMethod(rec(
  546     kind:="FUNCTION",
  547     name:="MultPol",
  548     sin:=[[alg^mult,"PAMult"],[list,"vars"],[list,"coeffs"],[elt-alg^mult,"denom"]],
  549     sou:=[[elt-alg^mult,"f"]],
  550     short:=
  551       "Returns an element of PAMult with variables vars, coefficients coeffs and denominator denom ",
  552     ex:=[],
  553     see:=[]), _MultPol2);
  554 
  555 InstallMethod(rec(
  556     kind:="FUNCTION",
  557     name:="MultPol",
  558     sin:=[[alg^mult,"PAMult"],[list,"vars"],[elt-fld^fin,"coeffs"],[elt-alg^mult,"denom"]],
  559     sou:=[[elt-alg^mult,"f"]],
  560     short:=
  561       "Returns an element of PAMult with variables vars and coefficients coeffs ",
  562     ex:=[],
  563     see:=[]), _MultPol2);
  564 
  565 InstallMethod(rec(
  566     kind:="FUNCTION",
  567     name:="MultPol",
  568     sin:=[[alg^mult,"PAMult"],[list,"vars"],[elt-ord^rat,"coeffs"],[elt-alg^mult,"denom"]],
  569     sou:=[[elt-alg^mult,"f"]],
  570     short:=
  571       "Returns an element of PAMult with variables vars and coefficients coeffs ",
  572     ex:=[],
  573     see:=[]), _MultPol2);
  574 
  575 InstallMethod(rec(
  576     kind:="FUNCTION",
  577     name:="MultPol",
  578     sin:=[[alg^mult,"PAMult"],[list,"vars"],[list,"coeffs"]],
  579     sou:=[[elt-alg^mult,"f"]],
  580     short:=
  581       "Returns an element of PAMult with variables vars and coefficients coeffs ",
  582     ex:=[],
  583     see:=[]), _MultPol);
  584 
  585 InstallMethod(rec(
  586     kind:="FUNCTION",
  587     name:="MultPol",
  588     sin:=[[alg^mult,"PAMult"],[list,"vars"],[elt-fld^fin,"coeffs"]],
  589     sou:=[[elt-alg^mult,"f"]],
  590     short:=
  591       "Returns an element of PAMult with variables vars and coefficients coeffs ",
  592     ex:=[],
  593     see:=[]), _MultPol);
  594 
  595 InstallMethod(rec(
  596     kind:="FUNCTION",
  597     name:="MultPol",
  598     sin:=[[alg^mult,"PAMult"],[list,"vars"],[elt-ord^rat,"coeffs"]],
  599     sou:=[[elt-alg^mult,"f"]],
  600     short:=
  601       "Returns an element of PAMult with variables vars and coefficients coeffs ",
  602     ex:=[],
  603     see:=[]), _MultPol);