source:git/Tst/BuchDL/Ex_L8.tst@9558c5f

spielwiese
Last change on this file since 9558c5f was 9558c5f, checked in by Hans Schönemann <hannes@…>, 18 years ago
*hannes/lossen: very long examples from DL-Book git-svn-id: file:///usr/local/Singular/svn/trunk@8762 2c84dea3-7e68-4137-9b89-c4e89433aadc
• Property mode set to `100644`
File size: 11.2 KB
Line
1LIB "tst.lib";
2tst_init();
3
4
5//======================  Example 8.6 =============================
6LIB "finvar.lib";
7ring R = (0,a), (x(0..4)), dp;
8minpoly = a4+a3+a2+a+1;  // need fifth roots of unity
9matrix Si[5][5] = 0,0,0,0,1,
10                  1,0,0,0,0,
11                  0,1,0,0,0,
12                  0,0,1,0,0,
13                  0,0,0,1,0;
14matrix Ta[5][5];
15Ta[1,1] = 1;  Ta[2,2] = a;  Ta[3,3] = a2;  Ta[4,4] = a3;
16Ta[5,5] = a4;
17int aa = timer;          // time in seconds
18matrix P,S,IS = invariant_ring(Si,Ta,intvec(0,0,1));
19
20//->   Generating the entire matrix group. Whenever a new group element is
21//->   found, the corresponding ring homomorphism of the Reynolds operator
22//->   and the corresponding term of the Molien series is generated.
23
24//->   Group element 3 has been found.
25//->                        [...]
26//->   Group element 125 has been found.
27
28//->   Now we are done calculating Molien series and Reynolds operator.
29
30//->   We can start looking for primary invariants...
31
32//->   Computing primary invariants in degree 5:
33//->   We find: x(0)*x(1)*x(2)*x(3)*x(4)
34//->   We find: x(0)^3*x(2)*x(3)+x(0)*x(1)*x(3)^3+x(0)*x(2)^3*x(4)+[...]
35//->   We find: x(0)*x(1)^3*x(2)+x(1)*x(2)^3*x(3)+x(0)^3*x(1)*x(4)+[...]
36//->   Computing primary invariants in degree 10:
37//->   We find: x(0)^10+x(1)^10+x(2)^10+x(3)^10+x(4)^10
38//->   We find: -x(0)^5*x(1)^5+x(0)^5*x(2)^5-[...]
39
40//->   We found all primary invariants.
41
42//->   Polynomial telling us where to look for secondary invariants:
43//->    x(0)^30+3*x(0)^25+24*x(0)^20+44*x(0)^15+24*x(0)^10+3*x(0)^5+1
44
45//->   In degree 0 we have: 1
46
47//->   Searching in degree 5, we need to find 3 invariant(s)...
48//->   We find: x(0)^2*x(1)^2*x(3)+x(0)*x(2)^2*x(3)^2+[...]
49//->   We find: x(0)^2*x(1)*x(2)^2+x(1)^2*x(2)*x(3)^2+[...]
50//->   We find: x(0)^5+x(1)^5+x(2)^5+x(3)^5+x(4)^5 .
51
52//->   Searching in degree 10, we need to find 24 invariant(s)...
53//->                        [...]
54//->   Searching in degree 15, we need to find 44 invariant(s)...
55//->                        [...]
56//->   Searching in degree 20, we need to find 24 invariant(s)...
57//->                        [...]
58//->   Searching in degree 25, we need to find 3 invariant(s)...
59//->                        [...]
60//->   Searching in degree 30, we need to find 1 invariant(s)...
61//->   We find: x(0)^15*x(1)^10*x(3)^5+x(0)^10*x(1)^15*x(3)^5+[...]
62
63//->   We're done!
64
65ideal HMQ = invariant_basis(5,Si,Ta);
66print(HMQ);
67//->   x(0)*x(1)*x(2)*x(3)*x(4),
68//->   x(0)^2*x(1)*x(2)^2+x(1)^2*x(2)*x(3)^2+x(0)^2*x(3)^2*x(4)+[...]
69//->   x(0)^2*x(1)^2*x(3)+x(0)*x(2)^2*x(3)^2+x(1)^2*x(2)^2*x(4)+[...]
70//->   x(0)^3*x(2)*x(3)+x(0)*x(1)*x(3)^3+x(0)*x(2)^3*x(4)+[...]
71//->   x(0)*x(1)^3*x(2)+x(1)*x(2)^3*x(3)+x(0)^3*x(1)*x(4)+[...]
72//->   x(0)^5+x(1)^5+x(2)^5+x(3)^5+x(4)^5
73
74
75kill R,aa;
76//==================  Remark 8.7(new Session) =========================
77if (not(defined(invariant_ring))){ LIB "finvar.lib"; }
78ring R = 101, (x(0..4)), dp;
79matrix Si[5][5] = 0,0,0,0,1,
80                  1,0,0,0,0,
81                  0,1,0,0,0,
82                  0,0,1,0,0,
83                  0,0,0,1,0;
84number a = 36;              // primitive fifth root of unity
85matrix Ta[5][5];
86Ta[1,1] = 1;  Ta[2,2] = a;  Ta[3,3] = a^2;  Ta[4,4] = a^3;
87Ta[5,5] = a^4;
88int aa = timer;             // time in seconds
89matrix P,S,IS = invariant_ring(Si,Ta,intvec(0,0,0));
90size(S);
91//->   100
92print(S[100]);
93//->   [x(0)^15*x(1)^10*x(3)^5+x(0)^10*x(1)^15*x(3)^5+[...]
94ideal HMQ = invariant_basis(5,Si,Ta);
95
96//==================  Example 8.8 (continued Session) ===================
97ring P4 = 101, (x(0..4)), dp;
98ideal HMQ = fetch(R,HMQ);
99ideal CI = HMQ[1], HMQ[2];
100list DEG = primdecGTZ(CI);
101size(DEG);
102//->   10
103ideal I6 = DEG[6][1]; I6;
104//->   I6[1]=x(1)^2*x(3)+x(2)*x(4)^2
105//->   I6[2]=x(0)
106degree(std(I6));
107//->   // dimension (proj.)  = 2
108//->   // degree (proj.)   = 3
109ideal I10 = DEG[10][1]; I10;
110//->   I10[1]=x(2)^2
111//->   I10[2]=x(2)*x(3)^2+x(0)*x(4)^2
112//->   I10[3]=x(0)*x(2)
113//->   I10[4]=x(0)^2
114degree(std(I10));
115//->   // dimension (proj.)  = 2
116//->   // degree (proj.)   = 2
117ideal IX = intersect(DEG[3][1],DEG[4][1],DEG[8][1],DEG[9][1],
118                     DEG[10][1]);
119degree(std(IX));
120//->   // dimension (proj.)  = 2
121//->   // degree (proj.)   = 10
122resolution FIX = mres(IX,0);
123print(betti(FIX),"betti");
124//->              0     1     2     3     4
125//->   ------------------------------------
126//->       0:     1     -     -     -     -
127//->       1:     -     -     -     -     -
128//->       2:     -     -     -     -     -
129//->       3:     -     -     -     -     -
130//->       4:     -     3     -     -     -
131//->       5:     -    15    35    20     -
132//->       6:     -     -     -     -     2
133//->   ------------------------------------
134//->   total:     1    18    35    20     2
135
136module N = transpose(FIX[3]);
137homog(N);
138//->   1
139intvec deg_N = attrib(N,"isHomog");
140attrib(N,"isHomog",deg_N-3);        // set degrees
141resolution FN = mres(N,0);
142print(betti(FN),"betti");
143//->              0     1     2     3     4     5
144//->   ------------------------------------------
145//->      -3:    20    35    15     -     -     -
146//->      -2:     -     -     4     -     -     -
147//->      -1:     -     -     -     -     -     -
148//->       0:     -     -     5    15    10     2
149//->   ------------------------------------------
150//->   total:    20    35    24    15    10     2
151matrix NN = FN[2];
152matrix PRESMHM[35][19] = NN[1..35,1..19];
153PRESMHM = transpose(PRESMHM);
154resolution FMHM = mres(PRESMHM,0);
155print(betti(FMHM),"betti");
156//->              0     1     2     3
157//->   ------------------------------
158//->       0:     4     -     -     -
159//->       1:    15    35    20     -
160//->       2:     -     -     -     2
161//->   ------------------------------
162//->   total:    19    35    20     2
163
164matrix zero[1][15];
165matrix ran = random(100,1,4);
166matrix psi = transpose(concat(zero,ran));
167matrix pres = PRESMHM + module(psi);
168module dir = transpose(pres);
169resolution fdir = mres(dir,2);
170print(betti(fdir),"betti");
171//->              0     1     2
172//->   ------------------------
173//->       0:    35    15     -
174//->       1:     -     3     -
175//->       2:     -     -     -
176//->       3:     -     -     -
177//->       4:     -     -     -
178//->       5:     -     -     1
179//->   ------------------------
180//->   total:    35    18     1
181ideal IA = groebner(flatten(fdir[2]));
182int codimIA = nvars(P4) - dim(IA);
183ideal sIA = minor(jacob(IA),codimIA)+IA;
184nvars(P4) - dim(groebner(sIA));
185//->   5
186
187matrix dummy[1][3] = IA[1..3];     // the 3 quintics in IA
188ideal CI2 = dummy*random(100,3,2);
189ideal IA' = sat(CI2,IA)[1];
190resolution FIA' = mres(IA',0);
191print(betti(FIA'),"betti");
192//->              0     1     2     3     4
193//->   ------------------------------------
194//->       0:     1     -     -     -     -
195//->       1:     -     -     -     -     -
196//->       2:     -     -     -     -     -
197//->       3:     -     -     -     -     -
198//->       4:     -     3     -     -     -
199//->       5:     -     -     -     -     -
200//->       6:     -     5    15    10     2
201//->   ------------------------------------
202//->   total:     1     8    15    10     2
203int codimIA' = nvars(P4) - dim(IA');
204ideal sIA' = minor(jacob(IA'),codimIA')+IA';
205nvars(P4) - dim(groebner(sIA'));
206//->   5
207
208ideal QA = IA[1..3];
209ideal HMlines = sat(QA,IA)[1];     // result is a Groebner basis
210degree(HMlines);
211//->   // dimension (proj.)  = 1
212//->   // degree (proj.)   = 25
213
214
215kill R,aa,P4,codimIA,codimIA',deg_N;
216//================== Example 8.9 (new Session) =========================
217ring R = 2, (x(1..4)), dp;
218matrix A[4][4];
219A[1,4] = 1;  A[2,1] = 1;  A[3,2] = 1;  A[4,3] = 1;
220print(A);
221//->   0,0,0,1,
222//->   1,0,0,0,
223//->   0,1,0,0,
224//->   0,0,1,0
225if (not(defined(invariant_ring))){ LIB "finvar.lib"; }
226matrix P,S = invariant_ring(A);
227P;
228//->   P[1,1]=x(1)+x(2)+x(3)+x(4)
229//->   P[1,2]=x(1)*x(3)+x(2)*x(4)
230//->   P[1,3]=x(1)*x(2)+x(2)*x(3)+x(1)*x(4)+x(3)*x(4)
231//->   P[1,4]=x(1)*x(2)*x(3)*x(4)
232size(S);
233//->   5
234
235
236kill R;
237//================== Example 8.10 (new Session) =========================
238if (not(defined(invariant_ring))){ LIB "finvar.lib"; }
239ring R = (0,a), (x(0..3)), dp;
240minpoly = a4+a3+a2+a+1;
241matrix A[4][4];
242A[1,1] = a;  A[2,2] = a2;  A[3,3] = a3;  A[4,4] = a4;
243matrix P,S,IS = invariant_ring(A,intvec(0,0,0));
244size(P);
245//->   4
246size(S);
247//->   12
248
249proc min_generating_set (matrix P,S)
250"USAGE:  min_generating_set(P,S);   P,S matrix
251ASSUME: The entries of P,S are homogeneous and ordered by ascending
252        degrees. The first entry of S equals 1. (As satisfied by
253        the first two output matrices of invariant_ring(G).)
254RETURN: ideal
255NOTE:   The given generators for the output ideal form a minimal
256        generating set for the ring generated by the entries of
257        P,S. The generators are homogeneous and ordered by
258        descending degrees.
259"
260{
261  if (defined(flatten)==0) { LIB "matrix.lib"; }
262  ideal I1,I2 = flatten(P),flatten(S);
263  int i1,i2 = size(I1),size(I2);
264  // We order the generators by descending degrees
265  // (the first generator 1 of I2 is omitted):
266  int i,j,s = i1,i2,i1+i2-1;
267  ideal I;
268  for (int k=1; k<=s; k++)
269  {
270    if (i==0) { I[k]=I2[j]; j--; }
271    else
272    {
273      if (j==0) { I[k]=I1[i]; i--; }
274      else
275      {
276        if (deg(I1[i])>deg(I2[j])) { I[k]=I1[i]; i--; }
277        else { I[k]=I2[j]; j--; }
278      }
279    }
280  }
281  intvec deg_I = deg(I[1..s]);
282  int n = nvars(basering);
283  def BR = basering;
284
285  // Create a new ring with elimination order:
286  //---------------------------------------------------------------
287  // ****    this part uses the command ringlist which is      ****
288  // ****     only available in SINGULAR-3-0-0 or newer        ****
289  //---------------------------------------------------------------
290  list rData = ringlist(BR);
291  intvec wDp;
292  for (k=1; k<=n; k++) {
293    rData[2][k] ="x("+string(k)+ ")";
294    wDp[k]=1;
295  }
296  for (k=1; k<=s; k++) { rData[2][n+k] ="y("+string(k)+ ")"; }
297  rData[3][1] = list("dp",wDp);
298  rData[3][2] = list("wp",deg_I);
299  def R_aux = ring(rData);
300  setring R_aux;
301  //---------------------------------------------------------------
302
303  ideal J;
304  map phi = BR, x(1..n);
305  ideal I = phi(I);
306  for (k=1; k<=s; k++) { J[k] = y(k)-I[k]; }
307  option(redSB);
308  J = std(J);
309
310  // Remove all generators that are depending on some x(i) from J:
311  int s_J = size(J);
312  for (k=1; k<=s_J; k++) { if (J[k]>=x(n)) {J[k]=0;} }
313
314  // The monomial order on K[y] is chosen such that linear leading
315  // terms in J are in 1-1 correspondence to superfluous generators
316  // in I :
318  intvec to_remove;
319  i=1;
320  for (k=1; k<=s; k++)
321  {
322    if (reduce(y(k),J_1jet)==0){ to_remove[i]=k; i++; }
323  }
324  setring BR;
325  if (to_remove == 0) { return(ideal(I)); }
326  for (i=1; i<=size(to_remove); i++)
327  {
328    I[to_remove[i]] = 0;
329  }
330  I = simplify(I,2);
331  return(I);
332}
333
334ideal FSI = min_generating_set(P,S);
335size(FSI);
336//->   14
337ring Rnew = 0, (x(0..3)), dp;  // coefficient field is now Q
338ideal FSI = fetch(R,FSI);
339ideal ZERO;
340ring R1 = 0, (y(0..13)), wp(5,5,5,5,4,4,4,4,3,3,3,3,2,2);
341ideal REL = preimage(Rnew,FSI,ZERO);
342homog(REL);                    // check that REL is homogeneous
343//->   1
344size(REL);
345//->   54
346setring Rnew;
347FSI[4];
348//->   x(0)^5+x(1)^5+x(2)^5+x(3)^5
349ideal F = FSI[4];
350setring R1;
351ideal GODEAUX = preimage(Rnew,FSI,F);
352size(GODEAUX);
353//->   55
354GODEAUX[1];
355//->   y(3)
356dim(std(GODEAUX));
357//->   3
358
359tst_status(1);\$
360
Note: See TracBrowser for help on using the repository browser.