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

spielwiese
Last change on this file since f84afa was f84afa, checked in by Michael Brickenstein <bricken@…>, 19 years ago
*bricken: qring for slimgb forbidden git-svn-id: file:///usr/local/Singular/svn/trunk@8449 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 27.4 KB
Line 
1//////////////////////////////////////////////////////////////////////////////
2version="$Id: standard.lib,v 1.72 2005-07-25 14:46:50 bricken Exp $";
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])     standard basis of ideal using the Hilbert function
10 groebner(ideal/module) standard basis using a heuristically chosen method
11 quot(any,any[,n])      quotient using heuristically chosen method
12 res(ideal/module,[i])  free resolution of ideal or module
13 sprintf(fmt,...)       returns fomatted string
14 fprintf(link,fmt,..)   writes formatted string to link
15 printf(fmt,...)        displays formatted string
16";
17
18//////////////////////////////////////////////////////////////////////////////
19
20proc stdfglm (ideal i, list #)
21"SYNTAX: @code{stdfglm (} ideal_expression @code{)} @*
22         @code{stdfglm (} ideal_expression@code{,} string_expression @code{)}
23TYPE:    ideal
24PURPOSE: computes the standard basis of the ideal in the basering
25         via @code{fglm} (from the ordering given as the second argument
26         to the ordering of the basering).@*
27         If no second argument is given, \"dp\" is used.
28SEE ALSO: fglm, groebner, std, stdhilb
29KEYWORDS: fglm
30EXAMPLE: example stdfglm; shows an example"
31{
32   string os;
33   def dr= basering;
34   if( (size(#)==0) or (typeof(#[1]) != "string") )
35   {
36     os = "dp(" + string( nvars(dr) ) + ")";
37     if ( (find( ordstr(dr), os ) != 0) and (find( ordstr(dr), "a") == 0) )
38     {
39       os= "Dp";
40     }
41     else
42     {
43       os= "dp";
44     }
45   }
46   else { os = #[1]; }
47   execute("ring sr=("+charstr(dr)+"),("+varstr(dr)+"),"+os+";");
48   ideal i= fetch(dr,i);
49   intvec opt= option(get);
50   option(redSB);
51   i=std(i);
52   option(set,opt);
53   setring dr;
54   return (fglm(sr,i));
55}
56example
57{ "EXAMPLE:"; echo = 2;
58   ring r=0,(x,y,z),lp;
59   ideal i=y3+x2,x2y+x2,x3-x2,z4-x2-y;
60   ideal i1=stdfglm(i);         //uses fglm from "dp" to "lp"
61   i1;
62   ideal i2=stdfglm(i,"Dp");    //uses fglm from "Dp" to "lp"
63   i2;
64}
65/////////////////////////////////////////////////////////////////////////////
66
67proc stdhilb(ideal i,list #)
68"SYNTAX: @code{stdhilb (} ideal_expression @code{)} @*
69         @code{stdhilb (} ideal_expression@code{,} intvec_expression @code{)}
70TYPE:    ideal
71PURPOSE: computes the standard basis of the homogeneous ideal in the basering,
72         via a Hilbert driven standard basis computation.@*
73         An optional second argument will be used as 1st Hilbert function.
74ASSUME:  The optional second argument is the first Hilbert series as computed
75         by @code{hilb}.
76SEE ALSO: stdfglm, std, groebner
77KEYWORDS: Hilbert function
78EXAMPLE: example stdhilb;  shows an example"
79{
80   def R=basering;
81
82   if((homog(i)==1)||(ordstr(basering)[1]=="d"))
83   {
84      if ((size(#)!=0)&&(homog(i)==1))
85      {
86         return(std(i,#[1]));
87      }
88      return(std(i));
89   }
90
91   execute("ring S = ("+charstr(R)+"),("+varstr(R)+",@t),dp;");
92   ideal i=homog(imap(R,i),@t);
93   intvec v=hilb(std(i),1);
94   execute("ring T = ("+charstr(R)+"),("+varstr(R)+",@t),("+ordstr(R)+");");
95   ideal i=fetch(S,i);
96   ideal a=std(i,v);
97   setring R;
98   map phi=T,maxideal(1),1;
99   ideal a=phi(a);
100
101   int k,j;
102   poly m;
103   int c=ncols(i);
104
105   for(j=1;j<c;j++)
106   {
107     if(deg(a[j])==0)
108     {
109       a=ideal(1);
110       attrib(a,"isSB",1);
111       return(a);
112     }
113     if(deg(a[j])>0)
114     {
115       m=lead(a[j]);
116       for(k=j+1;k<=c;k++)
117       {
118          if(size(lead(a[k])/m)>0)
119          {
120            a[k]=0;
121          }
122       }
123     }
124   }
125   a=simplify(a,2);
126   attrib(a,"isSB",1);
127   return(a);
128}
129example
130{ "EXAMPLE:"; echo = 2;
131   ring  r=0,(x,y,z),dp;
132   ideal i=y3+x2,x2y+x2,x3-x2,z4-x2-y;
133   ideal i1=stdhilb(i); i1;
134   // the latter computation is equivalent to:
135   intvec v=hilb(i,1);
136   ideal i2=stdhilb(i,v); i2;
137}
138//////////////////////////////////////////////////////////////////////////
139
140proc groebner(def i, list #)
141"SYNTAX: @code{groebner (} ideal_expression @code{)} @*
142         @code{groebner (} module_expression @code{)} @*
143         @code{groebner (} ideal_expression@code{,} int_expression @code{)} @*
144         @code{groebner (} module_expression@code{,} int_expression @code{)}
145TYPE:    type of the first argument
146PURPOSE: computes the standard basis of the first argument @code{I}
147         (ideal or module), by a heuristically chosen method: if the
148         ordering of the current ring is a local ordering, or if it is a
149         non-block ordering and the current ring has no parameters, then
150         @code{std(I)} is returned. Otherwise, @code{I} is mapped into a
151         ring with no parameters and ordering dp, where its Hilbert series
152         is computed. This is followed by a Hilbert-series based std
153         computation in the original ring.
154NOTE: If a 2nd argument @code{wait} is given, then the computation proceeds
155      at most @code{wait} seconds. That is, if no result could be computed in
156      @code{wait} seconds, then the computation is interrupted, 0 is returned,
157      a warning message is displayed, and the global variable
158      @code{groebner_error} is defined.
159SEE ALSO: stdhilb, stdfglm, std
160KEYWORDS: time limit on computations; MP, groebner basis computations
161EXAMPLE: example groebner;  shows an example"
162{
163  def P=basering;
164
165  // we have two arguments -- try to use MPfork links
166  if (size(#) > 0)
167  {
168    if (system("with", "MP"))
169    {
170      if (typeof(#[1]) == "int")
171      {
172        int wait = #[1];
173        int j = 10;
174
175        string bs = nameof(basering);
176        link l_fork = "MPtcp:fork";
177        open(l_fork);
178        write(l_fork, quote(system("pid")));
179        int pid = read(l_fork);
180        write(l_fork, quote(groebner(eval(i))));
181
182        // sleep in small intervalls for appr. one second
183        if (wait > 0)
184        {
185          while(j < 1000000)
186          {
187            if (status(l_fork, "read", "ready", j)) {break;}
188            j = j + j;
189          }
190        }
191
192        // sleep in intervalls of one second from now on
193        j = 1;
194        while (j < wait)
195        {
196          if (status(l_fork, "read", "ready", 1000000)) {break;}
197          j = j + 1;
198        }
199
200        if (status(l_fork, "read", "ready"))
201        {
202          def result = read(l_fork);
203          if (bs != nameof(basering))
204          {
205            def PP = basering;
206            setring P;
207            def result = imap(PP, result);
208            kill PP;
209          }
210          if (defined(groebner_error))
211          {
212            kill groebner_error;
213          }
214          kill l_fork;
215        }
216        else
217        {
218          ideal result;
219          if (! defined(groebner_error))
220          {
221            int groebner_error = 1;
222            export groebner_error;
223          }
224          "// ** groebner did not finish";
225          j = system("sh", "kill " + string(pid));
226        }
227        return (result);
228      }
229      else
230      {
231        "// ** groebner needs int as 2nd arg";
232      }
233    }
234    else
235    {
236      "// ** groebner with two args is not supported in this configuration";
237    }
238  }
239
240  // we are still here -- do the actual computation
241  string ordstr_P = ordstr(P);
242  if ((find(ordstr_P,"s") > 0)
243  ||(find(ordstr_P,"M") > 0)
244  ||(find(ordstr_P,"w") > 0)
245  ||(find(ordstr_P,"W") > 0))
246  {
247    //spaeter den lokalen fall ueber lp oder aehnlich behandeln
248    return(std(i));
249  }
250  //if ordering is global, there are parameters and minpoly is 0
251  if (((npars(basering)>0) &&(minpoly==0)) &&(typeof(basering)=="ring"))
252  {
253    return(slimgb(i));
254  }
255  int IsSimple_P;
256  if (system("nblocks") <= 2)
257  {
258    if (find(ordstr_P, "M") <= 0)
259    {
260      IsSimple_P = 1;
261    }
262  }
263  int npars_P = npars(P);
264
265  // return std if no parameters and (dp or wp)
266  if ((npars_P <= 1) && IsSimple_P)
267  {
268    if (find(ordstr_P, "d") > 0)
269    {
270      return (std(i));
271    }
272    if (find(ordstr_P,"w") > 0)
273    {
274      return (std(i));
275    }
276  }
277
278  // reset options
279  intvec opt=option(get);
280  int p_opt;
281  string s_opt = option();
282  option(none);
283  // turn on option(prot) and/or option(mem), if previously set
284  if (find(s_opt, "prot"))
285  {
286    option(prot);
287    p_opt = 1;
288  }
289  if (find(s_opt, "mem"))
290  { option(mem); }
291  if (find(s_opt, "intStrategy"))
292  { option(intStrategy); }
293
294  // construct ring in which first std computation is done
295  string varstr_P = varstr(P);
296  string parstr_P = parstr(P);
297  int is_homog = (homog(i) && (npars_P <= 1));
298  int add_vars = 0;
299  string ri = "ring Phelp =";
300
301  if (npars_P > 0)
302  {
303    for(int k=ncols(i); k>0; k--) { i[k]=cleardenom(i[k]); }
304  }
305  // more than one parameters are converted to ring variables
306  if (npars_P > 1)
307  {
308    ri = ri + string(char(P)) + ",(" + varstr_P + "," + parstr_P;
309    add_vars = npars_P;
310  }
311  else
312  {
313    ri = ri + "(" + charstr(P) + "),(" + varstr_P;
314  }
315
316  // a homogenizing variable is added, if necessary
317  if (! is_homog)
318  {
319    ri = ri + ",@t";
320    add_vars = add_vars + 1;
321  }
322  // ordering is set to (dp, C)
323  ri = ri + "),(dp,C);";
324
325  // change the ring
326  execute(ri);
327
328  // get ideal from previous ring
329  if (is_homog)
330  {
331    ideal qh = imap(P, i);
332  }
333  else
334  {
335    // and homogenize
336    ideal qh=homog(imap(P,i),@t);
337  }
338
339  // compute std and hilbert series
340  if (p_opt)
341  {
342    "std in " + ri[13, size(ri) - 13];
343  }
344  intvec hi=hilb(std(qh),1);
345
346  if (add_vars == 0)
347  {
348    // no additional variables were introduced
349    setring P; // can immediately change to original ring
350    // simply compute std with hilbert series in original ring
351    if (p_opt)
352    {
353      "std with hilb in basering";
354    }
355    i = std(i, hi);
356  }
357  else
358  {
359    // additional variables were introduced
360    // need another intermediate ring
361    ri = "ring Phelp1 = (" + charstr(Phelp)
362      + "),(" + varstr(Phelp) + "),(" + ordstr_P;
363
364    // for lp wit at most one parameter, we do not need a block ordering
365    if ( ! (IsSimple_P && (add_vars <2) && find(ordstr_P, "l")))
366    {
367      // need block ordering
368      ri = ri + ", dp(" + string(add_vars) + ")";
369    }
370    ri = ri + ");";
371
372    // change to intermediate ring
373    execute(ri);
374    ideal qh = imap(Phelp, qh);
375    kill Phelp;
376    if (p_opt)
377    {
378      "std with hilb in " + ri[14,size(ri)-14];
379    }
380    // compute std with Hilbert series
381    qh = std(qh, hi);
382    // subst 1 for homogenizing var
383    if (!is_homog)
384    {
385      if (p_opt)
386      {
387        "dehomogenization";
388      }
389      qh = subst(qh, @t, 1);
390    }
391
392    // go back to original ring
393    setring P;
394    // get ideal, delete zeros and clean SB
395    if (p_opt)
396    {
397      "imap to original ring";
398    }
399    i = imap(Phelp1,qh);
400    if (p_opt)
401    {
402      "simplification";
403    }
404    i = simplify(i, 34);
405    kill Phelp1;
406  }
407
408  // clean-up time
409  option(set, opt);
410  if (find(s_opt, "redSB") > 0)
411  {
412    if (p_opt)
413    {
414      "interreduction";
415    }
416    i=interred(i);
417  }
418  attrib(i, "isSB", 1);
419  return (i);
420}
421example
422{ "EXAMPLE: "; echo=2;
423  ring r=0,(a,b,c,d),lp;
424  option(prot);
425  ideal i=a+b+c+d,ab+ad+bc+cd,abc+abd+acd+bcd,abcd-1; // cyclic 4
426  groebner(i);
427  ring rp=(0,a,b),(c,d), lp;
428  ideal i=imap(r,i);
429  ideal j=groebner(i);
430  option(noprot);
431  j; simplify(j,1); std(i);
432  if (system("with","MP")) {groebner(i,0);}
433  defined(groebner_error);
434}
435//////////////////////////////////////////////////////////////////////////
436
437proc res(list #)
438"@c we do texinfo here:
439@cindex resolution, computation of
440@table @code
441@item @strong{Syntax:}
442@code{res (} ideal_expression@code{,} int_expression @code{[,} any_expression @code{])}
443@*@code{res (} module_expression@code{,} int_expression @code{[,} any_expression @code{])}
444@item @strong{Type:}
445resolution
446@item @strong{Purpose:}
447computes a (possibly minimal) free resolution of an ideal or module using
448a heuristically chosen method.
449@* The second (int) argument (say, @code{k}) specifies the length of
450the resolution. If it is not positive then @code{k} is assumed to be the
451number of variables of the basering.
452@* If a third argument is given, the returned resolution is minimized.
453
454Depending on the input, the returned resolution is computed using the
455following methods:
456@table @asis
457@item @strong{quotient rings:}
458@code{nres} (classical method using syzygies) , see @ref{nres}.
459
460@item @strong{homogeneous ideals and k=0:}
461@code{lres} (La'Scala's method), see @ref{lres}.
462
463@item @strong{not minimized resolution and (homogeneous input with k not 0, or local rings):}
464@code{sres} (Schreyer's method), see @ref{sres}.
465
466@item @strong{all other inputs:}
467@code{mres} (classical method), see @ref{mres}.
468@end table
469@item @strong{Note:}
470Accessing single elements of a resolution may require that some partial
471computations have to be finished and may therefore take some time.
472@end table
473@c ref
474See also
475@ref{betti};
476@ref{ideal};
477@ref{minres};
478@ref{module};
479@ref{mres};
480@ref{nres};
481@ref{lres};
482@ref{hres};
483@ref{sres}.
484@ref{resolution}
485@c ref
486"
487{
488   def P=basering;
489   if (size(#) < 2)
490   {
491     ERROR("res: need at least two arguments: ideal/module, int");
492   }
493
494   def m=#[1]; //the ideal or module
495   int i=#[2]; //the length of the resolution
496   if (i< 0) { i=0;}
497
498   string varstr_P = varstr(P);
499
500   int p_opt;
501   string s_opt = option();
502   // set p_opt, if option(prot) is set
503   if (find(s_opt, "prot"))
504   {
505     p_opt = 1;
506   }
507
508   if(size(ideal(basering)) > 0)
509   {
510     // the quick hack for qrings - seems to fit most needs
511     // (lres is not implemented for qrings, sres is not so efficient)
512     if (p_opt) { "using nres";}
513     return(nres(m,i));
514   }
515
516   if(homog(m)==1)
517   {
518      resolution re;
519      if (((i==0) or (i>=nvars(basering))) && typeof(m) != "module")
520      {
521        //LaScala for the homogeneous case and i == 0
522        if (p_opt) { "using lres";}
523        re=lres(m,i);
524        if(size(#)>2)
525        {
526           re=minres(re);
527        }
528      }
529      else
530      {
531        if(size(#)>2)
532        {
533          if (p_opt) { "using mres";}
534          re=mres(m,i);
535        }
536        else
537        {
538          if (p_opt) { "using sres";}
539          re=sres(std(m),i);
540        }
541      }
542      return(re);
543   }
544
545   //mres for the global non homogeneous case
546   if(find(ordstr(P),"s")==0)
547   {
548      string ri= "ring Phelp ="
549                  +string(char(P))+",("+varstr_P+"),(dp,C);";
550      execute(ri);
551      def m=imap(P,m);
552      if (p_opt) { "using mres in another ring";}
553      list re=mres(m,i);
554      setring P;
555      resolution result=imap(Phelp,re);
556      if (size(#) > 2) {result = minres(result);}
557      return(result);
558   }
559
560   //sres for the local case and not minimal resolution
561   if(size(#)<=2)
562   {
563      string ri= "ring Phelp ="
564                  +string(char(P))+",("+varstr_P+"),(ls,c);";
565      execute(ri);
566      def m=imap(P,m);
567      m=std(m);
568      if (p_opt) { "using sres in another ring";}
569      list re=sres(m,i);
570      setring P;
571      resolution result=imap(Phelp,re);
572      return(result);
573   }
574
575   //mres for the local case and minimal resolution
576   string ri= "ring Phelp ="
577                  +string(char(P))+",("+varstr_P+"),(ls,C);";
578   execute(ri);
579   def m=imap(P,m);
580    if (p_opt) { "using mres in another ring";}
581   list re=mres(m,i);
582   setring P;
583   resolution result=imap(Phelp,re);
584   result = minres(result);
585   return(result);
586}
587example
588{"EXAMPLE:"; echo = 2;
589  ring r=0,(x,y,z),dp;
590  ideal i=xz,yz,x3-y3;
591  def l=res(i,0); // homogeneous ideal: uses lres
592  l;
593  print(betti(l), "betti"); // input to betti may be of type resolution
594  l[2];         // element access may take some time
595  i=i,x+1;
596  l=res(i,0);   // inhomogeneous ideal: uses mres
597  l;
598  ring rs=0,(x,y,z),ds;
599  ideal i=imap(r,i);
600  def l=res(i,0); // local ring not minimized: uses sres
601  l;
602  res(i,0,0);     // local ring and minimized: uses mres
603}
604/////////////////////////////////////////////////////////////////////////
605
606proc quot (m1,m2,list #)
607"SYNTAX: @code{quot (} module_expression@code{,} module_expression @code{)}
608         @*@code{quot (} module_expression@code{,} module_expression@code{,}
609            int_expression @code{)}
610         @*@code{quot (} ideal_expression@code{,} ideal_expression @code{)}
611         @*@code{quot (} ideal_expression@code{,} ideal_expression@code{,}
612            int_expression @code{)}
613TYPE:    ideal
614SYNTAX:  @code{quot (} module_expression@code{,} ideal_expression @code{)}
615TYPE:    module
616PURPOSE: computes the quotient of the 1st and the 2nd argument.
617         If a 3rd argument 'n' is given the n-th method is used
618         (n=1...5).
619SEE ALSO: quotient
620EXAMPLE: example quot; shows an example"
621{
622  if (((typeof(m1)!="ideal") and (typeof(m1)!="module"))
623     or ((typeof(m2)!="ideal") and (typeof(m2)!="module")))
624  {
625    "USAGE:   quot(m1, m2[, n]); m1, m2 two submodules of k^s,";
626    "         n (optional) integer (1<= n <=5)";
627    "RETURN:  the quotient of m1 and m2";
628    "EXAMPLE: example quot; shows an example";
629    return();
630  }
631  if (typeof(m1)!=typeof(m2))
632  {
633    return(quotient(m1,m2));
634  }
635  if (size(#)>0)
636  {
637    if (typeof(#[1])=="int" )
638    {
639      return(quot1(m1,m2,#[1]));
640    }
641  }
642  else
643  {
644    return(quot1(m1,m2,2));
645  }
646}
647example
648{ "EXAMPLE:"; echo = 2;
649  ring r=181,(x,y,z),(c,ls);
650  ideal id1=maxideal(4);
651  ideal id2=x2+xyz,y2-z3y,z3+y5xz;
652  option(prot);
653  ideal id3=quotient(id1,id2);
654  id3;
655  ideal id4=quot(id1,id2,1);
656  id4;
657  ideal id5=quot(id1,id2,2);
658  id5;
659}
660
661static proc quot1 (module m1, module m2,int n)
662"USAGE:   quot1(m1, m2, n); m1, m2 two submodules of k^s,
663         n integer (1<= n <=5)
664RETURN:  the quotient of m1 and m2
665EXAMPLE: example quot1; shows an example"
666{
667  if (n==1)
668  {
669    return(quotient1(m1,m2));
670  }
671  else
672  {
673    if (n==2)
674    {
675      return(quotient2(m1,m2));
676    }
677    else
678    {
679      if (n==3)
680      {
681        return(quotient3(m1,m2));
682      }
683      else
684      {
685        if (n==4)
686        {
687          return(quotient4(m1,m2));
688        }
689        else
690        {
691          if (n==5)
692          {
693            return(quotient5(m1,m2));
694          }
695          else
696          {
697            return(quotient(m1,m2));
698          }
699        }
700      }
701    }
702  }
703}
704example
705{ "EXAMPLE:"; echo = 2;
706  ring r=181,(x,y,z),(c,ls);
707  ideal id1=maxideal(4);
708  ideal id2=x2+xyz,y2-z3y,z3+y5xz;
709  option(prot);
710  ideal id6=quotient(id1,id2);
711  id6;
712  ideal id7=quot1(id1,id2,1);
713  id7;
714  ideal id8=quot1(id1,id2,2);
715  id8;
716}
717
718static proc quotient0(module a,module b)
719{
720  module mm=b+a;
721  resolution rs=lres(mm,0);
722  list I=list(rs);
723  matrix M=I[2];
724  matrix A[1][nrows(M)]=M[1..nrows(M),1];
725  ideal i=A;
726  return (i);
727}
728proc quotient1(module a,module b)  //17sec
729"USAGE:   quotient1(m1, m2); m1, m2 two submodules of k^s,
730RETURN:  the quotient of m1 and m2"
731{
732  int i;
733  a=std(a);
734  module dummy;
735  module B=NF(b,a)+dummy;
736  ideal re=quotient(a,module(B[1]));
737  for(i=2;i<=ncols(B);i++)
738  {
739     re=intersect1(re,quotient(a,module(B[i])));
740  }
741  return(re);
742}
743proc quotient2(module a,module b)    //13sec
744"USAGE:   quotient2(m1, m2); m1, m2 two submodules of k^s,
745RETURN:  the quotient of m1 and m2"
746{
747  a=std(a);
748  module dummy;
749  module bb=NF(b,a)+dummy;
750  int i=ncols(bb);
751  ideal re=quotient(a,module(bb[i]));
752  bb[i]=0;
753  module temp;
754  module temp1;
755  module bbb;
756  int mx;
757  i=i-1;
758  while (1)
759  {
760    if (i==0) break;
761    temp = a+bb*re;
762    temp1 = lead(interred(temp));
763    mx=ncols(a);
764    if (ncols(temp1)>ncols(a))
765    {
766      mx=ncols(temp1);
767    }
768    temp1 = matrix(temp1,1,mx)-matrix(lead(a),1,mx);
769    temp1 = dummy+temp1;
770    if (deg(temp1[1])<0) break;
771    re=intersect1(re,quotient(a,module(bb[i])));
772    bb[i]=0;
773    i = i-1;
774  }
775  return(re);
776}
777proc quotient3(module a,module b)   //89sec
778"USAGE:   quotient3(m1, m2); m1, m2 two submodules of k^s,
779         only for global rings
780RETURN:  the quotient of m1 and m2"
781{
782  string s="ring @newr=("+charstr(basering)+
783           "),("+varstr(basering)+",@t,@w),dp;";
784  def @newP=basering;
785  execute(s);
786  module b=imap(@newP,b);
787  module a=imap(@newP,a);
788  int i;
789  int j=ncols(b);
790  vector @b;
791  for(i=1;i<=j;i++)
792  {
793     @b=@b+@t^(i-1)*@w^(j-i+1)*b[i];
794  }
795  ideal re=quotient(a,module(@b));
796  setring @newP;
797  ideal re=imap(@newr,re);
798  return(re);
799}
800proc quotient5(module a,module b)   //89sec
801"USAGE:   quotient5(m1, m2); m1, m2 two submodules of k^s,
802         only for global rings
803RETURN:  the quotient of m1 and m2"
804{
805  string s="ring @newr=("+charstr(basering)+
806           "),("+varstr(basering)+",@t),dp;";
807  def @newP=basering;
808  execute(s);
809  module b=imap(@newP,b);
810  module a=imap(@newP,a);
811  int i;
812  int j=ncols(b);
813  vector @b;
814  for(i=1;i<=j;i++)
815  {
816     @b=@b+@t^(i-1)*b[i];
817  }
818  @b=homog(@b,@w);
819  ideal re=quotient(a,module(@b));
820  setring @newP;
821  ideal re=imap(@newr,re);
822  return(re);
823}
824proc quotient4(module a,module b)   //95sec
825"USAGE:   quotient4(m1, m2); m1, m2 two submodules of k^s,
826         only for global rings
827RETURN:  the quotient of m1 and m2"
828{
829  string s="ring @newr=("+charstr(basering)+
830           "),("+varstr(basering)+",@t),dp;";
831  def @newP=basering;
832  execute(s);
833  module b=imap(@newP,b);
834  module a=imap(@newP,a);
835  int i;
836  vector @b=b[1];
837  for(i=2;i<=ncols(b);i++)
838  {
839     @b=@b+@t^(i-1)*b[i];
840  }
841  matrix sy=modulo(@b,a);
842  ideal re=sy;
843  setring @newP;
844  ideal re=imap(@newr,re);
845  return(re);
846}
847static proc intersect1(ideal i,ideal j)
848{
849  def R=basering;
850  execute("ring gnir = ("+charstr(basering)+"),
851                       ("+varstr(basering)+",@t),(C,dp);");
852  ideal i=var(nvars(basering))*imap(R,i)+(var(nvars(basering))-1)*imap(R,j);
853  ideal j=eliminate(i,var(nvars(basering)));
854  setring R;
855  map phi=gnir,maxideal(1);
856  return(phi(j));
857}
858
859//////////////////////////////////////////////////////////////////
860///
861/// sprintf, fprintf printf
862///
863proc sprintf(string fmt, list #)
864"SYNTAX:  @code{sprintf (} string_expression @code{[,} any_expressions
865               @code{] )}
866RETURN:   string
867PURPOSE:  @code{sprintf(fmt,...);} performs output formatting. The first
868          argument is a format control string. Additional arguments may be
869          required, depending on the content of the control string. A series
870          of output characters is generated as directed by the control string;
871          these characters are returned as a string. @*
872          The control string @code{fmt} is simply text to be copied,
873          except that the string may contain conversion specifications.@*
874          Do @code{help print;} for a listing of valid conversion
875          specifications. As an addition to the conversions of @code{print},
876          the @code{%n} and @code{%2} conversion specification does not
877          consume an additional argument, but simply generates a newline
878          character.
879NOTE:     If one of the additional arguments is a list, then it should be
880          enclosed once more into a @code{list()} command, since passing a list
881          as an argument flattens the list by one level.
882SEE ALSO: fprintf, printf, print, string
883EXAMPLE : example sprintf; shows an example
884"
885{
886  int sfmt = size(fmt);
887  if (sfmt  <= 1)
888  {
889    return (fmt);
890  }
891  int next, l, nnext;
892  string ret;
893  list formats = "%l", "%s", "%2l", "%2s", "%t", "%;", "%p", "%b", "%n", "%2";
894  while (1)
895  {
896    if (size(#) <= 0)
897    {
898      return (ret + fmt);
899    }
900    nnext = 0;
901    while (nnext < sfmt)
902    {
903      nnext = find(fmt, "%", nnext + 1);
904      if (nnext == 0)
905      {
906        next = 0;
907        break;
908      }
909      l = 1;
910      while (l <= size(formats))
911      {
912        next = find(fmt, formats[l], nnext);
913        if (next == nnext) break;
914        l++;
915      }
916      if (next == nnext) break;
917    }
918    if (next == 0)
919    {
920      return (ret + fmt);
921    }
922    if (formats[l] != "%2" && formats[l] != "%n")
923    {
924      ret = ret + fmt[1, next - 1] + print(#[1], formats[l]);
925      # = delete(#, 1);
926    }
927    else
928    {
929      ret = ret + fmt[1, next - 1] + print("", "%2s");
930    }
931    if (size(fmt) <= (next + size(formats[l]) - 1))
932    {
933      return (ret);
934    }
935    fmt = fmt[next + size(formats[l]), size(fmt)-next-size(formats[l]) + 1];
936  }
937}
938example
939{ "EXAMPLE:"; echo=2;
940  ring r=0,(x,y,z),dp;
941  module m=[1,y],[0,x+z];
942  intmat M=betti(mres(m,0));
943  list l = r, m, M;
944  string s = sprintf("s:%s,%n l:%l", 1, 2); s;
945  s = sprintf("s:%n%s", l); s;
946  s = sprintf("s:%2%s", list(l)); s;
947  s = sprintf("2l:%n%2l", list(l)); s;
948  s = sprintf("%p", list(l)); s;
949  s = sprintf("%;", list(l)); s;
950  s = sprintf("%b", M); s;
951}
952
953proc printf(string fmt, list #)
954"SYNTAX:  @code{printf (} string_expression @code{[,} any_expressions@code{] )}
955RETURN:   none
956PURPOSE:  @code{printf(fmt,...);} performs output formatting. The first
957          argument is a format control string. Additional arguments may be
958          required, depending on the content of the control string. A series
959          of output characters is generated as directed by the control string;
960          these characters are displayed (i.e., printed to standard out). @*
961          The control string @code{fmt} is simply text to be copied, except
962          that the string may contain conversion specifications. @*
963          Do @code{help print;} for a listing of valid conversion
964          specifications. As an addition to the conversions of @code{print},
965          the @code{%n} and @code{%2} conversion specification does not
966          consume an additional argument, but simply generates a newline
967          character.
968NOTE:     If one of the additional arguments is a list, then it should be
969          enclosed once more into a @code{list()} command, since passing a
970          list as an argument flattens the list by one level.
971SEE ALSO: sprintf, fprintf, print, string
972EXAMPLE : example printf; shows an example
973"
974{
975  write("", sprintf(fmt, #));
976}
977example
978{ "EXAMPLE:"; echo=2;
979  ring r=0,(x,y,z),dp;
980  module m=[1,y],[0,x+z];
981  intmat M=betti(mres(m,0));
982  list l=r,m,M;
983  printf("s:%s,l:%l",1,2);
984  printf("s:%s",l);
985  printf("s:%s",list(l));
986  printf("2l:%2l",list(l));
987  printf("%p",list(l));
988  printf("%;",list(l));
989  printf("%b",M);
990}
991
992
993proc fprintf(link l, string fmt, list #)
994"SYNTAX:  @code{fprintf (} link_expression@code{,} string_expression @code{[,}
995                any_expressions@code{] )}
996RETURN:   none
997PURPOSE:  @code{fprintf(l,fmt,...);} performs output formatting.
998          The second argument is a format control string. Additional
999          arguments may be required, depending on the content of the
1000          control string. A series of output characters is generated as
1001          directed by the control string; these characters are
1002          written to the link l.
1003          The control string @code{fmt} is simply text to be copied, except
1004          that the string may contain conversion specifications.@*
1005          Do @code{help print;} for a listing of valid conversion
1006          specifications. As an addition to the conversions of @code{print},
1007          the @code{%n} and @code{%2} conversion specification does not
1008          consume an additional argument, but simply generates a newline
1009          character.
1010NOTE:     If one of the additional arguments is a list, then it should be
1011          enclosed once more into a @code{list()} command, since passing
1012          a list as an argument flattens the list by one level.
1013SEE ALSO: sprintf, printf, print, string
1014EXAMPLE : example fprintf; shows an example
1015"
1016{
1017  write(l, sprintf(fmt, #));
1018}
1019example
1020{ "EXAMPLE:"; echo=2;
1021  ring r=0,(x,y,z),dp;
1022  module m=[1,y],[0,x+z];
1023  intmat M=betti(mres(m,0));
1024  list l=r,m,M;
1025  link li="";   // link to stdout
1026  fprintf(li,"s:%s,l:%l",1,2);
1027  fprintf(li,"s:%s",l);
1028  fprintf(li,"s:%s",list(l));
1029  fprintf(li,"2l:%2l",list(l));
1030  fprintf(li,"%p",list(l));
1031  fprintf(li,"%;",list(l));
1032  fprintf(li,"%b",M);
1033}
1034
1035//////////////////////////////////////////////////////////////////////////
1036
1037/*
1038proc minres(list #)
1039{
1040  if (size(#) == 2)
1041  {
1042    if (typeof(#[1]) == "ideal" || typeof(#[1]) == "module")
1043    {
1044      if (typeof(#[2] == "int"))
1045      {
1046        return (res(#[1],#[2],1));
1047      }
1048    }
1049  }
1050
1051  if (typeof(#[1]) == "resolution")
1052  {
1053    return minimizeres(#[1]);
1054  }
1055  else
1056  {
1057    return minimizeres(#);
1058  }
1059
1060}
1061*/
Note: See TracBrowser for help on using the repository browser.