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

spielwiese
Last change on this file since b79605 was b79605, checked in by Hans Schönemann <hannes@…>, 23 years ago
*hannes:fixed groebner(hilb) for w-orderings git-svn-id: file:///usr/local/Singular/svn/trunk@5622 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 27.1 KB
Line 
1//////////////////////////////////////////////////////////////////////////////
2version="$Id: standard.lib,v 1.60 2001-08-28 12:54:05 Singular 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=size(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
251  int IsSimple_P;
252  if (system("nblocks") <= 2)
253  {
254    if (find(ordstr_P, "M") <= 0)
255    {
256      IsSimple_P = 1;
257    }
258  }
259  int npars_P = npars(P);
260
261  // return std if no parameters and (dp or wp)
262  if ((npars_P <= 1) && IsSimple_P)
263  {
264    if (find(ordstr_P, "d") > 0)
265    {
266      return (std(i));
267    }
268    if (find(ordstr_P,"w") > 0)
269    {
270      return (std(i));
271    }
272  }
273
274  // reset options
275  intvec opt=option(get);
276  int p_opt;
277  string s_opt = option();
278  option(none);
279  // turn on option(prot) and/or option(mem), if previously set
280  if (find(s_opt, "prot"))
281  {
282    option(prot);
283    p_opt = 1;
284  }
285  if (find(s_opt, "mem"))
286  {
287    option(mem);
288  }
289
290  // construct ring in which first std computation is done
291  string varstr_P = varstr(P);
292  string parstr_P = parstr(P);
293  int is_homog = (homog(i) && (npars_P <= 1));
294  int add_vars = 0;
295  string ri = "ring Phelp =";
296
297  // more than one parameters are converted to ring variables
298  if (npars_P > 1)
299  {
300    ri = ri + string(char(P)) + ",(" + varstr_P + "," + parstr_P;
301    add_vars = npars_P;
302  }
303  else
304  {
305    ri = ri + "(" + charstr(P) + "),(" + varstr_P;
306  }
307
308  // a homogenizing variable is added, if necessary
309  if (! is_homog)
310  {
311    ri = ri + ",@t";
312    add_vars = add_vars + 1;
313  }
314  // ordering is set to (dp, C)
315  ri = ri + "),(dp,C);";
316
317  // change the ring
318  execute(ri);
319
320  // get ideal from previous ring
321  if (is_homog)
322  {
323    ideal qh = imap(P, i);
324  }
325  else
326  {
327    // and homogenize
328    ideal qh=homog(imap(P,i),@t);
329  }
330
331  // compute std and hilbert series
332  if (p_opt)
333  {
334    "std in " + ri[13, size(ri) - 13];
335  }
336  intvec hi=hilb(std(qh),1);
337
338  if (add_vars == 0)
339  {
340    // no additional variables were introduced
341    setring P; // can immediately change to original ring
342    // simply compute std with hilbert series in original ring
343    if (p_opt)
344    {
345      "std with hilb in basering";
346    }
347    i = std(i, hi);
348  }
349  else
350  {
351    // additional variables were introduced
352    // need another intermediate ring
353    ri = "ring Phelp1 = (" + charstr(Phelp)
354      + "),(" + varstr(Phelp) + "),(" + ordstr_P;
355
356    // for lp wit at most one parameter, we do not need a block ordering
357    if ( ! (IsSimple_P && (add_vars <2) && find(ordstr_P, "l")))
358    {
359      // need block ordering
360      ri = ri + ", dp(" + string(add_vars) + ")";
361    }
362    ri = ri + ");";
363
364    // change to intermediate ring
365    execute(ri);
366    ideal qh = imap(Phelp, qh);
367    kill Phelp;
368    if (p_opt)
369    {
370      "std with hilb in " + ri[14,size(ri)-14];
371    }
372    // compute std with Hilbert series
373    qh = std(qh, hi);
374    // subst 1 for homogenizing var
375    if (!is_homog)
376    {
377      if (p_opt)
378      {
379        "dehomogenization";
380      }
381      qh = subst(qh, @t, 1);
382    }
383
384    // go back to original ring
385    setring P;
386    // get ideal, delete zeros and clean SB
387    if (p_opt)
388    {
389      "imap to original ring";
390    }
391    i = imap(Phelp1,qh);
392    if (p_opt)
393    {
394      "simplification";
395    }
396    i = simplify(i, 34);
397    kill Phelp1;
398  }
399
400  // clean-up time
401  option(set, opt);
402  if (find(s_opt, "redSB") > 0)
403  {
404    if (p_opt)
405    {
406      "interreduction";
407    }
408    i=interred(i);
409  }
410  attrib(i, "isSB", 1);
411  return (i);
412}
413example
414{ "EXAMPLE: "; echo=2;
415  ring r=0,(a,b,c,d),lp;
416  option(prot);
417  ideal i=a+b+c+d,ab+ad+bc+cd,abc+abd+acd+bcd,abcd-1; // cyclic 4
418  groebner(i);
419  ring rp=(0,a,b),(c,d), lp;
420  ideal i=imap(r,i);
421  ideal j=groebner(i);
422  option(noprot);
423  j; simplify(j,1); std(i);
424  if (system("with","MP")) {groebner(i,0);}
425  defined(groebner_error);
426}
427//////////////////////////////////////////////////////////////////////////
428
429proc res(list #)
430"@c we do texinfo here:
431@cindex resolution, computation of
432@table @code
433@item @strong{Syntax:}
434@code{res (} ideal_expression@code{,} int_expression @code{[,} any_expression @code{])}
435@*@code{res (} module_expression@code{,} int_expression @code{[,} any_expression @code{])}
436@item @strong{Type:}
437resolution
438@item @strong{Purpose:}
439computes a (possibly minimal) free resolution of an ideal or module using
440a heuristically chosen method.
441@* The second (int) argument (say, @code{k}) specifies the length of
442the resolution. If it is not positive then @code{k} is assumed to be the
443number of variables of the basering.
444@* If a third argument is given, the returned resolution is minimized.
445
446Depending on the input, the returned resolution is computed using the
447following methods:
448@table @asis
449@item @strong{quotient rings:}
450@code{nres} (classical method using syzygies) , see @ref{nres}.
451
452@item @strong{homogeneous ideals and k=0:}
453@code{lres} (La'Scala's method), see @ref{lres}.
454
455@item @strong{not minimized resolution and (homogeneous input with k not 0, or local rings):}
456@code{sres} (Schreyer's method), see @ref{sres}.
457
458@item @strong{all other inputs:}
459@code{mres} (classical method), see @ref{mres}.
460@end table
461@item @strong{Note:}
462Accessing single elements of a resolution may require that some partial
463computations have to be finished and may therefore take some time.
464@end table
465@c ref
466See also
467@ref{betti};
468@ref{ideal};
469@ref{minres};
470@ref{module};
471@ref{mres};
472@ref{nres};
473@ref{lres};
474@ref{hres};
475@ref{sres}.
476@ref{resolution}
477@c ref
478"
479{
480   def P=basering;
481   if (size(#) < 2)
482   {
483     ERROR("res: need at least two arguments: ideal/module, int");
484   }
485
486   def m=#[1]; //the ideal or module
487   int i=#[2]; //the length of the resolution
488   if (i< 0) { i=0;}
489
490   string varstr_P = varstr(P);
491
492   int p_opt;
493   string s_opt = option();
494   // set p_opt, if option(prot) is set
495   if (find(s_opt, "prot"))
496   {
497     p_opt = 1;
498   }
499
500   if(size(ideal(basering)) > 0)
501   {
502     // the quick hack for qrings - seems to fit most needs
503     // (lres is not implemented for qrings, sres is not so efficient)
504     if (p_opt) { "using nres";}
505     return(nres(m,i));
506   }
507
508   if(homog(m)==1)
509   {
510      resolution re;
511      if (((i==0) or (i>=nvars(basering))) && typeof(m) != "module")
512      {
513        //LaScala for the homogeneous case and i == 0
514        if (p_opt) { "using lres";}
515        re=lres(m,i);
516        if(size(#)>2)
517        {
518           re=minres(re);
519        }
520      }
521      else
522      {
523        if(size(#)>2)
524        {
525          if (p_opt) { "using mres";}
526          re=mres(m,i);
527        }
528        else
529        {
530          if (p_opt) { "using sres";}
531          re=sres(std(m),i);
532        }
533      }
534      return(re);
535   }
536
537   //mres for the global non homogeneous case
538   if(find(ordstr(P),"s")==0)
539   {
540      string ri= "ring Phelp ="
541                  +string(char(P))+",("+varstr_P+"),(dp,C);";
542      execute(ri);
543      def m=imap(P,m);
544      if (p_opt) { "using mres in another ring";}
545      list re=mres(m,i);
546      setring P;
547      resolution result=imap(Phelp,re);
548      if (size(#) > 2) {result = minres(result);}
549      return(result);
550   }
551
552   //sres for the local case and not minimal resolution
553   if(size(#)<=2)
554   {
555      string ri= "ring Phelp ="
556                  +string(char(P))+",("+varstr_P+"),(ls,c);";
557      execute(ri);
558      def m=imap(P,m);
559      m=std(m);
560      if (p_opt) { "using sres in another ring";}
561      list re=sres(m,i);
562      setring P;
563      resolution result=imap(Phelp,re);
564      return(result);
565   }
566
567   //mres for the local case and minimal resolution
568   string ri= "ring Phelp ="
569                  +string(char(P))+",("+varstr_P+"),(ls,C);";
570   execute(ri);
571   def m=imap(P,m);
572    if (p_opt) { "using mres in another ring";}
573   list re=mres(m,i);
574   setring P;
575   resolution result=imap(Phelp,re);
576   result = minres(result);
577   return(result);
578}
579example
580{"EXAMPLE:"; echo = 2;
581  ring r=0,(x,y,z),dp;
582  ideal i=xz,yz,x3-y3;
583  def l=res(i,0); // homogeneous ideal: uses lres
584  l;
585  print(betti(l), "betti"); // input to betti may be of type resolution
586  l[2];         // element access may take some time
587  i=i,x+1;
588  l=res(i,0);   // inhomogeneous ideal: uses mres
589  l;
590  ring rs=0,(x,y,z),ds;
591  ideal i=imap(r,i);
592  def l=res(i,0); // local ring not minimized: uses sres
593  l;
594  res(i,0,0);     // local ring and minimized: uses mres
595}
596/////////////////////////////////////////////////////////////////////////
597
598proc quot (m1,m2,list #)
599"SYNTAX: @code{quot (} module_expression@code{,} module_expression @code{)}
600         @*@code{quot (} module_expression@code{,} module_expression@code{,}
601            int_expression @code{)}
602         @*@code{quot (} ideal_expression@code{,} ideal_expression @code{)}
603         @*@code{quot (} ideal_expression@code{,} ideal_expression@code{,}
604            int_expression @code{)}
605TYPE:    ideal
606SYNTAX:  @code{quot (} module_expression@code{,} ideal_expression @code{)}
607TYPE:    module
608PURPOSE: computes the quotient of the 1st and the 2nd argument.
609         If a 3rd argument 'n' is given the n-th method is used
610         (n=1...5).
611SEE ALSO: quotient
612EXAMPLE: example quot; shows an example"
613{
614  if (((typeof(m1)!="ideal") and (typeof(m1)!="module"))
615     or ((typeof(m2)!="ideal") and (typeof(m2)!="module")))
616  {
617    "USAGE:   quot(m1, m2[, n]); m1, m2 two submodules of k^s,";
618    "         n (optional) integer (1<= n <=5)";
619    "RETURN:  the quotient of m1 and m2";
620    "EXAMPLE: example quot; shows an example";
621    return();
622  }
623  if (typeof(m1)!=typeof(m2))
624  {
625    return(quotient(m1,m2));
626  }
627  if (size(#)>0)
628  {
629    if (typeof(#[1])=="int" )
630    {
631      return(quot1(m1,m2,#[1]));
632    }
633  }
634  else
635  {
636    return(quot1(m1,m2,2));
637  }
638}
639example
640{ "EXAMPLE:"; echo = 2;
641  ring r=181,(x,y,z),(c,ls);
642  ideal id1=maxideal(4);
643  ideal id2=x2+xyz,y2-z3y,z3+y5xz;
644  option(prot);
645  ideal id3=quotient(id1,id2);
646  id3;
647  ideal id4=quot(id1,id2,1);
648  id4;
649  ideal id5=quot(id1,id2,2);
650  id5;
651}
652
653static proc quot1 (module m1, module m2,int n)
654"USAGE:   quot1(m1, m2, n); m1, m2 two submodules of k^s,
655         n integer (1<= n <=5)
656RETURN:  the quotient of m1 and m2
657EXAMPLE: example quot1; shows an example"
658{
659  if (n==1)
660  {
661    return(quotient1(m1,m2));
662  }
663  else
664  {
665    if (n==2)
666    {
667      return(quotient2(m1,m2));
668    }
669    else
670    {
671      if (n==3)
672      {
673        return(quotient3(m1,m2));
674      }
675      else
676      {
677        if (n==4)
678        {
679          return(quotient4(m1,m2));
680        }
681        else
682        {
683          if (n==5)
684          {
685            return(quotient5(m1,m2));
686          }
687          else
688          {
689            return(quotient(m1,m2));
690          }
691        }
692      }
693    }
694  }
695}
696example
697{ "EXAMPLE:"; echo = 2;
698  ring r=181,(x,y,z),(c,ls);
699  ideal id1=maxideal(4);
700  ideal id2=x2+xyz,y2-z3y,z3+y5xz;
701  option(prot);
702  ideal id6=quotient(id1,id2);
703  id6;
704  ideal id7=quot1(id1,id2,1);
705  id7;
706  ideal id8=quot1(id1,id2,2);
707  id8;
708}
709
710static proc quotient0(module a,module b)
711{
712  module mm=b+a;
713  resolution rs=lres(mm,0);
714  list I=list(rs);
715  matrix M=I[2];
716  matrix A[1][nrows(M)]=M[1..nrows(M),1];
717  ideal i=A;
718  return (i);
719}
720proc quotient1(module a,module b)  //17sec
721"USAGE:   quotient1(m1, m2); m1, m2 two submodules of k^s,
722RETURN:  the quotient of m1 and m2"
723{
724  int i;
725  a=std(a);
726  module dummy;
727  module B=NF(b,a)+dummy;
728  ideal re=quotient(a,module(B[1]));
729  for(i=2;i<=size(B);i++)
730  {
731     re=intersect1(re,quotient(a,module(B[i])));
732  }
733  return(re);
734}
735proc quotient2(module a,module b)    //13sec
736"USAGE:   quotient2(m1, m2); m1, m2 two submodules of k^s,
737RETURN:  the quotient of m1 and m2"
738{
739  a=std(a);
740  module dummy;
741  module bb=NF(b,a)+dummy;
742  int i=size(bb);
743  ideal re=quotient(a,module(bb[i]));
744  bb[i]=0;
745  module temp;
746  module temp1;
747  module bbb;
748  int mx;
749  i=i-1;
750  while (1)
751  {
752    if (i==0) break;
753    temp = a+bb*re;
754    temp1 = lead(interred(temp));
755    mx=ncols(a);
756    if (ncols(temp1)>ncols(a))
757    {
758      mx=ncols(temp1);
759    }
760    temp1 = matrix(temp1,1,mx)-matrix(lead(a),1,mx);
761    temp1 = dummy+temp1;
762    if (deg(temp1[1])<0) break;
763    re=intersect1(re,quotient(a,module(bb[i])));
764    bb[i]=0;
765    i = i-1;
766  }
767  return(re);
768}
769proc quotient3(module a,module b)   //89sec
770"USAGE:   quotient3(m1, m2); m1, m2 two submodules of k^s,
771         only for global rings
772RETURN:  the quotient of m1 and m2"
773{
774  string s="ring @newr=("+charstr(basering)+
775           "),("+varstr(basering)+",@t,@w),dp;";
776  def @newP=basering;
777  execute(s);
778  module b=imap(@newP,b);
779  module a=imap(@newP,a);
780  int i;
781  int j=size(b);
782  vector @b;
783  for(i=1;i<=j;i++)
784  {
785     @b=@b+@t^(i-1)*@w^(j-i+1)*b[i];
786  }
787  ideal re=quotient(a,module(@b));
788  setring @newP;
789  ideal re=imap(@newr,re);
790  return(re);
791}
792proc quotient5(module a,module b)   //89sec
793"USAGE:   quotient5(m1, m2); m1, m2 two submodules of k^s,
794         only for global rings
795RETURN:  the quotient of m1 and m2"
796{
797  string s="ring @newr=("+charstr(basering)+
798           "),("+varstr(basering)+",@t),dp;";
799  def @newP=basering;
800  execute(s);
801  module b=imap(@newP,b);
802  module a=imap(@newP,a);
803  int i;
804  int j=size(b);
805  vector @b;
806  for(i=1;i<=j;i++)
807  {
808     @b=@b+@t^(i-1)*b[i];
809  }
810  @b=homog(@b,@w);
811  ideal re=quotient(a,module(@b));
812  setring @newP;
813  ideal re=imap(@newr,re);
814  return(re);
815}
816proc quotient4(module a,module b)   //95sec
817"USAGE:   quotient4(m1, m2); m1, m2 two submodules of k^s,
818         only for global rings
819RETURN:  the quotient of m1 and m2"
820{
821  string s="ring @newr=("+charstr(basering)+
822           "),("+varstr(basering)+",@t),dp;";
823  def @newP=basering;
824  execute(s);
825  module b=imap(@newP,b);
826  module a=imap(@newP,a);
827  int i;
828  vector @b=b[1];
829  for(i=2;i<=size(b);i++)
830  {
831     @b=@b+@t^(i-1)*b[i];
832  }
833  matrix sy=modulo(@b,a);
834  ideal re=sy;
835  setring @newP;
836  ideal re=imap(@newr,re);
837  return(re);
838}
839static proc intersect1(ideal i,ideal j)
840{
841  def R=basering;
842  execute("ring gnir = ("+charstr(basering)+"),
843                       ("+varstr(basering)+",@t),(C,dp);");
844  ideal i=var(nvars(basering))*imap(R,i)+(var(nvars(basering))-1)*imap(R,j);
845  ideal j=eliminate(i,var(nvars(basering)));
846  setring R;
847  map phi=gnir,maxideal(1);
848  return(phi(j));
849}
850
851//////////////////////////////////////////////////////////////////
852///
853/// sprintf, fprintf printf
854///
855proc sprintf(string fmt, list #)
856"SYNTAX:  @code{sprintf (} string_expression @code{[,} any_expressions
857               @code{] )}
858RETURN:   string
859PURPOSE:  @code{sprintf(fmt,...);} performs output formatting. The first
860          argument is a format control string. Additional arguments may be
861          required, depending on the content of the control string. A series
862          of output characters is generated as directed by the control string;
863          these characters are returned as a string. @*
864          The control string @code{fmt} is simply text to be copied,
865          except that the string may contain conversion specifications.@*
866          Do @code{help print;} for a listing of valid conversion
867          specifications. As an addition to the conversions of @code{print},
868          the @code{%n} and @code{%2} conversion specification does not
869          consume an additional argument, but simply generates a newline
870          character.
871NOTE:     If one of the additional arguments is a list, then it should be
872          enclosed once more into a @code{list()} command, since passing a list
873          as an argument flattens the list by one level.
874SEE ALSO: fprintf, printf, print, string
875EXAMPLE : example sprintf; shows an example
876"
877{
878  int sfmt = size(fmt);
879  if (sfmt  <= 1)
880  {
881    return (fmt);
882  }
883  int next, l, nnext;
884  string ret;
885  list formats = "%l", "%s", "%2l", "%2s", "%t", "%;", "%p", "%b", "%n", "%2";
886  while (1)
887  {
888    if (size(#) <= 0)
889    {
890      return (ret + fmt);
891    }
892    nnext = 0;
893    while (nnext < sfmt)
894    {
895      nnext = find(fmt, "%", nnext + 1);
896      if (nnext == 0)
897      {
898        next = 0;
899        break;
900      }
901      l = 1;
902      while (l <= size(formats))
903      {
904        next = find(fmt, formats[l], nnext);
905        if (next == nnext) break;
906        l++;
907      }
908      if (next == nnext) break;
909    }
910    if (next == 0)
911    {
912      return (ret + fmt);
913    }
914    if (formats[l] != "%2" && formats[l] != "%n")
915    {
916      ret = ret + fmt[1, next - 1] + print(#[1], formats[l]);
917      # = delete(#, 1);
918    }
919    else
920    {
921      ret = ret + fmt[1, next - 1] + print("", "%2s");
922    }
923    if (size(fmt) <= (next + size(formats[l]) - 1))
924    {
925      return (ret);
926    }
927    fmt = fmt[next + size(formats[l]), size(fmt)-next-size(formats[l]) + 1];
928  }
929}
930example
931{ "EXAMPLE:"; echo=2;
932  ring r=0,(x,y,z),dp;
933  module m=[1,y],[0,x+z];
934  intmat M=betti(mres(m,0));
935  list l = r, m, M;
936  string s = sprintf("s:%s,%n l:%l", 1, 2); s;
937  s = sprintf("s:%n%s", l); s;
938  s = sprintf("s:%2%s", list(l)); s;
939  s = sprintf("2l:%n%2l", list(l)); s;
940  s = sprintf("%p", list(l)); s;
941  s = sprintf("%;", list(l)); s;
942  s = sprintf("%b", M); s;
943}
944
945proc printf(string fmt, list #)
946"SYNTAX:  @code{printf (} string_expression @code{[,} any_expressions@code{] )}
947RETURN:   none
948PURPOSE:  @code{printf(fmt,...);} performs output formatting. The first
949          argument is a format control string. Additional arguments may be
950          required, depending on the content of the control string. A series
951          of output characters is generated as directed by the control string;
952          these characters are displayed (i.e., printed to standard out). @*
953          The control string @code{fmt} is simply text to be copied, except
954          that the string may contain conversion specifications. @*
955          Do @code{help print;} for a listing of valid conversion
956          specifications. As an addition to the conversions of @code{print},
957          the @code{%n} and @code{%2} conversion specification does not
958          consume an additional argument, but simply generates a newline
959          character.
960NOTE:     If one of the additional arguments is a list, then it should be
961          enclosed once more into a @code{list()} command, since passing a
962          list as an argument flattens the list by one level.
963SEE ALSO: sprintf, fprintf, print, string
964EXAMPLE : example printf; shows an example
965"
966{
967  write("", sprintf(fmt, #));
968}
969example
970{ "EXAMPLE:"; echo=2;
971  ring r=0,(x,y,z),dp;
972  module m=[1,y],[0,x+z];
973  intmat M=betti(mres(m,0));
974  list l=r,m,M;
975  printf("s:%s,l:%l",1,2);
976  printf("s:%s",l);
977  printf("s:%s",list(l));
978  printf("2l:%2l",list(l));
979  printf("%p",list(l));
980  printf("%;",list(l));
981  printf("%b",M);
982}
983
984
985proc fprintf(link l, string fmt, list #)
986"SYNTAX:  @code{fprintf (} link_expression@code{,} string_expression @code{[,}
987                any_expressions@code{] )}
988RETURN:   none
989PURPOSE:  @code{fprintf(l,fmt,...);} performs output formatting.
990          The second argument is a format control string. Additional
991          arguments may be required, depending on the content of the
992          control string. A series of output characters is generated as
993          directed by the control string; these characters are
994          written to the link l.
995          The control string @code{fmt} is simply text to be copied, except
996          that the string may contain conversion specifications.@*
997          Do @code{help print;} for a listing of valid conversion
998          specifications. As an addition to the conversions of @code{print},
999          the @code{%n} and @code{%2} conversion specification does not
1000          consume an additional argument, but simply generates a newline
1001          character.
1002NOTE:     If one of the additional arguments is a list, then it should be
1003          enclosed once more into a @code{list()} command, since passing
1004          a list as an argument flattens the list by one level.
1005SEE ALSO: sprintf, printf, print, string
1006EXAMPLE : example fprintf; shows an example
1007"
1008{
1009  write(l, sprintf(fmt, #));
1010}
1011example
1012{ "EXAMPLE:"; echo=2;
1013  ring r=0,(x,y,z),dp;
1014  module m=[1,y],[0,x+z];
1015  intmat M=betti(mres(m,0));
1016  list l=r,m,M;
1017  link li="";   // link to stdout
1018  fprintf(li,"s:%s,l:%l",1,2);
1019  fprintf(li,"s:%s",l);
1020  fprintf(li,"s:%s",list(l));
1021  fprintf(li,"2l:%2l",list(l));
1022  fprintf(li,"%p",list(l));
1023  fprintf(li,"%;",list(l));
1024  fprintf(li,"%b",M);
1025}
1026
1027/*
1028proc minres(list #)
1029{
1030  if (size(#) == 2)
1031  {
1032    if (typeof(#[1]) == "ideal" || typeof(#[1]) == "module")
1033    {
1034      if (typeof(#[2] == "int"))
1035      {
1036        return (res(#[1],#[2],1));
1037      }
1038    }
1039  }
1040
1041  if (typeof(#[1]) == "resolution")
1042  {
1043    return minimizeres(#[1]);
1044  }
1045  else
1046  {
1047    return minimizeres(#);
1048  }
1049
1050}
1051*/
Note: See TracBrowser for help on using the repository browser.