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

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