source: git/Singular/LIB/standard.lib @ 942c79

spielwiese
Last change on this file since 942c79 was 942c79, checked in by Hans Schönemann <hannes@…>, 20 years ago
*hannes: size -> nclos for indices git-svn-id: file:///usr/local/Singular/svn/trunk@7264 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 32.1 KB
Line 
1//////////////////////////////////////////////////////////////////////////////
2version="$Id: standard.lib,v 1.66 2004-07-15 16:33:46 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 timeStd(i,d)           std(i) if the standard basis computation finished after
17                        d-1 seconds and i otherwhise
18 timeFactorize(p,d)     factorize(p) if the factorization finished after d-1
19                        seconds otherwhise f is considered to be irreducible
20 factorH(p)             changes variables to become the last variable the
21                        principal one in the multivariate factorization and
22                        factorizes then the polynomial
23
24";
25
26//////////////////////////////////////////////////////////////////////////////
27
28proc stdfglm (ideal i, list #)
29"SYNTAX: @code{stdfglm (} ideal_expression @code{)} @*
30         @code{stdfglm (} ideal_expression@code{,} string_expression @code{)}
31TYPE:    ideal
32PURPOSE: computes the standard basis of the ideal in the basering
33         via @code{fglm} (from the ordering given as the second argument
34         to the ordering of the basering).@*
35         If no second argument is given, \"dp\" is used.
36SEE ALSO: fglm, groebner, std, stdhilb
37KEYWORDS: fglm
38EXAMPLE: example stdfglm; shows an example"
39{
40   string os;
41   def dr= basering;
42   if( (size(#)==0) or (typeof(#[1]) != "string") )
43   {
44     os = "dp(" + string( nvars(dr) ) + ")";
45     if ( (find( ordstr(dr), os ) != 0) and (find( ordstr(dr), "a") == 0) )
46     {
47       os= "Dp";
48     }
49     else
50     {
51       os= "dp";
52     }
53   }
54   else { os = #[1]; }
55   execute("ring sr=("+charstr(dr)+"),("+varstr(dr)+"),"+os+";");
56   ideal i= fetch(dr,i);
57   intvec opt= option(get);
58   option(redSB);
59   i=std(i);
60   option(set,opt);
61   setring dr;
62   return (fglm(sr,i));
63}
64example
65{ "EXAMPLE:"; echo = 2;
66   ring r=0,(x,y,z),lp;
67   ideal i=y3+x2,x2y+x2,x3-x2,z4-x2-y;
68   ideal i1=stdfglm(i);         //uses fglm from "dp" to "lp"
69   i1;
70   ideal i2=stdfglm(i,"Dp");    //uses fglm from "Dp" to "lp"
71   i2;
72}
73/////////////////////////////////////////////////////////////////////////////
74
75proc stdhilb(ideal i,list #)
76"SYNTAX: @code{stdhilb (} ideal_expression @code{)} @*
77         @code{stdhilb (} ideal_expression@code{,} intvec_expression @code{)}
78TYPE:    ideal
79PURPOSE: computes the standard basis of the homogeneous ideal in the basering,
80         via a Hilbert driven standard basis computation.@*
81         An optional second argument will be used as 1st Hilbert function.
82ASSUME:  The optional second argument is the first Hilbert series as computed
83         by @code{hilb}.
84SEE ALSO: stdfglm, std, groebner
85KEYWORDS: Hilbert function
86EXAMPLE: example stdhilb;  shows an example"
87{
88   def R=basering;
89
90   if((homog(i)==1)||(ordstr(basering)[1]=="d"))
91   {
92      if ((size(#)!=0)&&(homog(i)==1))
93      {
94         return(std(i,#[1]));
95      }
96      return(std(i));
97   }
98
99   execute("ring S = ("+charstr(R)+"),("+varstr(R)+",@t),dp;");
100   ideal i=homog(imap(R,i),@t);
101   intvec v=hilb(std(i),1);
102   execute("ring T = ("+charstr(R)+"),("+varstr(R)+",@t),("+ordstr(R)+");");
103   ideal i=fetch(S,i);
104   ideal a=std(i,v);
105   setring R;
106   map phi=T,maxideal(1),1;
107   ideal a=phi(a);
108
109   int k,j;
110   poly m;
111   int c=ncols(i);
112
113   for(j=1;j<c;j++)
114   {
115     if(deg(a[j])==0)
116     {
117       a=ideal(1);
118       attrib(a,"isSB",1);
119       return(a);
120     }
121     if(deg(a[j])>0)
122     {
123       m=lead(a[j]);
124       for(k=j+1;k<=c;k++)
125       {
126          if(size(lead(a[k])/m)>0)
127          {
128            a[k]=0;
129          }
130       }
131     }
132   }
133   a=simplify(a,2);
134   attrib(a,"isSB",1);
135   return(a);
136}
137example
138{ "EXAMPLE:"; echo = 2;
139   ring  r=0,(x,y,z),dp;
140   ideal i=y3+x2,x2y+x2,x3-x2,z4-x2-y;
141   ideal i1=stdhilb(i); i1;
142   // the latter computation is equivalent to:
143   intvec v=hilb(i,1);
144   ideal i2=stdhilb(i,v); i2;
145}
146//////////////////////////////////////////////////////////////////////////
147
148proc groebner(def i, list #)
149"SYNTAX: @code{groebner (} ideal_expression @code{)} @*
150         @code{groebner (} module_expression @code{)} @*
151         @code{groebner (} ideal_expression@code{,} int_expression @code{)} @*
152         @code{groebner (} module_expression@code{,} int_expression @code{)}
153TYPE:    type of the first argument
154PURPOSE: computes the standard basis of the first argument @code{I}
155         (ideal or module), by a heuristically chosen method: if the
156         ordering of the current ring is a local ordering, or if it is a
157         non-block ordering and the current ring has no parameters, then
158         @code{std(I)} is returned. Otherwise, @code{I} is mapped into a
159         ring with no parameters and ordering dp, where its Hilbert series
160         is computed. This is followed by a Hilbert-series based std
161         computation in the original ring.
162NOTE: If a 2nd argument @code{wait} is given, then the computation proceeds
163      at most @code{wait} seconds. That is, if no result could be computed in
164      @code{wait} seconds, then the computation is interrupted, 0 is returned,
165      a warning message is displayed, and the global variable
166      @code{groebner_error} is defined.
167SEE ALSO: stdhilb, stdfglm, std
168KEYWORDS: time limit on computations; MP, groebner basis computations
169EXAMPLE: example groebner;  shows an example"
170{
171  def P=basering;
172
173  // we have two arguments -- try to use MPfork links
174  if (size(#) > 0)
175  {
176    if (system("with", "MP"))
177    {
178      if (typeof(#[1]) == "int")
179      {
180        int wait = #[1];
181        int j = 10;
182
183        string bs = nameof(basering);
184        link l_fork = "MPtcp:fork";
185        open(l_fork);
186        write(l_fork, quote(system("pid")));
187        int pid = read(l_fork);
188        write(l_fork, quote(groebner(eval(i))));
189
190        // sleep in small intervalls for appr. one second
191        if (wait > 0)
192        {
193          while(j < 1000000)
194          {
195            if (status(l_fork, "read", "ready", j)) {break;}
196            j = j + j;
197          }
198        }
199
200        // sleep in intervalls of one second from now on
201        j = 1;
202        while (j < wait)
203        {
204          if (status(l_fork, "read", "ready", 1000000)) {break;}
205          j = j + 1;
206        }
207
208        if (status(l_fork, "read", "ready"))
209        {
210          def result = read(l_fork);
211          if (bs != nameof(basering))
212          {
213            def PP = basering;
214            setring P;
215            def result = imap(PP, result);
216            kill PP;
217          }
218          if (defined(groebner_error))
219          {
220            kill(groebner_error);
221          }
222          kill (l_fork);
223        }
224        else
225        {
226          ideal result;
227          if (! defined(groebner_error))
228          {
229            int groebner_error = 1;
230            export groebner_error;
231          }
232          "// ** groebner did not finish";
233          j = system("sh", "kill " + string(pid));
234        }
235        return (result);
236      }
237      else
238      {
239        "// ** groebner needs int as 2nd arg";
240      }
241    }
242    else
243    {
244      "// ** groebner with two args is not supported in this configuration";
245    }
246  }
247
248  // we are still here -- do the actual computation
249  string ordstr_P = ordstr(P);
250  if ((find(ordstr_P,"s") > 0)
251  ||(find(ordstr_P,"M") > 0)
252  ||(find(ordstr_P,"w") > 0)
253  ||(find(ordstr_P,"W") > 0))
254  {
255    //spaeter den lokalen fall ueber lp oder aehnlich behandeln
256    return(std(i));
257  }
258
259  int IsSimple_P;
260  if (system("nblocks") <= 2)
261  {
262    if (find(ordstr_P, "M") <= 0)
263    {
264      IsSimple_P = 1;
265    }
266  }
267  int npars_P = npars(P);
268
269  // return std if no parameters and (dp or wp)
270  if ((npars_P <= 1) && IsSimple_P)
271  {
272    if (find(ordstr_P, "d") > 0)
273    {
274      return (std(i));
275    }
276    if (find(ordstr_P,"w") > 0)
277    {
278      return (std(i));
279    }
280  }
281
282  // reset options
283  intvec opt=option(get);
284  int p_opt;
285  string s_opt = option();
286  option(none);
287  // turn on option(prot) and/or option(mem), if previously set
288  if (find(s_opt, "prot"))
289  {
290    option(prot);
291    p_opt = 1;
292  }
293  if (find(s_opt, "mem"))
294  { option(mem); }
295  if (find(s_opt, "intStrategy"))
296  { option(intStrategy); }
297
298  // construct ring in which first std computation is done
299  string varstr_P = varstr(P);
300  string parstr_P = parstr(P);
301  int is_homog = (homog(i) && (npars_P <= 1));
302  int add_vars = 0;
303  string ri = "ring Phelp =";
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
1035proc timeFactorize(poly i,list #)
1036"USAGE:  timeFactorize(p,d)  poly p , integer d
1037RETURN:  factorize(p) if the factorization finished after d-1
1038         seconds otherwhise f is considered to be irreducible
1039EXAMPLE: example timeFactorize; shows an example
1040"
1041{
1042  def P=basering;
1043  if (size(#) > 0)
1044  {
1045    if (system("with", "MP"))
1046    {
1047      if ((typeof(#[1]) == "int")&&(#[1]))
1048      {
1049        int wait = #[1];
1050        int j = 10;
1051
1052        string bs = nameof(basering);
1053        link l_fork = "MPtcp:fork";
1054        open(l_fork);
1055        write(l_fork, quote(system("pid")));
1056        int pid = read(l_fork);
1057        write(l_fork, quote(timeFactorize(eval(i))));
1058
1059        // sleep in small intervalls for appr. one second
1060        if (wait > 0)
1061        {
1062          while(j < 1000000)
1063          {
1064            if (status(l_fork, "read", "ready", j)) {break;}
1065            j = j + j;
1066          }
1067        }
1068
1069        // sleep in intervalls of one second from now on
1070        j = 1;
1071        while (j < wait)
1072        {
1073          if (status(l_fork, "read", "ready", 1000000)) {break;}
1074          j = j + 1;
1075        }
1076
1077        if (status(l_fork, "read", "ready"))
1078        {
1079          def result = read(l_fork);
1080          if (bs != nameof(basering))
1081          {
1082            def PP = basering;
1083            setring P;
1084            def result = imap(PP, result);
1085            kill PP;
1086          }
1087          kill (l_fork);
1088        }
1089        else
1090        {
1091          list result;
1092          intvec v=1,1;
1093          result[1]=list(1,i);
1094          result[2]=v;
1095          j = system("sh", "kill " + string(pid));
1096        }
1097        return (result);
1098      }
1099    }
1100  }
1101  return(factorH(i));
1102}
1103example
1104{ "EXAMPLE:"; echo = 2;
1105   ring r=0,(x,y),dp;
1106   poly p=((x2+y3)^2+xy6)*((x3+y2)^2+x10y);
1107   p=p^2;
1108   timeFactorize(p,2);
1109   "ERROR: will not run currently:";
1110   //timeFactorize(p,20);
1111}
1112
1113proc timeStd(ideal i,list #)
1114"USAGE:  timeStd(i,d), i ideal, d integer
1115RETURN:  std(i) if the standard basis computation finished after
1116         d-1 seconds and i otherwhise
1117EXAMPLE: example timeStd; shows an example
1118"
1119{
1120  def P=basering;
1121  if (size(#) > 0)
1122  {
1123    if (system("with", "MP"))
1124    {
1125      if ((typeof(#[1]) == "int")&&(#[1]))
1126      {
1127        int wait = #[1];
1128        int j = 10;
1129
1130        string bs = nameof(basering);
1131        link l_fork = "MPtcp:fork";
1132        open(l_fork);
1133        write(l_fork, quote(system("pid")));
1134        int pid = read(l_fork);
1135        write(l_fork, quote(timeStd(eval(i))));
1136
1137        // sleep in small intervalls for appr. one second
1138        if (wait > 0)
1139        {
1140          while(j < 1000000)
1141          {
1142            if (status(l_fork, "read", "ready", j)) {break;}
1143            j = j + j;
1144          }
1145        }
1146        j = 1;
1147        while (j < wait)
1148        {
1149          if (status(l_fork, "read", "ready", 1000000)) {break;}
1150          j = j + 1;
1151        }
1152        if (status(l_fork, "read", "ready"))
1153        {
1154          def result = read(l_fork);
1155          if (bs != nameof(basering))
1156          {
1157            def PP = basering;
1158            setring P;
1159            def result = imap(PP, result);
1160            kill PP;
1161          }
1162          kill (l_fork);
1163        }
1164        else
1165        {
1166          ideal result=i;
1167          j = system("sh", "kill " + string(pid));
1168        }
1169        return (result);
1170      }
1171    }
1172  }
1173  return(std(i));
1174}
1175example
1176{ "EXAMPLE:"; echo = 2;
1177   ring r=32003,(a,b,c,d,e),dp;
1178   int n=6;
1179   ideal i=
1180   a^n-b^n,
1181   b^n-c^n,
1182   c^n-d^n,
1183   d^n-e^n,
1184   a^(n-1)*b+b^(n-1)*c+c^(n-1)*d+d^(n-1)*e+e^(n-1)*a;
1185   timeStd(i,2);
1186   timeStd(i,20);
1187}
1188
1189proc factorH(poly p)
1190"USAGE:  factorH(p)  p poly
1191RETURN:  factorize(p)
1192NOTE:    changes variables to become the last variable the principal
1193         one in the multivariate factorization and factorizes then
1194         the polynomial
1195EXAMPLE: example factorH; shows an example
1196"
1197{
1198   def R=basering;
1199   int i,j;
1200   int n=1;
1201   int d=nrows(coeffs(p,var(1)));
1202   for(i=1;i<=nvars(R);i++)
1203   {
1204      j=nrows(coeffs(p,var(i)));
1205      if(d>j)
1206      {
1207         n=i;
1208         d=j;
1209      }
1210   }
1211   ideal ma=maxideal(1);    //die letzte Variable ist die Hauptvariable
1212   ma[nvars(R)]=var(n);
1213   ma[n]=var(nvars(R));
1214   map phi=R,ma;
1215   list fac=factorize(phi(p));
1216   list re=phi(fac);
1217   return(re);
1218}
1219example
1220{ "EXAMPLE:"; echo = 2;
1221  system("random",992851144);
1222  ring r=32003,(x,y,z,w,t),lp;
1223  poly p=y2w9+yz7t-yz5w4-z2w4t4-w8t3;
1224  factorize(p);  //fast
1225  system("random",992851262);
1226  //factorize(p);  //slow
1227  system("random",992851262);
1228  factorH(p);
1229}
1230//////////////////////////////////////////////////////////////////////////
1231
1232/*
1233proc minres(list #)
1234{
1235  if (size(#) == 2)
1236  {
1237    if (typeof(#[1]) == "ideal" || typeof(#[1]) == "module")
1238    {
1239      if (typeof(#[2] == "int"))
1240      {
1241        return (res(#[1],#[2],1));
1242      }
1243    }
1244  }
1245
1246  if (typeof(#[1]) == "resolution")
1247  {
1248    return minimizeres(#[1]);
1249  }
1250  else
1251  {
1252    return minimizeres(#);
1253  }
1254
1255}
1256*/
Note: See TracBrowser for help on using the repository browser.