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

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