// last change ML: 12.08.99
///////////////////////////////////////////////////////////////////////////////
// This library is for Singular 1.2 or newer

version="$Id: primitiv.lib,v 1.16 2001-08-27 14:47:59 Singular Exp $";
category="Commutative Algebra";
info="
LIBRARY:    primitiv.lib    Computing a Primitive Element
AUTHOR:     Martin Lamm,    email: lamm@mathematik.uni-kl.de

PROCEDURES:
 primitive(ideal i);   find minimal polynomial for a primitive element
 primitive_extra(i);   find primitive element for two generators
 splitring(f,R[,L]);   define ring extension with name R and switch to it
";

LIB "random.lib";
///////////////////////////////////////////////////////////////////////////////

proc primitive(ideal i)
"USAGE:   primitive(i); i ideal
ASSUME:   i is given by generators m[1],...,m[n] such that for j=1,...,n @*
          -  m[j] is a polynomial in k[x(1),...,x(j)] @*
          -  m[j](a[1],...,a[j-1],x(j)) is the minimal polynomial for a[j] over
             k(a[1],...,a[j-1]) @*
          (k the ground field of the current basering and x(1),...,x(n)
          the ring variables).
RETURN:   ideal j in k[x(n)] with
          - j[1] a minimal polynomial for a primitive element b of
                 k(a[1],...,a[n]) over k,
          - j[2],...,j[n+1] polynomials in k[x(n)] such that j[i+1](b)=a[i]
                 for i=1,...,n.
NOTE:     the number of variables in the basering has to be exactly n,
          the number of given generators (i.e., minimal polynomials).@*
          If the ground field k has only a few elements it may happen that no
          linear combination of a[1],...,a[n] is a primitive element. In this
          case @code{primitive(i)} returns the zero ideal, and one should use
          @code{primitive_extra(i)} instead.
SEE ALSO: primitive_extra
KEYWORDS: primitive element
EXAMPLE:  example primitive;  shows an example
"
{
 def altring=basering;
 execute("ring deglexring=("+charstr(altring)+"),("+varstr(altring)+"),dp;");
 ideal j;
 execute("ring lexring=("+charstr(altring)+"),("+varstr(altring)+"),lp;");
 ideal i=fetch(altring,i);

 int k,schlecht,Fehlversuche,maxtry;
 int nva = nvars(basering);
 int p=char(basering);
 if (p==0) {
   p=100000;
   if (nva<3) { maxtry= 100000000; }
   else       { maxtry=2147483647; }
 }
 else {
   if ((nva<4) || (p<60)) {
     maxtry=p^(nva-1); }
   else {
     maxtry=2147483647;          // int overflow(^)  vermeiden
   }
 }
 ideal jmap,j;
 map phi;
 option(redSB);

 //-------- Mache so lange Random-Koord.wechsel, bis letztes Poly -------------
 //--------------- das Minpoly eines primitiven Elements ist : ----------------
 for (Fehlversuche=0; Fehlversuche<maxtry; Fehlversuche++) {
   schlecht=0;
   if ((p<60) && (nva==2)) {  // systematische Suche statt random
      jmap=ideal(var(1),var(2)+Fehlversuche*var(1));
   }
   else {
    if (Fehlversuche==0) { jmap=maxideal(1);}
    else {
      if (Fehlversuche<5) { jmap=randomLast(10);}
      else {
       if (Fehlversuche<20) { jmap=randomLast(100);}
       else                 { jmap=randomLast(100000000);}
    }}                        // groessere Werte machen keinen Sinn
   }
   phi=lexring,jmap;
   j=phi(i);
   setring deglexring;
 //--------------- Berechne reduzierte Standardbasis mit fglm: ----------------
   j=std(fetch(lexring,j));
   setring lexring;
   j=fglm(deglexring,j);
 //-- teste, ob SB n Elemente enthaelt (falls ja, ob lead(Fi)=xi i=1... n-1): -
   if (size(j)==nva) {
     for (k=1; k<nva; k++) {
       j[k+1]=j[k+1]/leadcoef(j[k+1]);        // normiere die Erzeuger
       if (lead(j[k+1]) != var(nva-k)) { schlecht=1;}
     }
     if (schlecht==0) {
 //--- Random-Koord.wechsel war gut: Berechne das zurueckzugebende Ideal: -----
       ideal erg;
       for (k=1; k<nva; k++) { erg[k]=var(k)-j[nva-k+1]; }
                               // =g_k(x_n) mit a_k=g_k(a_n)
       erg[nva]=var(nva);
       map chi=lexring,erg;
       ideal extra=maxideal(1);extra=phi(extra);
                               // sonst: "argument of a map must have a name"
       erg=j[1],chi(extra);    // j[1] = Minimalpolynom
       setring altring;
       return(fetch(lexring,erg));
     }
   }
   dbprint("The random coordinate change was bad!");
 }
 if (voice==2) {
   "// ** Warning: No primitive element could be found.";
   "//    If the given ideal really describes the minimal polynomials of";
   "//    a series of algebraic elements (cf. `help primitive;') then";
   "//    try `primitive_extra'.";
 }
 setring altring;
 return(ideal(0));
}
example
{ "EXAMPLE:"; echo = 2;
 ring exring=0,(x,y),dp;
 ideal i=x2+1,y2-x;                  // compute Q(i,i^(1/2))=:L
 ideal j=primitive(i);
 j[1];                               // L=Q(a) with a=(-1)^(1/4)
 j[2];                               // i=a^2
 j[3];                               // i^(1/2)=a
 // the 2nd element was already primitive!
 j=primitive(ideal(x2-2,y2-3));      // compute Q(sqrt(2),sqrt(3))
 j[1];
 j[2];
 j[3];
 // no element was primitive -- the calculation of primitive elements
 // is based on a random choice.
}
///////////////////////////////////////////////////////////////////////////////

