source: git/Singular/LIB/toric.lib @ 380a17b

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