source: git/Singular/LIB/toric.lib @ 334c21f

spielwiese
Last change on this file since 334c21f was f79c9d, checked in by Anne Frühbis-Krüger <anne@…>, 23 years ago
*anne: corrected some typing mistakes git-svn-id: file:///usr/local/Singular/svn/trunk@5194 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 27.3 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="$Id: toric.lib,v 1.11 2001-02-06 11:43:23 anne Exp $";
3category="Commutative Algebra";
4info="
5LIBRARY: toric.lib   Standard Basis of Toric Ideals
6AUTHOR:  Christine Theis, email: ctheis@math.uni-sb.de
7
8PROCEDURES:
9toric_ideal(A,..);   computes the toric ideal of A
10toric_std(ideal I);  standard basis of I by a specialized Buchberger algorithm
11";
12
13///////////////////////////////////////////////////////////////////////////////
14
15static proc toric_ideal_1(intmat A, string alg)
16{
17  ideal I;
18  // to be returned
19
20  // check suitability of actual basering
21  if(nvars(basering)<ncols(A))
22  {
23    "ERROR: The number of matrix columns is smaller than the number of ring variables.";
24    return(I);
25  }
26
27  // check suitability of actual term ordering
28  // the "module ordering" c or C is ignored
29  string singular_ord=ordstr(basering);
30  string test_ord;
31  string external_ord="";
32  int i,j;
33  intvec weightvec;
34
35  if(find(singular_ord,"lp")==1)
36  {
37    external_ord="W_LEX";
38    for(i=1;i<=nvars(basering);i++)
39    {
40      weightvec[i]=0;
41    }
42    test_ord="lp("+string(nvars(basering))+"),";
43    if(singular_ord!=(test_ord+"C") && singular_ord!=(test_ord+"c"))
44    {
45      "Warning: Block orderings are not supported; lp used for computation.";
46    }
47  }
48  if(external_ord=="" && find(singular_ord,"lp")==3)
49  {
50    external_ord="W_LEX";
51    for(i=1;i<=nvars(basering);i++)
52    {
53      weightvec[i]=0;
54    }
55    test_ord="lp("+string(nvars(basering))+"),";
56    if(singular_ord!=("C"+test_ord) && singular_ord!=("c"+test_ord))
57    {
58      "Warning: Block orderings are not supported; lp used for computation.";
59    }
60  }
61
62  if(external_ord=="" && find(singular_ord,"dp")==1)
63  {
64    external_ord="W_REV_LEX";
65    for(i=1;i<=nvars(basering);i++)
66    {
67      weightvec[i]=1;
68    }
69    test_ord="dp("+string(nvars(basering))+"),";
70    if(singular_ord!=(test_ord+"C") && singular_ord!=(test_ord+"c"))
71    {
72      "Warning: Block orderings are not supported; dp used for computation.";
73    }
74  }
75  if(external_ord=="" && find(singular_ord,"dp")==3)
76  {
77    external_ord="W_REV_LEX";
78    for(i=1;i<=nvars(basering);i++)
79    {
80      weightvec[i]=1;
81    }
82    test_ord="dp("+string(nvars(basering))+"),";
83    if(singular_ord!=("C"+test_ord) && singular_ord!=("c"+test_ord))
84    {
85      "Warning: Block orderings are not supported; dp used for computation.";
86    }
87  }
88
89  if(external_ord=="" && find(singular_ord,"Dp")==1)
90  {
91    external_ord="W_LEX";
92    for(i=1;i<=nvars(basering);i++)
93    {
94      weightvec[i]=1;
95    }
96    test_ord="Dp("+string(nvars(basering))+"),";
97    if(singular_ord!=(test_ord+"C") && singular_ord!=(test_ord+"c"))
98    {
99      "Warning: Block orderings are not supported; Dp used for computation.";
100    }
101  }
102  if(external_ord=="" && find(singular_ord,"Dp")==3)
103  {
104    external_ord="W_LEX";
105    for(i=1;i<=nvars(basering);i++)
106    {
107      weightvec[i]=1;
108    }
109    test_ord="Dp("+string(nvars(basering))+"),";
110    if(singular_ord!=("C"+test_ord) && singular_ord!=("c"+test_ord))
111    {
112      "Warning: Block orderings are not supported; Dp used for computation.";
113    }
114  }
115
116  int pos;
117  string number_string;
118
119  if(external_ord=="" && find(singular_ord,"wp")==1)
120  {
121    external_ord="W_REV_LEX";
122    pos=3;
123    for(i=1;i<=nvars(basering);i++)
124    {
125      pos++;
126      number_string="";
127      while(singular_ord[pos]!=",")
128      {
129        number_string=number_string+singular_ord[pos];
130        pos++;
131      }
132      execute("weightvec[i]="+number_string);
133    }
134    test_ord="wp("+string(weightvec)+"),";
135    if(singular_ord!=(test_ord+"C") && singular_ord!=(test_ord+"c"))
136    {
137      "Warning: Block orderings are not supported; wp("+string(weightvec)+") used for computation.";
138    }
139  }
140  if(external_ord=="" && find(singular_ord,"wp")==3)
141  {
142    external_ord="W_REV_LEX";
143    pos=5;
144    for(i=1;i<=nvars(basering);i++)
145    {
146      pos++;
147      number_string="";
148      while(singular_ord[pos]!=",")
149      {
150        number_string=number_string+singular_ord[pos];
151        pos++;
152      }
153      execute("weightvec[i]="+number_string);
154    }
155    test_ord="wp("+string(weightvec)+"),";
156    if(singular_ord!=("C"+test_ord) && singular_ord!=("c"+test_ord))
157    {
158      "Warning: Block orderings are not supported; wp("+string(weightvec)+") used for computation.";
159    }
160  }
161
162  if(external_ord=="" && find(singular_ord,"Wp")==1)
163  {
164    external_ord="W_LEX";
165    pos=3;
166    for(i=1;i<=nvars(basering);i++)
167    {
168      pos++;
169      number_string="";
170      while(singular_ord[pos]!=",")
171      {
172        number_string=number_string+singular_ord[pos];
173        pos++;
174      }
175      execute("weightvec[i]="+number_string);
176    }
177    test_ord="Wp("+string(weightvec)+"),";
178    if(singular_ord!=(test_ord+"C") && singular_ord!=(test_ord+"c"))
179    {
180      "Warning: Block orderings are not supported; Wp("+string(weightvec)+") used for computation.";
181    }
182  }
183  if(external_ord=="" && find(singular_ord,"Wp")==3)
184  {
185    external_ord="W_LEX";
186    pos=5;
187    for(i=1;i<=nvars(basering);i++)
188    {
189      pos++;
190      number_string="";
191      while(singular_ord[pos]!=",")
192      {
193        number_string=number_string+singular_ord[pos];
194        pos++;
195      }
196      execute("weightvec[i]="+number_string);
197    }
198    test_ord="Wp("+string(weightvec)+"),";
199    if(singular_ord!=("C"+test_ord) && singular_ord!=("c"+test_ord))
200    {
201      "Warning: Block orderings are not supported; Wp("+string(weightvec)+") used for computation.";
202    }
203  }
204
205  if(external_ord=="")
206  {
207    "ERROR: The term ordering of the actual basering is not supported.";
208    return(I);
209  }
210
211  // check algorithm
212  if(alg=="ct" || alg=="pct")
213    // these algorithms will not cause an error in the external program;
214    // however, they do not compute the toric ideal of A, but of an
215    // extended matrix
216  {
217    "ERROR: The chosen algorithm is not suitable.";
218    return(I);
219  }
220
221  // create temporary file with which the external program is called
222
223  int dummy;
224  int process=system("pid");
225  string matrixfile="temp_MATRIX"+string(process);
226  link MATRIX=":w "+matrixfile;
227  open(MATRIX);
228
229  write(MATRIX,"MATRIX","columns:",ncols(A),"cost vector:");
230  for(j=1;j<=ncols(A);j++)
231  {
232    write(MATRIX,weightvec[j]);
233  }
234  write(MATRIX,"rows:",nrows(A),"matrix:");
235  for(i=1;i<=nrows(A);i++)
236  {
237    for(j=1;j<=ncols(A);j++)
238    {
239      write(MATRIX,A[i,j]);
240    }
241  }
242
243  // search for positive row space vector, if required by the
244  // algorithm
245  int found=0;
246  if((alg=="blr") || (alg=="hs"))
247  {
248    for(i=1;i<=nrows(A);i++)
249    {
250      found=i;
251      for(j=1;j<=ncols(A);j++)
252      {
253        if(A[i,j]<=0)
254        {
255          found=0;
256        }
257      }
258      if(found>0)
259      {
260        break;
261      }
262    }
263    if(found==0)
264    {
265      "ERROR: The chosen algorithm needs a positive vector in the row space of the matrix.";
266      close(MATRIX);
267      dummy=system("sh","rm -f "+matrixfile);
268      return(I);
269    }
270    write(MATRIX,"positive row space vector:");
271    for(j=1;j<=ncols(A);j++)
272    {
273      write(MATRIX,A[found,j]);
274    }
275  }
276  close(MATRIX);
277
278
279  // call external program
280  dummy=system("sh","toric_ideal -alg "+alg+" "+matrixfile);
281
282  // read toric ideal from created file
283  link TORIC_IDEAL=":r "+matrixfile+".GB."+alg;
284  string toric_id=read(TORIC_IDEAL);
285
286  int generators;
287  pos=find(toric_id,"size");
288  pos=find(toric_id,":",pos);
289  pos++;
290
291  while(toric_id[pos]==" " || toric_id[pos]==newline)
292  {
293    pos++;
294  }
295  number_string="";
296  while(toric_id[pos]!=" " && toric_id[pos]!=newline)
297  {
298    number_string=number_string+toric_id[pos];
299    pos++;
300  }
301  execute("generators="+number_string+";");
302
303  intvec v;
304  poly head;
305  poly tail;
306
307  pos=find(toric_id,"basis");
308  pos=find(toric_id,":",pos);
309  pos++;
310
311  for(i=1;i<=generators;i++)
312  {
313    head=1;
314    tail=1;
315
316    for(j=1;j<=ncols(A);j++)
317    {
318      while(toric_id[pos]==" " || toric_id[pos]==newline)
319      {
320        pos++;
321      }
322      number_string="";
323      while(toric_id[pos]!=" " && toric_id[pos]!=newline)
324      {
325        number_string=number_string+toric_id[pos];
326        pos++;
327      }
328      execute("v[j]="+number_string+";");
329      if(v[j]<0)
330      {
331        tail=tail*var(j)^(-v[j]);
332      }
333      if(v[j]>0)
334      {
335        head=head*var(j)^v[j];
336      }
337    }
338    I[i]=head-tail;
339  }
340
341  // delete all created files
342  dummy=system("sh","rm -f "+matrixfile);
343  dummy=system("sh","rm -f "+matrixfile+".GB."+alg);
344
345  return(I);
346}
347///////////////////////////////////////////////////////////////////////////////
348
349static proc toric_ideal_2(intmat A, string alg, intvec prsv)
350{
351  ideal I;
352  // to be returned
353
354  // check arguments
355  if(size(prsv)<ncols(A))
356  {
357    "ERROR: The number of matrix columns does not equal the size of the positive row space vector.";
358    return(I);
359  }
360
361  // check suitability of actual basering
362  if(nvars(basering)!=ncols(A))
363  {
364    "ERROR: The number of matrix columns is smaller than the number of ring variables.";
365    return(I);
366  }
367
368  // check suitability of actual term ordering
369  // the "module ordering" c or C is ignored
370  string singular_ord=ordstr(basering);
371  string test_ord;
372  string external_ord="";
373  int i,j;
374  intvec weightvec;
375
376  if(find(singular_ord,"lp")==1)
377  {
378    external_ord="W_LEX";
379    for(i=1;i<=nvars(basering);i++)
380    {
381      weightvec[i]=0;
382    }
383    test_ord="lp("+string(nvars(basering))+"),";
384    if(singular_ord!=(test_ord+"C") && singular_ord!=(test_ord+"c"))
385    {
386      "Warning: Block orderings are not supported; lp used for computation.";
387    }
388  }
389  if(external_ord=="" && find(singular_ord,"lp")==3)
390  {
391    external_ord="W_LEX";
392    for(i=1;i<=nvars(basering);i++)
393    {
394      weightvec[i]=0;
395    }
396    test_ord="lp("+string(nvars(basering))+"),";
397    if(singular_ord!=("C"+test_ord) && singular_ord!=("c"+test_ord))
398    {
399      "Warning: Block orderings are not supported; lp used for computation.";
400    }
401  }
402
403  if(external_ord=="" && find(singular_ord,"dp")==1)
404  {
405    external_ord="W_REV_LEX";
406    for(i=1;i<=nvars(basering);i++)
407    {
408      weightvec[i]=1;
409    }
410    test_ord="dp("+string(nvars(basering))+"),";
411    if(singular_ord!=(test_ord+"C") && singular_ord!=(test_ord+"c"))
412    {
413      "Warning: Block orderings are not supported; dp used for computation.";
414    }
415  }
416  if(external_ord=="" && find(singular_ord,"dp")==3)
417  {
418    external_ord="W_REV_LEX";
419    for(i=1;i<=nvars(basering);i++)
420    {
421      weightvec[i]=1;
422    }
423    test_ord="dp("+string(nvars(basering))+"),";
424    if(singular_ord!=("C"+test_ord) && singular_ord!=("c"+test_ord))
425    {
426      "Warning: Block orderings are not supported; dp used for computation.";
427    }
428  }
429
430  if(external_ord=="" && find(singular_ord,"Dp")==1)
431  {
432    external_ord="W_LEX";
433    for(i=1;i<=nvars(basering);i++)
434    {
435      weightvec[i]=1;
436    }
437    test_ord="Dp("+string(nvars(basering))+"),";
438    if(singular_ord!=(test_ord+"C") && singular_ord!=(test_ord+"c"))
439    {
440      "Warning: Block orderings are not supported; Dp used for computation.";
441    }
442  }
443  if(external_ord=="" && find(singular_ord,"Dp")==3)
444  {
445    external_ord="W_LEX";
446    for(i=1;i<=nvars(basering);i++)
447    {
448      weightvec[i]=1;
449    }
450    test_ord="Dp("+string(nvars(basering))+"),";
451    if(singular_ord!=("C"+test_ord) && singular_ord!=("c"+test_ord))
452    {
453      "Warning: Block orderings are not supported; Dp used for computation.";
454    }
455  }
456
457  int pos;
458  string number_string;
459
460  if(external_ord=="" && find(singular_ord,"wp")==1)
461  {
462    external_ord="W_REV_LEX";
463    pos=3;
464    for(i=1;i<=nvars(basering);i++)
465    {
466      pos++;
467      number_string="";
468      while(singular_ord[pos]!=",")
469      {
470        number_string=number_string+singular_ord[pos];
471        pos++;
472      }
473      execute("weightvec[i]="+number_string);
474    }
475    test_ord="wp("+string(weightvec)+"),";
476    if(singular_ord!=(test_ord+"C") && singular_ord!=(test_ord+"c"))
477    {
478      "Warning: Block orderings are not supported; wp("+string(weightvec)+") used for computation.";
479    }
480  }
481  if(external_ord=="" && find(singular_ord,"wp")==3)
482  {
483    external_ord="W_REV_LEX";
484    pos=5;
485    for(i=1;i<=nvars(basering);i++)
486    {
487      pos++;
488      number_string="";
489      while(singular_ord[pos]!=",")
490      {
491        number_string=number_string+singular_ord[pos];
492        pos++;
493      }
494      execute("weightvec[i]="+number_string);
495    }
496    test_ord="wp("+string(weightvec)+"),";
497    if(singular_ord!=("C"+test_ord) && singular_ord!=("c"+test_ord))
498    {
499      "Warning: Block orderings are not supported; wp("+string(weightvec)+") used for computation.";
500    }
501  }
502
503  if(external_ord=="" && find(singular_ord,"Wp")==1)
504  {
505    external_ord="W_LEX";
506    pos=3;
507    for(i=1;i<=nvars(basering);i++)
508    {
509      pos++;
510      number_string="";
511      while(singular_ord[pos]!=",")
512      {
513        number_string=number_string+singular_ord[pos];
514        pos++;
515      }
516      execute("weightvec[i]="+number_string);
517    }
518    test_ord="Wp("+string(weightvec)+"),";
519    if(singular_ord!=(test_ord+"C") && singular_ord!=(test_ord+"c"))
520    {
521      "Warning: Block orderings are not supported; Wp("+string(weightvec)+") used for computation.";
522    }
523  }
524  if(external_ord=="" && find(singular_ord,"Wp")==3)
525  {
526    external_ord="W_LEX";
527    pos=5;
528    for(i=1;i<=nvars(basering);i++)
529    {
530      pos++;
531      number_string="";
532      while(singular_ord[pos]!=",")
533      {
534        number_string=number_string+singular_ord[pos];
535        pos++;
536      }
537      execute("weightvec[i]="+number_string);
538    }
539    test_ord="Wp("+string(weightvec)+"),";
540    if(singular_ord!=("C"+test_ord) && singular_ord!=("c"+test_ord))
541    {
542      "Warning: Block orderings are not supported; Wp("+string(weightvec)+") used for computation.";
543    }
544  }
545
546  if(external_ord=="")
547  {
548    "ERROR: The term ordering of the actual basering is not supported.";
549    return(I);
550  }
551
552  // check algorithm
553  if(alg=="ct" || alg=="pct")
554    // these algorithms will not cause an error in the external program;
555    // however, they do not compute the toric ideal of A, but of an
556    // extended matrix
557  {
558    "ERROR: The chosen algorithm is not suitable.";
559    return(I);
560  }
561
562  // create temporary file with that the external program is called
563
564  int dummy;
565  int process=system("pid");
566  string matrixfile="temp_MATRIX"+string(process);
567  link MATRIX=":w "+matrixfile;
568  open(MATRIX);
569
570  write(MATRIX,"MATRIX","columns:",ncols(A),"cost vector:");
571  for(j=1;j<=ncols(A);j++)
572  {
573    write(MATRIX,weightvec[j]);
574  }
575  write(MATRIX,"rows:",nrows(A),"matrix:");
576  for(i=1;i<=nrows(A);i++)
577  {
578    for(j=1;j<=ncols(A);j++)
579    {
580      write(MATRIX,A[i,j]);
581    }
582  }
583
584  // enter positive row space vector, if required by the algorithm
585  if((alg=="blr") || (alg=="hs"))
586  {
587    write(MATRIX,"positive row space vector:");
588    for(j=1;j<=ncols(A);j++)
589    {
590      write(MATRIX,prsv[j]);
591    }
592  }
593  close(MATRIX);
594
595  // call external program
596  dummy=system("sh","toric_ideal -alg "+alg+" "+matrixfile);
597
598  // read toric ideal from created file
599  link TORIC_IDEAL=":r "+matrixfile+".GB."+alg;
600  string toric_id=read(TORIC_IDEAL);
601
602  int generators;
603  pos=find(toric_id,"size");
604  pos=find(toric_id,":",pos);
605  pos++;
606
607  while(toric_id[pos]==" " || toric_id[pos]==newline)
608  {
609    pos++;
610  }
611  number_string="";
612  while(toric_id[pos]!=" " && toric_id[pos]!=newline)
613  {
614    number_string=number_string+toric_id[pos];
615    pos++;
616  }
617  execute("generators="+number_string+";");
618
619  intvec v;
620  poly head;
621  poly tail;
622
623  pos=find(toric_id,"basis");
624  pos=find(toric_id,":",pos);
625  pos++;
626
627  for(i=1;i<=generators;i++)
628  {
629    head=1;
630    tail=1;
631
632    for(j=1;j<=ncols(A);j++)
633    {
634      while(toric_id[pos]==" " || toric_id[pos]==newline)
635      {
636        pos++;
637      }
638      number_string="";
639      while(toric_id[pos]!=" " && toric_id[pos]!=newline)
640      {
641        number_string=number_string+toric_id[pos];
642        pos++;
643      }
644      execute("v[j]="+number_string+";");
645      if(v[j]<0)
646      {
647        tail=tail*var(j)^(-v[j]);
648      }
649      if(v[j]>0)
650      {
651        head=head*var(j)^v[j];
652      }
653    }
654    I[i]=head-tail;
655  }
656
657  // delete all created files
658  dummy=system("sh","rm -f "+matrixfile);
659  dummy=system("sh","rm -f "+matrixfile+".GB."+alg);
660
661  return(I);
662}
663///////////////////////////////////////////////////////////////////////////////
664
665proc toric_ideal
666"USAGE: toric_ideal(A,alg);      A intmat, alg string
667        toric_ideal(A,alg,prsv); A intmat, alg string, prsv intvec
668RETURN: ideal: standard basis of the toric ideal of A
669NOTE:   These procedures return the standard basis of the toric ideal of A
670        with respect to the term ordering in the current basering. Not all
671        term orderings are supported: The usual global term orderings may be
672        used, but no block orderings combining them.
673        One may call the procedure with several different algorithms: @*
674        - the algorithm of Conti/Traverso using elimination (ect), @*
675        - the algorithm of Pottier (pt),
676        - an algorithm of Bigatti/La Scala/Robbiano (blr),
677        - the algorithm of Hosten/Sturmfels (hs),
678        - the algorithm of DiBiase/Urbanke (du).
679       The argument `alg' should be the abbreviation for an algorithm as
680       above: ect, pt, blr, hs or du.
681
682      If `alg' is chosen to be `blr' or `hs', the algorithm needs a vector
683      with positive coefficients in the row space of A.
684      If no row of A contains only positive entries, one has to use the
685      second version of toric_ideal which takes such a vector as its third
686      argument.
687      For the mathematical background, see
688  @texinfo
689  @ref{Toric ideals and integer programming}.
690  @end texinfo
691EXAMPLE:  example toric_ideal; shows an example
692SEE ALSO: toric_std, toric_lib, intprog_lib, Toric ideals
693"
694{
695  if(size(#)==2)
696  {
697    return(toric_ideal_1(#[1],#[2]));
698  }
699  else
700  {
701    return(toric_ideal_2(#[1],#[2],#[3]));
702  }
703}
704
705
706
707example
708{
709"EXAMPLE"; echo=2;
710
711ring r=0,(x,y,z),dp;
712
713// call with two arguments
714intmat A[2][3]=1,1,0,0,1,1;
715A;
716
717ideal I=toric_ideal(A,"du");
718I;
719
720I=toric_ideal(A,"blr");
721I;
722
723// call with three arguments
724intvec prsv=1,2,1;
725I=toric_ideal(A,"blr",prsv);
726I;
727
728}
729///////////////////////////////////////////////////////////////////////////////
730
731proc toric_std(ideal I)
732"USAGE:   toric_std(I);      I ideal
733RETURN:   ideal: standard basis of I
734NOTE:     This procedure computes the standard basis of I using a specialized
735          Buchberger algorithm. The generating system by which I is given has
736          to consist of binomials of the form x^u-x^v. There is no real check
737          if I is toric. If I is generated by binomials of the above form,
738          but not toric, toric_std computes an ideal `between' I and its
739          saturation with respect to all variables.
740          For the mathematical background, see
741   @texinfo
742   @ref{Toric ideals and integer programming}.
743   @end texinfo
744EXAMPLE:  example toric_std; shows an example
745SEE ALSO: toric_ideal, toric_lib, intprog_lib, Toric ideals
746"
747{
748  ideal J;
749  // to be returned
750
751  // check suitability of actual term ordering
752  // the "module ordering" c or C is ignored
753  string singular_ord=ordstr(basering);
754  string test_ord;
755  string external_ord="";
756  int i,j;
757  intvec weightvec;
758
759  if(find(singular_ord,"lp")==1)
760  {
761    external_ord="W_LEX";
762    for(i=1;i<=nvars(basering);i++)
763    {
764      weightvec[i]=0;
765    }
766    test_ord="lp("+string(nvars(basering))+"),";
767    if(singular_ord!=(test_ord+"C") && singular_ord!=(test_ord+"c"))
768    {
769      "Warning: Block orderings are not supported; lp used for computation.";
770    }
771  }
772  if(external_ord=="" && find(singular_ord,"lp")==3)
773  {
774    external_ord="W_LEX";
775    for(i=1;i<=nvars(basering);i++)
776    {
777      weightvec[i]=0;
778    }
779    test_ord="lp("+string(nvars(basering))+"),";
780    if(singular_ord!=("C"+test_ord) && singular_ord!=("c"+test_ord))
781    {
782      "Warning: Block orderings are not supported; lp used for computation.";
783    }
784  }
785
786  if(external_ord=="" && find(singular_ord,"dp")==1)
787  {
788    external_ord="W_REV_LEX";
789    for(i=1;i<=nvars(basering);i++)
790    {
791      weightvec[i]=1;
792    }
793    test_ord="dp("+string(nvars(basering))+"),";
794    if(singular_ord!=(test_ord+"C") && singular_ord!=(test_ord+"c"))
795    {
796      "Warning: Block orderings are not supported; dp used for computation.";
797    }
798  }
799  if(external_ord=="" && find(singular_ord,"dp")==3)
800  {
801    external_ord="W_REV_LEX";
802    for(i=1;i<=nvars(basering);i++)
803    {
804      weightvec[i]=1;
805    }
806    test_ord="dp("+string(nvars(basering))+"),";
807    if(singular_ord!=("C"+test_ord) && singular_ord!=("c"+test_ord))
808    {
809      "Warning: Block orderings are not supported; dp used for computation.";
810    }
811  }
812
813  if(external_ord=="" && find(singular_ord,"Dp")==1)
814  {
815    external_ord="W_LEX";
816    for(i=1;i<=nvars(basering);i++)
817    {
818      weightvec[i]=1;
819    }
820    test_ord="Dp("+string(nvars(basering))+"),";
821    if(singular_ord!=(test_ord+"C") && singular_ord!=(test_ord+"c"))
822    {
823      "Warning: Block orderings are not supported; Dp used for computation.";
824    }
825  }
826  if(external_ord=="" && find(singular_ord,"Dp")==3)
827  {
828    external_ord="W_LEX";
829    for(i=1;i<=nvars(basering);i++)
830    {
831      weightvec[i]=1;
832    }
833    test_ord="Dp("+string(nvars(basering))+"),";
834    if(singular_ord!=("C"+test_ord) && singular_ord!=("c"+test_ord))
835    {
836      "Warning: Block orderings are not supported; Dp used for computation.";
837    }
838  }
839
840  int pos;
841  string number_string;
842
843  if(external_ord=="" && find(singular_ord,"wp")==1)
844  {
845    external_ord="W_REV_LEX";
846    pos=3;
847    for(i=1;i<=nvars(basering);i++)
848    {
849      pos++;
850      number_string="";
851      while(singular_ord[pos]!="," && singular_ord[pos]!=")")
852      {
853        number_string=number_string+singular_ord[pos];
854        pos++;
855      }
856      execute("weightvec["+string(i)+"]="+number_string+";");
857    }
858    test_ord="wp("+string(weightvec)+"),";
859    if(singular_ord!=(test_ord+"C") && singular_ord!=(test_ord+"c"))
860    {
861      "Warning: Block orderings are not supported; wp("+string(weightvec)+") used for computation.";
862    }
863  }
864
865  if(external_ord=="" && find(singular_ord,"wp")==3)
866  {
867    external_ord="W_REV_LEX";
868    pos=5;
869    for(i=1;i<=nvars(basering);i++)
870    {
871      pos++;
872      number_string="";
873      while(singular_ord[pos]!=",")
874      {
875        number_string=number_string+singular_ord[pos];
876        pos++;
877      }
878      execute("weightvec[i]="+number_string);
879    }
880    test_ord="wp("+string(weightvec)+"),";
881    if(singular_ord!=("C"+test_ord) && singular_ord!=("c"+test_ord))
882    {
883      "Warning: Block orderings are not supported; wp("+string(weightvec)+") used for computation.";
884    }
885  }
886
887  if(external_ord=="" && find(singular_ord,"Wp")==1)
888  {
889    external_ord="W_LEX";
890    pos=3;
891    for(i=1;i<=nvars(basering);i++)
892    {
893      pos++;
894      number_string="";
895      while(singular_ord[pos]!=",")
896      {
897        number_string=number_string+singular_ord[pos];
898        pos++;
899      }
900      execute("weightvec[i]="+number_string);
901    }
902    test_ord="Wp("+string(weightvec)+"),";
903    if(singular_ord!=(test_ord+"C") && singular_ord!=(test_ord+"c"))
904    {
905      "Warning: Block orderings are not supported; Wp("+string(weightvec)+") used for computation.";
906    }
907  }
908  if(external_ord=="" && find(singular_ord,"Wp")==3)
909  {
910    external_ord="W_LEX";
911    pos=5;
912    for(i=1;i<=nvars(basering);i++)
913    {
914      pos++;
915      number_string="";
916      while(singular_ord[pos]!=",")
917      {
918        number_string=number_string+singular_ord[pos];
919        pos++;
920      }
921      execute("weightvec[i]="+number_string);
922    }
923    test_ord="Wp("+string(weightvec)+"),";
924    if(singular_ord!=("C"+test_ord) && singular_ord!=("c"+test_ord))
925    {
926      "Warning: Block orderings are not supported; Wp("+string(weightvec)+") used for computation.";
927    }
928  }
929
930  if(external_ord=="")
931  {
932    "ERROR: The term ordering of the actual basering is not supported.";
933    return(I);
934  }
935
936  // create first temporary file with which the external program is called
937
938  int dummy;
939  int process=system("pid");
940  string groebnerfile="temp_GROEBNER"+string(process);
941  link GROEBNER=":w "+groebnerfile;
942  open(GROEBNER);
943
944  write(GROEBNER,"GROEBNER","computed with algorithm:","pt","term ordering:","elimination block",0,"weighted block",nvars(basering),external_ord);
945  // algorithm is totally unimportant, only required by the external program
946
947  for(i=1;i<=nvars(basering);i++)
948  {
949    write(GROEBNER,weightvec[i]);
950  }
951
952  write(GROEBNER,"size:",size(I),"Groebner basis:");
953  poly head;
954  poly tail;
955  poly rest;
956  intvec v;
957
958  for(j=1;j<=size(I);j++)
959  {
960    // test suitability of generator j
961    rest=I[j];
962    head=lead(rest);
963    rest=rest-head;
964    tail=lead(rest);
965    rest=rest-tail;
966
967    if(head==0 && tail==0 && rest!=0)
968    {
969      "ERROR: Generator "+string(j)+" of the input ideal is no binomial.";
970      close(GROEBNER);
971      dummy=system("sh","rm -f "+groebnerfile);
972      return(J);
973    }
974
975    if(leadcoef(tail)!=-leadcoef(head))
976      // generator is no difference of monomials (or a constant multiple)
977    {
978      "ERROR: Generator "+string(j)+" of the input ideal is no difference of monomials.";
979      close(GROEBNER);
980      dummy=system("sh","rm -f "+groebnerfile);
981      return(J);
982    }
983
984    if(gcd(head,tail)!=1)
985    {
986      "Warning: The monomials of generator "+string(j)+" of the input ideal are not relatively prime.";
987    }
988
989    // write vector representation of generator j into the file
990    v=leadexp(head)-leadexp(tail);
991    for(i=1;i<=nvars(basering);i++)
992    {
993      write(GROEBNER,v[i]);
994    }
995  }
996  close(GROEBNER);
997
998  // create second temporary file
999
1000  string newcostfile="temp_NEW_COST"+string(process);
1001  link NEW_COST=":w "+newcostfile;
1002  open(NEW_COST);
1003
1004  write(NEW_COST,"NEW_COST","variables:",nvars(basering),"cost vector:");
1005  for(i=1;i<=nvars(basering);i++)
1006  {
1007    write(NEW_COST,weightvec[i]);
1008  }
1009
1010  // call external program
1011  dummy=system("sh","change_cost "+groebnerfile+" "+newcostfile);
1012
1013  // read toric standard basis from created file
1014  link TORIC_IDEAL=":r "+newcostfile+".GB.pt";
1015  string toric_id=read(TORIC_IDEAL);
1016
1017  int generators;
1018  pos=find(toric_id,"size");
1019  pos=find(toric_id,":",pos);
1020  pos++;
1021
1022  while(toric_id[pos]==" " || toric_id[pos]==newline)
1023  {
1024    pos++;
1025  }
1026  number_string="";
1027  while(toric_id[pos]!=" " && toric_id[pos]!=newline)
1028  {
1029    number_string=number_string+toric_id[pos];
1030    pos++;
1031  }
1032  execute("generators="+number_string+";");
1033
1034  pos=find(toric_id,"basis");
1035  pos=find(toric_id,":",pos);
1036  pos++;
1037
1038  for(j=1;j<=generators;j++)
1039  {
1040    head=1;
1041    tail=1;
1042
1043    for(i=1;i<=nvars(basering);i++)
1044    {
1045      while(toric_id[pos]==" " || toric_id[pos]==newline)
1046      {
1047        pos++;
1048      }
1049      number_string="";
1050      while(toric_id[pos]!=" " && toric_id[pos]!=newline)
1051      {
1052        number_string=number_string+toric_id[pos];
1053        pos++;
1054      }
1055      execute("v[i]="+number_string+";");
1056      if(v[i]<0)
1057      {
1058        tail=tail*var(i)^(-v[i]);
1059      }
1060      if(v[i]>0)
1061      {
1062        head=head*var(i)^v[i];
1063      }
1064    }
1065    J[j]=head-tail;
1066  }
1067
1068  // delete all created files
1069  dummy=system("sh","rm -f "+groebnerfile);
1070  dummy=system("sh","rm -f "+groebnerfile+".GB.pt");
1071  dummy=system("sh","rm -f "+newcostfile);
1072
1073  return(J);
1074}
1075
1076example
1077{
1078"EXAMPLE"; echo=2;
1079
1080ring r=0,(x,y,z),wp(3,2,1);
1081
1082// call with toric ideal (of the matrix A=(1,1,1))
1083ideal I=x-y,x-z;
1084ideal J=toric_std(I);
1085J;
1086
1087// call with the same ideal, but badly chosen generators:
1088// 1) not only binomials
1089I=x-y,2x-y-z;
1090J=toric_std(I);
1091// 2) binomials whose monomials are not relatively prime
1092I=x-y,xy-yz,y-z;
1093J=toric_std(I);
1094J;
1095
1096// call with a non-toric ideal that seems to be toric
1097I=x-yz,xy-z;
1098J=toric_std(I);
1099J;
1100// comparison with real standard basis and saturation
1101ideal H=std(I);
1102H;
1103LIB "elim.lib";
1104sat(H,xyz);
1105}
1106///////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.