# source:git/Singular/LIB/poly.lib@5480da

spielwiese
Last change on this file since 5480da was 5480da, checked in by Kai Krüger <krueger@…>, 26 years ago
• Property mode set to 100644
File size: 20.2 KB
Line
1// \$Id: poly.lib,v 1.10 1998-04-03 22:47:09 krueger Exp \$
2//system("random",787422842);
5///////////////////////////////////////////////////////////////////////////////
6
7version="\$Id: poly.lib,v 1.10 1998-04-03 22:47:09 krueger Exp \$";
8info="
9LIBRARY:  poly.lib      PROCEDURES FOR MANIPULATING POLYS, IDEALS, MODULES
10
11 cyclic(int);           ideal of cyclic n-roots
12 katsura([i]);          katsura [i] ideal
13 freerank(poly/...)     rank of coker(input) if coker is free else -1
14 is_homog(poly/...);    int, =1 resp. =0 if input is homogeneous resp. not
15 is_zero(poly/...);     int, =1 resp. =0 if coker(input) is 0 resp. not
16 lcm(ideal);            lcm of given generators of ideal
17 maxcoef(poly/...);     maximal length of coefficient occuring in poly/...
18 maxdeg(poly/...);      int/intmat = degree/s of terms of maximal order
19 maxdeg1(poly/...);     int = [weighted] maximal degree of input
20 mindeg(poly/...);      int/intmat = degree/s of terms of minimal order
21 mindeg1(poly/...);     int = [weighted] minimal degree of input
22 normalize(poly/...);   normalize poly/... such that leading coefficient is 1
24           (parameters in square brackets [] are optional)
25";
26
27LIB "general.lib";
28///////////////////////////////////////////////////////////////////////////////
29
30proc cyclic (int n)
31USAGE:   cyclic(n);  n integer
32RETURN:  ideal of cyclic n-roots from 1-st n variables of basering
33EXAMPLE: example cyclic; shows examples
34{
35//----------------------------- procedure body --------------------------------
36   ideal m = maxideal(1);
37   m = m[1..n],m[1..n];
38   int i,j;
39   ideal s; poly t;
40   for ( j=0; j<=n-2; j=j+1 )
41   {
42      t=0;
43      for( i=1;i<=n;i=i+1 ) { t=t+product(m,i..i+j); }
44      s=s+t;
45   }
46   s=s,product(m,1..n)-1;
47   return (s);
48}
49//-------------------------------- examples -----------------------------------
50example
51{ "EXAMPLE:"; echo = 2;
52   ring r=0,(u,v,w,x,y,z),lp;
53   cyclic(nvars(basering));
54   homog(cyclic(5),z);
55}
56///////////////////////////////////////////////////////////////////////////////
57
58proc katsura
59USAGE:  katsura([n]); n integer
60RETURN: katsura(n) : n-th katsura ideal of newly created and set ring
61                     (32003, x(0..n), dp)
62        katsura()  : katsura ideal of basering
63EXAMPLE: example katsura; shows examples
64{
65  if ( size(#) == 1 && typeof(#[1]) == "int")
66  {
67    ring katsura_ring = 32003, x(0..#[1]), dp;
68    keepring katsura_ring;
69  }
70  ideal s;
71  int i, j;
72  int n = nvars(basering) -1;
73  poly p;
74
75  p = -1;
76  for (i = -n; i <= n; i++)
77  {
78    p = p + kat_var(i, n);
79  }
80  s[1] = p;
81
82  for (i = 0; i < n; i++)
83  {
84    p = -1 * kat_var(i,n);
85    for (j = -n; j <= n; j++)
86    {
87      p = p + kat_var(j,n) * kat_var(i-j, n);
88    }
89    s = s,p;
90  }
91  return (s);
92}
93//-------------------------------- examples -----------------------------------
94example
95{
96  "EXAMPLE:"; echo = 2;
97  ring r;
98  katsura();
99  katsura(3);
100}
101
102proc kat_var(int i, int n)
103{
104  poly p;
105  if (i < 0)  { i = -i;}
106  if (i <= n) { p = var(i+1); }
107  return (p);
108}
109///////////////////////////////////////////////////////////////////////////////
110
111proc freerank
112USAGE:   freerank(M[,any]);  M=poly/ideal/vector/module/matrix
113COMPUTE: rank of module presented by M in case it is free. By definition this
114         is vdim(coker(M)/m*coker(M)) if coker(M) is free, where m = maximal
115         ideal of basering and M is considered as matrix (the 0-module is
116         free of rank 0)
117RETURN:  rank of coker(M) if coker(M) is free and -1 else;
118         in case of a second argument return a list:
119                L[1] = rank of coker(M) or -1
120                L[2] = minbase(M)
121NOTE:    freerank(syz(M)); computes the rank of M if M is free (and -1 else)
122         //* Zur Zeit noch ein Bug, da erste Bettizahl falsch berechnet wird:
123         //betti(0) ist -1 statt 0
124EXAMPLE: example freerank; shows examples
125{
126  int rk;
127  def M = simplify(#[1],10);
128  list mre = mres(M,2);
129  intmat B = betti(mre);
130  if ( ncols(B)>1 ) { rk = -1; }
131  else { rk = sum(B[1..nrows(B),1]); }
132  if (size(#) == 2) { list L=rk,mre[1]; return(L);}
133  return(rk);
134}
135example
136{"EXAMPLE";   echo=2;
137  ring r;
138  ideal i=x;
139  module M=[x,0,1],[-x,0,-1];
140  freerank(M);           // should be -1, coker(M) is not free
141                         // [1] should be 1, coker(syz(M))=M is free of rank 1
142  freerank(syz (M),"");  // [2] should be gen(2)+gen(1) (minimal relation of M)
143  freerank(i);
144  freerank(syz(i));      //* bug, should be 1, coker(syz(i))=i is free of rank 1
145}
146///////////////////////////////////////////////////////////////////////////////
147
148proc is_homog (id)
149USAGE:   is_homog(id);  id  poly/ideal/vector/module/matrix
150RETURN:  integer which is 1 if input is homogeneous (resp. weighted homogeneous
151         if the monomial ordering consists of one block of type ws,Ws,wp or Wp,
152         assuming that all weights are positive) and 0 otherwise
153NOTE:    A vector is homogeneous, if the components are homogeneous of same
154         degree, a module/matrix is homogeneous if all column vectors are
155         homogeneous
156         //*** ergaenzen, wenn Matrizen-Spalten Gewichte haben
157EXAMPLE: example is_homog; shows examples
158{
159//----------------------------- procedure body --------------------------------
160   module M = module(matrix(id));
161   M = simplify(M,2);                        // remove 0-columns
162   intvec v = ringweights(basering);
163   int i,j=1,1;
164   for (i=1; i<=ncols(M); i=i+1)
165   {
167      { return(0); }
168   }
169   return(1);
170}
171//-------------------------------- examples -----------------------------------
172example
173{ "EXAMPLE:"; echo = 2;
174   ring r = 0,(x,y,z),wp(1,2,3);
175   is_homog(x5-yz+y3);
176   ideal i = x6+y3+z2, x9-z3;
177   is_homog(i);
178   ring s = 0,(a,b,c),ds;
179   vector v = [a2,0,ac+bc];
180   vector w = [a3,b3,c4];
181   is_homog(v);
182   is_homog(w);
183}
184///////////////////////////////////////////////////////////////////////////////
185
186proc is_zero
187USAGE:   is_zero(M[,any]); M=poly/ideal/vector/module/matrix
188RETURN:  integer, 1 if coker(M)=0 resp. 0 if coker(M)!=0, where M is considered
189         as matrix
190         if a second argument is given, return a list:
191                L[1] = 1 if coker(M)=0 resp. 0 if coker(M)!=0
192                L[2] = dim(M)
193EXAMPLE: example is_zero; shows examples
194{
195  int d=dim(std(#[1]));
196  int a = ( d==-1 );
197  if( size(#) >1 ) { list L=a,d; return(L); }
198  return(a);
199}
200example
201{ "EXAMPLE:";   echo=2;
202  ring r;
203  module m = [x],[y],[1,z];
204  is_zero(m,1);
205  qring q = std(ideal(x2+y3+z2));
206  ideal j = x2+y3+z2-37;
207  is_zero(j);
208}
209////////////////////////////////////////////////////////////////////////////////
210
211proc maxcoef (f)
212USAGE:   maxcoef(f);  f  poly/ideal/vector/module/matrix
213RETURN:  maximal length of coefficient of f of type int (by counting the
214         length of the string of each coefficient)
215EXAMPLE: example maxcoef; shows examples
216{
217//----------------------------- procedure body --------------------------------
218   int max,s,ii,jj; string t;
219   ideal i = ideal(matrix(f));
220   i = simplify(i,6);            // delete 0's and keep first of equal elements
221   poly m = var(1); matrix C;
222   for (ii=2;ii<=nvars(basering);ii=ii+1) { m = m*var(ii); }
223   for (ii=1; ii<=size(i); ii=ii+1)
224   {
225      C = coef(i[ii],m);
226      for (jj=1; jj<=ncols(C); jj=jj+1)
227      {
228         t = string(C[2,jj]);  s = size(t);
229         if ( t[1] == "-" ) { s = s - 1; }
230         if ( s > max ) { max = s; }
231      }
232   }
233   return(max);
234}
235//-------------------------------- examples -----------------------------------
236example
237{ "EXAMPLE:"; echo = 2;
238   ring r= 0,(x,y,z),ds;
239   poly g = 345x2-1234567890y+7/4z;
240   maxcoef(g);
241   ideal i = g,10/1234567890;
242   maxcoef(i);
243   // since i[2]=1/123456789
244}
245///////////////////////////////////////////////////////////////////////////////
246
247proc maxdeg (id)
248USAGE:   maxdeg(id);  id  poly/ideal/vector/module/matrix
249RETURN:  int/intmat, each component equals maximal degree of monomials in the
250         corresponding component of id, independent of ring ordering
251         (maxdeg of each var is 1)
252         of type int if id is of type poly, of type intmat else
253NOTE:    proc maxdeg1 returns 1 integer, the absolut maximum; moreover, it has
254         an option for computing weighted degrees
255EXAMPLE: example maxdeg; shows examples
256{
257   //-------- subprocedure to find maximal degree of given component ----------
258   proc findmaxdeg
259   {
260      poly c = #[1];
261      if (c==0) { return(-1); }
262   //--- guess upper 'o' and lower 'u' bound, in case of negative weights -----
263      int d = (deg(c)>=0)*deg(c)-(deg(c)<0)*deg(c);
264      int i = d;
265      while ( c-jet(c,i) != 0 ) { i = 2*(i+1); }
266      int o = i-1;
267      int u = (d != i)*((i div  2)-1);
268   //----------------------- "quick search" for maxdeg ------------------------
269      while ( (c-jet(c,i)==0)*(c-jet(c,i-1)!=0) == 0)
270      {
271         i = (o+1+u) div  2;
272         if (c-jet(c,i)!=0) { u = i+1; }
273         else { o = i-1; }
274      }
275      return(i);
276   }
277//------------------------------ main program ---------------------------------
278   matrix M = matrix(id);
279   int r,c = nrows(M), ncols(M); int i,j;
280   intmat m[r][c];
281   for (i=r; i>0; i=i-1)
282   {
283      for (j=c; j>0; j=j-1) { m[i,j] = findmaxdeg(M[i,j]); }
284   }
285   if (typeof(id)=="poly") { return(m[1,1]); }
286   return(m);
287}
288//-------------------------------- examples -----------------------------------
289example
290{ "EXAMPLE:"; echo = 2;
291   ring r = 0,(x,y,z),wp(-1,-2,-3);
292   poly f = x+y2+z3;
293   deg(f);               //deg; returns weighted degree (in case of 1 block)!
294   maxdeg(f);
295   matrix m[2][2]=f+x10,1,0,f^2;
296   maxdeg(m);
297}
298///////////////////////////////////////////////////////////////////////////////
299
300proc maxdeg1 (id,list #)
301USAGE:   maxdeg1(id[,v]);  id=poly/ideal/vector/module/matrix, v=intvec
302RETURN:  integer, maximal [weighted] degree of monomials of id independent of
303         ring ordering, maxdeg1 of i-th variable is v[i] (default: v=1..1).
304NOTE:    This proc returns one integer while maxdeg returns, in general,
305         a matrix of integers. For one polynomial and if no intvec v is given
306         maxdeg is faster
307EXAMPLE: example maxdeg1; shows examples
308{
309   //-------- subprocedure to find maximal degree of given component ----------
310   proc findmaxdeg
311   {
312      poly c = #[1];
313      if (c==0) { return(-1); }
314      intvec v = #[2];
315   //--- guess upper 'o' and lower 'u' bound, in case of negative weights -----
316      int d = (deg(c)>=0)*deg(c)-(deg(c)<0)*deg(c);
317      int i = d;
318      if ( c == jet(c,-1,v))      //case: maxdeg is negative
319      {
320         i = -d;
321         while ( c == jet(c,i,v) ) { i = 2*(i-1); }
322         int o = (d != -i)*((i div  2)+2) - 1;
323         int u = i+1;
324         int e = -1;
325      }
326      else                        //case: maxdeg is nonnegative
327      {
328         while ( c != jet(c,i,v) ) { i = 2*(i+1); }
329         int o = i-1;
330         int u = (d != i)*((i div  2)-1);
331         int e = 1;
332      }
333   //----------------------- "quick search" for maxdeg ------------------------
334      while ( ( c==jet(c,i,v) )*( c!=jet(c,i-1,v) ) == 0 )
335      {
336         i = (o+e+u) div  2;
337         if ( c!=jet(c,i,v) ) { u = i+1; }
338         else { o = i-1; }
339      }
340      return(i);
341   }
342//------------------------------ main program ---------------------------------
343   ideal M = simplify(ideal(matrix(id)),8);   //delete scalar multiples from id
344   int c = ncols(M);
345   int i,n;
346   if( size(#)==0 )
347   {
348      int m = maxdeg(M[c]);
349      for (i=c-1; i>0; i=i-1)
350      {
351          n = maxdeg(M[i]);
352          m = (m>=n)*m + (m<n)*n;             //let m be the maximum of m and n
353      }
354   }
355   else
356   {
357      intvec v=#[1];                          //weight vector for the variables
358      int m = findmaxdeg(M[c],v);
359      for (i=c-1; i>0; i--)
360      {
361         n = findmaxdeg(M[i],v);
362         if( n>m ) { m=n; }
363      }
364   }
365   return(m);
366}
367//-------------------------------- examples -----------------------------------
368example
369{ "EXAMPLE:"; echo = 2;
370   ring r = 0,(x,y,z),wp(-1,-2,-3);
371   poly f = x+y2+z3;
372   deg(f);                  //deg returns weighted degree (in case of 1 block)!
373   maxdeg1(f);
374   intvec v = ringweights(r);
375   maxdeg1(f,v);                             //weighted maximal degree
376   matrix m[2][2]=f+x10,1,0,f^2;
377   maxdeg1(m,v);                             //absolut weighted maximal degree
378}
379///////////////////////////////////////////////////////////////////////////////
380
381proc mindeg (id)
382USAGE:   mindeg(id);  id  poly/ideal/vector/module/matrix
383RETURN:  minimal degree/s of monomials of id, independent of ring ordering
384         (mindeg of each variable is 1) of type int if id of type poly, else
385         of type intmat.
386NOTE:    proc mindeg1 returns one integer, the absolut minimum; moreover it
387         has an option for computing weighted degrees.
388EXAMPLE: example mindeg; shows examples
389{
390   //--------- subprocedure to find minimal degree of given component ---------
391   proc findmindeg
392   {
393      poly c = #[1];
394      if (c==0) { return(-1); }
395   //--- guess upper 'o' and lower 'u' bound, in case of negative weights -----
396      int d = (ord(c)>=0)*ord(c)-(ord(c)<0)*ord(c);
397      int i = d;
398      while ( jet(c,i) == 0 ) { i = 2*(i+1); }
399      int o = i-1;
400      int u = (d != i)*((i div  2)-1);
401   //----------------------- "quick search" for mindeg ------------------------
402      while ( (jet(c,u)==0)*(jet(c,o)!=0) )
403      {
404         i = (o+u) div  2;
405         if (jet(c,i)==0) { u = i+1; }
406         else { o = i-1; }
407      }
408      if (jet(c,u)!=0) { return(u); }
409      else { return(o+1); }
410   }
411//------------------------------ main program ---------------------------------
412   matrix M = matrix(id);
413   int r,c = nrows(M), ncols(M); int i,j;
414   intmat m[r][c];
415   for (i=r; i>0; i=i-1)
416   {
417      for (j=c; j>0; j=j-1) { m[i,j] = findmindeg(M[i,j]); }
418   }
419   if (typeof(id)=="poly") { return(m[1,1]); }
420   return(m);
421}
422//-------------------------------- examples -----------------------------------
423example
424{ "EXAMPLE:"; echo = 2;
425   ring r = 0,(x,y,z),ls;
426   poly f = x5+y2+z3;
427   ord(f);                      // ord returns weighted order of leading term!
428   mindeg(f);                   // computes minimal degree
429   matrix m[2][2]=x10,1,0,f^2;
430   mindeg(m);                   // computes matrix of minimum degrees
431}
432///////////////////////////////////////////////////////////////////////////////
433
434proc mindeg1 (id, list #)
435USAGE:   mindeg1(id[,v]);  id=poly/ideal/vector/module/matrix, v=intvec
436RETURN:  integer, minimal [weighted] degree of monomials of id independent of
437         ring ordering, mindeg1 of i-th variable is v[i] (default v=1..1).
438NOTE:    This proc returns one integer while mindeg returns, in general,
439         a matrix of integers. For one polynomial and if no intvec v is given
440         mindeg is faster.
441EXAMPLE: example mindeg1; shows examples
442{
443   //--------- subprocedure to find minimal degree of given component ---------
444   proc findmindeg
445   {
446      poly c = #[1];
447      intvec v = #[2];
448      if (c==0) { return(-1); }
449   //--- guess upper 'o' and lower 'u' bound, in case of negative weights -----
450      int d = (ord(c)>=0)*ord(c)-(ord(c)<0)*ord(c);
451      int i = d;
452      if ( jet(c,-1,v) !=0 )      //case: mindeg is negative
453      {
454         i = -d;
455         while ( jet(c,i,v) != 0 ) { i = 2*(i-1); }
456         int o = (d != -i)*((i div  2)+2) - 1;
457         int u = i+1;
458         int e = -1; i=u;
459      }
460      else                        //case: inded is nonnegative
461      {
462         while ( jet(c,i,v) == 0 ) { i = 2*(i+1); }
463         int o = i-1;
464         int u = (d != i)*((i div  2)-1);
465         int e = 1; i=u;
466      }
467   //----------------------- "quick search" for mindeg ------------------------
468      while ( (jet(c,i-1,v)==0)*(jet(c,i,v)!=0) == 0 )
469      {
470         i = (o+e+u) div  2;
471         if (jet(c,i,v)==0) { u = i+1; }
472         else { o = i-1; }
473      }
474      return(i);
475   }
476//------------------------------ main program ---------------------------------
477   ideal M = simplify(ideal(matrix(id)),8);   //delete scalar multiples from id
478   int c = ncols(M);
479   int i,n;
480   if( size(#)==0 )
481   {
482      int m = mindeg(M[c]);
483      for (i=c-1; i>0; i=i-1)
484      {
485          n = mindeg(M[i]);
486          m = (m<=n)*m + (m>n)*n;             //let m be the maximum of m and n
487      }
488   }
489   else
490   {
491      intvec v=#[1];                          //weight vector for the variables
492      int m = findmindeg(M[c],v);
493      for (i=c-1; i>0; i=i-1)
494      {
495         n = findmindeg(M[i],v);
496         m = (m<=n)*m + (m>n)*n;              //let m be the maximum of m and n
497      }
498   }
499   return(m);
500}
501//-------------------------------- examples -----------------------------------
502example
503{ "EXAMPLE:"; echo = 2;
504   ring r = 0,(x,y,z),ls;
505   poly f = x5+y2+z3;
506   ord(f);                      // ord returns weighted order of leading term!
507   intvec v = 1,-3,2;
508   mindeg1(f,v);                // computes minimal weighted degree
509   matrix m[2][2]=x10,1,0,f^2;
510   mindeg1(m,1..3);             // computes absolut minimum of weighted degrees
511}
512///////////////////////////////////////////////////////////////////////////////
513
514proc normalize (id)
515USAGE:   normalize(id);  id=poly/vector/ideal/module
516RETURN:  object of same type with leading coefficient equal to 1
517EXAMPLE: example normalize; shows an example
518{
519   return(simplify(id,1));
520}
521//-------------------------------- examples -----------------------------------
522example
523{  "EXAMPLE:"; echo = 2;
524   ring r = 0,(x,y,z),ls;
525   poly f = 2x5+3y2+4z3;
526   normalize(f);
527   module m=[9xy,0,3z3],[4z,6y,2x];
528   normalize(m);
529   ring s = 0,(x,y,z),(c,ls);
530   module m=[9xy,0,3z3],[4z,6y,2x];
531   normalize(m);
532   normalize(matrix(m));             // by automatic type conversion to module!
533}
534///////////////////////////////////////////////////////////////////////////////
535
536////////////////////////////////////////////////////////////////////////////////
537// Input: <ideal>=<f1,f2,...,fm> and <polynomial> g
538// Question: Does g lie in the radical of <ideal>?
539// Solution: Compute a standard basis G for <f1,f2,...,fm,gz-1> where z is a new
540//           variable. Then g is contained in the radical of <ideal> <=> 1 is
541//           generator in G.
542////////////////////////////////////////////////////////////////////////////////
545  RETURNS: 1 (TRUE) (type <int>) if <poly> is contained in the radical of
546           <ideal>, 0 (FALSE) (type <int>) otherwise
547  EXAMPLE: example rad_con; shows an example
548{ def br=basering;
549  int n=nvars(br);
550  int dB=degBound;
551  degBound=0;
552  string mp=string(minpoly);
553  execute "ring R=("+charstr(br)+"),(x(1..n),z),dp;";
554  execute "minpoly=number("+mp+");";
555  ideal irrel=x(1..n);
556  map f=br,irrel;
557  poly p=f(g);
558  ideal J=f(I)+ideal(p*z-1);
559  J=std(J);
560  degBound=dB;
561  if (J[1]==1)
562  { return(1);
563  }
564  else
565  { return(0);
566  }
567}
568example
569{ "  EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7.";
570  echo=2;
571           ring R=0,(x,y,z),dp;
572           ideal I=x2+y2,z2;
573           poly f=x4+y4;
575           ideal J=x2+y2,z2,x4+y4;
576           poly g=z;
578}
579
580///////////////////////////////////////////////////////////////////////////////
581
582proc lcm (ideal i)
583USAGE:   lcm(i); i ideal
584RETURN:  poly = lcm(i[1],...,i[size(i)])
585NOTE:
586EXAMPLE: example lcm; shows an example
587{
588  int k,j;
589   poly p,q;
590  i=simplify(i,10);
591  for(j=1;j<=size(i);j++)
592  {
593    if(deg(i[j])>0)
594    {
595      p=i[j];
596      break;
597    }
598  }
599  if(deg(p)==-1)
600  {
601    return(1);
602  }
603  for (k=j+1;k<=size(i);k++)
604  {
605     if(deg(i[k])!=0)
606     {
607        q=gcd(p,i[k]);
608        if(deg(q)==0)
609        {
610           p=p*i[k];
611        }
612        else
613        {
614           p=p/q;
615           p=p*i[k];
616        }
617     }
618   }
619  return(p);
620}
621example
622{ "EXAMPLE:"; echo = 2;
623   ring  r = 0,(x,y,z),lp;
624   poly  p = (x+y)*(y+z);
625   poly  q = (z4+2)*(y+z);
626   ideal l=p,q;
627   poly  pr= lcm(l);
628   pr;
629   l=1,-1,p,1,-1,q,1;
630   pr=lcm(l);
631   pr;
632}
633
634///////////////////////////////////////////////////////////////////////////////
635
636proc content(f)
637USAGE:   content(f); f polynomial/vector
638RETURN:  number, the content (greatest common factor of coefficients)
639         of the polynomial/vector f
640EXAMPLE: example content; shows an example
641{