Annotation of OpenXM_contrib/TiGERS_0.9/toric.c, Revision 1.1.1.1
1.1 maekawa 1: /*
2: ** toric.c Birk Huber, 4/99
3: ** -- support for computing first grobner bases of toric ideals.
4: ** -- includes functions gset_rlex_rgb() and gset_compute_colon()
5: ** which once tested and working will move back to the gset.c file.
6: **
7: **
8: ** TiGERS, Toric Groebner Basis Enumeration by Reverse Search
9: ** copyright (c) 1999 Birk Huber
10: **
11: */
12: #include<stdio.h>
13: #include<stdlib.h>
14: #include "utils.h"
15: #include "gset.h"
16: #include "matrices.h"
17: #include "Rsimp.h"
18:
19:
20:
21:
22: /*
23: ** void gset_compute_colon(gset g,int lv):
24: ** Compute deg rev lex rgb (lv least) for colon of g with variable lv.
25: ** - set ring_lv=lv.
26: ** - convert g to rgb wrt deg_rev_lex as implemented in binomial.h.
27: ** - divide each binomial by highest power of var lv possible.
28: ** - autoreduce.
29: **
30: */
31: #define min(m,n) (((m)<(n))?(m):(n))
32: void gset_compute_colon(gset g,int lv){
33: binomial ptr;
34: int lold,jk;
35:
36: /* progress report */
37: fprintf(stderr,"taking colon of J with %c:\n",'a'+lv);
38:
39: /* set lex last variable to be lv */
40: lold=ring_lv;
41: ring_lv=lv;
42:
43: /* compute grobner basis wrt to this deg rev lex term order */
44: gset_rgb(g,monomial_rlexcomp);
45:
46: /* now devide each binomial by highest power of lvar which divides it */
47: for(ptr=gset_first(g); ptr!=0; ptr=binomial_next(ptr)){
48: jk=min(binomial_lead(ptr)[lv],binomial_trail(ptr)[lv]);
49: binomial_lead(ptr)[lv]-=jk;
50: binomial_trail(ptr)[lv]-=jk;
51: }
52:
53: /* result is a grobner basis, reduce it to make it an rgb */
54: gset_autoreduce(g);
55:
56: /* restore original value of lvar*/
57: ring_lv=lold;
58: }
59:
60: /*
61: ** gset gset_toric_ideal(int **M,int m, int n):
62: ** Given an mxn integer matrix M -- compute an rgb for the toric ideal I_M.
63: ** Uses repeated colon computations.
64: */
65: gset gset_toric_ideal(int **M,int m, int n){
66: int i,j,jk,crk;
67: binomial b;
68: int **V,**S,**new_imatrix();
69: gset g;
70:
71: /* set ring dimension*/
72: if (ring_N!=n){
73: fprintf(stderr,"ERROR gset_toric_ideal():");
74: fprintf(stderr," matrix dimensions incompatible with ring\n");
75: return 0;
76: }
77: /* set wieghting of A */
78: for(i=0;i<n;i++){
79: ring_weight[i]=0;
80: for(j=0;j<m;j++) ring_weight[i]+=M[j][i];
81: }
82:
83: /* reserve space for S and V */
84: S=new_imatrix(n,m);
85: V=new_imatrix(n,n);
86:
87: /* copy transpose of M to S*/
88: for(i=0;i<n;i++)for(j=0;j<m;j++) IMref(S,i,j)=IMref(M,j,i);
89:
90: /* compute hermite normal form of M,S,V and corank of M*/
91: crk=ihermite(S,V,m,n);
92:
93: /* create new gset and read off equations corresponding to rows of V*/
94: g=gset_new();
95: for(i=crk;i<n;i++){
96: b=binomial_new();
97: for(j=0;j<n;j++){
98: jk=IMref(V,i,j);
99: if (jk>0){
100: binomial_lead(b)[j]=jk; binomial_trail(b)[j]=0;}
101: else {binomial_trail(b)[j]=-1*jk; binomial_lead(b)[j]=0;}
102: }
103: gset_insert(g,b);
104: }
105:
106: /* compute successive colon ideals */
107: for(i=0;i<n;i++)gset_compute_colon(g,i);
108:
109: /*unset weighting */
110: for(i=0;i<n;i++) ring_weight[i]=1;
111:
112: /* free space used by matrices */
113: free_imatrix(M);
114: free_imatrix(S);
115: free_imatrix(V);
116:
117: /* make sure result is the lex gbasis */
118: gset_rgb(g,monomial_lexcomp);
119:
120: return g;
121: }
122:
123:
124:
125:
126:
127:
128:
129:
130:
131:
132:
133:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>