// $Id: primitiv.lib,v 1.8 1999-04-14 12:11:53 Singular Exp $
// author:  Martin Lamm,  email: lamm@mathematik.uni-kl.de
// last change:            19.03.99
///////////////////////////////////////////////////////////////////////////////
// This library is for Singular 1.2 or newer

version="$Id: primitiv.lib,v 1.8 1999-04-14 12:11:53 Singular Exp $";
info="
LIBRARY:    primitiv.lib    PROCEDURES FOR FINDING A PRIMITIVE ELEMENT

 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
 randomLast(b);        random transformation of the last variable
";

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

proc randomLast(int b)
"USAGE:   randomLast(b);  b int
RETURN:  ideal = maxideal(1), but the last variable is exchanged by
         a sum of it with a linear random combination of the other
         variables. The coefficients are in the interval [-b,b].
EXAMPLE: example randomLast; shows an example
"
{
  ideal i=maxideal(1);
  int k=size(i);
  i[k]=0;
  i=randomid(i,size(i),b);
  ideal ires=maxideal(1);
  ires[k]=i[1]+var(k);
  return(ires);
}
example
{ "EXAMPLE:"; echo = 2;
   ring  r = 0,(x,y,z),lp;
   ideal i = randomLast(10);
   i;
}
///////////////////////////////////////////////////////////////////////////////

proc primitive(ideal i)
"USAGE:   primitive(i); i ideal of the following form:
 Let k be the ground field of your basering, a_1,...,a_n algebraic over k,
 m_1(x_1), m_2(x_1,x_2),...,m_n(x_1,...,x_n) polynomials in k such that
 m_j(a_1,...,a_(j-1),x_j) is minimal polynomial for a_j over k(a_1,...,a_(j-1))
                                                        for all j=1,...,n.
 Then i has to be generated by m_1,...,m_n.

RETURN:  ideal j in k[x_n] such that
 j[1] is minimal polynomial for a primitive element b of k(a_1,...,a_n)=k(b)
         over k
 j[2],...,j[n+1] polynomials in k[x_n] : j[i+1](b)=a_i for i=1,...,n
NOTE:    the number of variables in the basering has to be exactly the number n
         of given algebraic elements (and minimal polynomials).

         If k has few elements it can be that no one of the linear combinations
         of a_1,...,a_n is a primitive element. In this case `primitive'
         returns the zero ideal. If this occurs use primitive_extra instead.
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);               // -> we have L=Q(a):
 "minimal polynomial of a:",j[1];    // => a=(-1)^(1/4)
 "polynomial for i:       ",j[2];    // => i=a^2
 "polynomial for i^(1/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))
 "minimal polynomial:",j[1];
 "polynomial p s.t. p(a)=sqrt(2):",j[2];
 "polynomial r s.t. r(a)=sqrt(3):",j[3];
 // ==> no element was primitive -- the calculation of a is based on a random
 //     choice.
}
///////////////////////////////////////////////////////////////////////////////

proc primitive_extra(ideal i)
"USAGE:   primitive_extra(i);  ideal i=f,g;  with the following properties:
 Let k=Q or k=Z/pZ be the ground field of the basering, a,b algebraic over k,
 x the name of the first ring variable, y the name of the second, then:
 f is the minimal polynomial of a in k[x], g is a polynomial in k[x,y] s.t.
 g(a,y) is the minimal polynomial of b in k(a)[y]

RETURN:  ideal j in k[y] such that
 j[1] is minimal polynomial over k for a primitive element c of k(a,b)=k(c)
 j[2] is a polynomial s.t. j[2](c)=a

NOTE: - While `primitive' may fail for finite fields, this proc tries all
        elements of k(a,b) and hence finds by assurance 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.
      - primitive_extra assumes that 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, univariate, irreducible(!), R string,
                     L list of polys and/or ideals (optional)
CREATE: defines a ring with name R, in which f is reducible, and changes to it
        If the old ring has no parameter, the name 'a' is chosen for the
        parameter of R (if a is no variable; if it is, the proc takes 'b',
        etc.; if a,b,c,o are variables of the ring, produce an error message),
        otherwise the name of the parameter is kept and only the
        minimal polynomial is changed.
        The names of variables and orderings are not affected.

        It is also allowed to call splitring with R==\"\". Then the old basering
        will be REPLACED by the new ring (with the same name as the old ring).

RETURN: list L mapped into the new ring R, if L is given; else nothing
ASSUME: the active ring must allow an algebraic extension
         (e.g. it cannot be a transcendent ring extension of Q or Z/p)
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]);
  "// 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))
 splitring(x2-a,"r2",a); // change to Q(sqrt(2),sqrt(sqrt(2)))=Q(a)
                         // and return the transformed old parameter
 // the result is (a2) == (sqrt(sqrt(2)))^2
 nameof(basering);
 r2;
 kill r1; kill r2;
}