source: git/Singular/LIB/toric.lib @ 0ae4ce

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