# source:git/Singular/LIB/poly.lib@ff8d25

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