# source:git/IntegerProgramming/toric.lib@6ba162

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