source: git/Singular/LIB/standard.lib @ 1d5227

fieker-DuValspielwiese
Last change on this file since 1d5227 was 679bcd, checked in by Hans Schoenemann <hannes@…>, 14 years ago
datetime git-svn-id: file:///usr/local/Singular/svn/trunk@13235 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 70.3 KB
Line 
1//////////////////////////////////////////////////////////////////////////////
2//major revision Jan/Feb. 2007, GMG (groebner with several options)
3//Change of default methods in groebner (June 2008, GMG)
4//stdhilb can be called with std or slimgb (Jan 2008, GMG)
5//### Todo: im lokalen Fall die Hilbert-Samuel Funktion verwenden
6//////////////////////////////////////////////////////////////////////////////
7version="$Id$";
8category="Miscellaneous";
9info="
10LIBRARY: standard.lib   Procedures which are always loaded at Start-up
11
12PROCEDURES:
13 stdfglm(ideal[,ord])   standard basis of ideal via fglm [and ordering ord]
14 stdhilb(ideal[,h])     Hilbert driven Groebner basis of ideal
15 groebner(ideal,...)    standard basis using a heuristically chosen method
16 res(ideal/module,[i])  free resolution of ideal or module
17 sprintf(fmt,...)       returns fomatted string
18 fprintf(link,fmt,..)   writes formatted string to link
19 printf(fmt,...)        displays formatted string
20 weightKB(stc,dd,vl)    degree dd part of a kbase w.r.t. some weigths
21 qslimgb(i)             computes a standard basis with slimgb in a qring
22 par2varRing([i])       create a ring making pars to vars, together with i
23 datetime()             return date and time as a string
24
25";
26//AUXILIARY PROCEDURES:
27// hilbRing([i])          ring for computing the (weighted) hilbert series
28// quotientList(L,...)    ringlist for creating a correct quotient ring
29
30//////////////////////////////////////////////////////////////////////////////
31
32proc stdfglm (ideal i, list #)
33"SYNTAX: @code{stdfglm (} ideal_expression @code{)} @*
34         @code{stdfglm (} ideal_expression@code{,} string_expression @code{)}
35TYPE:    ideal
36PURPOSE: computes the standard basis of the ideal in the basering
37         via @code{fglm} from the ordering given as the second argument
38         to the ordering of the basering. If no second argument is given,
39         \"dp\" is used. The standard basis for the given ordering (resp. for
40         \"dp\") is computed via the command groebner except if a further
41         argument \"std\" or \"slimgb\" is given in which case std resp.
42         slimgb is used.
43SEE ALSO: fglm, groebner, std, slimgb, stdhilb
44KEYWORDS: fglm
45EXAMPLE: example stdfglm; shows an example"
46{
47  string os;
48  int s = size(#);
49  def P= basering;
50  string algorithm;
51  int ii;
52  for( ii=1; ii<=s; ii++)
53  {
54    if ( typeof(#[ii])== "string" )
55    {
56      if ( #[ii]=="std" || #[ii]=="slimgb" )
57      {
58        algorithm =  #[ii];
59        # = delete(#,ii);
60        s--;
61        ii--;
62      }
63    }
64  }
65
66  if((s > 0) && (typeof(#[1]) == "string"))
67  {
68    os = #[1];
69    ideal Qideal = ideal(P);
70    int sQ = size(Qideal);
71    int sM = size(minpoly);
72    if ( sM!=0 )
73    {
74      string mpoly = string(minpoly);
75    }
76    if (sQ!=0 )
77    {
78      execute("ring Rfglm=("+charstr(P)+"),("+varstr(P)+"),"+os+";");
79      ideal Qideal = fetch(P,Qideal);
80      qring Pfglm = groebner(Qideal,"std","slimgb");
81    }
82    else
83    {
84      execute("ring Pfglm=("+charstr(P)+"),("+varstr(P)+"),"+os+";");
85    }
86    if ( sM!=0 )
87    {
88      execute("minpoly="+mpoly+";");
89    }
90  }
91  else
92  {
93    list BRlist = ringlist(P);
94    int nvarP = nvars(P);
95    intvec w;                       //for ringweights of basering P
96    int k;
97    for(k=1;  k <= nvarP; k++)
98    {
99      w[k]=deg(var(k));
100    }
101
102    BRlist[3] = list();
103    if( s==0 or (typeof(#[1]) != "string") )
104    {
105      if( w==1 )
106      {
107        BRlist[3][1]=list("dp",w);
108      }
109      else
110      {
111        BRlist[3][1]=list("wp",w);
112      }
113      BRlist[3][2]=list("C",intvec(0));
114      def Pfglm = ring(quotientList(BRlist));
115      setring Pfglm;
116    }
117  }
118  ideal i = fetch(P,i);
119
120  intvec opt = option(get);            //save options
121  option(redSB);
122  if (size(algorithm) > 0)
123  {
124    i = groebner(i,algorithm);
125  }
126  else
127  {
128    i = groebner(i);
129  }
130  option(set,opt);
131  setring P;
132  return (fglm(Pfglm,i));
133}
134example
135{ "EXAMPLE:"; echo = 2;
136   ring r  = 0,(x,y,z),lp;
137   ideal i = y3+x2,x2y+x2,x3-x2,z4-x2-y;
138   stdfglm(i);                   //uses fglm from "dp" (with groebner) to "lp"
139   stdfglm(i,"std");             //uses fglm from "dp" (with std) to "lp"
140
141   ring s  = (0,x),(y,z,u,v),lp;
142   minpoly = x2+1;
143   ideal i = u5-v4,zv-u2,zu3-v3,z2u-v2,z3-uv,yv-zu,yu-z2,yz-v,y2-u,u-xy2;
144   weight(i);
145   stdfglm(i,"(a(2,3,4,5),dp)"); //uses fglm from "(a(2,3,4,5),dp)" to "lp"
146}
147
148/////////////////////////////////////////////////////////////////////////////
149
150proc stdhilb(i,list #)
151"SYNTAX: @code{stdhilb (} ideal_expression @code{)} @*
152         @code{stdhilb (} module_expression @code{)} @*
153         @code{stdhilb (} ideal_expression, intvec_expression @code{)}@*
154         @code{stdhilb (} module_expression, intvec_expression @code{)}@*
155         @code{stdhilb (} ideal_expression@code{,} list of string_expressions,
156               and intvec_expression @code{)} @*
157TYPE:    type of the first argument
158PURPOSE: Compute a Groebner basis of the ideal/module in the basering by
159         using the Hilbert driven Groebner basis algorithm.
160         If an argument of type string, stating @code{\"std\"} resp. @code{\"slimgb\"},
161         is given, the standard basis computation uses @code{std} or
162         @code{slimgb}, otherwise a heuristically chosen method (default)@*
163         If an optional second argument w of type intvec is given, w is used
164         as variable weights. If w is not given, it is computed as w[i] =
165         deg(var(i)). If the ideal is homogeneous w.r.t. w then the
166         Hilbert series is computed w.r.t. to these weights.
167THEORY:  If the ideal is not homogeneous compute first a Groebner basis
168         of the homogenization [w.r.t. the weights w] of the ideal/module,
169         then the Hilbert function and, finally, a Groebner basis in the
170         original ring by using the computed Hilbert function. If the given
171         w does not coincide with the variable weights of the basering, the
172         result may not be a groebner basis in the original ring.
173NOTE:    'Homogeneous' means weighted homogeneous with respect to the weights
174         w[i] of the variables var(i) of the basering. Parameters are not
175         converted to variables.
176SEE ALSO: stdfglm, std, slimgb, groebner
177KEYWORDS: Hilbert function
178EXAMPLE: example stdhilb;  shows an example"
179{
180
181//--------------------- save data from basering --------------------------
182  def P=basering;
183  int nr;
184  if (typeof(i)=="ideal") { nr=1;}
185  else                    { nr= nrows(i); }    //nr=1 if i is an ideal
186  ideal Qideal = ideal(P);      //defining the quotient ideal if P is a qring
187  int was_qring;                //remembers if basering was a qring
188  int is_homog =homog(i);       //check for homogeneity of i and Qideal
189  if (size(Qideal) > 0)
190  {
191     was_qring = 1;
192  }
193
194  // save ordering of basering P for later use
195  list ord_P = ringlist(P)[3];     //ordering of basering in ringlist
196  string ordstr_P = ordstr(P);     //ordering of basering as string
197  int nvarP = nvars(P);
198
199  //save options:
200  intvec gopt = option(get);
201  int p_opt;
202  string s_opt = option();
203  if (find(s_opt, "prot"))  { p_opt = 1; }
204
205//-------------------- check the given method and weights ---------------------
206//Note: stdhilb is used in elim where it is applied to an elimination ordering
207//a(1..1,0..0),wp(w). In such a ring deg(var(k)=0 for all vars corresponding to
208//0 in a(1..1,0..0), hence we cannot identify w via w[k] = deg(var(k));
209//Therefore hilbstd has the option to give ringweights.
210
211  int k;
212  string method;
213  for (k=1; k<=size(#); k++)
214  {
215    if (typeof(#[k]) == "intvec")
216    {
217        intvec w = #[k];         //given ringweights of basering P
218    }
219    if (typeof(#[k]) == "string")
220    {
221      method = method + "," + #[k];
222    }
223  }
224
225  if ( defined(w)!=voice )
226  {
227    intvec w;
228    for(k=nvarP;  k>=1; k--)
229    {
230       w[k] = deg(var(k));     //compute ring weights
231    }
232  }
233
234
235  if (npars(P) > 0)             //clear denominators of parameters
236  {
237    for( k=ncols(i); k>0; k-- )
238    {
239      i[k]=cleardenom(i[k]);
240    }
241  }
242
243//---------- exclude cases to which stdhilb should no be applied  ----------
244//Note that quotient ideal of qring must be homogeneous too
245
246  int neg=1-attrib (P,"global");
247
248  if( //find(ordstr_P,"s") ||// covered by neg
249     find(ordstr_P,"M") || neg )
250  {
251    // if( defined(hi) && is_homog )
252    // {
253    //  if (p_opt){"std with given Hilbert function in basering";}
254    //  return( std(i,hi,w) );
255    //### here we would need Hibert-Samuel function
256    // }
257
258    if (p_opt)
259    {"//-- stdhilb not implemented, we use std in ring:"; string(P);}
260    return( std(i) );
261  }
262
263//------------------------ change to hilbRing ----------------------------
264//The ground field of P and Philb coincide, Philb has an extra variable
265//@ or @(k). Philb is no qring and the predefined ideal/module Id(1) in
266//Philb is homogeneous (it is the homogenized i w.r.t. @ or @(k))
267//Parameters of P are not converted in Philb
268//Philb has only 1 block dp or wp(w)
269
270  list hiRi = hilbRing(i,w);
271  intvec W = hiRi[2];
272  def Philb = hiRi[1];
273  setring Philb;
274
275//-------- compute Hilbert series of homogenized ideal in Philb ---------
276//There are three cases
277
278  string algorithm;       //possibilities: std, slimgb, stdorslimgb
279  //define algorithm:
280  if( find(method,"std") && !find(method,"slimgb") )
281  {
282    algorithm = "std";
283  }
284  if( find(method,"slimgb") && !find(method,"std") )
285  {
286    algorithm = "slimgb";
287  }
288  if( find(method,"std") && find(method,"slimgb") ||
289    (!find(method,"std") && !find(method,"slimgb")) )
290  {
291    algorithm = "stdorslimgb";
292  }
293
294//### geaendert Dez08: es wird std(Id(1)) statt Id(1) aus Philb nach Phelp
295// weitergegeben fuer hilbertgetriebenen std
296
297  if (( algorithm=="std" || ( algorithm=="stdorslimgb" && char(P)>0 ) )
298  && (defined(hi)!=voice))
299  {
300    if (p_opt) {"compute hilbert series with std in ring " + string(Philb);
301                "weights used for hilbert series:",W;}
302    Id(1) = std(Id(1));
303    intvec hi = hilb( Id(1),1,W );
304  }
305  if (( algorithm=="slimgb" || ( algorithm=="stdorslimgb" && char(P)==0 ) )
306  && (defined(hi)!=voice))
307  {
308    if (p_opt) {"compute hilbert series with slimgb in ring " + string(Philb);
309                "weights used for hilbert series:",W;}
310    Id(1) = qslimgb(Id(1));
311    intvec hi = hilb( Id(1),1,W );
312  }
313
314  //-------------- we need another intermediate ring Phelp ----------------
315  //In Phelp we change only the ordering from Philb (otherwise it coincides
316  //with Philb). Phelp has in addition to P an extra homogenizing variable
317  //with name @ (resp. @(i) if @ and @(1), ..., @(i-1) are defined) with
318  //ordering an extra last block dp(1).
319  //Phelp has the same ordering as P on common variables. In Phelp
320  //a quotient ideal from P is added to the input
321
322  list BRlist = ringlist(Philb);
323  BRlist[3] = list();
324  int so = size(ord_P);
325  if( ord_P[so][1] =="c" || ord_P[so][1] =="C" )
326  {
327    list moduleord = ord_P[so];
328    so = so-1;
329  }
330  for (k=1; k<=so; k++)
331  {
332    BRlist[3][k] = ord_P[k];
333  }
334
335  BRlist[3][so+1] = list("dp",1);
336  w = w,1;
337
338  if( defined(moduleord)==voice )
339  {
340    BRlist[3][so+2] = moduleord;
341  }
342
343//--- change to extended ring Phelp and compute std with hilbert series ----
344  def Phelp = ring(quotientList(BRlist));
345  setring Phelp;
346  def i = imap(Philb, Id(1));
347  kill Philb;
348
349  // compute std with Hilbert series
350  option(redThrough);
351  if (w == 1)
352  {
353    if (p_opt){ "std with hilb in " + string(Phelp);}
354    i = std(i, hi);
355  }
356  else
357  {
358    if(p_opt){"std with weighted hilb in "+string(Phelp);}
359    i = std(i, hi, w);
360  }
361
362//-------------------- go back to original ring ---------------------------
363//The main computation is done. Do not forget to simplfy before maping.
364
365  // subst 1 for homogenizing var
366  if ( p_opt ) { "dehomogenization"; }
367  i = subst(i, var(nvars(basering)), 1);
368
369  if (p_opt) { "simplification"; }
370  i= simplify(i,34);
371
372  setring P;
373  if (p_opt) { "imap to ring "+string(P); }
374  i = imap(Phelp,i);
375  kill Phelp;
376  if( was_qring )
377  {
378    i = NF(i,std(0));
379  }
380  i = simplify(i,34);
381  // compute reduced SB
382  if (find(s_opt, "redSB") > 0)
383  {
384    if (p_opt) { "//interreduction"; }
385    i=interred(i);
386  }
387  attrib(i, "isSB", 1);
388  option(set,gopt);
389  return (i);
390}
391example
392{ "EXAMPLE:"; echo = 2;
393   ring  r = 0,(x,y,z),lp;
394   ideal i = y3+x2,x2y+x2z2,x3-z9,z4-y2-xz;
395   ideal j = stdhilb(i); j;
396
397   ring  r1 = 0,(x,y,z),wp(3,2,1);
398   ideal  i = y3+x2,x2y+x2z2,x3-z9,z4-y2-xz;  //ideal is homogeneous
399   ideal j = stdhilb(i,"std"); j;
400   //this is equivalent to:
401   intvec v = hilb(std(i),1);
402   ideal j1 = std(i,v,intvec(3,2,1)); j1;
403
404   size(NF(j,j1))+size(NF(j1,j));            //j and j1 define the same ideal
405}
406
407///////////////////////////////////////////////////////////////////////////////
408proc quotientList (list RL, list #)
409"SYNTAX: @code{quotientList (} list_expression @code{)} @*
410         @code{quotientList (} list_expression @code{,} string_expression@code{)}
411TYPE:    list
412PURPOSE: define a ringlist, say QL, of the first argument, say RL, which is
413         assumed to be the ringlist of a qring, but where the quotient ideal
414         RL[4] is not a standard basis with respect to the given monomial
415         order in RL[3]. Then QL will be obtained from RL just by replacing
416         RL[4] by a standard of it with respect to this order. RL itself
417         will be returnd if size(RL[4]) <= 1 (in which case it is known to be
418         a standard basis w.r.t. any ordering) or if a second argument
419         \"isSB\" of type string is given.
420NOTE:    the command ring(quotientList(RL)) defines a quotient ring correctly
421         and should be used instead of ring(RL) if the quotient ideal RL[4]
422         is not (or not known to be) a standard basis with respect to the
423         monomial ordering specified in RL[3].
424SEE ALSO: ringlist, ring
425EXAMPLE: example quotientList; shows an example"
426{
427  def P = basering;
428  if( size(#) > 0 )
429  {
430    if ( #[1] == "isSB")
431    {
432      return (RL);
433    }
434  }
435  ideal Qideal  = RL[4];  //##Achtung: falls basering Nullteiler hat, kann
436                           //die SB eines Elements mehrere Elemente enthalten
437  if( size(Qideal) <= 0)
438  {
439    return (RL);
440  }
441
442  RL[4] = ideal(0);
443  def Phelp = ring(RL);
444  setring Phelp;
445  ideal Qideal = groebner(fetch(P,Qideal));
446  setring P;
447  RL[4]=fetch(Phelp,Qideal);
448  return (RL);
449}
450example
451{ "EXAMPLE:"; echo = 2;
452   ring P = 0,(y,z,u,v),lp;
453   ideal i = y+u2+uv3, z+uv3;            //i is an lp-SB but not a dp_SB
454   qring Q = std(i);
455   list LQ = ringlist(Q);
456   LQ[3][1][1]="dp";
457   def Q1 = ring(quotientList(LQ));
458   setring Q1;
459   Q1;
460
461   setring Q;
462   ideal q1 = uv3+z, u2+y-z, yv3-zv3-zu; //q1 is a dp-standard basis
463   LQ[4] = q1;
464   def Q2 = ring(quotientList(LQ,"isSB"));
465   setring Q2;
466   Q2;
467}
468
469///////////////////////////////////////////////////////////////////////////////
470proc par2varRing (list #)
471"USAGE:   par2varRing([l]); l list of ideals/modules [default:l=empty list]
472RETURN:  list, say L, with L[1] a ring where the parameters of the
473         basering have been converted to an additional last block of
474         variables, all of weight 1, and ordering dp.
475         If a list l with l[i] an ideal/module is given, then
476         l[i] + minpoly*freemodule(nrows(l[i])) is mapped to an ideal/module
477         in L[1] with name Id(i).
478         If the basering has no parameters then L[1] is the basering.
479EXAMPLE: example par2varRing; shows an example"
480{
481  def P = basering;
482  int npar = npars(P);  //number of parameters
483  int s = size(#);
484  int ii;
485  if ( npar == 0)
486  {
487    dbprint(printlevel-voice+3,"// ** no parameters, ring was not changed");
488    for( ii = 1; ii <= s; ii++)
489    {
490      def Id(ii) = #[ii];
491      export (Id(ii));
492    }
493    return(list(P));
494  }
495
496  list rlist = ringlist(P);
497  list parlist = rlist[1];
498  rlist[1] = parlist[1];
499  poly Minpoly = minpoly;     //check for minpoly:
500  int sm = size(Minpoly);
501
502  //now create new ring
503  for( ii = 1; ii <= s; ii++)
504  {
505    def Id(ii) = #[ii];
506  }
507  int nvar = size(rlist[2]);
508  int nblock = size(rlist[3]);
509  int k;
510  for (k=1; k<=npar; k++)
511  {
512    rlist[2][nvar+k] = parlist[2][k];   //change variable list
513  }
514
515  //converted parameters get one block dp. If module ordering was in front
516  //it stays in front, otherwise it will be moved to the end
517  intvec OW = 1:npar;
518  if( rlist[3][nblock][1] =="c" || rlist[3][nblock][1] =="C" )
519  {
520    rlist[3][nblock+1] = rlist[3][nblock];
521    rlist[3][nblock] = list("dp",OW);
522  }
523  else
524  {
525    rlist[3][nblock+1] = list("dp",OW);
526  }
527
528  def Ppar2var = ring(quotientList(rlist));
529  setring Ppar2var;
530  if ( sm == 0 )
531  {
532    for( ii = 1; ii <= s; ii++)
533    {
534      def Id(ii) = imap(P,Id(ii));
535      export (Id(ii));
536    }
537  }
538  else
539  {
540    if( find(option(),"prot") ){"//add minpoly to input";}
541    poly Minpoly = imap(P,Minpoly);
542    for( ii = 1; ii <= s; ii++)
543    {
544      def Id(ii) = imap(P,Id(ii));
545      if (typeof(Id(ii))=="module")
546      {
547        Id(ii) = Id(ii),Minpoly*freemodule(nrows(Id(ii)));
548      }
549      else
550      {
551        Id(ii) = Id(ii),Minpoly;
552      }
553      export (Id(ii));
554    }
555  }
556  list Lpar2var = Ppar2var;
557  return(Lpar2var);
558}
559example
560{ "EXAMPLE:"; echo = 2;
561   ring R = (0,x),(y,z,u,v),lp;
562   minpoly = x2+1;
563   ideal i = x3,x2+y+z+u+v,xyzuv-1; i;
564   def P = par2varRing(i)[1]; P;
565   setring(P);
566   Id(1);
567
568   setring R;
569   module m = x3*[1,1,1], (xyzuv-1)*[1,0,1];
570   def Q = par2varRing(m)[1]; Q;
571   setring(Q);
572   print(Id(1));
573}
574
575//////////////////////////////////////////////////////////////////////////////
576proc hilbRing ( list # )
577"USAGE:   hilbRing([w,l]); w = intvec, l = list of ideals/modules
578RETURN:  list, say L: L[1] is a ring and L[2] an intvec
579         L[1] is a ring whith an extra homogenizing variable with name @,
580         resp. @(i) if @ and @(1), ..., @(i-1) are defined.
581         The monomial ordering of L[1] is consists of 1 block: dp if the
582         weights of the variables of the basering, say R, are all 1, resp.
583         wp(w,1) wehre w is either given or the intvec of weights of the
584         variables of R, i.e. w[k]=deg(var(k)).
585         If R is a quotient ring P/Q, then L[1] is not a quotient ring but
586         contains the ideal @Qidealhilb@, the homogenized ideal Q of P.
587         (Parameters of R are not touched).
588         If a list l is given with l[i] an ideal/module, then l[i] is mapped
589         to Id(i), the homogenized l[i]+Q*freemodule(nrows(l[i]) in L[1]
590         (Id(i) = l[i] if l[i] is already homogeneous).
591         L[2] is the intvec (w,1).
592PURPOSE: Prepare a ring for computing the (weighted) hilbert series of
593         an ideal/module with an easy monomial ordering.
594NOTE:    For this purpose we need w[k]=deg(var(k)). However, if the ordering
595         contains an extra weight vector a(v,0..0)) deg(var(k)) returns 0 for
596         k being an index which is 0 in a. Therefore we must compute w
597         beforehand and give it to hilbRing.
598EXAMPLE: example hilbRing; shows an example
599"
600{
601  def P = basering;
602  ideal Qideal = ideal(P);    //defining the quotient ideal if P is a qring
603  if( size(Qideal) != 0 )
604  {
605    int is_qring =1;
606  }
607  list BRlist = ringlist(P);
608  BRlist[4] = ideal(0);      //kill quotient ideal in BRlist
609
610  int nvarP = nvars(P);
611  int s = size(#);
612  int k;
613
614  for(k = 1; k <= s; k++)
615  {
616    if ( typeof(#[k]) == "intvec" )
617    {
618       intvec w = #[k];      //given weights for the variables
619       # = delete (#,k);
620    }
621  }
622
623  s = size(#);
624  for(k = 1; k <= s; k++)
625  {
626     def Id(k) = #[k];
627     int nr(k) = nrows(Id(k));
628  }
629
630  if ( defined(w)!=voice )
631  {
632    intvec w;                   //for ringweights of basering P
633    for(k=1;  k<=nvarP; k++)
634    {
635      w[k]=deg(var(k));        //degree of kth variable
636    }
637  }
638  //--------------------- a homogenizing variable is added ------------------
639  // call it @, resp. @(k) if @(1),...,@(k-1) are defined
640  string homvar;
641  if ( defined(@)==0 )
642  {
643    homvar = "@";
644  }
645  else
646  {
647    k=1;
648    while( defined(@(k)) != 0 )
649    {
650      k++;
651    }
652    homvar = "@("+string(k)+")";
653  }
654  BRlist[2][nvarP+1] = homvar;
655  w[nvarP +1]=1;
656
657  //ordering is set to (dp,C) if weights of all variables are 1
658  //resp. to (wp(w,1),C) where w are the ringweights of basering P
659  //homogenizing var gets weight 1:
660
661  BRlist[3] = list();
662  BRlist[3][2]=list("C",intvec(0));  //put module ordering always last
663  if(w==1)
664  {
665    BRlist[3][1]=list("dp",w);
666  }
667  else
668  {
669    BRlist[3][1]=list("wp",w);
670  }
671
672  //-------------- change ring and get ideal from previous ring ---------------
673  def Philb = ring(quotientList(BRlist));
674  kill BRlist;
675  setring Philb;
676  if( defined(is_qring)==voice )
677  {
678    ideal @Qidealhilb@ =  imap(P,Qideal);
679    if ( ! homog(@Qidealhilb@) )
680    {
681       @Qidealhilb@ =  homog( @Qidealhilb@, `homvar` );
682    }
683    export(@Qidealhilb@);
684
685    if( find(option(),"prot") ){"add quotient ideal to input";}
686
687    for(k = 1; k <= s; k++)
688    { //homogenize if necessary
689      def Id(k) =  imap(P,Id(k));
690      if ( ! homog(Id(k)) )
691      {
692         Id(k) =  homog( imap(P,Id(k)), `homvar` );
693      }
694      if (typeof(Id(k))=="module")
695      {
696        Id(k) =  Id(k),@Qidealhilb@*freemodule(nr(k)) ;
697      }
698      else
699      {
700        Id(k) =  Id(k),@Qidealhilb@ ;
701      }
702      export(Id(k));
703    }
704  }
705  else
706  {
707    for(k = 1; k <= s; k++)
708    { //homogenize if  necessary
709      def Id(k) =  imap(P,Id(k));
710      if ( ! homog(Id(k)) )
711      {
712         Id(k) =  homog( imap(P,Id(k)), `homvar` );
713      }
714      export(Id(k));
715    }
716  }
717  list Lhilb = Philb,w;
718  return(Lhilb);
719}
720example
721{ "EXAMPLE:"; echo = 2;
722   ring R = 0,(x,y,z,u,v),lp;
723   ideal i = x+y2+z3,xy+xv+yz+zu+uv,xyzuv-1;
724   intvec w = 6,3,2,1,1;
725   hilbRing(i,w);
726   def P = hilbRing(w,i)[1];
727   setring P;
728   Id(1);
729   hilb(std(Id(1)),1);
730
731   ring S =  0,(x,y,z,u,v),lp;
732   qring T = std(x+y2+z3);
733   ideal i = xy+xv+yz+zu+uv,xyzuv-v5;
734   module m = i*[0,1,1] + (xyzuv-v5)*[1,1,0];
735   def Q = hilbRing(m)[1];  Q;
736   setring Q;
737   print(Id(1));
738}
739
740//////////////////////////////////////////////////////////////////////////////
741proc qslimgb (i)
742"USAGE:   qslimgb(i); i ideal or module
743RETURN:  same type as input, a standard basis of i computed with slimgb
744NOTE:    Only as long as slimgb does not know qrings qslimgb should be used
745         in case the basering is (possibly) a quotient ring.
746         The quotient ideal is added to the input and slimgb is applied.
747EXAMPLE: example qslimgb; shows an example"
748{
749  def P = basering;
750  ideal Qideal = ideal(P);      //defining the quotient ideal if P is a qring
751  int p_opt;
752  if( find(option(),"prot") )
753  {
754    p_opt=1;
755  }
756  if (size(Qideal) == 0)
757  {
758    if (p_opt) { "slimgb in ring " + string(P); }
759    return(slimgb(i));
760  }
761
762  //case of a qring; since slimgb does not know qrings we
763  //delete the quotient ideal and add it to i
764
765  list BRlist = ringlist(P);
766  BRlist[4] = ideal(0);
767  def Phelp = ring(BRlist);
768  kill BRlist;
769  setring Phelp;
770  // module case:
771  def iq = imap(P,i);
772  iq = iq, imap(P,Qideal)*freemodule(nrows(iq));
773  if (p_opt)
774  {
775    "slimgb in ring " + string(Phelp);
776    "(with quotient ideal added to input)";
777  }
778  iq = slimgb(iq);
779
780  setring P;
781  if (p_opt) { "//imap to original ring"; }
782  i = imap(Phelp,iq);
783  kill Phelp;
784
785  if (find(option(),"redSB") > 0)
786  {
787    if (p_opt) { "//interreduction"; }
788    i=reduce(i,std(0));
789    i=interred(i);
790  }
791  attrib(i, "isSB", 1);
792  return (i);
793}
794example
795{ "EXAMPLE:"; echo = 2;
796   ring R  = (0,v),(x,y,z,u),dp;
797   qring Q = std(x2-y3);
798   ideal i = x+y2,xy+yz+zu+u*v,xyzu*v-1;
799   ideal j = qslimgb(i); j;
800
801   module m = [x+y2,1,0], [1,1,x2+y2+xyz];
802   print(qslimgb(m));
803}
804
805//////////////////////////////////////////////////////////////////////////////
806proc groebner(def i_par, list #)
807"SYNTAX: @code{groebner (} ideal_expression @code{)} @*
808         @code{groebner (} module_expression @code{)} @*
809         @code{groebner (} ideal_expression@code{,} int_expression @code{)} @*
810         @code{groebner (} module_expression@code{,} int_expression @code{)} @*
811         @code{groebner (} ideal_expression@code{,} list of string_expressions
812               @code{)} @*
813         @code{groebner (} ideal_expression@code{,} list of string_expressions
814               and int_expression @code{)} @*
815         @code{groebner (} ideal_expression@code{,} int_expression @code{)}
816TYPE:    type of the first argument
817PURPOSE: computes a standard basis of the first argument @code{I}
818         (ideal or module) by a heuristically chosen method (default)
819         or by a method specified by further arguments of type string.
820         Possible methods are:  @*
821         - the direct methods @code{\"std\"} or @code{\"slimgb\"} without
822           conversion, @*
823         - conversion methods @code{\"hilb\"} or @code{\"fglm\"} where
824           a Groebner basis is first computed with an \"easy\" ordering
825           and then converted to the ordering of the basering by the
826           Hilbert driven Groebner basis computation or by linear algebra.
827           The actual computation of the Groebner basis can be
828           specified by @code{\"std\"} or by @code{\"slimgb\"}
829           (not for all orderings implemented). @*
830         A further string @code{\"par2var\"} converts parameters to an extra
831         block of variables before a Groebner basis computation (and
832         afterwards back).
833         @code{option(prot)} informs about the chosen method.
834NOTE:    If an additional argument, say @code{wait}, of type int is given,
835         then the computation runs for at most @code{wait} seconds.
836         That is, if no result could be computed in @code{wait} seconds,
837         then the computation is interrupted, 0 is returned, a warning
838         message is displayed, and the global variable
839         @code{Standard::groebner_error} is defined.
840         This feature uses MP and hence it is available on UNIX platforms, only.
841HINT:    Since there exists no uniform best method for computing standard
842         bases, and since the difference in performance of a method on
843         different examples can be huge, it is recommended to test, for hard
844         examples, first various methods on a simplified example (e.g. use
845         characteristic 32003 instead of 0 or substitute a subset of
846         parameters/variables by integers, etc.). @*
847SEE ALSO: stdhilb, stdfglm, std, slimgb
848KEYWORDS: time limit on computations; MP, groebner basis computations
849EXAMPLE: example groebner;  shows an example"
850
851{
852//Vorgabe einer Teilmenge aus {hilb,fglm,par2var,std,slimgb}
853//V1: Erste Einstellungen (Jan 2007)
854//V2: Aktuelle Aenderungen (Juni 2008)
855//---------------------------------
856//0. Immer Aufruf von std unabhaengig von der Vorgabe:
857//   gemischte Ordnungen, extra Gewichtsvektor, Matrix Ordnungen
858//   ### Todo: extra Gewichtsvektor sollte nicht immer mit std wirken,
859//   sondern z.B. mit "hilb" arbeiten koennen
860//   ### Todo: es sollte ein Gewichtsvektor mitgegeben werden koennen (oder
861//   berechnet werden), z.B. groebner(I,"hilb",w) oder groebner(I,"withWeights")
862//   wie bei elim in elim.lib
863
864//1. Keine Vorgabe: es wirkt die aktuelle Heuristk:
865//   - Char = p: std
866//V1 - Char = 0: slimgb (im qring wird Quotientenideal zum Input addiert)
867//V2 - Char = 0: std
868//   - 1-Block-Ordnungen/non-commutative: direkt Aufruf von std oder slimgb
869//   - Komplizierte Ordnungen (lp oder > 1 Block): hilb
870//V1 - Parameter werden grundsaetzlich nicht in Variable umgewandelt
871//V2 - Mehr als ein Parmeter wird zu Variable konvertiert
872//   - fglm is keine Heuristik, da sonst vorher dim==0 peprueft werden muss
873
874//2. Vorgabe aus {std,slimgb}: es wird wo immer moeglich das Angegebene
875//   gewaehlt (da slimgb keine Hilbertfunktion kennt, wird std verwendet).
876//   Bei slimgb im qring, wird das Quotientenideal zum Ideal addiert.
877//   Bei Angabe von std zusammen mit slimgb (aequivalent zur Angabe von
878//   keinem von beidem) wirkt obige Heuristik.
879
880//3. Nichtleere Vorgabe aus {hilb,fglm,std,slimgb}:
881//   es wird nur das Angegebene und Moegliche sowie das Notwendige verwendet
882//   und bei Wahlmoeglickeit je nach Heuristik.
883//   Z.B. Vorgabe von {hilb} ist aequivalent zu {hilb,std,slimgb} und es wird
884//   hilb und nach Heuristik std oder slimgb verwendet,
885//   (V1: aber nicht par2var)
886//   bei Vorgabe von {hilb,slimgb} wird hilb und wo moeglich slimgb verwendet.
887
888//4. Bei Vorgabe von {par2var} wird par2var immer mit hilb und nach Heuristik
889//   std oder slimgb verwendet. Zu Variablen konvertierte Parameter haben
890//   extra letzten Block und Gewichte 1.
891
892  def P=basering;
893  if ((typeof(i_par)=="vector")||(typeof(i_par)=="module")||(typeof(i_par)=="matrix")) {module i=i_par;}
894  else {ideal i=i_par; } // int, poly, number, ideal
895  kill i_par;
896
897//----------------------- save the given method ---------------------------
898  string method;                //all given methods as a coma separated string
899  list Method;                  //all given methods as a list
900  int k;
901  for (k=1; k<=size(#); k++)
902  {
903     if (typeof(#[k]) == "int")
904     {
905       int wait = #[k];
906     }
907     if (typeof(#[k]) == "string")
908     {
909       method = method + "," + #[k];
910       Method = Method + list(#[k]);
911     }
912  }
913
914 //======= we have an argument of type int -- try to use MPfork links =======
915  if ( defined(wait) == voice )
916  {
917    if ( system("with", "MP") )
918    {
919        int j = 10;
920        string bs = nameof(basering);
921        link l_fork = "MPtcp:fork";
922        open(l_fork);
923        write(l_fork, quote(system("pid")));
924        int pid = read(l_fork);
925//        write(l_fork, quote(groebner(eval(i))));
926        write(l_fork, quote(groebner(eval(i),eval(Method))));
927//###Fehlermeldung:
928// ***dError: undef. ringorder used
929// occured at:
930
931        // sleep in small intervalls for appr. one second
932        if (wait > 0)
933        {
934          while(j < 1000000)
935          {
936            if (status(l_fork, "read", "ready", j)) {break;}
937            j = j + j;
938          }
939        }
940
941        // sleep in intervalls of one second from now on
942        j = 1;
943        while (j < wait)
944        {
945          if (status(l_fork, "read", "ready", 1000000)) {break;}
946          j = j + 1;
947        }
948
949        if (status(l_fork, "read", "ready"))
950        {
951          def result = read(l_fork);
952          if (bs != nameof(basering))
953          {
954            def PP = basering;
955            setring P;
956            def result = imap(PP, result);
957            kill PP;
958          }
959          if (defined(groebner_error)==1)
960          {
961            kill groebner_error;
962          }
963          kill l_fork;
964        }
965        else
966        {
967          ideal result;
968          if (! defined(groebner_error))
969          {
970            int groebner_error = 1;
971            export groebner_error;
972          }
973          "** groebner did not finish";
974          j = system("sh", "kill " + string(pid));
975        }
976        return (result);
977    }
978    else
979    {
980      "** groebner with a time limit on computation is not supported
981          in this configuration";
982    }
983  }
984
985 //=========== we are still here -- do the actual computation =============
986
987//--------------------- save data from basering ---------------------------
988  poly Minpoly = minpoly;      //minimal polynomial
989  int was_minpoly;             //remembers if there was a minpoly in P
990  if (size(Minpoly) > 0)
991  {
992     was_minpoly = 1;
993  }
994
995  ideal Qideal = ideal(P);      //defining the quotient ideal if P is a qring
996  int was_qring;                //remembers if basering was a qring
997  //int is_homog = 1;
998  if (size(Qideal) > 0)
999  {
1000     was_qring = 1;
1001     //is_homog = homog(Qideal); //remembers if Qideal was homog (homog(0)=1)
1002  }
1003  list BRlist = ringlist(P);     //ringlist of basering
1004
1005  // save ordering of basering P for later use
1006  list ord_P = BRlist[3];       //should be available in all rings
1007  string ordstr_P = ordstr(P);
1008  int nvars_P = nvars(P);
1009  int npars_P = npars(P);
1010  intvec w;                     //for ringweights of basering P
1011  for(k=1;  k<=nvars_P; k++)
1012  {
1013     w[k]=deg(var(k));
1014  }
1015  int neg=1-attrib (P,"global");
1016
1017  //save options:
1018  intvec opt=option(get);
1019  string s_opt = option();
1020  int p_opt;
1021  if (find(s_opt, "prot"))  { p_opt = 1; }
1022
1023//------------------ cases where std is always used ------------------------
1024//If other methods are not implemented or do not make sense, i.e. for
1025//local or mixed orderings, matrix orderings, extra weight vector
1026//### Todo: extra weight vector should be allowed for e.g. with "hilb"
1027
1028  if(  //( find(ordstr_P,"s") > 0 ) || // covered by neg
1029       ( find(ordstr_P,"M") > 0 )  || ( find(ordstr_P,"a") > 0 )  || neg )
1030  {
1031    if (p_opt) { "std in basering"; }
1032    return(std(i));
1033  }
1034
1035//now we have:
1036//ideal or module, global ordering, no matrix ordering, no extra weight vector
1037//The interesting cases start now.
1038
1039 //------------------ classify the possible settings ---------------------
1040  string algorithm;       //possibilities: std, slimgb, stdorslimgb
1041  string conversion;      //possibilities: hilb, fglm, hilborfglm, no
1042  string partovar;        //possibilities: yes, no
1043  string order;           //possibilities: simple, !simple
1044  string direct;          //possibilities: yes, no
1045
1046  //define algorithm:
1047  if( find(method,"std") && !find(method,"slimgb") )
1048  {
1049    algorithm = "std";
1050  }
1051  if( find(method,"slimgb") && !find(method,"std") )
1052  {
1053    algorithm = "slimgb";
1054  }
1055  if( find(method,"std") && find(method,"slimgb") ||
1056      (!find(method,"std") && !find(method,"slimgb")) )
1057  {
1058    algorithm = "stdorslimgb";
1059  }
1060
1061  //define conversion:
1062  if( find(method,"hilb") && !find(method,"fglm") )
1063  {
1064     conversion = "hilb";
1065  }
1066  if( find(method,"fglm") && !find(method,"hilb") )
1067  {
1068    conversion = "fglm";
1069  }
1070  if( find(method,"fglm") && find(method,"hilb") )
1071  {
1072    conversion = "hilborfglm";
1073  }
1074  if( !find(method,"fglm") && !find(method,"hilb") )
1075  {
1076    conversion = "no";
1077  }
1078
1079  //define partovar:
1080  //if( find(method,"par2var") && npars_P > 0 )   //V1
1081  if( find(method,"par2var") || npars_P > 1 )     //V2
1082  {
1083     partovar = "yes";
1084  }
1085  else
1086  {
1087     partovar = "no";
1088  }
1089
1090  //define order:
1091  if (system("nblocks") <= 2)
1092  {
1093    if ( find(ordstr_P,"M")+find(ordstr_P,"lp")+find(ordstr_P,"rp") <= 0 )
1094    {
1095      order = "simple";
1096    }
1097  }
1098
1099  //define direct:
1100  if ( (order=="simple" && (size(method)==0)) ||
1101       (size(BRlist)>4) ||
1102        (order=="simple" && (method==",par2var" && npars_P==0 )) ||
1103         (conversion=="no" && partovar=="no" &&
1104           (algorithm=="std" || algorithm=="slimgb" ||
1105            (find(method,"std") && find(method,"slimgb")) ) ) )
1106  {
1107    direct = "yes";
1108  }
1109  else
1110  {
1111    direct = "no";
1112  }
1113
1114  //order=="simple" means that the ordering of the variables consists of one
1115  //block which is not a matrix ordering and not a lexicographical ordering.
1116  //(Note:Singular counts always least 2 blocks, one is for module component):
1117  //Call a method "direct" if conversion=="no" && partovar="no" which means
1118  //that we apply std or slimgb dircet in the basering (exception
1119  //as long as slimgb does not know qrings: in a qring of a ring P
1120  //the ideal Qideal is added to the ideal and slimgb is applied in P).
1121  //We apply a direct method if we have a simple monomial ordering, if no
1122  //conversion (fglm or hilb) is specified and if the parameters shall
1123  //not be made to variables
1124  //BRlist (=ringlist of basering) > 4 if the basering is non-commutative
1125//---------------------------- direct methods -----------------------------
1126  if ( direct == "yes" )
1127  {
1128  //if ( algorithm=="std" || (algorithm=="stdorslimgb" && char(P)>0) )   //V1
1129    if ( algorithm=="std" || (algorithm=="stdorslimgb") )                //V2
1130    {
1131      if (p_opt) { "std in " + string(P); }
1132      return(std(i));
1133    }
1134  //if( algorithm=="slimgb" || (algorithm=="stdorslimgb" && char(P)==0)) //V1
1135    if ( algorithm=="slimgb" )                                           //V2
1136    {
1137      return(qslimgb(i));
1138    }
1139  }
1140
1141//--------------------------- indirect methods -----------------------------
1142//indirect methods are methods where a conversion is used with a ring change
1143//We are in the following situation:
1144//direct=="no" (i.e. "hilb" or "fglm" or "par2var" is given)
1145//or no method is given and we have a complicated monomial ordering
1146//V1: "par2var" is not a default strategy, it must be explicitely
1147//given in order to be performed.
1148//V2: "par2var" is a default strategy if there are more than 1 parameters
1149
1150//------------ case where no parameters are made to variables  -------------
1151  if (  partovar == "no" && conversion == "hilb"
1152    || (partovar == "no" && conversion == "fglm" )
1153    || (partovar == "no" && conversion == "hilborfglm" )
1154    || (partovar == "no" && conversion == "no" && direct == "no") )
1155  //last case: heuristic
1156  {
1157    if ( conversion=="fglm" )
1158    {
1159    //if ( algorithm=="std" || (algorithm=="stdorslimgb" && char(P)>0) ) //V1
1160      if ( algorithm=="std" || (algorithm=="stdorslimgb") )              //V2
1161      {
1162        return (stdfglm(i,"std"));
1163      }
1164    //if(algorithm=="slimgb" || (algorithm=="stdorslimgb" && char(P)==0))//V1
1165      if( algorithm=="slimgb" )                                          //V2
1166      {
1167        return (stdfglm(i,"slimgb"));
1168      }
1169    }
1170    else
1171    {
1172    //if ( algorithm=="std" || (algorithm=="stdorslimgb" && char(P)>0) )//V1
1173      if ( algorithm=="std" || (algorithm=="stdorslimgb" ) )            //V2
1174      {
1175        return (stdhilb(i,"std"));
1176      }
1177    //if(algorithm=="slimgb" || (algorithm=="stdorslimgb" && char(P)==0))//V1
1178      if ( algorithm=="slimgb" )                                         //V2
1179      {
1180        return (stdhilb(i,"slimgb"));
1181      }
1182    }
1183  }
1184
1185//------------ case where parameters are made to variables  ----------------
1186//define a ring Phelp via par2varRing in which the parameters are variables
1187
1188  else
1189  {
1190    // reset options
1191    option(none);
1192    // turn on options prot, mem, redSB, intStrategy if previously set
1193    if ( find(s_opt, "prot") )
1194      { option(prot); }
1195    if ( find(s_opt, "mem") )
1196      { option(mem); }
1197    if ( find(s_opt, "redSB") )
1198      { option(redSB); }
1199    if ( find(s_opt, "intStrategy") )
1200      { option(intStrategy); }
1201
1202    //first clear denominators of parameters
1203    if (npars_P > 0)
1204    {
1205      for( k=ncols(i); k>0; k-- )
1206      { i[k]=cleardenom(i[k]); }
1207    }
1208
1209    def Phelp = par2varRing(i)[1];   //minpoly is mapped with i
1210    setring Phelp;
1211    def i = Id(1);
1212    //is_homog = homog(i);
1213
1214    //If parameters are converted to ring variables, they appear in an extra
1215    //block. Therefore we use always hilb for this block ordering:
1216    if ( conversion=="fglm" )
1217    {
1218      i = (stdfglm(i));       //only uesful for 1 parameter with minpoly
1219    }
1220    else
1221    {
1222    //if ( algorithm=="std" || (algorithm=="stdorslimgb" && char(P)>0) )//V1
1223      if ( algorithm=="std" || (algorithm=="stdorslimgb" ))             //V2
1224      {
1225        i = stdhilb(i,"std");
1226      }
1227    //if(algorithm=="slimgb" || (algorithm=="stdorslimgb" && char(P)==0))//V1
1228      if ( algorithm=="slimgb" )                                         //V2
1229      {
1230        i = stdhilb(i,"slimgb");
1231      }
1232    }
1233  }
1234
1235//-------------------- go back to original ring ---------------------------
1236//The main computation is done. However, the SB coming from a ring with
1237//extra variables is in general too big. We simplify it before mapping it
1238//to the basering.
1239
1240  if (p_opt) { "//simplification"; }
1241
1242  if (was_minpoly)
1243  {
1244    ideal Minpoly = imap(P,Minpoly);
1245    attrib(Minpoly,"isSB",1);
1246    i = simplify(NF(i,Minpoly),2);
1247  }
1248
1249  def Li = lead(i);
1250  setring P;
1251  def Li = imap(Phelp,Li);
1252  Li = simplify(Li,32);
1253  intvec vi;
1254  for (k=1; k<=ncols(Li); k++)
1255  {
1256    vi[k] = Li[k]==0;
1257  }
1258
1259  setring Phelp;
1260  for (k=1;  k<=size(i) ;k++)
1261  {
1262    if(vi[k]==1)
1263    {
1264      i[k]=0;
1265    }
1266  }
1267  i = simplify(i,2);
1268
1269  setring P;
1270  if (p_opt) { "//imap to original ring"; }
1271  i = imap(Phelp,i);
1272  kill Phelp;
1273  i = simplify(i,34);
1274
1275  // clean-up time
1276  option(set, opt);
1277  if (find(s_opt, "redSB") > 0)
1278  {
1279    if (p_opt) { "//interreduction"; }
1280    i=interred(i);
1281  }
1282  attrib(i, "isSB", 1);
1283  return (i);
1284}
1285example
1286{ "EXAMPLE: "; echo=2;
1287  intvec opt = option(get);
1288  option(prot);
1289  ring r  = 0,(a,b,c,d),dp;
1290  ideal i = a+b+c+d,ab+ad+bc+cd,abc+abd+acd+bcd,abcd-1;
1291  groebner(i);
1292
1293  ring s  = 0,(a,b,c,d),lp;
1294  ideal i = imap(r,i);
1295  groebner(i,"hilb");
1296
1297  ring R  = (0,a),(b,c,d),lp;
1298  minpoly = a2+1;
1299  ideal i = a+b+c+d,ab+ad+bc+cd,abc+abd+acd+bcd,d2-c2b2;
1300  groebner(i,"par2var","slimgb");
1301
1302  groebner(i,"fglm");          //computes a reduced standard basis
1303
1304  option(set,opt);
1305}
1306
1307//////////////////////////////////////////////////////////////////////////
1308
1309proc res(list #)
1310"@c we do texinfo here:
1311@cindex resolution, computation of
1312@table @code
1313@item @strong{Syntax:}
1314@code{res (} ideal_expression@code{,} int_expression @code{[,} any_expression @code{])}
1315@*@code{res (} module_expression@code{,} int_expression @code{[,} any_expression @code{])}
1316@item @strong{Type:}
1317resolution
1318@item @strong{Purpose:}
1319computes a (possibly minimal) free resolution of an ideal or module using
1320a heuristically chosen method.
1321@* The second (int) argument (say @code{k}) specifies the length of
1322the resolution. If it is not positive then @code{k} is assumed to be the
1323number of variables of the basering.
1324@* If a third argument is given, the returned resolution is minimized.
1325
1326Depending on the input, the returned resolution is computed using the
1327following methods:
1328@table @asis
1329@item @strong{quotient rings:}
1330@code{nres} (classical method using syzygies) , see @ref{nres}.
1331
1332@item @strong{homogeneous ideals and k=0:}
1333@code{lres} (La'Scala's method), see @ref{lres}.
1334
1335@item @strong{not minimized resolution and (homogeneous input with k not 0, or local rings):}
1336@code{sres} (Schreyer's method), see @ref{sres}.
1337
1338@item @strong{all other inputs:}
1339@code{mres} (classical method), see @ref{mres}.
1340@end table
1341@item @strong{Note:}
1342Accessing single elements of a resolution may require some partial
1343computations to be finished and may therefore take some time.
1344@end table
1345@c ref
1346See also
1347@ref{betti};
1348@ref{ideal};
1349@ref{minres};
1350@ref{module};
1351@ref{mres};
1352@ref{nres};
1353@ref{lres};
1354@ref{hres};
1355@ref{sres};
1356@ref{resolution}.
1357@c ref
1358"
1359{
1360   def P=basering;
1361   if (size(#) < 2)
1362   {
1363     ERROR("res: need at least two arguments: ideal/module, int");
1364   }
1365
1366   def m=#[1]; //the ideal or module
1367   int i=#[2]; //the length of the resolution
1368   if (i< 0) { i=0;}
1369
1370   string varstr_P = varstr(P);
1371
1372   int p_opt;
1373   string s_opt = option();
1374   // set p_opt, if option(prot) is set
1375   if (find(s_opt, "prot"))
1376   {
1377     p_opt = 1;
1378   }
1379
1380   if( (size(ideal(basering)) > 0) || (size(ringlist(P)) > 4) )
1381   {
1382     // the quick hack for qrings - seems to fit most needs
1383     // (lres is not implemented for qrings, sres is not so efficient)
1384     // || non-commutative, since only n/m-res are implemented for NC rings
1385     if (p_opt) { "using nres";}
1386     return(nres(m,i));
1387   }
1388
1389   if(homog(m)==1)
1390   {
1391      resolution re;
1392      if (((i==0) or (i>=nvars(basering))) && typeof(m) != "module")
1393      {
1394        //LaScala for the homogeneous case and i == 0
1395        if (p_opt) { "using lres";}
1396        re=lres(m,i);
1397        if(size(#)>2)
1398        {
1399           re=minres(re);
1400        }
1401      }
1402      else
1403      {
1404        if(size(#)>2)
1405        {
1406          if (p_opt) { "using mres";}
1407          re=mres(m,i);
1408        }
1409        else
1410        {
1411          if (p_opt) { "using sres";}
1412          re=sres(std(m),i);
1413        }
1414      }
1415      return(re);
1416   }
1417
1418   //mres for the global non homogeneous case
1419   if(find(ordstr(P),"s")==0)
1420   {
1421      string ri= "ring Phelp ="
1422                  +string(char(P))+",("+varstr_P+"),(dp,C);";
1423      ri = ri + "minpoly = "+string(minpoly) + ";";
1424      execute(ri);
1425      def m=imap(P,m);
1426      if (p_opt) { "using mres in another ring";}
1427      list re=mres(m,i);
1428      setring P;
1429      resolution result=imap(Phelp,re);
1430      if (size(#) > 2) {result = minres(result);}
1431      return(result);
1432   }
1433
1434   //sres for the local case and not minimal resolution
1435   if(size(#)<=2)
1436   {
1437      string ri= "ring Phelp ="
1438                  +string(char(P))+",("+varstr_P+"),(ls,c);";
1439      ri = ri + "minpoly = "+string(minpoly) + ";";
1440      execute(ri);
1441      def m=imap(P,m);
1442      m=std(m);
1443      if (p_opt) { "using sres in another ring";}
1444      list re=sres(m,i);
1445      setring P;
1446      resolution result=imap(Phelp,re);
1447      return(result);
1448   }
1449
1450   //mres for the local case and minimal resolution
1451   string ri= "ring Phelp ="
1452                  +string(char(P))+",("+varstr_P+"),(ls,C);";
1453   ri = ri + "minpoly = "+string(minpoly) + ";";
1454   execute(ri);
1455   def m=imap(P,m);
1456    if (p_opt) { "using mres in another ring";}
1457   list re=mres(m,i);
1458   setring P;
1459   resolution result=imap(Phelp,re);
1460   result = minres(result);
1461   return(result);
1462}
1463example
1464{"EXAMPLE:"; echo = 2;
1465  ring r=0,(x,y,z),dp;
1466  ideal i=xz,yz,x3-y3;
1467  def l=res(i,0); // homogeneous ideal: uses lres
1468  l;
1469  print(betti(l), "betti"); // input to betti may be of type resolution
1470  l[2];         // element access may take some time
1471  i=i,x+1;
1472  l=res(i,0);   // inhomogeneous ideal: uses mres
1473  l;
1474  ring rs=0,(x,y,z),ds;
1475  ideal i=imap(r,i);
1476  def l=res(i,0); // local ring not minimized: uses sres
1477  l;
1478  res(i,0,0);     // local ring and minimized: uses mres
1479}
1480/////////////////////////////////////////////////////////////////////////
1481
1482proc quot (m1,m2,list #)
1483"SYNTAX: @code{quot (} module_expression@code{,} module_expression @code{)}
1484         @*@code{quot (} module_expression@code{,} module_expression@code{,}
1485            int_expression @code{)}
1486         @*@code{quot (} ideal_expression@code{,} ideal_expression @code{)}
1487         @*@code{quot (} ideal_expression@code{,} ideal_expression@code{,}
1488            int_expression @code{)}
1489TYPE:    ideal
1490SYNTAX:  @code{quot (} module_expression@code{,} ideal_expression @code{)}
1491TYPE:    module
1492PURPOSE: computes the quotient of the 1st and the 2nd argument.
1493         If a 3rd argument @code{n} is given the @code{n}-th method is used
1494         (@code{n}=1...5).
1495SEE ALSO: quotient
1496EXAMPLE: example quot; shows an example"
1497{
1498  if (((typeof(m1)!="ideal") and (typeof(m1)!="module"))
1499     or ((typeof(m2)!="ideal") and (typeof(m2)!="module")))
1500  {
1501    "USAGE:   quot(m1, m2[, n]); m1, m2 two submodules of k^s,";
1502    "         n (optional) integer (1<= n <=5)";
1503    "RETURN:  the quotient of m1 and m2";
1504    "EXAMPLE: example quot; shows an example";
1505    return();
1506  }
1507  if (typeof(m1)!=typeof(m2))
1508  {
1509    return(quotient(m1,m2));
1510  }
1511  if (size(#)>0)
1512  {
1513    if (typeof(#[1])=="int" )
1514    {
1515      return(quot1(m1,m2,#[1]));
1516    }
1517  }
1518  else
1519  {
1520    return(quot1(m1,m2,2));
1521  }
1522}
1523example
1524{ "EXAMPLE:"; echo = 2;
1525  ring r=181,(x,y,z),(c,ls);
1526  ideal id1=maxideal(4);
1527  ideal id2=x2+xyz,y2-z3y,z3+y5xz;
1528  option(prot);
1529  ideal id3=quotient(id1,id2);
1530  id3;
1531  ideal id4=quot(id1,id2,1);
1532  id4;
1533  ideal id5=quot(id1,id2,2);
1534  id5;
1535}
1536
1537static proc quot1 (module m1, module m2,int n)
1538"USAGE:   quot1(m1, m2, n); m1, m2 two submodules of k^s,
1539         n integer (1<= n <=5)
1540RETURN:  the quotient of m1 and m2
1541EXAMPLE: example quot1; shows an example"
1542{
1543  if (n==1)
1544  {
1545    return(quotient1(m1,m2));
1546  }
1547  else
1548  {
1549    if (n==2)
1550    {
1551      return(quotient2(m1,m2));
1552    }
1553    else
1554    {
1555      if (n==3)
1556      {
1557        return(quotient3(m1,m2));
1558      }
1559      else
1560      {
1561        if (n==4)
1562        {
1563          return(quotient4(m1,m2));
1564        }
1565        else
1566        {
1567          if (n==5)
1568          {
1569            return(quotient5(m1,m2));
1570          }
1571          else
1572          {
1573            return(quotient(m1,m2));
1574          }
1575        }
1576      }
1577    }
1578  }
1579}
1580example
1581{ "EXAMPLE:"; echo = 2;
1582  ring r=181,(x,y,z),(c,ls);
1583  ideal id1=maxideal(4);
1584  ideal id2=x2+xyz,y2-z3y,z3+y5xz;
1585  option(prot);
1586  ideal id6=quotient(id1,id2);
1587  id6;
1588  ideal id7=quot1(id1,id2,1);
1589  id7;
1590  ideal id8=quot1(id1,id2,2);
1591  id8;
1592}
1593
1594static proc quotient0(module a,module b)
1595{
1596  module mm=b+a;
1597  resolution rs=lres(mm,0);
1598  list I=list(rs);
1599  matrix M=I[2];
1600  matrix A[1][nrows(M)]=M[1..nrows(M),1];
1601  ideal i=A;
1602  return (i);
1603}
1604proc quotient1(module a,module b)  //17sec
1605"USAGE:   quotient1(m1, m2); m1, m2 two submodules of k^s,
1606RETURN:  the quotient of m1 and m2"
1607{
1608  int i;
1609  a=std(a);
1610  module dummy;
1611  module B=NF(b,a)+dummy;
1612  ideal re=quotient(a,module(B[1]));
1613  for(i=2;i<=ncols(B);i++)
1614  {
1615     re=intersect1(re,quotient(a,module(B[i])));
1616  }
1617  return(re);
1618}
1619proc quotient2(module a,module b)    //13sec
1620"USAGE:   quotient2(m1, m2); m1, m2 two submodules of k^s,
1621RETURN:  the quotient of m1 and m2"
1622{
1623  a=std(a);
1624  module dummy;
1625  module bb=NF(b,a)+dummy;
1626  int i=ncols(bb);
1627  ideal re=quotient(a,module(bb[i]));
1628  bb[i]=0;
1629  module temp;
1630  module temp1;
1631  module bbb;
1632  int mx;
1633  i=i-1;
1634  while (1)
1635  {
1636    if (i==0) break;
1637    temp = a+bb*re;
1638    temp1 = lead(interred(temp));
1639    mx=ncols(a);
1640    if (ncols(temp1)>ncols(a))
1641    {
1642      mx=ncols(temp1);
1643    }
1644    temp1 = matrix(temp1,1,mx)-matrix(lead(a),1,mx);
1645    temp1 = dummy+temp1;
1646    if (deg(temp1[1])<0) break;
1647    re=intersect1(re,quotient(a,module(bb[i])));
1648    bb[i]=0;
1649    i = i-1;
1650  }
1651  return(re);
1652}
1653proc quotient3(module a,module b)   //89sec
1654"USAGE:   quotient3(m1, m2); m1, m2 two submodules of k^s,
1655         only for global rings
1656RETURN:  the quotient of m1 and m2"
1657{
1658  string s="ring @newr=("+charstr(basering)+
1659           "),("+varstr(basering)+",@t,@w),dp;";
1660  def @newP=basering;
1661  execute(s);
1662  module b=imap(@newP,b);
1663  module a=imap(@newP,a);
1664  int i;
1665  int j=ncols(b);
1666  vector @b;
1667  for(i=1;i<=j;i++)
1668  {
1669     @b=@b+@t^(i-1)*@w^(j-i+1)*b[i];
1670  }
1671  ideal re=quotient(a,module(@b));
1672  setring @newP;
1673  ideal re=imap(@newr,re);
1674  return(re);
1675}
1676proc quotient5(module a,module b)   //89sec
1677"USAGE:   quotient5(m1, m2); m1, m2 two submodules of k^s,
1678         only for global rings
1679RETURN:  the quotient of m1 and m2"
1680{
1681  string s="ring @newr=("+charstr(basering)+
1682           "),("+varstr(basering)+",@t),dp;";
1683  def @newP=basering;
1684  execute(s);
1685  module b=imap(@newP,b);
1686  module a=imap(@newP,a);
1687  int i;
1688  int j=ncols(b);
1689  vector @b;
1690  for(i=1;i<=j;i++)
1691  {
1692     @b=@b+@t^(i-1)*b[i];
1693  }
1694  @b=homog(@b,@w);
1695  ideal re=quotient(a,module(@b));
1696  setring @newP;
1697  ideal re=imap(@newr,re);
1698  return(re);
1699}
1700proc quotient4(module a,module b)   //95sec
1701"USAGE:   quotient4(m1, m2); m1, m2 two submodules of k^s,
1702         only for global rings
1703RETURN:  the quotient of m1 and m2"
1704{
1705  string s="ring @newr=("+charstr(basering)+
1706           "),("+varstr(basering)+",@t),dp;";
1707  def @newP=basering;
1708  execute(s);
1709  module b=imap(@newP,b);
1710  module a=imap(@newP,a);
1711  int i;
1712  vector @b=b[1];
1713  for(i=2;i<=ncols(b);i++)
1714  {
1715     @b=@b+@t^(i-1)*b[i];
1716  }
1717  matrix sy=modulo(@b,a);
1718  ideal re=sy;
1719  setring @newP;
1720  ideal re=imap(@newr,re);
1721  return(re);
1722}
1723static proc intersect1(ideal i,ideal j)
1724{
1725  def R=basering;
1726  execute("ring gnir = ("+charstr(basering)+"),
1727                       ("+varstr(basering)+",@t),(C,dp);");
1728  ideal i=var(nvars(basering))*imap(R,i)+(var(nvars(basering))-1)*imap(R,j);
1729  ideal j=eliminate(i,var(nvars(basering)));
1730  setring R;
1731  map phi=gnir,maxideal(1);
1732  return(phi(j));
1733}
1734
1735//////////////////////////////////////////////////////////////////
1736///
1737/// sprintf, fprintf printf
1738///
1739proc sprintf(string fmt, list #)
1740"SYNTAX:  @code{sprintf (} string_expression @code{[,} any_expressions
1741               @code{] )}
1742RETURN:   string
1743PURPOSE:  @code{sprintf(fmt,...);} performs output formatting. The first
1744          argument is a format control string. Additional arguments may be
1745          required, depending on the content of the control string. A series
1746          of output characters is generated as directed by the control string;
1747          these characters are returned as a string. @*
1748          The control string @code{fmt} is simply text to be copied,
1749          except that the string may contain conversion specifications.@*
1750          Type @code{help print;} for a listing of valid conversion
1751          specifications. As an addition to the conversions of @code{print},
1752          the @code{%n} and @code{%2} conversion specification does not
1753          consume an additional argument, but simply generates a newline
1754          character.
1755NOTE:     If one of the additional arguments is a list, then it should be
1756          wrapped in an additional @code{list()} command, since passing a list
1757          as an argument flattens the list by one level.
1758SEE ALSO: fprintf, printf, print, string
1759EXAMPLE : example sprintf; shows an example
1760"
1761{
1762  int sfmt = size(fmt);
1763  if (sfmt  <= 1)
1764  {
1765    return (fmt);
1766  }
1767  int next, l, nnext;
1768  string ret;
1769  list formats = "%l", "%s", "%2l", "%2s", "%t", "%;", "%p", "%b", "%n", "%2";
1770  while (1)
1771  {
1772    if (size(#) <= 0)
1773    {
1774      return (ret + fmt);
1775    }
1776    nnext = 0;
1777    while (nnext < sfmt)
1778    {
1779      nnext = find(fmt, "%", nnext + 1);
1780      if (nnext == 0)
1781      {
1782        next = 0;
1783        break;
1784      }
1785      l = 1;
1786      while (l <= size(formats))
1787      {
1788        next = find(fmt, formats[l], nnext);
1789        if (next == nnext) break;
1790        l++;
1791      }
1792      if (next == nnext) break;
1793    }
1794    if (next == 0)
1795    {
1796      return (ret + fmt);
1797    }
1798    if (formats[l] != "%2" && formats[l] != "%n")
1799    {
1800      ret = ret + fmt[1, next - 1] + print(#[1], formats[l]);
1801      # = delete(#, 1);
1802    }
1803    else
1804    {
1805      ret = ret + fmt[1, next - 1] + print("", "%2s");
1806    }
1807    if (size(fmt) <= (next + size(formats[l]) - 1))
1808    {
1809      return (ret);
1810    }
1811    fmt = fmt[next + size(formats[l]), size(fmt)-next-size(formats[l]) + 1];
1812  }
1813}
1814example
1815{ "EXAMPLE:"; echo=2;
1816  ring r=0,(x,y,z),dp;
1817  module m=[1,y],[0,x+z];
1818  intmat M=betti(mres(m,0));
1819  list l = r, m, M;
1820  string s = sprintf("s:%s,%n l:%l", 1, 2); s;
1821  s = sprintf("s:%n%s", l); s;
1822  s = sprintf("s:%2%s", list(l)); s;
1823  s = sprintf("2l:%n%2l", list(l)); s;
1824  s = sprintf("%p", list(l)); s;
1825  s = sprintf("%;", list(l)); s;
1826  s = sprintf("%b", M); s;
1827}
1828
1829proc printf(string fmt, list #)
1830"SYNTAX:  @code{printf (} string_expression @code{[,} any_expressions@code{] )}
1831RETURN:   none
1832PURPOSE:  @code{printf(fmt,...);} performs output formatting. The first
1833          argument is a format control string. Additional arguments may be
1834          required, depending on the content of the control string. A series
1835          of output characters is generated as directed by the control string;
1836          these characters are displayed (i.e., printed to standard out). @*
1837          The control string @code{fmt} is simply text to be copied, except
1838          that the string may contain conversion specifications. @*
1839          Type @code{help print;} for a listing of valid conversion
1840          specifications. As an addition to the conversions of @code{print},
1841          the @code{%n} and @code{%2} conversion specification does not
1842          consume an additional argument, but simply generates a newline
1843          character.
1844NOTE:     If one of the additional arguments is a list, then it should be
1845          enclosed once more into a @code{list()} command, since passing a
1846          list as an argument flattens the list by one level.
1847SEE ALSO: sprintf, fprintf, print, string
1848EXAMPLE : example printf; shows an example
1849"
1850{
1851  write("", sprintf(fmt, #));
1852}
1853example
1854{ "EXAMPLE:"; echo=2;
1855  ring r=0,(x,y,z),dp;
1856  module m=[1,y],[0,x+z];
1857  intmat M=betti(mres(m,0));
1858  list l=r,m,matrix(M);
1859  printf("s:%s,l:%l",1,2);
1860  printf("s:%s",l);
1861  printf("s:%s",list(l));
1862  printf("2l:%2l",list(l));
1863  printf("%p",matrix(M));
1864  printf("%;",matrix(M));
1865  printf("%b",M);
1866}
1867
1868
1869proc fprintf(link l, string fmt, list #)
1870"SYNTAX:  @code{fprintf (} link_expression@code{,} string_expression @code{[,}
1871                any_expressions@code{] )}
1872RETURN:   none
1873PURPOSE:  @code{fprintf(l,fmt,...);} performs output formatting.
1874          The second argument is a format control string. Additional
1875          arguments may be required, depending on the content of the
1876          control string. A series of output characters is generated as
1877          directed by the control string; these characters are
1878          written to the link l.
1879          The control string @code{fmt} is simply text to be copied, except
1880          that the string may contain conversion specifications.@*
1881          Type @code{help print;} for a listing of valid conversion
1882          specifications. As an addition to the conversions of @code{print},
1883          the @code{%n} and @code{%2} conversion specification does not
1884          consume an additional argument, but simply generates a newline
1885          character.
1886NOTE:     If one of the additional arguments is a list, then it should be
1887          enclosed once more into a @code{list()} command, since passing
1888          a list as an argument flattens the list by one level.
1889SEE ALSO: sprintf, printf, print, string
1890EXAMPLE : example fprintf; shows an example
1891"
1892{
1893  write(l, sprintf(fmt, #));
1894}
1895example
1896{ "EXAMPLE:"; echo=2;
1897  ring r=0,(x,y,z),dp;
1898  module m=[1,y],[0,x+z];
1899  intmat M=betti(mres(m,0));
1900  list l=r,m,M;
1901  link li="";   // link to stdout
1902  fprintf(li,"s:%s,l:%l",1,2);
1903  fprintf(li,"s:%s",l);
1904  fprintf(li,"s:%s",list(l));
1905  fprintf(li,"2l:%2l",list(l));
1906  fprintf(li,"%p",list(l));
1907  fprintf(li,"%;",list(l));
1908  fprintf(li,"%b",M);
1909}
1910
1911//////////////////////////////////////////////////////////////////////////
1912
1913/*
1914proc minres(list #)
1915{
1916  if (size(#) == 2)
1917  {
1918    if (typeof(#[1]) == "ideal" || typeof(#[1]) == "module")
1919    {
1920      if (typeof(#[2] == "int"))
1921      {
1922        return (res(#[1],#[2],1));
1923      }
1924    }
1925  }
1926
1927  if (typeof(#[1]) == "resolution")
1928  {
1929    return minimizeres(#[1]);
1930  }
1931  else
1932  {
1933    return minimizeres(#);
1934  }
1935
1936}
1937*/
1938///////////////////////////////////////////////////////////////////////////////
1939
1940proc weightKB(def stc, int dd, list wim)
1941"SYNTAX: @code{weightKB (} module_expression@code{,} int_expression @code{,}
1942            list_expression @code{)}@*
1943         @code{weightKB (} ideal_expression@code{,} int_expression@code{,}
1944            list_expression @code{)}
1945RETURN:  the same as the input type of the first argument
1946PURPOSE: If @code{I,d,wim} denotes the three arguments then weightKB
1947         computes the weighted degree- @code{d} part of a vector space basis
1948         (consisting of monomials) of the quotient ring, resp. of the
1949         quotient module, modulo @code{I} w.r.t. weights given by @code{wim}
1950         The information about the weights is given as a list of two intvec:
1951            @code{wim[1]} weights for all variables (positive),
1952            @code{wim[2]} weights for the module generators.
1953NOTE:    This is a generalization of the command @code{kbase} with the same
1954         first two arguments.
1955SEE ALSO: kbase
1956EXAMPLE: example weightKB; shows an example
1957"
1958{
1959  if(checkww(wim)){ERROR("wrong weights";);}
1960  kbclass();
1961  wwtop=wim[1];
1962  stc=interred(lead(stc));
1963  if(typeof(stc)=="ideal")
1964  {
1965    stdtop=stc;
1966    ideal out=widkbase(dd);
1967    delkbclass();
1968    out=simplify(out,2); // delete 0
1969    return(out);
1970  }
1971  list mbase=kbprepare(stc);
1972  module mout;
1973  int im,ii;
1974  if(size(wim)>1){mmtop=wim[2];}
1975  else{mmtop=0;}
1976  for(im=size(mbase);im>0;im--)
1977  {
1978    stdtop=mbase[im];
1979    if(im>size(mmtop)){ii=dd;}
1980    else{ii=dd-mmtop[im];}
1981    mout=mout+widkbase(ii)*gen(im);
1982  }
1983  delkbclass();
1984  mout=simplify(mout,2); // delete 0
1985  return(mout);
1986}
1987example
1988{ "EXAMPLE:"; echo=2;
1989  ring R=0, (x,y), wp(1,2);
1990  weightKB(ideal(0),3,intvec(1,2));
1991}
1992
1993///////////////////////////////////////////////////////////////////////////////
1994
1995proc datetime()
1996"SYNTAX: @code{datetime ()}
1997RETURN:  string
1998PURPOSE: return the curent date and time as a string
1999EXAMPLE: example datetime; shows an example
2000"
2001{
2002  return(read("|: date"));
2003}
2004example
2005{ "EXAMPLE:"; echo=2;
2006  datetime();
2007}
2008
2009///////////////////////////////////////////////////////////////////////////////
2010// construct global values
2011static proc kbclass()
2012{
2013  intvec wwtop,mmtop;
2014  export (wwtop,mmtop);
2015  ideal stdtop,kbtop;
2016  export (stdtop,kbtop);
2017}
2018// delete global values
2019static proc delkbclass()
2020{
2021  kill wwtop,mmtop;
2022  kill stdtop,kbtop;
2023}
2024//  select parts of the modul
2025static proc kbprepare(module mstc)
2026{
2027  list rr;
2028  ideal kk;
2029  int i1,i2;
2030  mstc=transpose(mstc);
2031  for(i1=ncols(mstc);i1>0;i1--)
2032  {
2033    kk=0;
2034    for(i2=nrows(mstc[i1]);i2>0;i2--)
2035    {
2036      kk=kk+mstc[i1][i2];
2037    }
2038    rr[i1]=kk;
2039  }
2040  return(rr);
2041}
2042//  check for weights
2043static proc checkww(list vv)
2044{
2045  if(typeof(vv[1])!="intvec"){return(1);}
2046  intvec ww=vv[1];
2047  int mv=nvars(basering);
2048  if(size(ww)<mv){return(1);}
2049  while(mv>0)
2050  {
2051    if(ww[mv]<=0){return(1);}
2052    mv--;
2053  }
2054  if(size(vv)>1)
2055  {
2056    if(typeof(vv[2])!="intvec"){return(1);}
2057  }
2058  return(0);
2059}
2060///////////////////////////////////////////////////////
2061// The "Caller" for ideals
2062//    dd   - the degree of the result
2063static proc widkbase(int dd)
2064{
2065  if((size(stdtop)==1)&&(deg(stdtop[1])==0)){return(0);}
2066  if(dd<=0)
2067  {
2068    if(dd<0){return(0);}
2069    else{return(1);}
2070  }
2071  int m1,m2;
2072  m1=nvars(basering);
2073  while(wwtop[m1]>dd)
2074  {
2075    m1--;
2076    if(m1==0){return(0);}
2077  }
2078  attrib(stdtop,"isSB",1);
2079  poly mo=1;
2080  if(m1==1)
2081  {
2082    m2=dd/wwtop[1];
2083    if((m2*wwtop[1])==dd)
2084    {
2085      mo=var(1)^m2;
2086      if(reduce(mo,stdtop)==mo){return(mo);}
2087      else{return(0);}
2088    }
2089  }
2090  kbtop=0;
2091  m2=dd;
2092  weightmon(m1-1,m2,mo);
2093  while(m2>=wwtop[m1])
2094  {
2095    m2=m2-wwtop[m1];
2096    mo=var(m1)*mo;
2097    if(m2==0)
2098    {
2099      if((mo!=0) and (reduce(mo,stdtop)==mo))
2100      {
2101        kbtop[ncols(kbtop)+1]=mo;
2102        return(kbtop);
2103      }
2104    }
2105    weightmon(m1-1,m2,mo);
2106  }
2107  return(kbtop);
2108}
2109/////////////////////////////////////////////////////////
2110// the recursive procedure
2111//    va     - number of the variable
2112//    drest  - rest of the degree
2113//    mm     - the candidate
2114static proc weightmon(int va, int drest, poly mm)
2115{
2116  while(wwtop[va]>drest)
2117  {
2118    va--;
2119    if(va==0){return();}
2120  }
2121  int m2;
2122  if(va==1)
2123  {
2124    m2=drest/wwtop[1];
2125    if((m2*wwtop[1])==drest)
2126    {
2127      mm=var(1)^m2*mm;
2128      if ((mm!=0) and (reduce(mm,stdtop)==mm))
2129      {
2130        kbtop[ncols(kbtop)+1]=mm;
2131      }
2132    }
2133    return();
2134  }
2135  m2=drest;
2136  if ((mm!=0) and (reduce(mm,stdtop)==mm))
2137  {
2138    weightmon(va-1,m2,mm);
2139  }
2140  while(m2>=wwtop[va])
2141  {
2142    m2=m2-wwtop[va];
2143    mm=var(va)*mm;
2144    if(m2==0)
2145    {
2146      if ((mm!=0) and (reduce(mm,stdtop)==mm))
2147      {
2148        kbtop[ncols(kbtop)+1]=mm;
2149        return();
2150      }
2151    }
2152     if ((mm!=0) and (reduce(mm,stdtop)==mm))
2153     {
2154       weightmon(va-1,m2,mm);
2155     }
2156  }
2157  return();
2158}
2159example
2160{ "EXAMPLE:"; echo=2;
2161  ring r=0,(x,y,z),dp;
2162  ideal i = x6,y4,xyz;
2163  intvec w = 2,3,6;
2164  weightKB(i, 12, list(w));
2165}
2166
2167///////////////////////////////////////////////////////////////////////////////
2168/*
2169                                Versuche:
2170///////////////////////////////////////////////////////////////////////////////
2171proc downsizeSB (I, list #)
2172"USAGE:   downsizeSB(I [,l]); I ideal, l list of integers [default: l=0]
2173RETURN:  intvec, say v, with v[j] either 1 or 0. We have v[j]=1 if
2174         leadmonom(I[j]) is divisible by some leadmonom(I[k]) or if
2175         leadmonom(i[j]) == leadmonom(i[k]) and l[j] >= l[k], with k!=j.
2176PURPOSE: The procedure is applied in a situation where the standard basis
2177         computation in the basering R is done via a conversion through an
2178         overring Phelp with additional variables and where a direct
2179         imap from Phelp to R is too expensive.
2180         Assume Phelp is created by the procedure @code{par2varRing} or
2181         @code{hilbRing} and IPhelp is a SB in Phelp [ with l[j]=
2182         length(IPhelp(j)) or any other integer reflecting the complexity
2183         of a IPhelp[j] ]. Let I = lead(IPhelp) mapped to R and compute
2184         v = downsizeSB(imap(Phelp,I),l) in R. Then, if Ihelp[j] is deleted
2185         for all j with v[j]=1, we can apply imap to the remaining generators
2186         of Ihelp and still get SB in R  (in general not minimal).
2187EXAMPLE: example downsizeSB; shows an example"
2188{
2189   int k,j;
2190   intvec v,l;
2191   poly M,N,W;
2192   int c=size(I);
2193   if( size(#) != 0 )
2194   {
2195     if ( typeof(#[1]) == "intvec" )
2196     {
2197       l = #[1];
2198     }
2199     else
2200     {
2201        ERROR("// 2nd argument must be an intvec");
2202     }
2203   }
2204
2205   l[c+1]=0;
2206   v[c]=0;
2207
2208   j=0;
2209   while(j<c-1)
2210   {
2211     j++;
2212     M = leadmonom(I[j]);
2213     if( M != 0 )
2214     {
2215        for( k=j+1; k<=c; k++ )
2216        {
2217          N = leadmonom(I[k]);
2218          if( N != 0 )
2219          {
2220             if( (M==N) && (l[j]>l[k]) )
2221             {
2222               I[j]=0;
2223               v[j]=1;
2224               break;
2225             }
2226             if( (M==N) && (l[j]<=l[k]) || N/M != 0 )
2227             {
2228               I[k]=0;
2229               v[k]=1;
2230             }
2231          }
2232        }
2233      }
2234   }
2235   return(v);
2236}
2237example
2238{ "EXAMPLE:"; echo = 2;
2239   ring  r = 0,(x,y,z,t),(dp(3),dp);
2240   ideal i = x+y+z+t,xy+yz+xt+zt,xyz+xyt+xzt+yzt,xyzt-t4;
2241   ideal Id = std(i);
2242   ideal I = lead(Id);  I;
2243   ring S = (0,t),(x,y,z),dp;
2244   downsizeSB(imap(r,I));
2245   //Id[5] can be deleted, we still have a SB of i in the ring S
2246
2247   ring R = (0,x),(y,z,u),lp;
2248   ideal i = x+y+z+u,xy+xu+yz+zu,xyz+xyu+xzu+yzu,xyzu-1;
2249   def Phelp = par2varRing()[1];
2250   setring Phelp;
2251   ideal IPhelp = std(imap(R,i));
2252   ideal I = lead(IPhelp);
2253   setring R;
2254   ideal I = imap(Phelp,I); I;
2255   intvec v = downsizeSB(I); v;
2256}
2257///////////////////////////////////////////////////////////////////////////
2258// PROBLEM: Die Prozedur funktioniert nur fuer Ringe die global bekannt
2259//          sind, also interaktiv, aber nicht aus einer Prozedur.
2260//          Z.B. funktioniert example imapDownsize; nicht
2261
2262proc imapDownsize (string R, string I)
2263"SYNTAX: @code{imapDownsize (} string @code{,} string @code{)} *@
2264         First string must be the string of the name of a ring, second
2265         string must be the string of the name of an object in the ring.
2266TYPE:    same type as the object with name the second string
2267PURPOSE: maps the object given by the second string to the basering.
2268         If R resp. I are the first resp. second string, then
2269         imapDownsize(R,I) is equivalent to simplify(imap(`R`,`I`),34).
2270NOTE:    imapDownsize is usually faster than imap if `I` is large and if
2271         simplify has a great effect, since the procedure maps only those
2272         generators from `I` which are not killed by simplify( - ,34).
2273         This is useful if `I` is a standard bases for a block ordering of
2274         `R` and if some variables from the last block in `R` are mapped
2275         to parameters. Then the returned result is a standard basis in
2276         the basering.
2277SEE ALSO: imap, fetch, map
2278EXAMPLE: example imapDownsize; shows an example"
2279{
2280       def BR = basering;
2281       int k;
2282
2283       setring `R`;
2284       def @leadI@ = lead(`I`);
2285       int s = ncols(@leadI@);
2286       setring BR;
2287       ideal @leadI@ = simplify(imap(`R`,@leadI@),32);
2288       intvec vi;
2289       for (k=1; k<=s; k++)
2290       {
2291         vi[k] = @leadI@[k]==0;
2292       }
2293       kill @leadI@;
2294
2295       setring `R`;
2296       kill @leadI@;
2297       for (k=1;  k<=s; k++)
2298       {
2299           if( vi[k]==1 )
2300           {
2301              `I`[k]=0;
2302           }
2303       }
2304       `I` = simplify(`I`,2);
2305
2306       setring BR;
2307       return(imap(`R`,`I`));
2308}
2309example
2310{ "EXAMPLE:"; echo = 2;
2311   ring  r = 0,(x,y,z,t),(dp(3),dp);
2312   ideal i = x+y+z+t,xy+yz+xt+zt,xyz+xyt+xzt+yzt,xyzt-1;
2313   i = std(i); i;
2314
2315   ring s = (0,t),(x,y,z),dp;
2316   imapDownsize("r","i");     //i[5] is omitted since lead(i[2]) | lead(i[5])
2317}
2318///////////////////////////////////////////////////////////////////////////////
2319//die folgende proc war fuer groebner mit fglm vorgesehen, ist aber zu teuer.
2320//Um die projektive Dimension korrekt zu berechnen, muss man aber teuer
2321//voerher ein SB bzgl. einer Gradordnung berechnen und dann homogenisieren.
2322//Sonst koennen hoeherdimensionale Komponenten in Unendlich entstehen
2323
2324proc projInvariants(ideal i,list #)
2325"SYNTAX: @code{projInvariants (} ideal_expression @code{)} @*
2326         @code{projInvariants (} ideal_expression@code{,} list of string_expres          sions@code{)}
2327TYPE:    list, say L, with L[1] and L[2] of type int and L[3] of type intvec
2328PURPOSE: Computes the (projective) dimension (L[1]), degree (L[2]) and the
2329         first Hilbert series (L[3], as intvec) of the homogenized ideal
2330         in the ring given by the procedure @code{hilbRing} with global
2331         ordering dp (resp. wp if the variables have weights >1)
2332         If an argument of type string @code{\"std\"} resp. @code{\"slimgb\"}
2333         is given, the standard basis computatuion uses @code{std} or
2334         @code{slimgb}, otherwise a heuristically chosen method (default)
2335NOTE:    Homogenized means weighted homogenized with respect to the weights
2336         w[i] of the variables var(i) of the basering. The returned dimension,
2337         degree and Hilbertseries are the respective invariants of the
2338         projective variety defined by the homogenized ideal. The dimension
2339         is equal to the (affine) dimension of the ideal in the basering
2340         (degree and Hilbert series make only sense for homogeneous ideals).
2341SEE ALSO: dim, dmult, hilb
2342KEYWORDS: dimension, degree, Hilbert function
2343EXAMPLE: example projInvariants;  shows an example"
2344{
2345  def P = basering;
2346  int p_opt;
2347  string s_opt = option();
2348  if (find(option(), "prot"))  { p_opt = 1; }
2349
2350//---------------- check method and clear denomintors --------------------
2351  int k;
2352  string method;
2353  for (k=1; k<=size(#); k++)
2354  {
2355     if (typeof(#[k]) == "string")
2356     {
2357       method = method + "," + #[k];
2358     }
2359  }
2360
2361  if (npars(P) > 0)             //clear denominators of parameters
2362  {
2363    for( k=ncols(i); k>0; k-- )
2364    {
2365       i[k]=cleardenom(i[k]);
2366    }
2367  }
2368
2369//------------------------ change to hilbRing ----------------------------
2370     list hiRi = hilbRing(i);
2371     intvec W = hiRi[2];
2372     def Philb = hiRi[1];      //note: Philb is no qring and the predefined
2373     setring Philb;            //ideal Id(1) in Philb is homogeneous
2374     int di, de;               //for dimension, degree
2375     intvec hi;                //for hilbert series
2376
2377//-------- compute Hilbert function of homogenized ideal in Philb ---------
2378//Philb has only 1 block. There are three cases
2379
2380     string algorithm;       //possibilities: std, slimgb, stdorslimgb
2381     //define algorithm:
2382     if( find(method,"std") && !find(method,"slimgb") )
2383     {
2384        algorithm = "std";
2385     }
2386     if( find(method,"slimgb") && !find(method,"std") )
2387     {
2388         algorithm = "slimgb";
2389     }
2390     if( find(method,"std") && find(method,"slimgb") ||
2391         (!find(method,"std") && !find(method,"slimgb")) )
2392     {
2393         algorithm = "stdorslimgb";
2394     }
2395
2396     if ( algorithm=="std" || ( algorithm=="stdorslimgb" && char(P)>0 ) )
2397     {
2398        if (p_opt) {"std in ring " + string(Philb);}
2399        Id(1) = std(Id(1));
2400        di = dim(Id(1))-1;
2401        de = mult(Id(1));
2402        hi = hilb( Id(1),1,W );
2403     }
2404     if ( algorithm=="slimgb" || ( algorithm=="stdorslimgb" && char(P)==0 ) )
2405     {
2406        if (p_opt) {"slimgb in ring " + string(Philb);}
2407        Id(1) = slimgb(Id(1));
2408        di = dim( Id(1) );
2409        if (di > -1)
2410        {
2411           di = di-1;
2412        }
2413        de = mult( Id(1) );
2414        hi = hilb( Id(1),1,W );
2415     }
2416     kill Philb;
2417     list L = di,de,hi;
2418     return(L);
2419}
2420example
2421{ "EXAMPLE:"; echo = 2;
2422   ring r = 32003,(x,y,z),lp;
2423   ideal i = y2-xz,x2-z;
2424   projInvariants(i);
2425
2426   ring R = (0),(x,y,z,u,v),lp;
2427   //minpoly = x2+1;
2428   ideal i = x2+1,x2+y+z+u+v,xyzuv-1;
2429   projInvariants(i);
2430   qring S =std(x2+1);
2431   ideal i = imap(R,i);
2432   projInvariants(i);
2433}
2434
2435*/
2436///////////////////////////////////////////////////////////////////////////////
2437//                           EXAMPLES
2438///////////////////////////////////////////////////////////////////////////////
2439/*
2440example stdfglm;
2441example stdhilb;
2442example groebner;
2443example res;
2444example sprintf;
2445example fprintf;
2446example printf;
2447example weightKB;
2448example qslimgb;
2449example par2varRing;
2450*/
2451static proc mod_init()
2452{
2453  //int pagelength=24;
2454  //exportto(Top,pagelength);
2455}
Note: See TracBrowser for help on using the repository browser.