source: git/Singular/LIB/standard.lib @ 3b77465

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