source:git/Singular/LIB/invar.lib@e801fe

spielwiese
Last change on this file since e801fe was e801fe, checked in by Gerhard Pfister <pfister@…>, 26 years ago
* updated invar.lib normal.lib, merged prim_dec.lib with primdec.lib git-svn-id: file:///usr/local/Singular/svn/trunk@1145 2c84dea3-7e68-4137-9b89-c4e89433aadc
• Property mode set to `100644`
File size: 14.9 KB
Line
1// \$Id: invar.lib,v 1.4 1998-02-16 12:07:01 pfister Exp \$
2///////////////////////////////////////////////////////
3// invar.lib
4// algorithm for computing the ring of invariants under
5// the action of the additive group (C,+)
6// written by Gerhard Pfister
7//////////////////////////////////////////////////////
8
9LIBRARY: invar.lib PROCEDURES FOR COMPUTING INVARIANTS OF (C,+)-ACTIONS
10
11  invariantRing(matrix m,poly p,poly q,int choose)
12  // ring of invariants of the action of the additive group
13  // defined by the vectorfield corresponding to the matrix m
14  // (m[i,1] are the coefficients of d/dx(i))
15  // the polys p and q are assumed to be variables x(i) and x(j)
16  // such that m[j,1]=0 and m[i,1]=x(j)
17  // if choose=0 the computation stops if generators of the ring
18  // of invariants are computed (to be used only if you know that
19  // the ring of invariants is finitey generated)
20  // if choose<>0 it computes invariants up to degree choose
21
22  actionIsProper(matrix m)
23  // returns 1 if the action of the additive group defined by the
24  // matrix m as above i proper and 0 if not.
25
26///////////////////////////////////////////////////////////////////////////////
27
28LIB "inout.lib";
29LIB "general.lib";
30
31///////////////////////////////////////////////////////////////////////////////
32
33
34proc sortier(ideal id)
35{
36  if(size(id)==0)
37  {
38     return(id);
39  }
40  intvec i=sortvec(id);
41  int j;
42  ideal m;
43  for (j=1;j<=size(i);j++)
44  {
45    m[j] = id[i[j]];
46  }
47  return(m);
48}
49example
50{ "EXAMPLE:"; echo = 2;
51   ring q=0,(x,y,z,u,v,w),dp;
52   ideal i=w,x,z,y,v;
53   ideal j=sortier(i);
54   j;
55}
56
57
58///////////////////////////////////////////////////////////////////////////////
59
60
61proc der (matrix m, poly f)
62USAGE:   der(m,f);  m matrix, f poly
63RETURN:  poly= application of the vectorfield m befined by the matrix m
64         (m[i,1] are the coefficients of d/dx(i)) to f
65NOTE:
66EXAMPLE: example der; shows an example
67{
68  matrix mh=matrix(jacob(f))*m;
69  return(mh[1,1]);
70}
71example
72{ "EXAMPLE:"; echo = 2;
73   ring q=0,(x,y,z,u,v,w),dp;
74   poly f=2xz-y2;
75   matrix m[6][1];
76   m[2,1]=x;
77   m[3,1]=y;
78   m[5,1]=u;
79   m[6,1]=v;
80   der(m,f);
81}
82
83///////////////////////////////////////////////////////////////////////////////
84
85
86proc actionIsProper(matrix m)
87USAGE:   actionIsProper(m); m matrix
88RETURN:  int= 1 if is proper, 0 else
89NOTE:
90EXAMPLE: example actionIsProper; shows an example
91{
92  int i;
93  ideal id=maxideal(1);
94  def bsr=basering;
95
96  //changes the basering bsr to bsr[@t]
97  execute "ring s="+charstr(basering)+",("+varstr(basering)+",@t),dp;";
98  poly inv,delta,tee,j;
99  ideal id=imap(bsr,id);
100  matrix @m[size(id)+1][1];
101  @m=imap(bsr,m),0;
102
103  //computes the exp(@t*m)(var(i)) for all i
104  for(i=1;i<=nvars(basering)-1;i++)
105  {
106     inv=var(i);
107     delta=der(@m,inv);
108     j=1;
109     tee=@t;
110     while(delta!=0)
111     {
112        inv=inv+1/j*delta*tee;
113        j=j*(j+1);
114        tee=tee*@t;
115        delta=der(@m,delta);
116     }
117     id=id+ideal(inv);
118  }
119  i=inSubring(@t,id)[1];
120  setring(bsr);
121  return(i);
122}
123example
124{ "EXAMPLE:"; echo = 2;
125
126  ring rf=0,(x(1..7)),dp;
127  matrix m[7][1];
128  m[4,1]=x(1)^3;
129  m[5,1]=x(2)^3;
130  m[6,1]=x(3)^3;
131  m[7,1]=(x(1)*x(2)*x(3))^2;
132  actionIsProper(m);
133
134  ring rd=0,(x(1..5)),dp;
135  matrix m[5][1];
136  m[3,1]=x(1);
137  m[4,1]=x(2);
138  m[5,1]=1+x(1)*x(4)^2
139  actionIsProper(m);
140}
141///////////////////////////////////////////////////////////////////////////////
142
143
144proc reduction(poly p, ideal dom, list #)
145USAGE:   reduction(p,dom,q); p poly, dom ideal, q (optional) monomial
146RETURN:  poly= (p-H(f1,...,fr))/q^a, if Lt(p)=H(Lt(f1),...,Lt(fr)) for
147               some polynomial H in r variables over the base field,
148               a maximal such that q^a devides p-H(f1,...,fr),
149               dom =(f1,...,fr)
150NOTE:
151EXAMPLE: example reduction; shows an example
152{
153  int i;
154  int z=size(dom);
155  def bsr=basering;
156
157  //arranges the monomial v for elimination
158  poly v=var(1);
159  for (i=2;i<=nvars(basering);i=i+1)
160  {
161    v=v*var(i);
162  }
163
164  //changes the basering bsr to bsr[@(0),...,@(z)]
165  execute "ring s="+charstr(basering)+",("+varstr(basering)+",@(0..z)),dp;";
166
167  //costructes the ideal dom=(p-@(0),dom[1]-@(1),...,dom[z]-@(z))
168  ideal dom=imap(bsr,dom);
169  for (i=1;i<=z;i++)
170  {
172  }
174
175  //eliminates the variables of the basering bsr
176  //i.e. computes dom intersected with K[@(0),...,@(z)]
177  ideal kern=eliminate(dom,imap(bsr,v));
178
179  // test wether @(0)-h(@(1),...,@(z)) is in ker for some poly h
180  poly h;
181  z=size(kern);
182  for (i=1;i<=z;i++)
183  {
184     h=kern[i]/@(0);
185     if (deg(h)==0)
186     {
187        h=(1/h)*kern[i];
188        // defines the map psi : s ---> bsr defined by @(i) ---> p,dom[i]
189        setring bsr;
190        map psi=s,maxideal(1),p,dom;
191        poly re=psi(h);
192
193        // devides by the maximal power of #[1]
194        if (size(#)>0)
195        {
196           while ((re!=0) && (re!=#[1]) &&(subst(re,#[1],0)==0))
197           {
198             re=re/#[1];
199           }
200        }
201
202        return(re);
203     }
204  }
205  setring bsr;
206  return(p);
207}
208example
209{ "EXAMPLE:"; echo = 2;
210   ring q=0,(x,y,z,u,v,w),dp;
211   poly p=x2yz-x2v;
212   ideal dom =x-w,u2w+1,yz-v;
213   reduction(p,dom);
214   reduction(p,dom,w);
215}
216
217///////////////////////////////////////////////////////////////////////////////
218
219proc completeReduction(poly p, ideal dom, list #)
220USAGE:   completeReduction(p,dom,q); p poly, dom ideal,
221                                     q (optional) monomial
222RETURN:  poly= the polynomial p reduced with dom via the procedure
223               reduction as long as possible
224NOTE:
225EXAMPLE: example completeReduction; shows an example
226{
227  poly p1=p;
228  poly p2=reduction(p,dom,#);
229  while (p1!=p2)
230  {
231    p1=p2;
232    p2=reduction(p1,dom,#);
233  }
234  return(p2);
235}
236example
237{ "EXAMPLE:"; echo = 2;
238   ring q=0,(x,y,z,u,v,w),dp;
239   poly p=x2yz-x2v;
240   ideal dom =x-w,u2w+1,yz-v;
241   completeReduction(p,dom);
242   completeReduction(p,dom,w);
243}
244///////////////////////////////////////////////////////////////////////////////
245
246proc inSubring(poly p, ideal dom)
247USAGE:   inSubring(p,dom); p poly, dom ideal
248RETURN:  list= 1,string(@(0)-h(@(1),...,@(size(dom)))) :if p = h(dom[1],...,dom[size(dom)])
249              0,string(h(@(0),...,@(size(dom)))) :if there is only a nonlinear relation
250              h(p,dom[1],...,dom[size(dom)])=0.
251NOTE:
252EXAMPLE: example inSubring; shows an example
253{
254  int z=size(dom);
255  int i;
256  def gnir=basering;
257  list l;
258  poly mile=var(1);
259  for (i=2;i<=nvars(basering);i++)
260  {
261    mile=mile*var(i);
262  }
263  string eli=string(mile);
264  // the intersection of ideal nett=(p-@(0),dom[1]-@(1),...)
265  // with the ring k[@(0),...,@(n)] is computed, the result is ker
266  execute "ring r1="+charstr(basering)+",("+varstr(basering)+",@(0..z)),dp;";
267  ideal nett=imap(gnir,dom);
268  poly p;
269  for (i=1;i<=z;i++)
270  {
271    execute "p=@("+string(i)+");";
272    nett[i]=nett[i]-p;
273  }
274  nett=imap(gnir,p)-@(0),nett;
275  execute "ideal ker=eliminate(nett,"+eli+");";
276  // test wether @(0)-h(@(1),...,@(z)) is in ker
277  l[1]=0;
278  l[2]="";
279  for (i=1;i<=size(ker);i++)
280  {
281     if (deg(ker[i]/@(0))==0)
282     {
283        string str=string(ker[i]);
284        setring gnir;
285        l[1]=1;
286        l[2]=str;
287        return(l);
288     }
289     if (deg(ker[i]/@(0))>0)
290     {
291        l[2]=l[2]+string(ker[i]);
292     }
293  }
294  setring gnir;
295  return(l);
296}
297example
298{ "EXAMPLE:"; echo = 2;
299   ring q=0,(x,y,z,u,v,w),dp;
300   poly p=xyzu2w-1yzu2w2+u4w2-1xu2vw+u2vw2+xyz-1yzw+2u2w-1xv+vw+2;
301   ideal dom =x-w,u2w+1,yz-v;
302   inSubring(p,dom);
303}
304
305///////////////////////////////////////////////////////////////////////////////
306
307proc localInvar(matrix m, poly p, poly q, poly h)
308USAGE:   localInvar(m,p,q,h); m matrix, p,q,h poly
309RETURN:  poly= the invariant of the vectorfield m=Sum m[i,1]*d/dx(i) with respect
310               to p,q,h, i.e.
311               Sum (-1)^v*(1/v!)*m^v(p)*(q/m(q))^v)*m(q)^N, m^N(q)=0, m^(N-1)(q)<>0
312               it is assumed that m(q) and h are invariant
313               the sum above is divided by h as much as possible
314NOTE:
315EXAMPLE: example localInvar; shows an example
316{
317  if ((der(m,h) !=0) || (der(m,der(m,q)) !=0))
318  {
319    "the last variable defines not an invariant function ";
320    return(q);
321  }
322  poly inv=p;
323  poly dif= der(m,inv);
324  poly a=der(m,q);
325  poly sgn=-1;
326  poly coeff=sgn*q;
327  int k=1;
328  if (dif==0)
329  {
330    return(inv);
331  }
332  while (dif!=0)
333  {
334    inv=(a*inv)+(coeff*dif);
335    dif=der(m,dif);
336    k=k+1;
337    coeff=q*coeff*sgn/k;
338  }
339  while ((inv!=0) && (inv!=h) &&(subst(inv,h,0)==0))
340 {
341   inv=inv/h;
342  }
343  return(inv);
344}
345example
346{ "EXAMPLE:"; echo = 2;
347   ring q=0,(x,y,z),dp;
348   matrix m[3][1];
349   m[2,1]=x;
350   m[3,1]=y;
351   poly in=localInvar(m,z,y,x);
352   in;
353}
354///////////////////////////////////////////////////////////////////////////////
355
356
357
358proc furtherInvar(matrix m, ideal id, ideal karl, poly q)
359USAGE:   furtherInvar(m,id,karl,q); m matrix, id,karl ideal,q poly
360RETURN:  ideal= further invariants of the vectorfield m=Sum m[i,1]*d/dx(i) with respect
361               to id,p,q, i.e.
362               the ideal id contains invariants of m and we are looking for elements
363               in the subring generated by id which are divisible by q
364               it is assumed that m(p) and q are invariant
365               the elements mentioned  above are computed and divided by q
366               as much as possible
367               the ideal karl contains all invariants computed yet
368NOTE:
369EXAMPLE: example furtherInvar; shows an example
370{
371  int i;
372  ideal null;
373  int z=size(id);
374  intvec v;
375  def @r=basering;
376  ideal su;
377  for (i=1;i<=z;i++)
378  {
379    su[i]=subst(id[i],q,0);
380  }
381  // defines the map phi : r1 ---> @r defined by
382  // y(i) ---> id[i](q=0)
383  execute "ring r1="+charstr(basering)+",(y(1..z)),dp";
384  setring @r;
385  map phi=r1,su;
386  setring r1;
387  // computes the kernel of phi
388  execute "ideal ker=preimage(@r,phi,null)";
389  // defines the map psi : r1 ---> @r defined by y(i) ---> id[i]
390  setring @r;
391  map psi=r1,id;
392  // computes psi(ker(phi))
393  ideal rel=psi(ker);
394  // devides by the maximal power of q
395  // and tests wether we really obtain invariants
396  for (i=1;i<=size(rel);i++)
397  {
398    while ((rel[i]!=0) && (rel[i]!=q) &&(subst(rel[i],q,0)==0))
399    {
400      rel[i]=rel[i]/q;
401      if (der(m,rel[i])!=0)
402      {
403         "error in furtherInvar, function not invariant";
404         rel[i];
405      }
406    }
407    rel[i]=simplify(rel[i],1);
408  }
409  // test whether some variables occur linearly
410  // and deletes the corresponding invariant function
411  setring r1;
412  int j;
413  for (i=1;i<=size(ker);i=i+1)
414  {
415     for (j=1;j<=z;j++)
416     {
417        if (deg(ker[i]/y(j))==0)
418        {
419           setring @r;
420           rel[i]= completeReduction(rel[i],karl,q);
421           if(rel[i]!=0)
422           {
423              karl[j+1]=rel[i];
424              rel[i]=0;
425           }
426           setring r1;
427        }
428     }
429
430  }
431  setring @r;
432  list l=rel+null,karl;
433  return(l);
434}
435example
436{ "EXAMPLE:"; echo = 2;
437   ring r=0,(x,y,z,u),dp;
438   matrix m[4][1];
439   m[2,1]=x;
440   m[3,1]=y;
441   m[4,1]=z;
442   ideal id=localInvar(m,z,y,x),localInvar(m,u,y,x);
443   ideal karl=id,x;
444   list in=furtherInvar(m,id,karl,x);
445   in;
446}
447///////////////////////////////////////////////////////////////////////////////
448
449
450
451proc invariantRing(matrix m, poly p, poly q,list #)
452USAGE:   invariantRing(m,p,q); m matrix, p,q poly
453RETURN:  ideal= the invariants of the vectorfield m=Sum m[i,1]*d/dx(i)
454                p,q variables with m(p)=q invariant
455NOTE:
456EXAMPLE: example furtherInvar; shows an example
457{
458  ideal j;
459  int i,it;
460  int bou=-1;
461  if(size(#)>0)
462  {
463     bou=#[1];
464  }
465  int z;
466  ideal karl;
467  ideal k1=1;
468  list k2;
469  // computation of local invariants
470  for (i=1;i<=nvars(basering);i++)
471  {
472    karl=karl+localInvar(m,var(i),p,q);
473  }
474  if(bou==0)
475  {
476     "                     ";
477     "the local invariants:";
478     "                     ";
479     karl;
480    // pause;
481     "                     ";
482  }
483  // computation of further invariants
484  it=0;
485  while (size(k1)!=0)
486 {
487    // test if the new invariants are already in the ring generated
488    // by the invariants we constructed already
489    it++;
490    karl=sortier(karl);
491    j=q;
492    for (i=1;i<=size(karl);i++)
493    {
494       j=j+ simplify(completeReduction(karl[i],j,q),1);
495    }
496    karl=j;
497    j[1]=0;
498    j=simplify(j,2);
499    k2=furtherInvar(m,j,karl,q);
500    k1=k2[1];
501    karl=k2[2];
502    k1=sortier(k1);
503    z=size(k1);
504    for (i=1;i<=z;i++)
505    {
506      k1[i]= completeReduction(k1[i],karl,q);
507      if (k1[i]!=0)
508      {
509        karl=karl+simplify(k1[i],1);
510      }
511    }
512    if(bou==0)
513    {
514       "                                  ";
515       "the invariants after the iteration";
516       it;
517       "                                  ";
518       karl;
519      // pause;
520       "                                  ";
521    }
522    if((bou>0) && (size(k1)>0))
523    {
524      if(deg(k1[size(k1)])>bou)
525      {
526         return(karl);
527      }
528    }
529  }
530  return(karl);
531}
532example
533{ "EXAMPLE:"; echo = 2;
534
535  //Winkelmann: free action but Spec k[x(1),...,x(5)]---> Spec In-
536  //variantring is not surjective
537
538  ring rw=0,(x(1..5)),dp;
539  matrix m[5][1];
540  m[3,1]=x(1);
541  m[4,1]=x(2);
542  m[5,1]=1+x(1)*x(4)+x(2)*x(3);
543  ideal in=invariantRing(m,x(3),x(1),0);
544  in;
545
546  //Deveney/Finston: The ring of invariants is not finitely generated
547
548  ring rf=0,(x(1..7)),dp;
549  matrix m[7][1];
550  m[4,1]=x(1)^3;
551  m[5,1]=x(2)^3;
552  m[6,1]=x(3)^3;
553  m[7,1]=(x(1)*x(2)*x(3))^2;
554  ideal in=invariantRing(m,x(4),x(1),6);
555  in;
556
557
558  //Deveney/Finston:Proper Ga-action which is not locally trivial
559  //r[x(1),...,x(5)] is not flat over the ring of invariants
560
561  ring rd=0,(x(1..5)),dp;
562  matrix m[5][1];
563  m[3,1]=x(1);
564  m[4,1]=x(2);
565  m[5,1]=1+x(1)*x(4)^2;
566  ideal in=invariantRing(m,x(3),x(1));
567  in;
568
569  actionIsProper(m);
570
571  //computes the relations between the invariants
572  int z=size(in);
573  ideal null;
574  ring r1=0,(y(1..z)),dp;
575  setring rd;
576  map phi=r1,in;
577  setring r1;
578  ideal ker=preimage(rd,phi,null);
579  ker;
580
581  //the discriminant
582
583  ring r=0,(x(1..2),y(1..2),z,t),dp;
584  poly p=z+(1+x(1)*y(2)^2)*t+x(1)*y(1)*y(2)*t^2+(1/3)*x(1)*y(1)^2*t^3;
585
586  matrix m[5][5];
587  m[1,1]=z;
588  m[1,2]=x(1)*y(2)^2+1;
589  m[1,3]=x(1)*y(1)*y(2);
590  m[1,4]=1/3*x(1)*y(1)^2;
591  m[1,5]=0;
592  m[2,1]=0;
593  m[2,2]=z;
594  m[2,3]=x(1)*y(2)^2+1;
595  m[2,4]=x(1)*y(1)*y(2);
596  m[2,5]=1/3*x(1)*y(1)^2;
597  m[3,1]=x(1)*y(2)^2+1;
598  m[3,2]=2*x(1)*y(1)*y(2);
599  m[3,3]=x(1)*y(1)^2;
600  m[3,4]=0;
601  m[3,5]=0;
602  m[4,1]=0;
603  m[4,2]=x(1)*y(2)^2+1;
604  m[4,3]=2*x(1)*y(1)*y(2);
605  m[4,4]=x(1)*y(1)^2;
606  m[4,5]=0;
607  m[5,1]=0;
608  m[5,2]=0;
609  m[5,3]=x(1)*y(2)^2+1;
610  m[5,4]=2*x(1)*y(1)*y(2);
611  m[5,5]=x(1)*y(1)^2;
612
613  poly disc=9*det(m)/(x(1)^2*y(1)^4);
614
615  LIB "invar.lib";
616  matrix n[6][1];
617  n[2,1]=x(1);
618  n[4,1]=y(1);
619  n[5,1]=1+x(1)*y(2)^2;
620
621  der(n,disc);
622
623  //x(1)^3*y(2)^6-6*x(1)^2*y(1)*y(2)^3*z+6*x(1)^2*y(2)^4+9*x(1)*y(1)^2*z^2-18*x(1)*y(1)*y(2)*z+9*x(1)*y(2)^2+4
624
625
626  //constructive approach to Weizenbcks theorem
627
628  int n=5;
629
630  ring w=0,(x(1..n)),wp(1..n);
631
632  // definition of the vectorfield m=sum m[i]*d/dx(i)
633  matrix m[n][1];
634  int i;
635  for (i=1;i<=n-1;i=i+1)
636  {
637    m[i+1,1]=x(i);
638  }
639  ideal in=invariantRing(m,x(2),x(1),0);
640  in;
641
642
643
644}
Note: See TracBrowser for help on using the repository browser.