source:git/Tst/Buch/Proc_4_6.tst@891438c

spielwiese
Last change on this file since 891438c was 891438c, checked in by Gerhard Pfister <pfister@…>, 22 years ago
*GP: initial release git-svn-id: file:///usr/local/Singular/svn/trunk@5579 2c84dea3-7e68-4137-9b89-c4e89433aadc
• Property mode set to `100644`
File size: 9.0 KB
Line
1LIB "tst.lib";
2tst_init();
3
4LIB"ring.lib";
5proc primaryTest(ideal i, poly p)
6"USAGE:   primaryTest(i,p); i standard basis with respect to
7          lp, p irreducible polynomial in K[var(n)],
8          p^a=i[1] for some a;
9RETURN:   an ideal, the radical of i if i is primary and in
10          general position with respect to lp,
11          the zero ideal else.
12
13{
14   int m,e;
15   int n=nvars(basering);
16   poly t;
17   ideal prm=p;
18
19   for(m=2;m<=size(i);m++)
20   {
22     {
23       n--;
24//----------------i[m] has a power of var(n) as leading term
25       attrib(prm,"isSB",1);
26//--- ?? i[m]=(c*var(n)+h)^e modulo prm for h
27//     in K[var(n+1),...], c in K ??
30       /var(n)^(e-1);
32//---if not (0) is returned, else c*var(n)+h is added to prm
33       if (reduce(i[m]-t^e,prm,1) !=0)
34       {
35         return(ideal(0));
36       }
37       prm = prm,cleardenom(simplify(t,1));
38     }
39   }
40   return(prm);
41}
42
43   ring s=(0,x),(d,e,f,g),lp;
44   ideal i=g^5,(x*f-g)^3,5*e-g^2,x*d^3;
45   primaryTest(i,g);
46   kill s;
47
48proc zeroDecomp(ideal i)
49"USAGE:  zeroDecomp(i); i zerodimensional ideal
50RETURN:  list l of lists of two ideals such that the
51         intersection(l[j][1], j=1..)=i, the l[i][1] are
52         primary and the l[i][2] their radicals
53NOTE:    algorithm of Gianni/Trager/Zacharias
54
55{
56   def  BAS = basering;
57//----the ordering is changed to the lexicographical one
58   changeord("R","lp");
59   ideal i=fetch(BAS,i);
60   int n=nvars(R);
61   int k;
62   list result,rest;
63   ideal primary,prim;
64   option(redSB);
65
66//------the random coordinatechange and its inverse
67   ideal m=maxideal(1);
68   m[n]=0;
69   poly p=(random(100,1,n)*transpose(m))[1,1]+var(n);
70   m[n]=p;
71   map phi=R,m;
72   m[n]=2*var(n)-p;
73   map invphi=R,m;
74   ideal j=groebner(phi(i));
75//-------------factorization of the first element in i
76   list fac=factorize(j[1],2);
77//-------------computation of the primaries and primes
78   for(k=1;k<=size(fac[1]);k++)
79   {
80     p=fac[1][k]^fac[2][k];
81     primary=groebner(j+p);
82     prim=primaryTest(primary,fac[1][k]);
83//---test whether all ideals were primary and in general
84//   position
85     if(prim==0)
86     {
87       rest[size(rest)+1]=i+invphi(p);
88     }
89     else
90     {
91       result[size(result)+1]=
92         list(std(i+invphi(p)),std(invphi(prim)));
93     }
94   }
95//-------treat the bad cases collected in the rest again
96   for(k=1;k<=size(rest);k++)
97   {
98     result=result+zeroDecomp(rest[k]);
99   }
100   option(noredSB);
101   setring BAS;
102   list result=imap(R,result);
103   kill R;
104   return(result);
105}
106
107   ring  r = 32003,(x,y,z),dp;
108   poly  p = z2+1;
109   poly  q = z4+2;
110   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
111   list pr= zeroDecomp(i);
112   pr;
113   kill r;
114
115proc prepareQuotientring(ideal i)
116"USAGE:   prepareQuotientring(i); i standard basis
117RETURN:   a list l of two strings:
118          l[1] to define K[x\u,u ], u a maximal independent
119          set for i
120          l[2] to define K(u)[x\u ], u a maximal independent
121          set for i
122          both rings with lexicographical ordering
123
124{
125  string va,pa;
126//v describes the independent set u: var(j) is in
127//u iff v[j]!=0
128  intvec v=indepSet(i);
129  int k;
130
131  for(k=1;k<=size(v);k++)
132  {
133    if(v[k]!=0)
134    {
135      pa=pa+"var("+string(k)+"),";
136    }
137    else
138    {
139      va=va+"var("+string(k)+"),";
140    }
141  }
142
143  pa=pa[1..size(pa)-1];
144  va=va[1..size(va)-1];
145
146  string newring
147  ="ring nring=("+charstr(basering)+"),("+va+","+pa+"),lp;";
148  string quotring
149  ="ring quring=("+charstr(basering)+","+pa+"),("+va+"),lp;";
150  return(newring,quotring);
151}
152
153   ring s=(0,x),(a,b,c,d,e,f,g),dp;
154   ideal i=x*b*c,d^2,f-g;
155   i=std(i);
156   def Q=basering;
157   list l= prepareQuotientring(i);
158   l;
159   execute (l[1]);
160   basering;
161   execute (l[2]);
162   basering;
163   setring Q;
164
165proc prepareSat(ideal i)
166{
167   int k;
169   for(k=2;k<=size(i);k++)
170   {
172   }
173   return(p);
174}
175
176LIB"elim.lib";
177proc decomp(ideal i)
178"USAGE:  decomp(i); i ideal
179RETURN:  list l of lists of two ideals such that the
180         intersection(l[j][1], j=1..)=i, the l[i][1] are
181         primary and the l[i][2] their radicals
182NOTE:    algorithm of Gianni/Trager/Zacharias
183
184{
185   if(i==0)
186   {
187     return(list(i,i));
188   }
189   def  BAS = basering;
190   ideal j;
191   int n=nvars(BAS);
192   int k;
193
194   ideal SBi=std(i);
195   int d=dim(SBi);
196//---the trivial case and the zero-dimensional case
197   if ((d==0)||(d==-1))
198   {
199      return(zeroDecomp(i));
200   }
201//---prepare the quotientring with respect to a maximal
202//   independent set
203   list quotring=prepareQuotientring(SBi);
204   execute (quotring[1]);
205//---we use this ring to compute a standardbasis of i*quring
206//   which is in i
207   ideal i=std(imap(BAS,i));
208//---pass to the quotientring with respect to a maximal
209//   independent set
210   execute(quotring[2]);
211   ideal i=imap(nring,i);
212   kill nring;
213//---computation of the zerodimensional decomposition
214   list ra=zeroDecomp(i);
215//---preparation for saturation
216   list p;
217   for(k=1;k<=size(ra);k++)
218   {
219     p[k]=list(prepareSat(ra[k][1]),prepareSat(ra[k][2]));
220   }
221   poly q=prepareSat(i);
222//---back to the original ring
223   setring BAS;
224   list p=imap(quring,p);
225   list ra=imap(quring,ra);
226
227   poly q=imap(quring,q);
228
229   kill quring;
230//---compute the intersection of ra with BAS
231   for(k=1;k<=size(ra);k++)
232   {
233     ra[k]=list(sat(ra[k][1],p[k][1])[1],
234                sat(ra[k][2],p[k][2])[1]);
235   }
236   q=q^sat(i,q)[2];
237
238//---i=intersection((i:q),(i,q)) and ra is the primary
239//   decomposition of i:q
240   if(deg(q)>0)
241   {
242     ra=ra+decomp(i+q);
243   }
244   return(ra);
245}
246
247   ring  r = 0,(x,y,z),dp;
248   ideal i = intersect(ideal(x,y,z)^3,ideal(x-y-z)^2,
249             ideal(x-y,x-z)^2);
250   list pr= decomp(i);
251   pr;
252   kill r;
253
254proc equidimensional(ideal i)
255"USAGE:  equidimensional(i); i ideal
256RETURN:  list l of two ideals such that intersetion(l[1],
257         l[2])=i
258         l[1] is equidimensional and dim(l[1])>dim(l[2])
259"
260{
261   def  BAS = basering;
262
263   ideal SBi=std(i);
264   int d=dim(SBi);
265   int n=nvars(BAS);
266   int k;
267   list result;
268
269//-----the trivial cases
270   if ((d==-1)||(n==d)||(n==1)||(d==0))
271   {
272      result=i,ideal(1);
273      return(result);
274   }
275//-----prepare the quotientring with respect to a maximal
276//     independent set
277   list quotring=prepareQuotientring(SBi);
278   execute (quotring[1]);
279//-----we use this ring to compute a standardbasis of
280//     i*quring which is in i
281   ideal eq=std(imap(BAS,i));
282//-----pass to the quotientring with respect to a maximal
283//     independent set
284   execute (quotring[2]);
285   ideal eq=imap(nring,eq);
286   kill nring;
287//-----preparation for saturation
288   poly p=prepareSat(eq);
289//-----back to the original ring
290   setring BAS;
291   poly p=imap(quring,p);
292   ideal eq=imap(quring,eq);
293   kill quring;
294//-----compute the intersection of eq with BAS
295   eq=sat(eq,p)[1];
296
297   SBi=std(quotient(i,eq));
298
299   if(d>dim(SBi))
300   {
301     result=eq,SBi;
302     return(result);
303   }
304   result=equidimensional(i);
305   result=intersect(result[1],eq),result[2];
306   return(result);
307}
308
309   ring  r = 0,(x,y,z),dp;
310   ideal i = intersect(ideal(x,y,z)^3,ideal(x-y-z)^2,
311             ideal(x-y,x-z)^2);
312   list pr= equidimensional(i);
313   pr;
314   dim(std(pr[1]));
315   dim(std(pr[2]));
316   option(redSB);
317   std(i);
318   std(intersect(pr[1],pr[2]));
319   kill r;
320
321proc squarefree(poly f, int i)
322{
323   poly h=gcd(f,diff(f,var(i)));
324   poly g=lift(h,f)[1][1];    //  f/h
325   return(g);
326}
327
330RETURN:  ideal = the radical of i
331NOTE:    algorithm of Krick/Logar
332
333{
334   def  BAS = basering;
335   ideal j;
336   int n=nvars(BAS);
337   int k;
338
339   option(redSB);
340   ideal SBi=std(i);
341   option(noredSB);
342   int d=dim(SBi);
343
344//-----the trivial cases
345   if ((d==-1)||(n==d)||(n==1))
346   {
347      return(ideal(squarefree(SBi[1],1)));
348   }
349//-----the zero-dimensional case
350   if (d==0)
351   {
352      j=finduni(SBi);
353      for(k=1;k<=size(j);k++)
354      {
355         i=i,squarefree(cleardenom(j[k]),k);
356      }
357      return(std(i));
358   }
359//-----prepare the quotientring with respect to a maximal
360//     independent set
361   list quotring=prepareQuotientring(SBi);
362   execute (quotring[1]);
363//-----we use this ring to compute a standardbasis of
364//     i*quring which is in i
365   ideal i=std(imap(BAS,i));
366//-----pass to the quotientring with respect to a maximal
367//     independent set
368   execute( quotring[2]);
369   ideal i=imap(nring,i);
370   kill nring;
373//-----preparation for saturation
374   poly p=prepareSat(ra);
375   poly q=prepareSat(i);
376//-----back to the original ring
377   setring BAS;
378   poly p=imap(quring,p);
379   poly q=imap(quring,q);
380   ideal ra=imap(quring,ra);
381   kill quring;
382//-----compute the intersection of ra with BAS
383   ra=sat(ra,p)[1];