source: git/Singular/LIB/standard.lib @ 7de961

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