# source:git/Singular/LIB/standard.lib@27d4bb

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