source: git/Singular/LIB/poly.lib @ a57b65

spielwiese
Last change on this file since a57b65 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: 33.5 KB
Line 
1////////////////////////////////////////////////////////////////////////////
2version="version poly.lib 4.0.0.0 Jun_2013 "; // $Id$
3category="General purpose";
4info="
5LIBRARY:  poly.lib      Procedures for Manipulating Polys, Ideals, Modules
6AUTHORS:  O. Bachmann, G.-M. Greuel, A. Fruehbis
7
8PROCEDURES:
9 cyclic(int);           ideal of cyclic n-roots
10 elemSymmId(int);       ideal of elementary symmetric polynomials
11 katsura([i]);          katsura [i] ideal
12 freerank(poly/...)     rank of coker(input) if coker is free else -1
13 is_zero(poly/...);     int, =1 resp. =0 if coker(input) is 0 resp. not
14 lcm(ideal);            lcm of given generators of ideal
15 maxcoef(poly/...);     maximal length of coefficient occurring in poly/...
16 maxdeg(poly/...);      int/intmat = degree/s of terms of maximal order
17 maxdeg1(poly/...);     int = [weighted] maximal degree of input
18 mindeg(poly/...);      int/intmat = degree/s of terms of minimal order
19 mindeg1(poly/...);     int = [weighted] minimal degree of input
20 normalize(poly/...);   normalize poly/... such that leading coefficient is 1
21 rad_con(p,I);          check radical containment of polynomial p in ideal I
22 content(f);            content of polynomial/vector f
23 mod2id(M,iv);          conversion of a module M to an ideal
24 id2mod(i,iv);          conversion inverse to mod2id
25 substitute(I,...)      substitute in I variables by polynomials
26 subrInterred(i1,i2,iv);interred w.r.t. a subset of variables
27 newtonDiag(f);         Newton diagram of a polynomial
28 hilbPoly(I);           Hilbert polynomial of basering/I
29          (parameters in square brackets [] are optional)
30
31";
32
33LIB "general.lib";
34LIB "ring.lib";
35///////////////////////////////////////////////////////////////////////////////
36static proc bino(int a, int b)
37{
38//computes binomial var(1)+a over b
39   int i;
40   if(b==0){return(1);}
41   poly p=(var(1)+a)/b;
42   for(i=1;i<=b-1;i++)
43   {
44      p=p*(var(1)+a-i)/i;
45   }
46   return(p);
47}
48
49proc hilbPoly(ideal I)
50"USAGE: hilbPoly(I); I a homogeneous ideal
51RETURN: the Hilbert polynomial of basering/I as an intvec v=v_0,...,v_r
52        such that the Hilbert polynomial is (v_0+v_1*t+...v_r*t^r)/r!
53EXAMPLE: example hilbPoly; shows an example
54"
55{
56   def R=basering;
57   if(!attrib(I,"isSB")){I=std(I);}
58   intvec v=hilb(I,2);
59   int s=dim(I);
60   intvec hp;
61   if(s==0){return(hp);}
62   int d=size(v)-2;
63   ring S=0,t,dp;
64   poly p=v[1+d]*bino(s-1-d,s-1);
65   int i;
66   for(i=1;i<=d;i++)
67   {
68      p=p+v[d-i+1]*bino(s-1-d+i,s-1);
69   }
70   int n=1;
71   for(i=2;i<=s-1;i++){n=n*i;}
72   p=n*p;
73   hp[s]=int(leadcoef(p));
74   for(i=2;i<=size(p);i++)
75   {
76      hp[leadexp(p[i])+1]=int(leadcoef(p[i]));
77   }
78   setring R;
79   return(hp);
80}
81example
82{ "EXAMPLE:"; echo = 2;
83   ring r = 0,(b,c,t,h),dp;
84   ideal I=
85   bct-t2h+2th2+h3,
86   bt3-ct3-t4+b2th+c2th-2bt2h+2ct2h+2t3h-bch2-2bth2+2cth2+2th3,
87   b2c2+bt2h-ct2h-t3h+b2h2+2bch2+c2h2-2bth2+2cth2+t2h2-2bh3+2ch3+2th3+3h4,
88   c2t3+ct4-c3th-2c2t2h-2ct3h-t4h+bc2h2-2c2th2-bt2h2+4t3h2+2bth3-2cth3-t2h3
89   +bh4-6th4-2h5;
90   hilbPoly(I);
91}
92
93///////////////////////////////////////////////////////////////////////////////
94proc substitute (I,list #)
95"USAGE:  - case 1: typeof(#[1])==poly:
96           substitute (I,v,f[,v1,f1,v2,f2,...]); I object of basering which
97           can be mapped, v,v1,v2,.. ring variables, f,f1,f2,... poly
98@*       - case 2: typeof(#[1])==ideal:
99           substitute (I,v,f); I object of basering which can be mapped,
100           v ideal of ring variables, f ideal
101RETURN:  object of same type as I,
102@*       - case 1: ring variable v,v1,v2,... substituted by polynomials
103           f,f1,f2,..., in this order
104@*       - case 2: ring variables in v substituted by polynomials in f:
105           v[i] is substituted by f[i], i=1,...,i=min(size(v),ncols(f))
106NOTE:    this procedure extends the built-in command subst via maps
107EXAMPLE: example substitute; shows an example
108"
109{
110  def bas = basering;
111  ideal m = maxideal(1);
112  int i,ii;
113  if(typeof(#[1])=="poly")
114  {
115    poly v = #[1];
116    poly f = #[2];
117    map phi = bas,m;
118    def J = I;
119    for (ii=1; ii<=size(#) - 1; ii=ii+2)
120    {
121      m = maxideal(1);
122      i=rvar(#[ii]);
123      m[i] = #[ii+1];
124      phi = bas,m;
125      J = phi(J);
126    }
127    return(J);
128  }
129  if(typeof(#[1])=="ideal")
130  {
131    ideal v = #[1];
132    ideal f = #[2];
133    int mi = ncols(v);
134    if(ncols(f)<mi)
135    {
136      mi = ncols(f);
137    }
138    def J = I;
139    for (ii=1; ii<=mi; ii++)
140    {
141      m[rvar(v[ii])]=f[ii];
142    }
143    map phi = bas,m;
144    J = phi(I);
145    return(J);
146  }
147}
148example
149{ "EXAMPLE:"; echo = 2;
150   ring r = 0,(b,c,t),dp;
151   ideal I = -bc+4b2c2t,bc2t-5b2c;
152   substitute(I,c,b+c,t,0,b,b-1);
153   ideal v = c,t,b;
154   ideal f = b+c,0,b-1;
155   substitute(I,v,f);
156}
157
158///////////////////////////////////////////////////////////////////////////////
159proc cyclic (int n)
160"USAGE:   cyclic(n);  n integer
161RETURN:  ideal of cyclic n-roots from 1-st n variables of basering
162EXAMPLE: example cyclic; shows examples
163"
164{
165//----------------------------- procedure body --------------------------------
166   ideal m = maxideal(1);
167   m = m[1..n],m[1..n];
168   int i,j;
169   ideal s; poly t;
170   for ( j=0; j<=n-2; j++ )
171   {
172      t=0;
173      for( i=1;i<=n;i++ ) { t=t+product(m,i..i+j); }
174      s=s+t;
175   }
176   s=s,product(m,1..n)-1;
177   return (s);
178}
179//-------------------------------- examples -----------------------------------
180example
181{ "EXAMPLE:"; echo = 2;
182   ring r=0,(u,v,w,x,y,z),lp;
183   cyclic(nvars(basering));
184   homog(cyclic(5),z);
185}
186
187///////////////////////////////////////////////////////////////////////////////
188proc elemSymmPoly(int d, int lindex, int hindex)
189"USAGE:  elemSymmPoly(d,lindex,hindex); d,lindex,hindex integers
190RETURN:  elementary symmetric polynomial of degree d for variables
191@*       var(lindex),...,var(hindex) of basering
192EXAMPLE: example elemSymmPoly; shows an example
193"
194{
195 if(hindex - lindex + 1 < d)
196 {
197   int i = hindex - lindex + 1;
198   "========================================================";
199   "There is no elementary symmetric polynomial of degree "+string(d);
200   "for just "+string(i)+" variables.";
201   "========================================================";
202   return(poly(0));
203 }
204 if(d == 0)
205 {
206   return(poly(1));
207 }
208 else
209 {
210   int i;
211   poly p;
212   for (i = lindex; i <=  hindex - d + 1; i++)
213   {
214     p = p + var(i) * elemSymmPoly(d - 1, i + 1, hindex);
215   }
216   return(p);
217 }
218}
219example
220{ "EXAMPLE:"; echo = 2;
221  ring R = 0, (u,v,w,x,y,z), lp;
222  elemSymmPoly(3,2,5);
223  elemSymmPoly(6,1,6);
224}
225
226//////////////////////////////////////////////////////////////////////////
227proc elemSymmId(int n)
228"USAGE:  elemSymmId(n); n integer
229RETURN:  ideal of elementary symmetric polynomials for 1-st n
230@*       variables of basering
231EXAMPLE: example elemSymmId; shows an example
232{
233 int i;
234 ideal symm;
235 for(i = 1; i <= n; i++)
236 {
237   symm = symm + elemSymmPoly(i, 1, n);
238 }
239 return(symm);
240}
241example
242{ "EXAMPLE:"; echo = 2;
243  ring R = 0, (v,w,x,y,z), lp;
244  elemSymmId(3);
245  elemSymmId(nvars(basering));
246}
247
248//////////////////////////////////////////////////////////////////////////
249proc katsura
250"USAGE: katsura([n]); n integer
251RETURN: katsura(n) : n-th katsura ideal of
252         (1) newly created and set ring (32003, x(0..n), dp), if
253             nvars(basering) < n
254         (2) basering, if nvars(basering) >= n
255        katsura()  : katsura ideal of basering
256EXAMPLE: example katsura; shows examples
257"
258{
259  int n;
260  if ( size(#) == 1 && typeof(#[1]) == "int")
261  {
262    n = #[1] - 1;
263    while (1)
264    {
265      if (defined(basering))
266      {
267        if (nvars(basering) >= #[1]) {break;}
268      }
269      ring katsura_ring = 32003, x(0..#[1]), dp;
270      keepring katsura_ring;
271      break;
272    }
273  }
274  else
275  {
276    n = nvars(basering) -1;
277  }
278
279  ideal s;
280  int i, j;
281  poly p;
282
283  p = -1;
284  for (i = -n; i <= n; i++)
285  {
286    p = p + kat_var(i, n);
287  }
288  s[1] = p;
289
290  for (i = 0; i < n; i++)
291  {
292    p = -1 * kat_var(i,n);
293    for (j = -n; j <= n; j++)
294    {
295      p = p + kat_var(j,n) * kat_var(i-j, n);
296    }
297    s = s,p;
298  }
299  return (s);
300}
301//-------------------------------- examples -----------------------------------
302example
303{
304  "EXAMPLE:"; echo = 2;
305  ring r; basering;
306  katsura();
307  katsura(4); basering;
308}
309
310proc kat_var(int i, int n)
311{
312  poly p;
313  if (i < 0)  { i = -i;}
314  if (i <= n) { p = var(i+1); }
315  return (p);
316}
317///////////////////////////////////////////////////////////////////////////////
318
319proc freerank
320"USAGE:   freerank(M[,any]);  M=poly/ideal/vector/module/matrix
321COMPUTE: rank of module presented by M in case it is free.
322         By definition this is vdim(coker(M)/m*coker(M)) if coker(M)
323         is free, where m is the maximal ideal of the variables of the
324         basering and M is considered to be a matrix.
325         (the 0-module is free of rank 0)
326RETURN:  rank of coker(M) if coker(M) is free and -1 else;
327         in case of a second argument return a list:
328                L[1] = rank of coker(M) or -1
329                L[2] = minbase(M)
330NOTE:    freerank(syz(M)); computes the rank of M if M is free (and -1 else)
331EXAMPLE: example freerank; shows examples
332"
333{
334  int rk;
335  def M = simplify(#[1],10);
336  resolution mre = res(M,2);
337  intmat B = betti(mre);
338  if ( ncols(B)>1 ) { rk = -1; }
339  else { rk = sum(B[1..nrows(B),1]); }
340  if (size(#) == 2) { list L=rk,mre[1]; return(L);}
341  return(rk);
342}
343example
344{"EXAMPLE";   echo=2;
345  ring r;
346  ideal i=x;
347  module M=[x,0,1],[-x,0,-1];
348  freerank(M);          // should be 2, coker(M) is not free
349  freerank(syz (M),"");
350  // [1] should be 1, coker(syz(M))=M is free of rank 1
351  // [2] should be gen(2)+gen(1) (minimal relation of M)
352  freerank(i);
353  freerank(syz(i));     // should be 1, coker(syz(i))=i is free of rank 1
354}
355///////////////////////////////////////////////////////////////////////////////
356
357proc is_zero
358"USAGE:   is_zero(M[,any]); M=poly/ideal/vector/module/matrix
359RETURN:  integer, 1 if coker(M)=0 resp. 0 if coker(M)!=0, where M is
360         considered as matrix.
361         If a second argument is given, return a list:
362                L[1] = 1 if coker(M)=0 resp. 0 if coker(M)!=0
363                L[2] = dim(M)
364EXAMPLE: example is_zero; shows examples
365"
366{
367  int d=dim(std(#[1]));
368  int a = ( d==-1 );
369  if( size(#) >1 ) { return(list(a,d)); }
370  return(a);
371}
372example
373{ "EXAMPLE:";   echo=2;
374  ring r;
375  module m = [x],[y],[1,z];
376  is_zero(m,1);
377  qring q = std(ideal(x2+y3+z2));
378  ideal j = x2+y3+z2-37;
379  is_zero(j);
380}
381///////////////////////////////////////////////////////////////////////////////
382
383proc maxcoef (f)
384"USAGE:   maxcoef(f);  f  poly/ideal/vector/module/matrix
385RETURN:  maximal length of coefficient of f of type int (by measuring the
386         length of the string of each coefficient)
387EXAMPLE: example maxcoef; shows examples
388"
389{
390//----------------------------- procedure body --------------------------------
391   int max,s,ii,jj; string t;
392   ideal i = ideal(matrix(f));
393   i = simplify(i,6);            // delete 0's and keep first of equal elements
394   poly m = var(1); matrix C;
395   for (ii=2;ii<=nvars(basering);ii++) { m = m*var(ii); }
396   for (ii=1; ii<=size(i); ii++)
397   {
398      C = coef(i[ii],m);
399      for (jj=1; jj<=ncols(C); jj++)
400      {
401         t = string(C[2,jj]);  s = size(t);
402         if ( t[1] == "-" ) { s = s - 1; }
403         if ( s > max ) { max = s; }
404      }
405   }
406   return(max);
407}
408//-------------------------------- examples -----------------------------------
409example
410{ "EXAMPLE:"; echo = 2;
411   ring r= 0,(x,y,z),ds;
412   poly g = 345x2-1234567890y+7/4z;
413   maxcoef(g);
414   ideal i = g,10/1234567890;
415   maxcoef(i);
416   // since i[2]=1/123456789
417}
418///////////////////////////////////////////////////////////////////////////////
419
420proc maxdeg (id)
421"USAGE:   maxdeg(id);  id  poly/ideal/vector/module/matrix
422RETURN:  int/intmat, each component equals maximal degree of monomials in the
423         corresponding component of id, independent of ring ordering
424         (maxdeg of each var is 1).
425         Of type int, if id is of type poly; of type intmat otherwise
426SEE ALSO: maxdeg1
427EXAMPLE: example maxdeg; shows examples
428"
429{
430   //-------- subprocedure to find maximal degree of given component ----------
431   proc findmaxdeg
432   {
433      poly c = #[1];
434      if (c==0) { return(-1); }
435   //--- guess upper 'o' and lower 'u' bound, in case of negative weights -----
436      int d = (deg(c)>=0)*deg(c)-(deg(c)<0)*deg(c);
437      int i = d;
438      while ( c-jet(c,i) != 0 ) { i = 2*(i+1); }
439      int o = i-1;
440      int u = (d != i)*((i div  2)-1);
441   //----------------------- "quick search" for maxdeg ------------------------
442      while ( (c-jet(c,i)==0)*(c-jet(c,i-1)!=0) == 0)
443      {
444         i = (o+1+u) div  2;
445         if (c-jet(c,i)!=0) { u = i+1; }
446         else { o = i-1; }
447      }
448      return(i);
449   }
450//------------------------------ main program ---------------------------------
451   matrix M = matrix(id);
452   int r,c = nrows(M), ncols(M); int i,j;
453   intmat m[r][c];
454   for (i=r; i>0; i--)
455   {
456      for (j=c; j>0; j--) { m[i,j] = findmaxdeg(M[i,j]); }
457   }
458   if (typeof(id)=="poly") { return(m[1,1]); }
459   return(m);
460}
461//-------------------------------- examples -----------------------------------
462example
463{ "EXAMPLE:"; echo = 2;
464   ring r = 0,(x,y,z),wp(1,2,3);
465   poly f = x+y2+z3;
466   deg(f);             //deg; returns weighted degree (in case of 1 block)!
467   maxdeg(f);
468   matrix m[2][2]=f+x10,1,0,f^2;
469   maxdeg(m);
470}
471///////////////////////////////////////////////////////////////////////////////
472
473proc maxdeg1 (id,list #)
474"USAGE:   maxdeg1(id[,v]);  id=poly/ideal/vector/module/matrix, v=intvec
475RETURN:  integer, maximal [weighted] degree of monomials of id independent of
476         ring ordering, maxdeg1 of i-th variable is v[i] (default: v=1..1).
477NOTE:    This proc returns one integer while maxdeg returns, in general,
478         a matrix of integers. For one polynomial and if no intvec v is given
479         maxdeg is faster
480EXAMPLE: example maxdeg1; shows examples
481"
482{
483   //-------- subprocedure to find maximal degree of given component ----------
484   proc findmaxdeg
485   {
486      poly c = #[1];
487      if (c==0) { return(-1); }
488      intvec v = #[2];
489   //--- guess upper 'o' and lower 'u' bound, in case of negative weights -----
490      int d = (deg(c)>=0)*deg(c)-(deg(c)<0)*deg(c);
491      int i = d;
492      if ( c == jet(c,-1,v))      //case: maxdeg is negative
493      {
494         i = -d;
495         while ( c == jet(c,i,v) ) { i = 2*(i-1); }
496         int o = (d != -i)*((i div  2)+2) - 1;
497         int u = i+1;
498         int e = -1;
499      }
500      else                        //case: maxdeg is nonnegative
501      {
502         while ( c != jet(c,i,v) ) { i = 2*(i+1); }
503         int o = i-1;
504         int u = (d != i)*((i div  2)-1);
505         int e = 1;
506      }
507   //----------------------- "quick search" for maxdeg ------------------------
508      while ( ( c==jet(c,i,v) )*( c!=jet(c,i-1,v) ) == 0 )
509      {
510         i = (o+e+u) div  2;
511         if ( c!=jet(c,i,v) ) { u = i+1; }
512         else { o = i-1; }
513      }
514      return(i);
515   }
516//------------------------------ main program ---------------------------------
517   ideal M = simplify(ideal(matrix(id)),8);   //delete scalar multiples from id
518   int c = ncols(M);
519   int i,n;
520   if( size(#)==0 )
521   {
522      int m = maxdeg(M[c]);
523      for (i=c-1; i>0; i--)
524      {
525          n = maxdeg(M[i]);
526          m = (m>=n)*m + (m<n)*n;             //let m be the maximum of m and n
527      }
528   }
529   else
530   {
531      intvec v=#[1];                          //weight vector for the variables
532      int m = findmaxdeg(M[c],v);
533      for (i=c-1; i>0; i--)
534      {
535         n = findmaxdeg(M[i],v);
536         if( n>m ) { m=n; }
537      }
538   }
539   return(m);
540}
541//-------------------------------- examples -----------------------------------
542example
543{ "EXAMPLE:"; echo = 2;
544   ring r = 0,(x,y,z),wp(1,2,3);
545   poly f = x+y2+z3;
546   deg(f);            //deg returns weighted degree (in case of 1 block)!
547   maxdeg1(f);
548   intvec v = ringweights(r);
549   maxdeg1(f,v);                        //weighted maximal degree
550   matrix m[2][2]=f+x10,1,0,f^2;
551   maxdeg1(m,v);                        //absolute weighted maximal degree
552}
553///////////////////////////////////////////////////////////////////////////////
554
555proc mindeg (id)
556"USAGE:   mindeg(id);  id  poly/ideal/vector/module/matrix
557RETURN:  minimal degree/s of monomials of id, independent of ring ordering
558         (mindeg of each variable is 1) of type int if id of type poly, else
559         of type intmat.
560SEE ALSO: mindeg1
561EXAMPLE: example mindeg; shows examples
562"
563{
564   //--------- subprocedure to find minimal degree of given component ---------
565   proc findmindeg
566   {
567      poly c = #[1];
568      if (c==0) { return(-1); }
569   //--- guess upper 'o' and lower 'u' bound, in case of negative weights -----
570      int d = (ord(c)>=0)*ord(c)-(ord(c)<0)*ord(c);
571      int i = d;
572      while ( jet(c,i) == 0 ) { i = 2*(i+1); }
573      int o = i-1;
574      int u = (d != i)*((i div  2)-1);
575   //----------------------- "quick search" for mindeg ------------------------
576      while ( (jet(c,u)==0)*(jet(c,o)!=0) )
577      {
578         i = (o+u) div 2;
579         if (jet(c,i)==0) { u = i+1; }
580         else { o = i-1; }
581      }
582      if (jet(c,u)!=0) { return(u); }
583      else { return(o+1); }
584   }
585//------------------------------ main program ---------------------------------
586   matrix M = matrix(id);
587   int r,c = nrows(M), ncols(M); int i,j;
588   intmat m[r][c];
589   for (i=r; i>0; i--)
590   {
591      for (j=c; j>0; j--) { m[i,j] = findmindeg(M[i,j]); }
592   }
593   if (typeof(id)=="poly") { return(m[1,1]); }
594   return(m);
595}
596//-------------------------------- examples -----------------------------------
597example
598{ "EXAMPLE:"; echo = 2;
599   ring r = 0,(x,y,z),ls;
600   poly f = x5+y2+z3;
601   ord(f);                  // ord returns weighted order of leading term!
602   mindeg(f);               // computes minimal degree
603   matrix m[2][2]=x10,1,0,f^2;
604   mindeg(m);               // computes matrix of minimum degrees
605}
606///////////////////////////////////////////////////////////////////////////////
607
608proc mindeg1 (id, list #)
609"USAGE:   mindeg1(id[,v]);  id=poly/ideal/vector/module/matrix, v=intvec
610RETURN:  integer, minimal [weighted] degree of monomials of id independent of
611         ring ordering, mindeg1 of i-th variable is v[i] (default v=1..1).
612NOTE:    This proc returns one integer while mindeg returns, in general,
613         a matrix of integers. For one polynomial and if no intvec v is given
614         mindeg is faster.
615EXAMPLE: example mindeg1; shows examples
616"
617{
618   //--------- subprocedure to find minimal degree of given component ---------
619   proc findmindeg
620   {
621      poly c = #[1];
622      intvec v = #[2];
623      if (c==0) { return(-1); }
624   //--- guess upper 'o' and lower 'u' bound, in case of negative weights -----
625      int d = (ord(c)>=0)*ord(c)-(ord(c)<0)*ord(c);
626      int i = d;
627      if ( jet(c,-1,v) !=0 )      //case: mindeg is negative
628      {
629         i = -d;
630         while ( jet(c,i,v) != 0 ) { i = 2*(i-1); }
631         int o = (d != -i)*((i div  2)+2) - 1;
632         int u = i+1;
633         int e = -1; i=u;
634      }
635      else                        //case: inded is nonnegative
636      {
637         while ( jet(c,i,v) == 0 ) { i = 2*(i+1); }
638         int o = i-1;
639         int u = (d != i)*((i div  2)-1);
640         int e = 1; i=u;
641      }
642   //----------------------- "quick search" for mindeg ------------------------
643      while ( (jet(c,i-1,v)==0)*(jet(c,i,v)!=0) == 0 )
644      {
645         i = (o+e+u) div  2;
646         if (jet(c,i,v)==0) { u = i+1; }
647         else { o = i-1; }
648      }
649      return(i);
650   }
651//------------------------------ main program ---------------------------------
652   ideal M = simplify(ideal(matrix(id)),8);   //delete scalar multiples from id
653   int c = ncols(M);
654   int i,n;
655   if( size(#)==0 )
656   {
657      int m = mindeg(M[c]);
658      for (i=c-1; i>0; i--)
659      {
660          n = mindeg(M[i]);
661          m = (m<=n)*m + (m>n)*n;             //let m be the maximum of m and n
662      }
663   }
664   else
665   {
666      intvec v=#[1];                          //weight vector for the variables
667      int m = findmindeg(M[c],v);
668      for (i=c-1; i>0; i--)
669      {
670         n = findmindeg(M[i],v);
671         m = (m<=n)*m + (m>n)*n;              //let m be the maximum of m and n
672      }
673   }
674   return(m);
675}
676//-------------------------------- examples -----------------------------------
677example
678{ "EXAMPLE:"; echo = 2;
679   ring r = 0,(x,y,z),ls;
680   poly f = x5+y2+z3;
681   ord(f);                  // ord returns weighted order of leading term!
682   intvec v = 1,-3,2;
683   mindeg1(f,v);            // computes minimal weighted degree
684   matrix m[2][2]=x10,1,0,f^2;
685   mindeg1(m,1..3);         // computes absolute minimum of weighted degrees
686}
687///////////////////////////////////////////////////////////////////////////////
688
689proc normalize (id)
690"USAGE:   normalize(id);  id=poly/vector/ideal/module
691RETURN:  object of same type
692         each element is normalized with leading coefficient equal to 1
693EXAMPLE: example normalize; shows an example
694"
695{
696   return(simplify(id,1));
697}
698//-------------------------------- examples -----------------------------------
699example
700{ "EXAMPLE:"; echo = 2;
701   ring r = 0,(x,y,z),ls;
702   poly f = 2x5+3y2+4z3;
703   normalize(f);
704   module m=[9xy,0,3z3],[4z,6y,2x];
705   normalize(m);
706   ring s = 0,(x,y,z),(c,ls);
707   module m=[9xy,0,3z3],[4z,6y,2x];
708   normalize(m);
709}
710///////////////////////////////////////////////////////////////////////////////
711
712///////////////////////////////////////////////////////////////////////////////
713// Input: <ideal>=<f1,f2,...,fm> and <polynomial> g
714// Question: Does g lie in the radical of <ideal>?
715// Solution: Compute a standard basis G for <f1,f2,...,fm,gz-1> where z is a
716//           new variable. Then g is contained in the radical of <ideal> <=>
717//           1 is generator in G.
718///////////////////////////////////////////////////////////////////////////////
719proc rad_con (poly g,ideal I)
720"USAGE:   rad_con(g,I); g polynomial, I ideal
721RETURN:  1 (TRUE) (type int) if g is contained in the radical of I
722@*       0 (FALSE) (type int) otherwise
723EXAMPLE: example rad_con; shows an example
724"
725{ def br=basering;
726  int n=nvars(br);
727  int dB=degBound;
728  degBound=0;
729  string mp=string(minpoly);
730  if (attrib(br,"global")==1)
731  {
732    execute("ring R=("+charstr(br)+"),(@x(1..n),@z),dp;");
733  }
734  else
735  {
736    execute("ring R=("+charstr(br)+"),(@z,@x(1..n)),(dp(1),"+ordstr(br)+");");
737  }
738  if (mp!="0")
739  {
740    execute("minpoly=number("+mp+");");
741  }
742  ideal irrel=@x(1..n);
743  map f=br,irrel;
744  poly p=f(g);
745  ideal J=f(I),ideal(p*@z-1);
746  J=std(J);
747  degBound=dB;
748  if (J[1]==1)
749  {
750     setring br;
751     return(1);
752  }
753  else
754  {
755     setring br;
756     return(0);
757  }
758}
759example
760{ "EXAMPLE:"; echo=2;
761   ring R=0,(x,y,z),dp;
762   ideal I=x2+y2,z2;
763   poly f=x4+y4;
764   rad_con(f,I);
765   ideal J=x2+y2,z2,x4+y4;
766   poly g=z;
767   rad_con(g,I);
768}
769///////////////////////////////////////////////////////////////////////////////
770
771proc lcm (id, list #)
772"USAGE:   lcm(p[,q]); p int/intvec q a list of integers or
773          p poly/ideal q a list of polynomials
774RETURN:  the least common multiple of p and q:
775@*         - of type int if p is an int/intvec
776@*         - of type poly if p is a poly/ideal
777EXAMPLE: example lcm; shows an example
778"
779{
780  int k,j;
781//------------------------------ integer case --------------------------------
782  if( typeof(id) == "int" or typeof(id) == "intvec" )
783  {
784     intvec i = id;
785     if (size(#)!=0)
786     {
787        for (j = 1; j<=size(#); j++)
788        {
789           if (typeof(#[j]) !="int" and typeof(#[j]) !="intvec")
790           { ERROR("// ** list element must be an integer");}
791           else
792           { i = i,#[j]; }
793        }
794     }
795     int p,q;
796     if( i == 0 )
797     {
798        return(0);
799     }
800     for(j=1;j<=size(i);j++)
801     {
802       if( i[j] != 0 )
803       {
804         p=i[j];
805         break;
806       }
807     }
808     for (k=j+1;k<=size(i);k++)
809     {
810        if( i[k] !=0)
811        {
812           q=gcd(p,i[k]);
813           p=p div q;
814           p=p*i[k];
815        }
816     }
817     if(p <0 )
818     {return(-p);}
819     return(p);
820   }
821
822//----------------------------- polynomial case ------------------------------
823  if( typeof(id) == "poly" or typeof(id) == "ideal" )
824  {
825     ideal i = id;
826     if (size(#)!=0)
827     {
828        for (j = 1; j<=size(#); j++)
829        {
830           if (typeof(#[j]) !="poly" and typeof(#[j]) !="ideal"
831              and typeof(#[j]) !="int" and typeof(#[j]) !="intvec")
832           { ERROR("// ** list element must be a polynomial");}
833           else
834           { i = i,#[j]; }
835        }
836     }
837     poly p,q;
838     i=simplify(i,10);
839     if(size(i) == 0)
840     {
841        return(0);
842     }
843     for(j=1;j<=size(i);j++)
844     {
845       if(maxdeg(i[j])!= 0)
846       {
847         p=i[j];
848         break;
849       }
850     }
851     if(p==0)
852     {
853       return(1);
854     }
855     for (k=j+1;k<=size(i);k++)
856     {
857        if(maxdeg(i[k])!=0)
858        {
859           q=gcd(p,i[k]);
860           if(maxdeg(q)==0)
861           {
862             p=p*i[k];
863           }
864           else
865           {
866             p=p/q;
867             p=p*i[k];
868           }
869        }
870     }
871     return(p);
872   }
873}
874example
875{ "EXAMPLE:"; echo = 2;
876   ring  r = 0,(x,y,z),lp;
877   poly  p = (x+y)*(y+z);
878   poly  q = (z4+2)*(y+z);
879   lcm(p,q);
880   ideal i=p,q,y+z;
881   lcm(i,p);
882   lcm(2,3,6);
883   lcm(2..6);
884}
885
886///////////////////////////////////////////////////////////////////////////////
887
888proc content(f)
889"USAGE:   content(f); f polynomial/vector
890RETURN:  number, the content (greatest common factor of coefficients)
891         of the polynomial/vector f
892SEE ALSO: cleardenom
893EXAMPLE: example content; shows an example
894"
895{
896  if (f==0) { return(number(1)); }
897  return(leadcoef(f)/leadcoef(cleardenom(f)));
898}
899example
900{ "EXAMPLE:"; echo = 2;
901   ring r=0,(x,y,z),(c,lp);
902   content(3x2+18xy-27xyz);
903   vector v=[3x2+18xy-27xyz,15x2+12y4,3];
904   content(v);
905}
906///////////////////////////////////////////////////////////////////////////////
907
908////////////////////////////////////////////////////////////////////////
909// The idea of the procedures mod2id, id2mod and subrInterred is, to
910// perform standard basis computation or interreduction of a submodule
911// of a free module with generators gen(1),...,gen(n) over a ring R
912// in a ring R[t1,...,tn]/<ti*tj> with gen(i) maped to ti
913////////////////////////////////////////////////////////////////////////
914
915proc mod2id(matrix M,intvec vpos)
916"USAGE:   mod2id(M,vpos); M matrix, vpos intvec
917ASSUME:  vpos is an integer vector such that gen(i) corresponds
918         to var(vpos[i]).
919         The basering contains variables var(vpos[i]) which do not occur
920         in M.
921RETURN:  ideal I in which each gen(i) from the module is replaced by
922         var(vpos[i]) and all monomials var(vpos[i])*var(vpos[j]) have
923         been added to the generating set of I.
924NOTE:    This procedure should be used in the following situation:
925         one wants to pass to a ring with new variables, say e(1),..,e(s),
926         which correspond to the components gen(1),..,gen(s) of the
927         module M such that e(i)*e(j)=0 for all i,j.
928         The new ring should already exist and be the current ring
929EXAMPLE: example mod2id; shows an example"
930{
931  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
932//----------------------------------------------------------------------
933// define the ideal generated by the var(vpos[i]) and set up the matrix
934// for the multiplication
935//----------------------------------------------------------------------
936  ideal vars=var(vpos[1]);
937  for(int i=2;i<=size(vpos);i++)
938  {
939    vars=vars,var(vpos[i]);
940  }
941  matrix varm[1][size(vpos)]=vars;
942  if (size(vpos) > nrows(M))
943  {
944    matrix Mt[size(vpos)][ncols(M)];
945    Mt[1..nrows(M),1..ncols(M)]=M;
946    kill M;
947    matrix M=Mt;
948  }
949//----------------------------------------------------------------------
950// define the desired ideal
951//----------------------------------------------------------------------
952  ideal ret=vars^2,varm*M;
953  return(ret);
954}
955example
956{ "EXAMPLE:"; echo=2;
957  ring r=0,(e(1),e(2),x,y,z),(dp(2),ds(3));
958  module mo=x*gen(1)+y*gen(2);
959  intvec iv=2,1;
960  mod2id(mo,iv);
961}
962////////////////////////////////////////////////////////////////////////
963
964proc id2mod(ideal i,intvec vpos)
965"USAGE:   id2mod(I,vpos); I ideal, vpos intvec
966RETURN:  module corresponding to the ideal by replacing var(vpos[i]) by
967         gen(i) and omitting all generators var(vpos[i])*var(vpos[j])
968NOTE:    * This procedure only makes sense if the ideal contains
969           all var(vpos[i])*var(vpos[j]) as monomial generators and
970           all other generators of I are linear combinations of the
971           var(vpos[i]) over the ring in the other variables.
972         * This is the inverse procedure to mod2id and should be applied
973           only to ideals created by mod2id using the same intvec vpos
974           (possibly after a standard basis computation)
975EXAMPLE: example id2mod; shows an example"
976{
977  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
978//---------------------------------------------------------------------
979// Initialization
980//---------------------------------------------------------------------
981  int n=size(i);
982  int v=size(vpos);
983  matrix tempmat;
984  matrix mm[v][n];
985//---------------------------------------------------------------------
986// Conversion
987//---------------------------------------------------------------------
988  for(int j=1;j<=v;j++)
989  {
990    tempmat=coeffs(i,var(vpos[j]));
991    mm[j,1..n]=tempmat[2,1..n];
992  }
993  for(j=1;j<=v;j++)
994  {
995    mm=subst(mm,var(vpos[j]),0);
996  }
997  module ret=simplify(mm,10);
998  return(ret);
999}
1000example
1001{ "EXAMPLE:"; echo=2;
1002  ring r=0,(e(1),e(2),x,y,z),(dp(2),ds(3));
1003  ideal i=e(2)^2,e(1)*e(2),e(1)^2,e(1)*y+e(2)*x;
1004  intvec iv=2,1;
1005  id2mod(i,iv);
1006}
1007///////////////////////////////////////////////////////////////////////
1008
1009proc subrInterred(ideal mon, ideal sm, intvec iv)
1010"USAGE:   subrInterred(mon,sm,iv);
1011         sm:   ideal in a ring r with n + s variables,
1012               e.g. x_1,..,x_n and t_1,..,t_s
1013         mon:  ideal with monomial generators (not divisible by
1014               any of the t_i) such that sm is contained in the module
1015               k[t_1,..,t_s]*mon[1]+..+k[t_1,..,t_s]*mon[size(mon)]
1016         iv:   intvec listing the variables which are supposed to be used
1017               as x_i
1018RETURN:  list l:
1019          l[1]=the monomials from mon in the order used
1020          l[2]=their coefficients after interreduction
1021          l[3]=l[1]*l[2]
1022PURPOSE:  Do interred only w.r.t. a subset of variables.
1023         The procedure returns an interreduced system of generators of
1024         sm considered as a k[t_1,..,t_s]-submodule of the free module
1025         k[t_1,..,t_s]*mon[1]+..+k[t_1,..,t_s]*mon[size(mon)]).
1026EXAMPLE: example subrInterred; shows an example
1027"
1028{
1029  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
1030//-----------------------------------------------------------------------
1031// check that mon is really generated by monomials
1032// and sort its generators with respect to the monomial ordering
1033//-----------------------------------------------------------------------
1034  int err;
1035  for(int i=1;i<=ncols(mon);i++)
1036  {
1037    if ( size(mon[i]) > 1 )
1038    {
1039      err=1;
1040    }
1041  }
1042  if (err==1)
1043  {
1044    ERROR("mon has to be generated by monomials");
1045  }
1046  intvec sv=sortvec(mon);
1047  ideal mons;
1048  for(i=1;i<=size(sv);i++)
1049  {
1050    mons[i]=mon[sv[i]];
1051  }
1052  ideal itemp=maxideal(1);
1053  for(i=1;i<=size(iv);i++)
1054  {
1055    itemp=subst(itemp,var(iv[i]),0);
1056  }
1057  itemp=simplify(itemp,10);
1058  def r=basering;
1059  string tempstr="ring rtemp=" + charstr(basering) + ",(" + string(itemp)
1060                               + "),(C,dp);";
1061//-----------------------------------------------------------------------
1062// express m in terms of the generators of mon and check whether m
1063// can be considered as a submodule of k[t_1,..,t_n]*mon
1064//-----------------------------------------------------------------------
1065  module motemp;
1066  motemp[ncols(sm)]=0;
1067  poly ptemp;
1068  int j;
1069  for(i=1;i<=ncols(mons);i++)
1070  {
1071    for(j=1;j<=ncols(sm);j++)
1072    {
1073      ptemp=sm[j]/mons[i];
1074      motemp[j]=motemp[j]+ptemp*gen(i);
1075    }
1076  }
1077  for(i=1;i<=size(iv);i++)
1078  {
1079    motemp=subst(motemp,var(iv[i]),0);
1080  }
1081  matrix monmat[1][ncols(mons)]=mons;
1082  ideal dummy=monmat*motemp;
1083  for(i=1;i<=size(sm);i++)
1084  {
1085    if(sm[i]-dummy[i]!=0)
1086    {
1087      ERROR("the second argument is not a submodule of the assumed structure");
1088    }
1089  }
1090//----------------------------------------------------------------------
1091// do the interreduction and convert back
1092//----------------------------------------------------------------------
1093  execute(tempstr);
1094  def motemp=imap(r,motemp);
1095  intvec save=option(get);
1096  option(redSB);
1097  motemp=interred(motemp);
1098  option(set,save);
1099  setring r;
1100  kill motemp;
1101  def motemp=imap(rtemp,motemp);
1102  //list ret=monmat,motemp,monmat*motemp;
1103  module motemp2=motemp;
1104  for(i=1;i<=ncols(motemp2);i++)
1105  {
1106    motemp2[i]=cleardenom(motemp2[i]);
1107  }
1108  module motemp3=monmat*motemp;
1109  for(i=1;i<=ncols(motemp3);i++)
1110  {
1111    motemp3[i]=cleardenom(motemp3[i]);
1112  }
1113  list ret=monmat,motemp2,matrix(motemp3);
1114  return(ret);
1115}
1116example
1117{ "EXAMPLE:"; echo=2;
1118  ring r=0,(x,y,z),dp;
1119  ideal i=x^2+x*y^2,x*y+x^2*y,z;
1120  ideal j=x^2+x*y^2,x*y,z;
1121  ideal mon=x^2,z,x*y;
1122  intvec iv=1,3;
1123  subrInterred(mon,i,iv);
1124  subrInterred(mon,j,iv);
1125}
1126
1127///////////////////////////////////////////////////////////////////////////////
1128// moved here from nctools.lib
1129// This procedure calculates the Newton diagram of the polynomial f
1130// The output is a intmat M, each row of M is the exp of a monomial in f
1131////////////////////////////////////////////////////////////////////////
1132proc newtonDiag(poly f)
1133"USAGE:  newtonDiag(f); f a poly
1134RETURN:  intmat
1135PURPOSE: compute the Newton diagram of f
1136NOTE:    each row is the exponent of a monomial of f
1137EXAMPLE: example newtonDiag; shows examples
1138"{
1139  int n = nvars(basering);
1140  intvec N=0;
1141  if ( f != 0 )
1142  {
1143    while ( f != 0 )
1144    {
1145      N = N, leadexp(f);
1146      f = f-lead(f);
1147    }
1148  }
1149  else
1150  {
1151    N=N, leadexp(f);
1152  }
1153  N = N[2..size(N)]; // Deletes the zero added in the definition of T
1154  intmat M=intmat(N,(size(N) div n),n); // Conversion from vector to matrix
1155  return (M);
1156}
1157example
1158{
1159  "EXAMPLE:";echo=2;
1160  ring r = 0,(x,y,z),lp;
1161  poly f = x2y+3xz-5y+3;
1162  newtonDiag(f);
1163}
1164
1165////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.