proc primitive_extra(ideal i)
"USAGE:   primitive_extra(i); i ideal
ASSUME:  The ground field of the basering is k=Q or k=Z/pZ and the ideal
         i is given by 2 generators f,g with the following properties:
@format
   f is the minimal polynomial of a in k[x],
   g is a polynomial in k[x,y] s.th. g(a,y) is the minpoly of b in k(a)[y].
@end format
          Here, x is the name of the first ring variable, y the name of the
          second.
RETURN:  ideal j in k[y] such that
@format
   j[1] is the minimal polynomial for a primitive element c of k(a,b) over k,
   j[2] is a polynomial s.th. j[2](c)=a.
@end format
NOTE:    While @code{primitive(i)} may fail for finite fields,
         @code{primitive_extra(i)} tries all elements of k(a,b) and, hence,
         always finds a primitive element. @*
         In order to do this (try all elements), field extensions like Z/pZ(a)
         are not allowed for the ground field k. @*
         @code{primitive_extra(i)} assumes that the second generator, g, is
         monic as polynomial in (k[x])[y].
EXAMPLE: example primitive_extra;  shows an example
"
{
 def altring=basering;
 int grad1,grad2=deg(i[1]),deg(jet(i[2],0,intvec(1,0)));
 int countx,countz;
 ring deglexring=char(altring),(x,y,z),dp;
 map transfer=altring,x,z;
 ideal i=transfer(i);
 if (size(i)!=2) {
   "//** Error -- either wrong number of given minimal polynomials";
   "//**          or wrong choice of ring variables (must use the first two)";
   setring altring;
   return(ideal(0));
 }
 matrix mat;
 ring lexring=char(altring),(x,y),lp;
 ideal j;
 ring deglex2ring=char(altring),(x,y),dp;
 ideal j;
 setring deglexring;
 ideal j;
 option(redSB);
 poly g=z;
 int found=0;

 //---------------- Schleife zum Finden des primitiven Elements ---------------
 //--- Schleife ist so angordnet, dass g in Charakteristik 0 linear bleibt ----
 while (found==0) {
   j=eliminate(i+ideal(g-y),z);
   setring deglex2ring;
   j=std(imap(deglexring,j));
   setring lexring;
   j=fglm(deglex2ring,j);
   if (size(j)==2) {
     if (deg(j[1])==grad1*grad2) {
       j[2]=j[2]/leadcoef(j[2]);    // Normierung
       if (lead(j[2])==x) {         // Alles ok
          found=1;
       }
     }
   }
   setring deglexring;
   if (found==0) {
 //------------------ waehle ein neues Polynom g ------------------------------
     dbprint("Still searching for primitive element...");
     countx=0;
     countz=0;
     while (found==0) {
      countx++;
      if (countx>=grad1) {
        countx=0;
        countz++;
        if (countz>=grad2) {
         "//** Error: No primitive element found!! This should NEVER happen!";
         setring altring;
         return(ideal(0));
        }
      }
      g = g +x^countx *z^countz;
      mat=coeffs(g,z);
      if (size(mat)>countz) {
        mat=coeffs(mat[countz+1,1],x);
        if (size(mat)>countx) {
          if (mat[countx+1,1] != 0) {
            found=1;         // d.h. hier: neues g gefunden
      }}}
     }
     found=0;
   }
 }
 //------------------- primitives Element gefunden; Rueckgabe -----------------
 setring lexring;
 j[2]=x-j[2];
 setring altring;
 map transfer=lexring,var(1),var(2);
 return(transfer(j));
}
example
{ "EXAMPLE:"; echo = 2;
 ring exring=3,(x,y),dp;
 ideal i=x2+1,y3+y2-1;
 primitive_extra(i);
 ring extension=(3,y),x,dp;
 minpoly=y6-y5+y4-y3-y-1;
 number a=y5+y4+y2+y+1;
 a^2;
 factorize(x2+1);
 factorize(x3+x2-1);
}
///////////////////////////////////////////////////////////////////////////////

