# source:git/Singular/LIB/poly.lib@1f9a84

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