source: git/Singular/LIB/standard.lib @ 179432

spielwiese
Last change on this file since 179432 was 179432, checked in by Hans Schönemann <hannes@…>, 7 years ago
version number for standrard.lib
  • Property mode set to 100644
File size: 70.2 KB
Line 
1///////////////////////////////////////////////////////////////////////////
2version="version standard.lib 4.1.0.0 Dec_2016 "; // $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 (attrib(basering,"ring_cf")==1) // 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, mathicgb
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( (was_minpoly == 0) && (npars_P == 0) && (was_qring == 0) && (attrib (P,"global") == 1) && (char(P) > 0) && (size(BRlist)<=4) )
964  {
965    if( defined(Singmathic) )
966    {
967      algorithm = "mathicgb"; // make it default for any appropriate setting... if mathicgb is available...
968    } else
969    {
970      if( p_opt && find(method,"mathicgb")  ) { "Sorry Singmathic::mathicgb is not available!"; }
971    }
972  }
973  if( find(method,"std") && !find(method,"slimgb") )
974  {
975    algorithm = "std";
976  }
977  if( find(method,"slimgb") && !find(method,"std") )
978  {
979    algorithm = "slimgb";
980  }
981  if( find(method,"std") && find(method,"slimgb") ||
982      (!find(method,"std") && !find(method,"slimgb")) )
983  {
984    algorithm = "stdorslimgb";
985  }
986
987  //define conversion:
988  if( find(method,"hilb") && !find(method,"fglm") )
989  {
990     conversion = "hilb";
991  }
992  if( find(method,"fglm") && !find(method,"hilb") )
993  {
994    conversion = "fglm";
995  }
996  if( find(method,"fglm") && find(method,"hilb") )
997  {
998    conversion = "hilborfglm";
999  }
1000  if( !find(method,"fglm") && !find(method,"hilb") )
1001  {
1002    conversion = "no";
1003  }
1004
1005  //define partovar:
1006  //if( find(method,"par2var") && npars_P > 0 )   //V1
1007  if( find(method,"par2var") || npars_P > 1 )     //V2
1008  {
1009     partovar = "yes";
1010  }
1011  else
1012  {
1013     partovar = "no";
1014  }
1015
1016  //define order:
1017  if (system("nblocks") <= 2)
1018  {
1019    if ( find(ordstr_P,"M")+find(ordstr_P,"lp")+find(ordstr_P,"rp") <= 0 )
1020    {
1021      order = "simple";
1022    }
1023  }
1024
1025  //define direct:
1026  if ( (order=="simple" && (size(method)==0)) ||
1027       (size(BRlist)>4) ||
1028        (order=="simple" && (method==",par2var" && npars_P==0 )) ||
1029         (conversion=="no" && partovar=="no" &&
1030           (algorithm=="std" || algorithm=="slimgb" || algorithm=="mathicgb" ||
1031            (find(method,"std") && find(method,"slimgb"))
1032           )
1033         )
1034     )
1035  {
1036    direct = "yes";
1037  }
1038  else
1039  {
1040    direct = "no";
1041  }
1042
1043  //order=="simple" means that the ordering of the variables consists of one
1044  //block which is not a matrix ordering and not a lexicographical ordering.
1045  //(Note:Singular counts always least 2 blocks, one is for module component):
1046  //Call a method "direct" if conversion=="no" && partovar="no" which means
1047  //that we apply std or slimgb dircet in the basering (exception
1048  //as long as slimgb does not know qrings: in a qring of a ring P
1049  //the ideal Qideal is added to the ideal and slimgb is applied in P).
1050  //We apply a direct method if we have a simple monomial ordering, if no
1051  //conversion (fglm or hilb) is specified and if the parameters shall
1052  //not be made to variables
1053  //BRlist (=ringlist of basering) > 4 if the basering is non-commutative
1054//---------------------------- direct methods -----------------------------
1055  if ( algorithm=="mathicgb" )
1056  {
1057    if (p_opt) { algorithm + " in " + string(P); }
1058    return( mathicgb(i) );
1059  }
1060  if ( direct == "yes" )
1061  {
1062  //if ( algorithm=="std" || (algorithm=="stdorslimgb" && char(P)>0) )   //V1
1063    if ( algorithm=="std" || (algorithm=="stdorslimgb") )                //V2
1064    {
1065      if (p_opt) { "std in " + string(P); }
1066      return(std(i));
1067    }
1068  //if( algorithm=="slimgb" || (algorithm=="stdorslimgb" && char(P)==0)) //V1
1069    if ( algorithm=="slimgb" )                                           //V2
1070    {
1071      return(qslimgb(i));
1072    }
1073  }
1074
1075//--------------------------- indirect methods -----------------------------
1076//indirect methods are methods where a conversion is used with a ring change
1077//We are in the following situation:
1078//direct=="no" (i.e. "hilb" or "fglm" or "par2var" is given)
1079//or no method is given and we have a complicated monomial ordering
1080//V1: "par2var" is not a default strategy, it must be explicitely
1081//given in order to be performed.
1082//V2: "par2var" is a default strategy if there are more than 1 parameters
1083
1084//------------ case where no parameters are made to variables  -------------
1085  if (  partovar == "no" && conversion == "hilb"
1086    || (partovar == "no" && conversion == "fglm" )
1087    || (partovar == "no" && conversion == "hilborfglm" )
1088    || (partovar == "no" && conversion == "no" && direct == "no") )
1089  //last case: heuristic
1090  {
1091    if ( conversion=="fglm" )
1092    {
1093    //if ( algorithm=="std" || (algorithm=="stdorslimgb" && char(P)>0) ) //V1
1094      if ( algorithm=="std" || (algorithm=="stdorslimgb") )              //V2
1095      {
1096        return (stdfglm(i,"std"));
1097      }
1098    //if(algorithm=="slimgb" || (algorithm=="stdorslimgb" && char(P)==0))//V1
1099      if( algorithm=="slimgb" )                                          //V2
1100      {
1101        return (stdfglm(i,"slimgb"));
1102      }
1103    }
1104    else
1105    {
1106    //if ( algorithm=="std" || (algorithm=="stdorslimgb" && char(P)>0) )//V1
1107      if ( algorithm=="std" || (algorithm=="stdorslimgb" ) )            //V2
1108      {
1109        return (stdhilb(i,"std"));
1110      }
1111    //if(algorithm=="slimgb" || (algorithm=="stdorslimgb" && char(P)==0))//V1
1112      if ( algorithm=="slimgb" )                                         //V2
1113      {
1114        return (stdhilb(i,"slimgb"));
1115      }
1116    }
1117  }
1118
1119//------------ case where parameters are made to variables  ----------------
1120//define a ring Phelp via par2varRing in which the parameters are variables
1121
1122  else
1123  {
1124    // reset options
1125    option(none);
1126    // turn on options prot, mem, redSB, intStrategy if previously set
1127    if ( find(s_opt, "prot") )
1128      { option(prot); }
1129    if ( find(s_opt, "mem") )
1130      { option(mem); }
1131    if ( find(s_opt, "redSB") )
1132      { option(redSB); }
1133    if ( find(s_opt, "intStrategy") )
1134      { option(intStrategy); }
1135
1136    //first clear denominators of parameters
1137    if (npars_P > 0)
1138    {
1139      for( k=ncols(i); k>0; k-- )
1140      { i[k]=cleardenom(i[k]); }
1141    }
1142
1143    def Phelp = par2varRing(i)[1];   //minpoly is mapped with i
1144    setring Phelp;
1145    def i = Id(1);
1146    //is_homog = homog(i);
1147
1148    //If parameters are converted to ring variables, they appear in an extra
1149    //block. Therefore we use always hilb for this block ordering:
1150    if ( conversion=="fglm" )
1151    {
1152      i = (stdfglm(i));       //only uesful for 1 parameter with minpoly
1153    }
1154    else
1155    {
1156    //if ( algorithm=="std" || (algorithm=="stdorslimgb" && char(P)>0) )//V1
1157      if ( algorithm=="std" || (algorithm=="stdorslimgb" ))             //V2
1158      {
1159        i = stdhilb(i,"std");
1160      }
1161    //if(algorithm=="slimgb" || (algorithm=="stdorslimgb" && char(P)==0))//V1
1162      if ( algorithm=="slimgb" )                                         //V2
1163      {
1164        i = stdhilb(i,"slimgb");
1165      }
1166    }
1167  }
1168
1169//-------------------- go back to original ring ---------------------------
1170//The main computation is done. However, the SB coming from a ring with
1171//extra variables is in general too big. We simplify it before mapping it
1172//to the basering.
1173
1174  if (p_opt) { "//simplification"; }
1175
1176  if (was_minpoly)
1177  {
1178    execute("ideal Minpoly = " + @Minpoly + ";");
1179    attrib(Minpoly,"isSB",1);
1180    i = simplify(NF(i,Minpoly),2);
1181  }
1182
1183  def Li = lead(i);
1184  setring P;
1185  def Li = imap(Phelp,Li);
1186  Li = simplify(Li,32);
1187  intvec vi;
1188  for (k=1; k<=ncols(Li); k++)
1189  {
1190    vi[k] = Li[k]==0;
1191  }
1192
1193  setring Phelp;
1194  for (k=1;  k<=size(i) ;k++)
1195  {
1196    if(vi[k]==1)
1197    {
1198      i[k]=0;
1199    }
1200  }
1201  i = simplify(i,2);
1202
1203  setring P;
1204  if (p_opt) { "//imap to original ring"; }
1205  i = imap(Phelp,i);
1206  kill Phelp;
1207  i = simplify(i,34);
1208
1209  // clean-up time
1210  option(set, opt);
1211  if (find(s_opt, "redSB") > 0)
1212  {
1213    if (p_opt) { "//interreduction"; }
1214    i=interred(i);
1215  }
1216  attrib(i, "isSB", 1);
1217  return (i);
1218}
1219example
1220{ "EXAMPLE: "; echo=2;
1221  intvec opt = option(get);
1222  option(prot);
1223  ring r  = 0,(a,b,c,d),dp;
1224  ideal i = a+b+c+d,ab+ad+bc+cd,abc+abd+acd+bcd,abcd-1;
1225  groebner(i);
1226
1227  ring s  = 0,(a,b,c,d),lp;
1228  ideal i = imap(r,i);
1229  groebner(i,"hilb");
1230
1231  ring R  = (0,a),(b,c,d),lp;
1232  minpoly = a2+1;
1233  ideal i = a+b+c+d,ab+ad+bc+cd,abc+abd+acd+bcd,d2-c2b2;
1234  groebner(i,"par2var","slimgb");
1235
1236  groebner(i,"fglm");          //computes a reduced standard basis
1237
1238  option(set,opt);
1239}
1240
1241//////////////////////////////////////////////////////////////////////////
1242
1243proc res(list #)
1244"@c we do texinfo here:
1245@cindex resolution, computation of
1246@table @code
1247@item @strong{Syntax:}
1248@code{res (} ideal_expression@code{,} int_expression @code{[,} any_expression @code{])}
1249@*@code{res (} module_expression@code{,} int_expression @code{[,} any_expression @code{])}
1250@item @strong{Type:}
1251resolution
1252@item @strong{Purpose:}
1253computes a (possibly minimal) free resolution of an ideal or module using
1254a heuristically chosen method.
1255@* The second (int) argument (say @code{k}) specifies the length of
1256the resolution. If it is not positive then @code{k} is assumed to be the
1257number of variables of the basering.
1258@* If a third argument is given, the returned resolution is minimized.
1259
1260Depending on the input, the returned resolution is computed using the
1261following methods:
1262@table @asis
1263@item @strong{quotient rings:}
1264@code{nres} (classical method using syzygies) , see @ref{nres}.
1265
1266@item @strong{homogeneous ideals and k=0:}
1267@code{lres} (La'Scala's method), see @ref{lres}.
1268
1269@item @strong{not minimized resolution and (homogeneous input with k not 0, or local rings):}
1270@code{sres} (Schreyer's method), see @ref{sres}.
1271
1272@item @strong{all other inputs:}
1273@code{mres} (classical method), see @ref{mres}.
1274@end table
1275@item @strong{Note:}
1276Accessing single elements of a resolution may require some partial
1277computations to be finished and may therefore take some time.
1278@end table
1279@c ref
1280See also
1281@ref{betti};
1282@ref{ideal};
1283@ref{minres};
1284@ref{module};
1285@ref{mres};
1286@ref{nres};
1287@ref{lres};
1288@ref{hres};
1289@ref{sres};
1290@ref{resolution}.
1291@c ref
1292"
1293{
1294   def P=basering;
1295   if (size(#) < 2)
1296   {
1297     ERROR("res: need at least two arguments: ideal/module, int");
1298   }
1299
1300   def m=#[1]; //the ideal or module
1301   int i=#[2]; //the length of the resolution
1302   if (i< 0) { i=0;}
1303
1304   string varstr_P = varstr(P);
1305
1306   int p_opt;
1307   string s_opt = option();
1308   // set p_opt, if option(prot) is set
1309   if (find(s_opt, "prot"))
1310   {
1311     p_opt = 1;
1312   }
1313
1314   if( (size(ideal(basering)) > 0) || (size(ringlist(P)) > 4) )
1315   {
1316     // the quick hack for qrings - seems to fit most needs
1317     // (lres is not implemented for qrings, sres is not so efficient)
1318     // || non-commutative, since only n/m-res are implemented for NC rings
1319     if (p_opt) { "using nres";}
1320     return(nres(m,i));
1321   }
1322
1323/* if( attrib(basering, "global") == 1 ) // preparations for s_res usage. in testing!
1324   {
1325     if (p_opt) { "using s_res";}
1326     if( !defined(s_res) )
1327     {
1328       def @@v=option(get); option(noloadLib); option(noloadProc); LIB( "schreyer.lib" ); // for s_res
1329       option(set, @@v); kill @@v;
1330     }
1331     resolution re = s_res(m,i);
1332     if(size(#)>2)
1333     {
1334       re=minres(re);
1335     }
1336     return(re);
1337   }*/
1338
1339   if(homog(m)==1)
1340   {
1341      resolution re;
1342      if (((i==0) or (i>=nvars(basering))) && (typeof(m) != "module") && (nvars(basering)>1))
1343      {
1344        //LaScala for the homogeneous case and i == 0
1345        if (p_opt) { "using lres";}
1346        re=lres(m,i);
1347        if(size(#)>2)
1348        {
1349           re=minres(re);
1350        }
1351      }
1352      else
1353      {
1354        if(size(#)>2)
1355        {
1356          if (p_opt) { "using mres";}
1357          re=mres(m,i);
1358        }
1359        else
1360        {
1361          if (p_opt) { "using sres";}
1362          re=sres(std(m),i);
1363        }
1364      }
1365      return(re);
1366   }
1367
1368   //mres for the global non homogeneous case
1369   if(find(ordstr(P),"s")==0)
1370   {
1371      string ri= "ring Phelp ="
1372                  +string(char(P))+",("+varstr_P+"),(dp,C);";
1373      ri = ri + "minpoly = "+string(minpoly) + ";";
1374      execute(ri);
1375      def m=imap(P,m);
1376      if (p_opt) { "using mres in another ring";}
1377      list re=mres(m,i);
1378      setring P;
1379      resolution result=imap(Phelp,re);
1380      if (size(#) > 2) {result = minres(result);}
1381      return(result);
1382   }
1383
1384   //sres for the local case and not minimal resolution
1385   if(size(#)<=2)
1386   {
1387      string ri= "ring Phelp ="
1388                  +string(char(P))+",("+varstr_P+"),(ls,c);";
1389      ri = ri + "minpoly = "+string(minpoly) + ";";
1390      execute(ri);
1391      def m=imap(P,m);
1392      m=std(m);
1393      if (p_opt) { "using sres in another ring";}
1394      list re=sres(m,i);
1395      setring P;
1396      resolution result=imap(Phelp,re);
1397      return(result);
1398   }
1399
1400   //mres for the local case and minimal resolution
1401   string ri= "ring Phelp ="
1402                  +string(char(P))+",("+varstr_P+"),(ls,C);";
1403   ri = ri + "minpoly = "+string(minpoly) + ";";
1404   execute(ri);
1405   def m=imap(P,m);
1406    if (p_opt) { "using mres in another ring";}
1407   list re=mres(m,i);
1408   setring P;
1409   resolution result=imap(Phelp,re);
1410   result = minres(result);
1411   return(result);
1412}
1413example
1414{"EXAMPLE:"; echo = 2;
1415  ring r=0,(x,y,z),dp;
1416  ideal i=xz,yz,x3-y3;
1417  def l=res(i,0); // homogeneous ideal: uses lres
1418  l;
1419  print(betti(l), "betti"); // input to betti may be of type resolution
1420  l[2];         // element access may take some time
1421  i=i,x+1;
1422  l=res(i,0);   // inhomogeneous ideal: uses mres
1423  l;
1424  ring rs=0,(x,y,z),ds;
1425  ideal i=imap(r,i);
1426  def l=res(i,0); // local ring not minimized: uses sres
1427  l;
1428  res(i,0,0);     // local ring and minimized: uses mres
1429}
1430/////////////////////////////////////////////////////////////////////////
1431
1432proc quot (def m1,def m2,list #)
1433"SYNTAX: @code{quot (} module_expression@code{,} module_expression @code{)}
1434         @*@code{quot (} module_expression@code{,} module_expression@code{,}
1435            int_expression @code{)}
1436         @*@code{quot (} ideal_expression@code{,} ideal_expression @code{)}
1437         @*@code{quot (} ideal_expression@code{,} ideal_expression@code{,}
1438            int_expression @code{)}
1439TYPE:    ideal
1440SYNTAX:  @code{quot (} module_expression@code{,} ideal_expression @code{)}
1441TYPE:    module
1442PURPOSE: computes the quotient of the 1st and the 2nd argument.
1443         If a 3rd argument @code{n} is given the @code{n}-th method is used
1444         (@code{n}=1...5).
1445SEE ALSO: quotient
1446EXAMPLE: example quot; shows an example"
1447{
1448  if (((typeof(m1)!="ideal") and (typeof(m1)!="module"))
1449     or ((typeof(m2)!="ideal") and (typeof(m2)!="module")))
1450  {
1451    "USAGE:   quot(m1, m2[, n]); m1, m2 two submodules of k^s,";
1452    "         n (optional) integer (1<= n <=5)";
1453    "RETURN:  the quotient of m1 and m2";
1454    "EXAMPLE: example quot; shows an example";
1455    return();
1456  }
1457  if (typeof(m1)!=typeof(m2))
1458  {
1459    return(quotient(m1,m2));
1460  }
1461  if (size(#)>0)
1462  {
1463    if (typeof(#[1])=="int" )
1464    {
1465      return(quot1(m1,m2,#[1]));
1466    }
1467  }
1468  else
1469  {
1470    return(quot1(m1,m2,2));
1471  }
1472}
1473example
1474{ "EXAMPLE:"; echo = 2;
1475  ring r=181,(x,y,z),(c,ls);
1476  ideal id1=maxideal(4);
1477  ideal id2=x2+xyz,y2-z3y,z3+y5xz;
1478  option(prot);
1479  ideal id3=quotient(id1,id2);
1480  id3;
1481  ideal id4=quot(id1,id2,1);
1482  id4;
1483  ideal id5=quot(id1,id2,2);
1484  id5;
1485}
1486
1487static proc quot1 (module m1, module m2,int n)
1488"USAGE:   quot1(m1, m2, n); m1, m2 two submodules of k^s,
1489         n integer (1<= n <=5)
1490RETURN:  the quotient of m1 and m2
1491EXAMPLE: example quot1; shows an example"
1492{
1493  if (n==1)
1494  {
1495    return(quotient1(m1,m2));
1496  }
1497  else
1498  {
1499    if (n==2)
1500    {
1501      return(quotient2(m1,m2));
1502    }
1503    else
1504    {
1505      if (n==3)
1506      {
1507        return(quotient3(m1,m2));
1508      }
1509      else
1510      {
1511        if (n==4)
1512        {
1513          return(quotient4(m1,m2));
1514        }
1515        else
1516        {
1517          if (n==5)
1518          {
1519            return(quotient5(m1,m2));
1520          }
1521          else
1522          {
1523            return(quotient(m1,m2));
1524          }
1525        }
1526      }
1527    }
1528  }
1529}
1530example
1531{ "EXAMPLE:"; echo = 2;
1532  ring r=181,(x,y,z),(c,ls);
1533  ideal id1=maxideal(4);
1534  ideal id2=x2+xyz,y2-z3y,z3+y5xz;
1535  option(prot);
1536  ideal id6=quotient(id1,id2);
1537  id6;
1538  ideal id7=quot1(id1,id2,1);
1539  id7;
1540  ideal id8=quot1(id1,id2,2);
1541  id8;
1542}
1543
1544static proc quotient0(module a,module b)
1545{
1546  module mm=b+a;
1547  resolution rs=lres(mm,0);
1548  list I=list(rs);
1549  matrix M=I[2];
1550  matrix A[1][nrows(M)]=M[1..nrows(M),1];
1551  ideal i=A;
1552  return (i);
1553}
1554proc quotient1(module a,module b)  //17sec
1555"USAGE:   quotient1(m1, m2); m1, m2 two submodules of k^s,
1556RETURN:  the quotient of m1 and m2"
1557{
1558  int i;
1559  a=std(a);
1560  module dummy;
1561  module B=NF(b,a)+dummy;
1562  ideal re=quotient(a,module(B[1]));
1563  for(i=2;i<=ncols(B);i++)
1564  {
1565     re=intersect1(re,quotient(a,module(B[i])));
1566  }
1567  return(re);
1568}
1569proc quotient2(module a,module b)    //13sec
1570"USAGE:   quotient2(m1, m2); m1, m2 two submodules of k^s,
1571RETURN:  the quotient of m1 and m2"
1572{
1573  a=std(a);
1574  module dummy;
1575  module bb=NF(b,a)+dummy;
1576  int i=ncols(bb);
1577  ideal re=quotient(a,module(bb[i]));
1578  bb[i]=0;
1579  module temp;
1580  module temp1;
1581  module bbb;
1582  int mx;
1583  i=i-1;
1584  while (1)
1585  {
1586    if (i==0) break;
1587    temp = a+bb*re;
1588    temp1 = lead(interred(temp));
1589    mx=ncols(a);
1590    if (ncols(temp1)>ncols(a))
1591    {
1592      mx=ncols(temp1);
1593    }
1594    temp1 = matrix(temp1,1,mx)-matrix(lead(a),1,mx);
1595    temp1 = dummy+temp1;
1596    if (deg(temp1[1])<0) break;
1597    re=intersect1(re,quotient(a,module(bb[i])));
1598    bb[i]=0;
1599    i = i-1;
1600  }
1601  return(re);
1602}
1603proc quotient3(module a,module b)   //89sec
1604"USAGE:   quotient3(m1, m2); m1, m2 two submodules of k^s,
1605         only for global rings
1606RETURN:  the quotient of m1 and m2"
1607{
1608  string s="ring @newr=("+charstr(basering)+
1609           "),("+varstr(basering)+",@t,@w),dp;";
1610  def @newP=basering;
1611  execute(s);
1612  module b=imap(@newP,b);
1613  module a=imap(@newP,a);
1614  int i;
1615  int j=ncols(b);
1616  vector @b;
1617  for(i=1;i<=j;i++)
1618  {
1619     @b=@b+@t^(i-1)*@w^(j-i+1)*b[i];
1620  }
1621  ideal re=quotient(a,module(@b));
1622  setring @newP;
1623  ideal re=imap(@newr,re);
1624  return(re);
1625}
1626proc quotient5(module a,module b)   //89sec
1627"USAGE:   quotient5(m1, m2); m1, m2 two submodules of k^s,
1628         only for global rings
1629RETURN:  the quotient of m1 and m2"
1630{
1631  string s="ring @newr=("+charstr(basering)+
1632           "),("+varstr(basering)+",@t),dp;";
1633  def @newP=basering;
1634  execute(s);
1635  module b=imap(@newP,b);
1636  module a=imap(@newP,a);
1637  int i;
1638  int j=ncols(b);
1639  vector @b;
1640  for(i=1;i<=j;i++)
1641  {
1642     @b=@b+@t^(i-1)*b[i];
1643  }
1644  @b=homog(@b,@w);
1645  ideal re=quotient(a,module(@b));
1646  setring @newP;
1647  ideal re=imap(@newr,re);
1648  return(re);
1649}
1650proc quotient4(module a,module b)   //95sec
1651"USAGE:   quotient4(m1, m2); m1, m2 two submodules of k^s,
1652         only for global rings
1653RETURN:  the quotient of m1 and m2"
1654{
1655  string s="ring @newr=("+charstr(basering)+
1656           "),("+varstr(basering)+",@t),dp;";
1657  def @newP=basering;
1658  execute(s);
1659  module b=imap(@newP,b);
1660  module a=imap(@newP,a);
1661  int i;
1662  vector @b=b[1];
1663  for(i=2;i<=ncols(b);i++)
1664  {
1665     @b=@b+@t^(i-1)*b[i];
1666  }
1667  matrix sy=modulo(@b,a);
1668  ideal re=sy;
1669  setring @newP;
1670  ideal re=imap(@newr,re);
1671  return(re);
1672}
1673static proc intersect1(ideal i,ideal j)
1674{
1675  def R=basering;
1676  execute("ring gnir = ("+charstr(basering)+"),
1677                       ("+varstr(basering)+",@t),(C,dp);");
1678  ideal i=var(nvars(basering))*imap(R,i)+(var(nvars(basering))-1)*imap(R,j);
1679  ideal j=eliminate(i,var(nvars(basering)));
1680  setring R;
1681  map phi=gnir,maxideal(1);
1682  return(phi(j));
1683}
1684
1685//////////////////////////////////////////////////////////////////
1686///
1687/// sprintf, fprintf printf
1688///
1689proc sprintf(string fmt, list #)
1690"SYNTAX:  @code{sprintf (} string_expression @code{[,} any_expressions
1691               @code{] )}
1692RETURN:   string
1693PURPOSE:  @code{sprintf(fmt,...);} performs output formatting. The first
1694          argument is a format control string. Additional arguments may be
1695          required, depending on the content of the control string. A series
1696          of output characters is generated as directed by the control string;
1697          these characters are returned as a string. @*
1698          The control string @code{fmt} is simply text to be copied,
1699          except that the string may contain conversion specifications.@*
1700          Type @code{help print;} for a listing of valid conversion
1701          specifications. As an addition to the conversions of @code{print},
1702          the @code{%n} and @code{%2} conversion specification does not
1703          consume an additional argument, but simply generates a newline
1704          character.
1705NOTE:     If one of the additional arguments is a list, then it should be
1706          wrapped in an additional @code{list()} command, since passing a list
1707          as an argument flattens the list by one level.
1708SEE ALSO: fprintf, printf, print, string
1709EXAMPLE : example sprintf; shows an example
1710"
1711{
1712  int sfmt = size(fmt);
1713  if (sfmt  <= 1)
1714  {
1715    return (fmt);
1716  }
1717  int next, l, nnext;
1718  string ret;
1719  list formats = "%l", "%s", "%2l", "%2s", "%t", "%;", "%p", "%b", "%n", "%2";
1720  while (1)
1721  {
1722    if (size(#) <= 0)
1723    {
1724      return (ret + fmt);
1725    }
1726    nnext = 0;
1727    while (nnext < sfmt)
1728    {
1729      nnext = find(fmt, "%", nnext + 1);
1730      if (nnext == 0)
1731      {
1732        next = 0;
1733        break;
1734      }
1735      l = 1;
1736      while (l <= size(formats))
1737      {
1738        next = find(fmt, formats[l], nnext);
1739        if (next == nnext) break;
1740        l++;
1741      }
1742      if (next == nnext) break;
1743    }
1744    if (next == 0)
1745    {
1746      return (ret + fmt);
1747    }
1748    if (formats[l] != "%2" && formats[l] != "%n")
1749    {
1750      ret = ret + fmt[1, next - 1] + print(#[1], formats[l]);
1751      # = delete(#, 1);
1752    }
1753    else
1754    {
1755      ret = ret + fmt[1, next - 1] + print("", "%2s");
1756    }
1757    if (size(fmt) <= (next + size(formats[l]) - 1))
1758    {
1759      return (ret);
1760    }
1761    fmt = fmt[next + size(formats[l]), size(fmt)-next-size(formats[l]) + 1];
1762  }
1763}
1764example
1765{ "EXAMPLE:"; echo=2;
1766  ring r=0,(x,y,z),dp;
1767  module m=[1,y],[0,x+z];
1768  intmat M=betti(mres(m,0));
1769  list l = r, m, M;
1770  string s = sprintf("s:%s,%n l:%l", 1, 2); s;
1771  s = sprintf("s:%n%s", l); s;
1772  s = sprintf("s:%2%s", list(l)); s;
1773  s = sprintf("2l:%n%2l", list(l)); s;
1774  s = sprintf("%p", list(l)); s;
1775  s = sprintf("%;", list(l)); s;
1776  s = sprintf("%b", M); s;
1777}
1778
1779proc printf(string fmt, list #)
1780"SYNTAX:  @code{printf (} string_expression @code{[,} any_expressions@code{] )}
1781RETURN:   none
1782PURPOSE:  @code{printf(fmt,...);} performs output formatting. The first
1783          argument is a format control string. Additional arguments may be
1784          required, depending on the content of the control string. A series
1785          of output characters is generated as directed by the control string;
1786          these characters are displayed (i.e., printed to standard out). @*
1787          The control string @code{fmt} is simply text to be copied, except
1788          that the string may contain conversion specifications. @*
1789          Type @code{help print;} for a listing of valid conversion
1790          specifications. As an addition to the conversions of @code{print},
1791          the @code{%n} and @code{%2} conversion specification does not
1792          consume an additional argument, but simply generates a newline
1793          character.
1794NOTE:     If one of the additional arguments is a list, then it should be
1795          enclosed once more into a @code{list()} command, since passing a
1796          list as an argument flattens the list by one level.
1797SEE ALSO: sprintf, fprintf, print, string
1798EXAMPLE : example printf; shows an example
1799"
1800{
1801  write("", sprintf(fmt, #));
1802}
1803example
1804{ "EXAMPLE:"; echo=2;
1805  ring r=0,(x,y,z),dp;
1806  module m=[1,y],[0,x+z];
1807  intmat M=betti(mres(m,0));
1808  list l=r,m,matrix(M);
1809  printf("s:%s,l:%l",1,2);
1810  printf("s:%s",l);
1811  printf("s:%s",list(l));
1812  printf("2l:%2l",list(l));
1813  printf("%p",matrix(M));
1814  printf("%;",matrix(M));
1815  printf("%b",M);
1816}
1817
1818
1819proc fprintf(link l, string fmt, list #)
1820"SYNTAX:  @code{fprintf (} link_expression@code{,} string_expression @code{[,}
1821                any_expressions@code{] )}
1822RETURN:   none
1823PURPOSE:  @code{fprintf(l,fmt,...);} performs output formatting.
1824          The second argument is a format control string. Additional
1825          arguments may be required, depending on the content of the
1826          control string. A series of output characters is generated as
1827          directed by the control string; these characters are
1828          written to the link l.
1829          The control string @code{fmt} is simply text to be copied, except
1830          that the string may contain conversion specifications.@*
1831          Type @code{help print;} for a listing of valid conversion
1832          specifications. As an addition to the conversions of @code{print},
1833          the @code{%n} and @code{%2} conversion specification does not
1834          consume an additional argument, but simply generates a newline
1835          character.
1836NOTE:     If one of the additional arguments is a list, then it should be
1837          enclosed once more into a @code{list()} command, since passing
1838          a list as an argument flattens the list by one level.
1839SEE ALSO: sprintf, printf, print, string
1840EXAMPLE : example fprintf; shows an example
1841"
1842{
1843  write(l, sprintf(fmt, #));
1844}
1845example
1846{ "EXAMPLE:"; echo=2;
1847  ring r=0,(x,y,z),dp;
1848  module m=[1,y],[0,x+z];
1849  intmat M=betti(mres(m,0));
1850  list l=r,m,M;
1851  link li="";   // link to stdout
1852  fprintf(li,"s:%s,l:%l",1,2);
1853  fprintf(li,"s:%s",l);
1854  fprintf(li,"s:%s",list(l));
1855  fprintf(li,"2l:%2l",list(l));
1856  fprintf(li,"%p",list(l));
1857  fprintf(li,"%;",list(l));
1858  fprintf(li,"%b",M);
1859}
1860
1861//////////////////////////////////////////////////////////////////////////
1862
1863/*
1864proc minres(list #)
1865{
1866  if (size(#) == 2)
1867  {
1868    if (typeof(#[1]) == "ideal" || typeof(#[1]) == "module")
1869    {
1870      if (typeof(#[2] == "int"))
1871      {
1872        return (res(#[1],#[2],1));
1873      }
1874    }
1875  }
1876
1877  if (typeof(#[1]) == "resolution")
1878  {
1879    return minimizeres(#[1]);
1880  }
1881  else
1882  {
1883    return minimizeres(#);
1884  }
1885
1886}
1887*/
1888///////////////////////////////////////////////////////////////////////////////
1889
1890proc weightKB(def stc, int dd, list wim)
1891"SYNTAX: @code{weightKB (} module_expression@code{,} int_expression @code{,}
1892            list_expression @code{)}@*
1893         @code{weightKB (} ideal_expression@code{,} int_expression@code{,}
1894            list_expression @code{)}
1895RETURN:  the same as the input type of the first argument
1896PURPOSE: If @code{I,d,wim} denotes the three arguments then weightKB
1897         computes the weighted degree- @code{d} part of a vector space basis
1898         (consisting of monomials) of the quotient ring, resp. of the
1899         quotient module, modulo @code{I} w.r.t. weights given by @code{wim}
1900         The information about the weights is given as a list of two intvec:
1901            @code{wim[1]} weights for all variables (positive),
1902            @code{wim[2]} weights for the module generators.
1903NOTE:    This is a generalization of the command @code{kbase} with the same
1904         first two arguments.
1905SEE ALSO: kbase
1906EXAMPLE: example weightKB; shows an example
1907"
1908{
1909  if(checkww(wim)){ERROR("wrong weights";);}
1910  kbclass();
1911  wwtop=wim[1];
1912  stc=interred(lead(stc));
1913  if(typeof(stc)=="ideal")
1914  {
1915    stdtop=stc;
1916    ideal out=widkbase(dd);
1917    delkbclass();
1918    out=simplify(out,2); // delete 0
1919    return(out);
1920  }
1921  list mbase=kbprepare(stc);
1922  module mout;
1923  int im,ii;
1924  if(size(wim)>1){mmtop=wim[2];}
1925  else{mmtop=0;}
1926  for(im=size(mbase);im>0;im--)
1927  {
1928    stdtop=mbase[im];
1929    if(im>size(mmtop)){ii=dd;}
1930    else{ii=dd-mmtop[im];}
1931    mout=mout+widkbase(ii)*gen(im);
1932  }
1933  delkbclass();
1934  mout=simplify(mout,2); // delete 0
1935  return(mout);
1936}
1937example
1938{ "EXAMPLE:"; echo=2;
1939  ring R=0, (x,y), wp(1,2);
1940  weightKB(ideal(0),3,intvec(1,2));
1941}
1942
1943///////////////////////////////////////////////////////////////////////////////
1944
1945proc datetime()
1946"SYNTAX: @code{datetime ()}
1947RETURN:  string
1948PURPOSE: return the curent date and time as a string
1949EXAMPLE: example datetime; shows an example
1950"
1951{
1952  return(read("|: date"));
1953}
1954example
1955{ "EXAMPLE:"; echo=2;
1956  datetime();
1957}
1958
1959///////////////////////////////////////////////////////////////////////////////
1960// construct global values
1961static proc kbclass()
1962{
1963  intvec wwtop,mmtop;
1964  export (wwtop,mmtop);
1965  ideal stdtop,kbtop;
1966  export (stdtop,kbtop);
1967}
1968// delete global values
1969static proc delkbclass()
1970{
1971  kill wwtop,mmtop;
1972  kill stdtop,kbtop;
1973}
1974//  select parts of the modul
1975static proc kbprepare(module mstc)
1976{
1977  list rr;
1978  ideal kk;
1979  int i1,i2;
1980  mstc=transpose(mstc);
1981  for(i1=ncols(mstc);i1>0;i1--)
1982  {
1983    kk=0;
1984    for(i2=nrows(mstc[i1]);i2>0;i2--)
1985    {
1986      kk=kk+mstc[i1][i2];
1987    }
1988    rr[i1]=kk;
1989  }
1990  return(rr);
1991}
1992//  check for weights
1993static proc checkww(list vv)
1994{
1995  if(typeof(vv[1])!="intvec"){return(1);}
1996  intvec ww=vv[1];
1997  int mv=nvars(basering);
1998  if(size(ww)<mv){return(1);}
1999  while(mv>0)
2000  {
2001    if(ww[mv]<=0){return(1);}
2002    mv--;
2003  }
2004  if(size(vv)>1)
2005  {
2006    if(typeof(vv[2])!="intvec"){return(1);}
2007  }
2008  return(0);
2009}
2010///////////////////////////////////////////////////////
2011// The "Caller" for ideals
2012//    dd   - the degree of the result
2013static proc widkbase(int dd)
2014{
2015  if((size(stdtop)==1)&&(deg(stdtop[1])==0)){return(0);}
2016  if(dd<=0)
2017  {
2018    if(dd<0){return(0);}
2019    else{return(1);}
2020  }
2021  int m1,m2;
2022  m1=nvars(basering);
2023  while(wwtop[m1]>dd)
2024  {
2025    m1--;
2026    if(m1==0){return(0);}
2027  }
2028  attrib(stdtop,"isSB",1);
2029  poly mo=1;
2030  if(m1==1)
2031  {
2032    m2=dd div wwtop[1];
2033    if((m2*wwtop[1])==dd)
2034    {
2035      mo=var(1)^m2;
2036      if(reduce(mo,stdtop)==mo){return(mo);}
2037      else{return(0);}
2038    }
2039  }
2040  kbtop=0;
2041  m2=dd;
2042  weightmon(m1-1,m2,mo);
2043  while(m2>=wwtop[m1])
2044  {
2045    m2=m2-wwtop[m1];
2046    mo=var(m1)*mo;
2047    if(m2==0)
2048    {
2049      if((mo!=0) and (reduce(mo,stdtop)==mo))
2050      {
2051        kbtop[ncols(kbtop)+1]=mo;
2052        return(kbtop);
2053      }
2054    }
2055    weightmon(m1-1,m2,mo);
2056  }
2057  return(kbtop);
2058}
2059/////////////////////////////////////////////////////////
2060// the recursive procedure
2061//    va     - number of the variable
2062//    drest  - rest of the degree
2063//    mm     - the candidate
2064static proc weightmon(int va, int drest, poly mm)
2065{
2066  if(va==0){return();}
2067  while(wwtop[va]>drest)
2068  {
2069    va--;
2070    if(va==0){return();}
2071  }
2072  int m2;
2073  if(va==1)
2074  {
2075    m2=drest div wwtop[1];
2076    if((m2*wwtop[1])==drest)
2077    {
2078      mm=var(1)^m2*mm;
2079      if ((mm!=0) and (reduce(mm,stdtop)==mm))
2080      {
2081        kbtop[ncols(kbtop)+1]=mm;
2082      }
2083    }
2084    return();
2085  }
2086  m2=drest;
2087  if ((mm!=0) and (reduce(mm,stdtop)==mm))
2088  {
2089    weightmon(va-1,m2,mm);
2090  }
2091  while(m2>=wwtop[va])
2092  {
2093    m2=m2-wwtop[va];
2094    mm=var(va)*mm;
2095    if(m2==0)
2096    {
2097      if ((mm!=0) and (reduce(mm,stdtop)==mm))
2098      {
2099        kbtop[ncols(kbtop)+1]=mm;
2100        return();
2101      }
2102    }
2103     if ((mm!=0) and (reduce(mm,stdtop)==mm))
2104     {
2105       weightmon(va-1,m2,mm);
2106     }
2107  }
2108  return();
2109}
2110example
2111{ "EXAMPLE:"; echo=2;
2112  ring r=0,(x,y,z),dp;
2113  ideal i = x6,y4,xyz;
2114  intvec w = 2,3,6;
2115  weightKB(i, 12, list(w));
2116}
2117
2118///////////////////////////////////////////////////////////////////////////////
2119proc max(def i,list #)
2120"SYNTAX: max (i_1, ..., i_k)
2121TYPE:    same as type of i_1, ..., i_k resp.
2122PURPOSE: returns the maximum for any arguments of a type
2123         for which '>' is defined
2124SEE ALSO: min
2125EXAMPLE: example max; shows an example"
2126{
2127  def maximum = i;
2128  for (int j=1; j<=size(#); j++)
2129  {
2130    if(#[j]>maximum)
2131    {
2132      maximum = #[j];
2133    }
2134  }
2135  return(maximum);
2136}
2137example
2138{ "EXAMPLE:"; echo=2;
2139  // biggest int
2140  max(2,3);
2141  max(1,4,3);
2142  // lexicographically biggest intvec
2143  max(intvec(1,2),intvec(0,1),intvec(1,1));
2144  // polynopmial with biggest leading monomial
2145  ring r = 0,x,dp;
2146  max(x+1,x2+x);
2147}
2148///////////////////////////////////////////////////////////////////////////////
2149proc min(def i,list #)
2150"SYNTAX: min (i_1, ..., i_k)
2151TYPE:    same as type of i_1, ..., i_k resp.
2152PURPOSE: returns the minimum for any arguments of a type
2153         for which '>' is defined
2154SEE ALSO: max
2155EXAMPLE: example min; shows an example"
2156{
2157  def minimum = i;
2158  for (int j=1; j<=size(#); j++)
2159  {
2160    if(#[j]<minimum)
2161    {
2162      minimum = #[j];
2163    }
2164  }
2165  return(minimum);
2166}
2167example
2168{ "EXAMPLE:"; echo=2;
2169  // smallest int
2170  min(2,3);
2171  min(1,4,3);
2172  // lexicographically smallest intvec
2173  min(intvec(1,2),intvec(0,1),intvec(1,1));
2174  // polynopmial with smallest leading monomial
2175  ring r = 0,x,dp;
2176  min(x+1,x2+x);
2177}
2178
2179
2180///////////////////////////////////////////////////////////////////////////////
2181/*
2182                                Versuche:
2183///////////////////////////////////////////////////////////////////////////////
2184proc downsizeSB (def I, list #)
2185"USAGE:   downsizeSB(I [,l]); I ideal, l list of integers [default: l=0]
2186RETURN:  intvec, say v, with v[j] either 1 or 0. We have v[j]=1 if
2187         leadmonom(I[j]) is divisible by some leadmonom(I[k]) or if
2188         leadmonom(i[j]) == leadmonom(i[k]) and l[j] >= l[k], with k!=j.
2189PURPOSE: The procedure is applied in a situation where the standard basis
2190         computation in the basering R is done via a conversion through an
2191         overring Phelp with additional variables and where a direct
2192         imap from Phelp to R is too expensive.
2193         Assume Phelp is created by the procedure @code{par2varRing} or
2194         @code{hilbRing} and IPhelp is a SB in Phelp [ with l[j]=
2195         length(IPhelp(j)) or any other integer reflecting the complexity
2196         of a IPhelp[j] ]. Let I = lead(IPhelp) mapped to R and compute
2197         v = downsizeSB(imap(Phelp,I),l) in R. Then, if Ihelp[j] is deleted
2198         for all j with v[j]=1, we can apply imap to the remaining generators
2199         of Ihelp and still get SB in R  (in general not minimal).
2200EXAMPLE: example downsizeSB; shows an example"
2201{
2202   int k,j;
2203   intvec v,l;
2204   poly M,N,W;
2205   int c=size(I);
2206   if( size(#) != 0 )
2207   {
2208     if ( typeof(#[1]) == "intvec" )
2209     {
2210       l = #[1];
2211     }
2212     else
2213     {
2214        ERROR("// 2nd argument must be an intvec");
2215     }
2216   }
2217
2218   l[c+1]=0;
2219   v[c]=0;
2220
2221   j=0;
2222   while(j<c-1)
2223   {
2224     j++;
2225     M = leadmonom(I[j]);
2226     if( M != 0 )
2227     {
2228        for( k=j+1; k<=c; k++ )
2229        {
2230          N = leadmonom(I[k]);
2231          if( N != 0 )
2232          {
2233             if( (M==N) && (l[j]>l[k]) )
2234             {
2235               I[j]=0;
2236               v[j]=1;
2237               break;
2238             }
2239             if( (M==N) && (l[j]<=l[k]) || N/M != 0 )
2240             {
2241               I[k]=0;
2242               v[k]=1;
2243             }
2244          }
2245        }
2246      }
2247   }
2248   return(v);
2249}
2250example
2251{ "EXAMPLE:"; echo = 2;
2252   ring  r = 0,(x,y,z,t),(dp(3),dp);
2253   ideal i = x+y+z+t,xy+yz+xt+zt,xyz+xyt+xzt+yzt,xyzt-t4;
2254   ideal Id = std(i);
2255   ideal I = lead(Id);  I;
2256   ring S = (0,t),(x,y,z),dp;
2257   downsizeSB(imap(r,I));
2258   //Id[5] can be deleted, we still have a SB of i in the ring S
2259
2260   ring R = (0,x),(y,z,u),lp;
2261   ideal i = x+y+z+u,xy+xu+yz+zu,xyz+xyu+xzu+yzu,xyzu-1;
2262   def Phelp = par2varRing()[1];
2263   setring Phelp;
2264   ideal IPhelp = std(imap(R,i));
2265   ideal I = lead(IPhelp);
2266   setring R;
2267   ideal I = imap(Phelp,I); I;
2268   intvec v = downsizeSB(I); v;
2269}
2270///////////////////////////////////////////////////////////////////////////
2271// PROBLEM: Die Prozedur funktioniert nur fuer Ringe die global bekannt
2272//          sind, also interaktiv, aber nicht aus einer Prozedur.
2273//          Z.B. funktioniert example imapDownsize; nicht
2274
2275proc imapDownsize (string R, string I)
2276"SYNTAX: @code{imapDownsize (} string @code{,} string @code{)} *@
2277         First string must be the string of the name of a ring, second
2278         string must be the string of the name of an object in the ring.
2279TYPE:    same type as the object with name the second string
2280PURPOSE: maps the object given by the second string to the basering.
2281         If R resp. I are the first resp. second string, then
2282         imapDownsize(R,I) is equivalent to simplify(imap(`R`,`I`),34).
2283NOTE:    imapDownsize is usually faster than imap if `I` is large and if
2284         simplify has a great effect, since the procedure maps only those
2285         generators from `I` which are not killed by simplify( - ,34).
2286         This is useful if `I` is a standard bases for a block ordering of
2287         `R` and if some variables from the last block in `R` are mapped
2288         to parameters. Then the returned result is a standard basis in
2289         the basering.
2290SEE ALSO: imap, fetch, map
2291EXAMPLE: example imapDownsize; shows an example"
2292{
2293       def BR = basering;
2294       int k;
2295
2296       setring `R`;
2297       def @leadI@ = lead(`I`);
2298       int s = ncols(@leadI@);
2299       setring BR;
2300       ideal @leadI@ = simplify(imap(`R`,@leadI@),32);
2301       intvec vi;
2302       for (k=1; k<=s; k++)
2303       {
2304         vi[k] = @leadI@[k]==0;
2305       }
2306       kill @leadI@;
2307
2308       setring `R`;
2309       kill @leadI@;
2310       for (k=1;  k<=s; k++)
2311       {
2312           if( vi[k]==1 )
2313           {
2314              `I`[k]=0;
2315           }
2316       }
2317       `I` = simplify(`I`,2);
2318
2319       setring BR;
2320       return(imap(`R`,`I`));
2321}
2322example
2323{ "EXAMPLE:"; echo = 2;
2324   ring  r = 0,(x,y,z,t),(dp(3),dp);
2325   ideal i = x+y+z+t,xy+yz+xt+zt,xyz+xyt+xzt+yzt,xyzt-1;
2326   i = std(i); i;
2327
2328   ring s = (0,t),(x,y,z),dp;
2329   imapDownsize("r","i");     //i[5] is omitted since lead(i[2]) | lead(i[5])
2330}
2331///////////////////////////////////////////////////////////////////////////////
2332//die folgende proc war fuer groebner mit fglm vorgesehen, ist aber zu teuer.
2333//Um die projektive Dimension korrekt zu berechnen, muss man aber teuer
2334//voerher ein SB bzgl. einer Gradordnung berechnen und dann homogenisieren.
2335//Sonst koennen hoeherdimensionale Komponenten in Unendlich entstehen
2336
2337proc projInvariants(ideal i,list #)
2338"SYNTAX: @code{projInvariants (} ideal_expression @code{)} @*
2339         @code{projInvariants (} ideal_expression@code{,} list of string_expres          sions@code{)}
2340TYPE:    list, say L, with L[1] and L[2] of type int and L[3] of type intvec
2341PURPOSE: Computes the (projective) dimension (L[1]), degree (L[2]) and the
2342         first Hilbert series (L[3], as intvec) of the homogenized ideal
2343         in the ring given by the procedure @code{hilbRing} with global
2344         ordering dp (resp. wp if the variables have weights >1)
2345         If an argument of type string @code{\"std\"} resp. @code{\"slimgb\"}
2346         is given, the standard basis computatuion uses @code{std} or
2347         @code{slimgb}, otherwise a heuristically chosen method (default)
2348NOTE:    Homogenized means weighted homogenized with respect to the weights
2349         w[i] of the variables var(i) of the basering. The returned dimension,
2350         degree and Hilbertseries are the respective invariants of the
2351         projective variety defined by the homogenized ideal. The dimension
2352         is equal to the (affine) dimension of the ideal in the basering
2353         (degree and Hilbert series make only sense for homogeneous ideals).
2354SEE ALSO: dim, dmult, hilb
2355KEYWORDS: dimension, degree, Hilbert function
2356EXAMPLE: example projInvariants;  shows an example"
2357{
2358  def P = basering;
2359  int p_opt;
2360  string s_opt = option();
2361  if (find(option(), "prot"))  { p_opt = 1; }
2362
2363//---------------- check method and clear denomintors --------------------
2364  int k;
2365  string method;
2366  for (k=1; k<=size(#); k++)
2367  {
2368     if (typeof(#[k]) == "string")
2369     {
2370       method = method + "," + #[k];
2371     }
2372  }
2373
2374  if (npars(P) > 0)             //clear denominators of parameters
2375  {
2376    for( k=ncols(i); k>0; k-- )
2377    {
2378       i[k]=cleardenom(i[k]);
2379    }
2380  }
2381
2382//------------------------ change to hilbRing ----------------------------
2383     list hiRi = hilbRing(i);
2384     intvec W = hiRi[2];
2385     def Philb = hiRi[1];      //note: Philb is no qring and the predefined
2386     setring Philb;            //ideal Id(1) in Philb is homogeneous
2387     int di, de;               //for dimension, degree
2388     intvec hi;                //for hilbert series
2389
2390//-------- compute Hilbert function of homogenized ideal in Philb ---------
2391//Philb has only 1 block. There are three cases
2392
2393     string algorithm;       //possibilities: std, slimgb, stdorslimgb
2394     //define algorithm:
2395     if( find(method,"std") && !find(method,"slimgb") )
2396     {
2397        algorithm = "std";
2398     }
2399     if( find(method,"slimgb") && !find(method,"std") )
2400     {
2401         algorithm = "slimgb";
2402     }
2403     if( find(method,"std") && find(method,"slimgb") ||
2404         (!find(method,"std") && !find(method,"slimgb")) )
2405     {
2406         algorithm = "stdorslimgb";
2407     }
2408
2409     if ( algorithm=="std" || ( algorithm=="stdorslimgb" && char(P)>0 ) )
2410     {
2411        if (p_opt) {"std in ring " + string(Philb);}
2412        Id(1) = std(Id(1));
2413        di = dim(Id(1))-1;
2414        de = mult(Id(1));
2415        hi = hilb( Id(1),1,W );
2416     }
2417     if ( algorithm=="slimgb" || ( algorithm=="stdorslimgb" && char(P)==0 ) )
2418     {
2419        if (p_opt) {"slimgb in ring " + string(Philb);}
2420        Id(1) = slimgb(Id(1));
2421        di = dim( Id(1) );
2422        if (di > -1)
2423        {
2424           di = di-1;
2425        }
2426        de = mult( Id(1) );
2427        hi = hilb( Id(1),1,W );
2428     }
2429     kill Philb;
2430     list L = di,de,hi;
2431     return(L);
2432}
2433example
2434{ "EXAMPLE:"; echo = 2;
2435   ring r = 32003,(x,y,z),lp;
2436   ideal i = y2-xz,x2-z;
2437   projInvariants(i);
2438
2439   ring R = (0),(x,y,z,u,v),lp;
2440   //minpoly = x2+1;
2441   ideal i = x2+1,x2+y+z+u+v,xyzuv-1;
2442   projInvariants(i);
2443   qring S =std(x2+1);
2444   ideal i = imap(R,i);
2445   projInvariants(i);
2446}
2447
2448*/
2449///////////////////////////////////////////////////////////////////////////////
2450//                           EXAMPLES
2451///////////////////////////////////////////////////////////////////////////////
2452/*
2453example stdfglm;
2454example stdhilb;
2455example groebner;
2456example res;
2457example sprintf;
2458example fprintf;
2459example printf;
2460example weightKB;
2461example qslimgb;
2462example par2varRing;
2463*/
2464static proc mod_init()
2465{
2466  if(!defined(Singmathic))
2467  {
2468    load("singmathic.so","try");
2469  }
2470  //int pagelength=24;
2471  //exportto(Top,pagelength);
2472}
Note: See TracBrowser for help on using the repository browser.