proc splitring
"USAGE:   splitring(f,R[,L]); f poly, R string, L list of polys and/or ideals
         (optional)
ASSUME:  f is univariate and irreducible over the active basering. @*
         The active ring must allow an algebraic extension (e.g., it cannot
         be a transcendent ring extension of Q or Z/p).
CREATE:  a ring with name R, in which f is reducible, and CHANGE to it.
RETURN:  list L mapped into the new ring R, if L is given; else nothing
NOTE:    If the old ring has no parameter, the name @code{a} is chosen for the
         parameter of R (if @code{a} is no ring variable; if it is, @code{b} is
         chosen, etc.; if @code{a,b,c,o} are ring variables,
         @code{splitring(f,R[,L])} produces an error message), otherwise the
         name of the parameter is kept and only the minimal polynomial is
         changed. @*
         The names of the ring variables and the orderings are not affected. @*
         It is also allowed to call @code{splitring} with R=\"\".
         Then the old basering will be REPLACED by the new ring (with the
         same name as the old ring).
KEYWORDS: algebraic field extension; extension of rings
EXAMPLE: example splitring;  shows an example
"
{
 //----------------- split ist bereits eine proc in 'inout.lib' ! -------------
 poly f=#[1]; string @R=#[2];
 if (size(#)>2) {
    list L=#[3];
    int L_groesse=size(L);
 }
 else { int L_groesse=-1; }
 //-------------- ermittle das Minimalpolynom des aktuellen Rings: ------------
 string minp=string(minpoly);

 if (@R=="") {
  string altrname=nameof(basering);
  @R="splt_temp";
 }

 def altring=basering;
 string charakt=string(char(altring));
 string varnames=varstr(altring);
 string algname;
 int i;
 int anzvar=size(maxideal(1));
 //--------------- Fall 1: Bisheriger Ring hatte kein Minimalpolynom ----------
 if (minp=="0") {
  if (find(varnames,"a")==0)        { algname="a";}
  else { if (find(varnames,"b")==0) { algname="b";}
         else { if (find(varnames,"c")==0)
                                    { algname="c";}
         else { if (find(varnames,"o")==0)
                                    { algname="o";}
         else {
           "** Sorry -- could not find a free name for the primitive element.";
           "** Try e.g. a ring without 'a' or 'b' as variable.";
           return();
         }}
       }
  }
 //-- erzeuge einen String, der das Minimalpolynom des neuen Rings enthaelt: --
  execute("ring splt1="+charakt+","+algname+",dp;");
  ideal abbnach=var(1);
  for (i=1; i<anzvar; i++) { abbnach=abbnach,var(1); }
  map nach_splt1=altring,abbnach;
  execute("poly mipol="+string(nach_splt1(f))+";");
  string Rminp=string(mipol);
 //--------------------- definiere den neuen Ring: ----------------------------
  execute("ring "+@R+" = ("+charakt+","+algname+"),("+varnames+"),("
           +ordstr(altring)+");");
  execute("minpoly="+Rminp+";");
  execute("export "+@R+";");
  def neuring=basering;
 //---------------------- Berechne die zurueckzugebende Liste: ----------------
  list erg;
  if (L_groesse>0) {
 // L ist ja nicht in 'neuring' def., daher merke man sich die Groesse als int
   map take=altring,maxideal(1);
   erg=take(L);
  }            // take(empty list) gibt nicht empty list, sondern Fehlermeldung
 }
 else {
 //------------- Fall 2: Bisheriger Ring hatte ein Minimalpolynom: ------------
  algname=parstr(altring);           // Name des algebraischen Elements
  if (size(algname)>1) {"only one Parameter is allowed!!"; return();}
 //---------------- Minimalpolynom in ein Polynom umwandeln: ------------------
  execute("ring splt2="+charakt+","+algname+",dp;");
  execute("poly mipol="+minp+";");
 // f ist Polynom in algname und einer weiteren Variablen --> mache f bivariat:
  execute("ring splt3="+charakt+",("+algname+","+varnames+"),dp;");
  poly f=imap(altring,f);
 //-------------- Vorbereitung des Aufrufes von primitive: --------------------
  execute("ring splt1="+charakt+",(x,y),dp;");
  ideal abbnach=x;
  for (i=1; i<=anzvar; i++) { abbnach=abbnach,y; }
  map nach_splt1_3=splt3,abbnach;
  map nach_splt1_2=splt2,x;
  ideal maxid=nach_splt1_2(mipol),nach_splt1_3(f);
  ideal primit=primitive(maxid);
  if (size(primit)==0) {             // Suche mit 1. Proc erfolglos
    primit=primitive_extra(maxid);
  }
 //-- erzeuge einen String, der das Minimalpolynom des neuen Rings enthaelt: --
  setring splt2;
  map nach_splt2=splt1,0,var(1);     // x->0, y->a
  minp=string(nach_splt2(primit)[1]);
  if (printlevel > -1) { "// new minimal polynomial:",minp; }
 //--------------------- definiere den neuen Ring: ----------------------------
  execute("ring "+@R+" = ("+charakt+","+algname+"),("+varnames+"),("
          +ordstr(altring)+");");
  execute("minpoly="+minp+";");
  execute("export "+@R+";");
  def neuring=basering;

 //--------------- Uebersicht: wenn altring=(p,a),(x,y),dp; dann: -------------
 //------------ splt1=p,(x,y),dp;  splt2=p,a,dp;  splt3=p,(a,x,y),dp; ---------

  list erg;
  if (L_groesse>0) {
 //---------------------- Berechne die zurueckzugebende Liste: ----------------
    setring splt3;
    list zwi=imap(altring,L);
    map nach_splt3_1=splt1,0,var(1);  // x->0, y->a
 //----- rechne das primitive Element von altring in das von neuring um: ------
    ideal convid=maxideal(1);
    convid[1]=nach_splt3_1(primit)[2];
    map convert=splt3,convid;
    zwi=convert(zwi);
    setring neuring;
    erg=imap(splt3,zwi);
  }
 }
 if (defined(altrname)) {
   if(system("with","Namespaces"))
   { kill Top::`altrname`; kill Top::splt_temp; }
   execute("kill "+altrname+";");
   execute("def "+altrname+" = splt_temp;");
   @R=altrname;
   execute("export "+altrname+";");
   kill splt_temp;
 }

 execute("keepring "+@R+";");
 if (L_groesse >= 0) { return(erg); }
}
example
{ "EXAMPLE:"; echo = 2;
 ring r=0,(x,y),dp;
 splitring(x2-2,"r1");   // change to Q(sqrt(2))
 // change to Q(sqrt(2),sqrt(sqrt(2)))=Q(a) and return the transformed
 // old parameter:
 splitring(x2-a,"r2",a);
 // the result is (a)^2 = (sqrt(sqrt(2)))^2
 nameof(basering);
 r2;
 kill r1; kill r2;
}
///////////////////////////////////////////////////////////////////////////////