source: git/Singular/LIB/primitiv.lib @ 3686937

spielwiese
Last change on this file since 3686937 was 3686937, checked in by Oleksandr Motsak <motsak@…>, 11 years ago
Added '$Id$' as a comment to all libs (LIB/*.lib)
  • Property mode set to 100644
File size: 14.0 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="version primitiv.lib 4.0.0.0 Jun_2013 "; // $Id$
3category="Commutative Algebra";
4info="
5LIBRARY:    primitiv.lib    Computing a Primitive Element
6AUTHOR:     Martin Lamm,    email: lamm@mathematik.uni-kl.de
7
8PROCEDURES:
9 primitive(ideal i);   find minimal polynomial for a primitive element
10 primitive_extra(i);   find primitive element for two generators
11 splitring(f,R[,L]);   define ring extension with name R and switch to it
12";
13
14LIB "random.lib";
15///////////////////////////////////////////////////////////////////////////////
16
17proc primitive(ideal i)
18"USAGE:   primitive(i); i ideal
19ASSUME:   i is given by generators m[1],...,m[n] such that for j=1,...,n @*
20          -  m[j] is a polynomial in k[x(1),...,x(j)] @*
21          -  m[j](a[1],...,a[j-1],x(j)) is the minimal polynomial for a[j] over
22             k(a[1],...,a[j-1]) @*
23          (k the ground field of the current basering and x(1),...,x(n)
24          the ring variables).
25RETURN:   ideal j in k[x(n)] with
26          - j[1] a minimal polynomial for a primitive element b of
27                 k(a[1],...,a[n]) over k,
28          - j[2],...,j[n+1] polynomials in k[x(n)] such that j[i+1](b)=a[i]
29                 for i=1,...,n.
30NOTE:     the number of variables in the basering has to be exactly n,
31          the number of given generators (i.e., minimal polynomials).@*
32          If the ground field k has only a few elements it may happen that no
33          linear combination of a[1],...,a[n] is a primitive element. In this
34          case @code{primitive(i)} returns the zero ideal, and one should use
35          @code{primitive_extra(i)} instead.
36SEE ALSO: primitive_extra
37KEYWORDS: primitive element
38EXAMPLE:  example primitive;  shows an example
39"
40{
41 def altring=basering;
42 execute("ring deglexring=("+charstr(altring)+"),("+varstr(altring)+"),dp;");
43 ideal j;
44 execute("ring lexring=("+charstr(altring)+"),("+varstr(altring)+"),lp;");
45 ideal i=fetch(altring,i);
46
47 int k,schlecht,Fehlversuche,maxtry;
48 int nva = nvars(basering);
49 int p=char(basering);
50 if (p==0) {
51   p=100000;
52   if (nva<3) { maxtry= 100000000; }
53   else       { maxtry=2147483647; }
54 }
55 else {
56   if ((nva<4) || (p<60)) {
57     maxtry=p^(nva-1); }
58   else {
59     maxtry=2147483647;          // int overflow(^)  vermeiden
60   }
61 }
62 ideal jmap,j;
63 map phi;
64 option(redSB);
65
66 //-------- Mache so lange Random-Koord.wechsel, bis letztes Polynom -------------
67 //--------------- das Minpoly eines primitiven Elements ist : ----------------
68 for (Fehlversuche=0; Fehlversuche<maxtry; Fehlversuche++) {
69   schlecht=0;
70   if ((p<60) && (nva==2)) {  // systematische Suche statt random
71      jmap=ideal(var(1),var(2)+Fehlversuche*var(1));
72   }
73   else {
74    if (Fehlversuche==0) { jmap=maxideal(1);}
75    else {
76      if (Fehlversuche<5) { jmap=randomLast(10);}
77      else {
78       if (Fehlversuche<20) { jmap=randomLast(100);}
79       else                 { jmap=randomLast(100000000);}
80    }}                        // groessere Werte machen keinen Sinn
81   }
82   phi=lexring,jmap;
83   j=phi(i);
84   setring deglexring;
85 //--------------- Berechne reduzierte Standardbasis mit fglm: ----------------
86   j=std(fetch(lexring,j));
87   setring lexring;
88   j=fglm(deglexring,j);
89 //-- teste, ob SB n Elemente enthaelt (falls ja, ob lead(Fi)=xi i=1... n-1): -
90   if (size(j)==nva) {
91     for (k=1; k<nva; k++) {
92       j[k+1]=j[k+1]/leadcoef(j[k+1]);        // normiere die Erzeuger
93       if (lead(j[k+1]) != var(nva-k)) { schlecht=1;}
94     }
95     if (schlecht==0) {
96 //--- Random-Koord.wechsel war gut: Berechne das zurueckzugebende Ideal: -----
97       ideal erg;
98       for (k=1; k<nva; k++) { erg[k]=var(k)-j[nva-k+1]; }
99                               // =g_k(x_n) mit a_k=g_k(a_n)
100       erg[nva]=var(nva);
101       map chi=lexring,erg;
102       ideal extra=maxideal(1);extra=phi(extra);
103                               // sonst: "argument of a map must have a name"
104       erg=j[1],chi(extra);    // j[1] = Minimalpolynom
105       setring altring;
106       return(fetch(lexring,erg));
107     }
108   }
109   dbprint("The random coordinate change was bad!");
110 }
111 if (voice==2) {
112   "// ** Warning: No primitive element could be found.";
113   "//    If the given ideal really describes the minimal polynomials of";
114   "//    a series of algebraic elements (cf. `help primitive;') then";
115   "//    try `primitive_extra'.";
116 }
117 setring altring;
118 return(ideal(0));
119}
120example
121{ "EXAMPLE:"; echo = 2;
122 ring exring=0,(x,y),dp;
123 ideal i=x2+1,y2-x;                  // compute Q(i,i^(1/2))=:L
124 ideal j=primitive(i);
125 j[1];                               // L=Q(a) with a=(-1)^(1/4)
126 j[2];                               // i=a^2
127 j[3];                               // i^(1/2)=a
128 // the 2nd element was already primitive!
129 j=primitive(ideal(x2-2,y2-3));      // compute Q(sqrt(2),sqrt(3))
130 j[1];
131 j[2];
132 j[3];
133 // no element was primitive -- the calculation of primitive elements
134 // is based on a random choice.
135}
136///////////////////////////////////////////////////////////////////////////////
137
138proc primitive_extra(ideal i)
139"USAGE:   primitive_extra(i); i ideal
140ASSUME:  The ground field of the basering is k=Q or k=Z/pZ and the ideal
141         i is given by 2 generators f,g with the following properties:
142@format
143   f is the minimal polynomial of a in k[x],
144   g is a polynomial in k[x,y] s.th. g(a,y) is the minpoly of b in k(a)[y].
145@end format
146          Here, x is the name of the first ring variable, y the name of the
147          second.
148RETURN:  ideal j in k[y] such that
149@format
150   j[1] is the minimal polynomial for a primitive element c of k(a,b) over k,
151   j[2] is a polynomial s.th. j[2](c)=a.
152@end format
153NOTE:    While @code{primitive(i)} may fail for finite fields,
154         @code{primitive_extra(i)} tries all elements of k(a,b) and, hence,
155         always finds a primitive element. @*
156         In order to do this (try all elements), field extensions like Z/pZ(a)
157         are not allowed for the ground field k. @*
158         @code{primitive_extra(i)} assumes that the second generator, g, is
159         monic as polynomial in (k[x])[y].
160EXAMPLE: example primitive_extra;  shows an example
161"
162{
163 def altring=basering;
164 int grad1=deg(i[1]);
165 int grad2=deg(jet(i[2],0,intvec(1,0)));
166 if (grad2==0) { ERROR("i[2] is not monic"); }
167 int countx,countz;
168  if (size(variables(i[1]))!=1) { ERROR("i[1] must be poly in x"); }
169  if (size(variables(i[2]))>2) { ERROR("i[2] must be poly in x,a"); }
170  //if (variables(i[2])[2]!=a) { ERROR("i[2] must be poly in x,a"); }
171 ring deglexring=char(altring),(x,y,z),dp;
172 map transfer=altring,x,z;
173 ideal i=transfer(i);
174 if (size(i)!=2)
175 {
176   ERROR("either wrong number of given minimal polynomials"+newline+
177   "or wrong choice of ring variables (must use the first two)");
178 }
179 matrix mat;
180 ring lexring=char(altring),(x,y),lp;
181 ideal j;
182 ring deglex2ring=char(altring),(x,y),dp;
183 ideal j;
184 setring deglexring;
185 ideal j;
186 option(redSB);
187 poly g=z;
188 int found=0;
189
190 //---------------- Schleife zum Finden des primitiven Elements ---------------
191 //--- Schleife ist so angordnet, dass g in Charakteristik 0 linear bleibt ----
192 while (found==0)
193 {
194   j=eliminate(i+ideal(g-y),z);
195   setring deglex2ring;
196   j=std(imap(deglexring,j));
197   setring lexring;
198   j=fglm(deglex2ring,j);
199   if (size(j)==2)
200   {
201     if (deg(j[1])==grad1*grad2)
202     {
203       j[2]=j[2]/leadcoef(j[2]);    // Normierung
204       if (lead(j[2])==x)
205       {         // Alles ok
206          found=1;
207       }
208     }
209   }
210   setring deglexring;
211   if (found==0)
212   {
213 //------------------ waehle ein neues Polynom g ------------------------------
214     dbprint("Still searching for primitive element...");
215     countx=0;
216     countz=0;
217     while (found==0)
218     {
219      countx++;
220      if (countx>=grad1)
221      {
222        countx=0;
223        countz++;
224        if (countz>=grad2)
225        { ERROR("No primitive element found!! This should NEVER happen!"); }
226      }
227      g = g +x^countx *z^countz;
228      mat=coeffs(g,z);
229      if (size(mat)>countz)
230      {
231        mat=coeffs(mat[countz+1,1],x);
232        if (size(mat)>countx)
233        {
234          if (mat[countx+1,1] != 0)
235          {
236            found=1;         // d.h. hier: neues g gefunden
237      }}}
238     }
239     found=0;
240   }
241 }
242 //------------------- primitives Element gefunden; Rueckgabe -----------------
243 setring lexring;
244 j[2]=x-j[2];
245 setring altring;
246 map transfer=lexring,var(1),var(2);
247 return(transfer(j));
248}
249example
250{ "EXAMPLE:"; echo = 2;
251 ring exring=3,(x,y),dp;
252 ideal i=x2+1,y3+y2-1;
253 primitive_extra(i);
254 ring extension=(3,y),x,dp;
255 minpoly=y6-y5+y4-y3-y-1;
256 number a=y5+y4+y2+y+1;
257 a^2;
258 factorize(x2+1);
259 factorize(x3+x2-1);
260}
261///////////////////////////////////////////////////////////////////////////////
262
263proc splitring(poly f,list #)
264"USAGE:   splitring(f[,L]); f poly, L list of polys and/or ideals
265         (optional)
266ASSUME:  f is univariate and irreducible over the active ring. @*
267         The active ring must allow an algebraic extension (e.g., it cannot
268         be a transcendent ring extension of Q or Z/p).
269RETURN:  ring; @*
270           if called with a nonempty second parameter L, then in the output
271           ring there is defined a list erg ( =L mapped to the new ring);
272           if the minpoly of the active ring is non-zero, then the image of
273           the primitive root of f in the output ring is appended as last
274           entry of the list erg.
275NOTE:    If the old ring has no parameter, the name @code{a} is chosen for the
276         parameter of R (if @code{a} is no ring variable; if it is, @code{b} is
277         chosen, etc.; if @code{a,b,c,o} are ring variables,
278         @code{splitring(f[,L])} produces an error message), otherwise the
279         name of the parameter is kept and only the minimal polynomial is
280         changed. @*
281         The names of the ring variables and the orderings are not affected. @*
282KEYWORDS: algebraic field extension; extension of rings
283EXAMPLE: example splitring;  shows an example
284"
285{
286 //----------------- split ist bereits eine proc in 'inout.lib' ! -------------
287 if (size(#)>=1) {
288    list L=#;
289    int L_groesse=size(L);
290 }
291 else { int L_groesse=-1; }
292 //-------------- ermittle das Minimalpolynom des aktuellen Rings: ------------
293 string minp=string(minpoly);
294
295 def altring=basering;
296 string charakt=string(char(altring));
297 string varnames=varstr(altring);
298 string algname;
299 int i;
300 int anzvar=size(maxideal(1));
301 //--------------- Fall 1: Bisheriger Ring hatte kein Minimalpolynom ----------
302 if (minp=="0") { // only possible without parameters (by assumption)
303  if (find(varnames,"a")==0)        { algname="a";}
304  else { if (find(varnames,"b")==0) { algname="b";}
305         else { if (find(varnames,"c")==0)
306                                    { algname="c";}
307         else { if (find(varnames,"o")==0)
308                                    { algname="o";}
309         else {
310           "** Sorry -- could not find a free name for the primitive element.";
311           "** Try e.g. a ring without 'a' or 'b' as variable.";
312           return();
313         }}
314       }
315  }
316  //-- erzeuge einen String, der das Minimalpolynom des neuen Rings enthaelt: -
317  execute("ring splt1="+charakt+","+algname+",dp;");
318  ideal abbnach=var(1);
319  for (i=1; i<anzvar; i++) { abbnach=abbnach,var(1); }
320  map nach_splt1=altring,abbnach;
321  execute("poly mipol="+string(nach_splt1(f))+";");
322  string Rminp=string(mipol);
323
324  //--------------------- definiere den neuen Ring: ---------------------------
325  execute("ring neuring = ("+charakt+","+algname+"),("+varnames+"),("
326           +ordstr(altring)+");");
327  execute("minpoly="+Rminp+";");
328
329  //---------------------- Berechne die zurueckzugebende Liste: ---------------
330  if (L_groesse>0) {
331   list erg;
332   map take=altring,maxideal(1);
333   erg=take(L);
334  }
335 }
336 else {
337
338  //------------- Fall 2: Bisheriger Ring hatte ein Minimalpolynom: -----------
339  algname=parstr(altring);           // Name des algebraischen Elements
340  if (npars(altring)>1) {"only one Parameter is allowed!!"; return(altring);}
341
342  //---------------- Minimalpolynom in ein Polynom umwandeln: -----------------
343  execute("ring splt2="+charakt+","+algname+",dp;");
344  execute("poly mipol="+minp+";");
345  // f ist Polynom in algname und einer weiteren Variablen -> mache f bivariat:
346  execute("ring splt3="+charakt+",("+algname+","+varnames+"),dp;");
347  poly f=imap(altring,f);
348
349  //-------------- Vorbereitung des Aufrufes von primitive: -------------------
350  execute("ring splt1="+charakt+",(x,y),dp;");
351  ideal abbnach=x;
352  for (i=1; i<=anzvar; i++) { abbnach=abbnach,y; }
353  map nach_splt1_3=splt3,abbnach;
354  map nach_splt1_2=splt2,x;
355  ideal maxid=nach_splt1_2(mipol),nach_splt1_3(f);
356  ideal primit=primitive(maxid);
357  if (size(primit)==0) {             // Suche mit 1. Proc erfolglos
358    primit=primitive_extra(maxid);
359  }
360  //-- erzeuge einen String, der das Minimalpolynom des neuen Rings enthaelt: -
361  setring splt2;
362  map nach_splt2=splt1,0,var(1);     // x->0, y->a
363  minp=string(nach_splt2(primit)[1]);
364  if (printlevel > -1) { "// new minimal polynomial:",minp; }
365  //--------------------- definiere den neuen Ring: ---------------------------
366  execute("ring neuring = ("+charakt+","+algname+"),("+varnames+"),("
367          +ordstr(altring)+");");
368  execute("minpoly="+minp+";");
369
370  if (L_groesse>0) {
371    //---------------------- Berechne die zurueckzugebende Liste: -------------
372    list erg;
373    setring splt3;
374    list zwi=imap(altring,L);
375    map nach_splt3_1=splt1,0,var(1);  // x->0, y->a
376    //----- rechne das primitive Element von altring in das von neuring um: ---
377    ideal convid=maxideal(1);
378    convid[1]=nach_splt3_1(primit)[2];
379    poly new_b=nach_splt3_1(primit)[3];
380    map convert=splt3,convid;
381    zwi=convert(zwi);
382    setring neuring;
383    erg=imap(splt3,zwi);
384    erg[size(erg)+1]=imap(splt3,new_b);
385  }
386 }
387 if (defined(erg)){export erg;}
388 return(neuring);
389}
390example
391{ "EXAMPLE:"; echo = 2;
392 ring r=0,(x,y),dp;
393 def r1=splitring(x2-2);
394 setring r1; basering;    // change to Q(sqrt(2))
395 // change to Q(sqrt(2),sqrt(sqrt(2)))=Q(a) and return the transformed
396 // old parameter:
397 def r2=splitring(x2-a,a);
398 setring r2; basering; erg;
399 // the result is (a)^2 = (sqrt(sqrt(2)))^2
400 kill r1; kill r2;
401}
402///////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.