# source:git/IntegerProgramming/toric.lib@27e935

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