Annotation of OpenXM_contrib/TiGERS_0.9/binomial.c, Revision 1.1.1.1
1.1 maekawa 1: /*
2: ** binomial.c Birk Huber, 4/99
3: ** -- define binomal type and main operations on binomials
4: **
5: **
6: ** TiGERS, Toric Groebner Basis Enumeration by Reverse Search
7: ** copyright (c) 1999 Birk Huber
8: **
9: */
10: #define BINOMIAL_H 1
11: #include<stdio.h>
12: #include<stdlib.h>
13: #include "utils.h"
14: #include "binomial.h"
15:
16:
17: /*
18: **-------------ring defs-----------------------------------------------------
19: **
20: ** Monomials in the ambiant polynomial ring will be represented by their
21: ** exponent vectors. Before any computations can be done using monomials the
22: ** dimension of this ring (and the exponent vectors) must be set and stored
23: ** in the global variable "ring_N". The mapping between indeces and variable
24: ** names for i.o. purposes is currently done by just using the first ring_N
25: ** letters of the alphabet.
26: **
27: ** Data:
28: ** int ring_N: global variable holding dimension of ambiant ring
29: ** int ring_lv: global variable used by deg-rev-lex term order, i.e. specify
30: ** which variable to make cheepest(others ordered alphabetically)
31: **
32: ** Functions (or macros):
33: ** int ring_set(int n): Set ring_N to n if ring_N has not yet been set.
34: ** int ring_read(File *infile): Read in a description of the ring from
35: ** infile and set ring values accordingly.
36: ** ring_getvar(FILE *): convert next variable name on FILE to variable index
37: ** ring_putvar(File *of,int v): convert index to variable name.
38: **
39: ** Note: for now a ring description is just given the number of variables.
40: ** the names are just the first ring_N letters of the alphabet:
41: ** i.e. ** name <====> index
42: ** 'a' <====> 1
43: ** 'b' <====> 1
44: ** eventually this scheme may be adjusted to allow names to be set
45: ** at run time.
46: */
47: int ring_N=0;
48: int ring_lv=-1;
49: int *ring_weight=0;
50:
51: int ring_set(int n){
52: int i;
53: if (ring_N!=0) {
54: fprintf(stderr,"ring_set(): can't reset dimension\n");
55: return FALSE;
56: }
57: if ((ring_weight=(int *)(malloc(n*sizeof(int))))==0){
58: fprintf(stderr,"ring_set(): malloc failed for weight vector\n");
59: return FALSE;
60: }
61: for(i=0;i<n;i++) ring_weight[i]=1;
62: return ring_N=n;
63: }
64:
65: int ring_read(FILE *infile){
66: int n;
67: eatwhite(infile);
68: if (fscanf(infile,"%d",&n)!=1){
69: fprintf(stderr,"ERROR: read_ring() unable to parse ring description\n");
70: return FALSE;
71: }
72: ring_set(n);
73: return TRUE;
74: }
75:
76: int ring_getvar(FILE *ifile){
77: char c;
78: c=fgetc(ifile);
79: if ('a' > c || c >'a'+ring_N-1){
80: fprintf(stderr,"ring_getvar(): invalid variable read\n");
81: ungetc(c,ifile);
82: return -1;
83: }
84: return c-'a';
85: }
86:
87: /*
88: **------------monomial defs-----------------------------------------------
89: ** monomials are represented by their exponent vectors:
90: ** C integer vectors indexed from 0 to ring_dim-1
91: **
92: */
93:
94: /*
95: ** void print_monomial(FILE *,int *):
96: ** print exponent vector as monomial
97: **
98: */
99: void print_monomial(FILE *of, int *exps){
100: int i,jk=0;
101: for(i=0;i<ring_N;i++){
102: if (exps[i]!=0){
103: if (jk!=0) fprintf(of,"*");
104: else jk=1;
105: ring_putvar(of,i);
106: if (exps[i]>1) fprintf(of,"^%d",exps[i]);
107: }
108: }
109: }
110:
111: /*
112: ** void get_monomial(FILE *,int *):
113: ** Read ascii representation of monomial as power product -- convert
114: ** to exponent vector stored in *exps.
115: **
116: */
117: void get_monomial(FILE *is,int *exps){
118: char c='*';
119: int v,e;
120:
121: for(v=0;v<ring_N;v++)exps[v]=0;
122: while(c=='*'){
123: eatwhite(is);
124: v=ring_getvar(is);
125: eatwhite(is);
126: c=fgetc(is);
127: if (c=='^'){
128: fscanf(is,"%d",&e);
129: c=fgetc(is);
130: }
131: else {
132: e=1;
133: }
134: exps[v]=e;
135: }
136: ungetc(c,is);
137: }
138:
139: /*
140: ** int monomial_divides(monomial m1, monomial m2):
141: ** return TRUE if m1 divides m2
142: **
143: */
144: int monomial_divides(monomial m1, monomial m2){
145: int i;
146: for(i=0;i<ring_N;i++)if (m1[i]>m2[i]) return FALSE;
147: return TRUE;
148: }
149:
150: /*
151: ** int monomial_rel_prime(monomial m1, monomial m2):
152: ** return TRUE if m1 is reletively prime to m2.
153: **
154: */
155: int monomial_rel_prime(monomial m1, monomial m2){
156: int i;
157: for(i=0;i<ring_N;i++) if (m1[i]!=0 && m2[i]!=0) return FALSE;
158: return TRUE;
159: }
160:
161: /*
162: ** int monomial_equal(monomial m1, monomial m2):
163: ** return TRUE if monomials m1 and m2 have same exponents
164: **
165: */
166: int monomial_equal(monomial m1, monomial m2){
167: int i;
168: for(i=0;i<ring_N;i++) if (m1[i]-m2[i]!=0) return FALSE;
169: return TRUE;
170: }
171:
172: /*
173: **int monomial_lcm(monomial m1, monomial m2, monomial m3):
174: ** copy least common multiple of m1 and m2 into m3.
175: **
176: */
177: int monomial_lcm(monomial m1, monomial m2, monomial m3){
178: int i;
179: for(i=0;i<ring_N;i++){
180: if (m1[i]>=m2[i]) m3[i]=m1[i];
181: else m3[i]=m2[i];
182: }
183: }
184:
185: /*
186: **
187: **
188: */
189: int monomial_stddegree(monomial m){
190: int i,d=0;
191: for(i=0;i<ring_N;i++) d+=m[i];
192: return d;
193: }
194:
195:
196: /*
197: ** int monomial_lexcomp(monomial m1, monomial m2):
198: ** return >0 if m1 is lexicographically greater then m2
199: **
200: */
201: int monomial_lexcomp(monomial m1, monomial m2){
202: int i,jk;
203: for(i=0;i<ring_N;i++) if ((jk=m1[i]-m2[i])!=0) break;
204: return jk;
205: }
206:
207: /*
208: ** int monomial_rlexcomp(monomial m1, monomial m2):
209: ** return >0 if m1 is degree reverse lexicographically greater than m2
210: ** with variables taken in order a>b...., except that variable ring_lv is
211: ** taken last (if it is >=0).
212: **
213: */
214: int monomial_rlexcomp(monomial m1, monomial m2){
215: int i,jk,d=0;
216: /* compute diference of total degrees */
217: for(i=0;i<ring_N;i++) d+=ring_weight[i]*(m1[i]-m2[i]);
218: if (d!=0) return d;
219: /* test ring_lv entree first */
220: if (ring_lv>=0 && ((jk=m1[ring_lv]-m2[ring_lv])!=0)) return -1*jk;
221: /* now test entrees in reverse order */
222: for(i=ring_N-1;i>=0;i--) if ((jk=m1[i]-m2[i])!=0) return -1*jk;
223: return 0;
224: }
225:
226: /*
227: **----------binomial defs------------------------------------------------
228: **
229: **
230: */
231:
232: /*
233: ** binomial binomial_new(): allocate storage for a new binomial--
234: ** initialize flags and set exponents to zero;
235: */
236: binomial binomial_new(){
237: binomial m;
238: int i,*E=0,*E2=0;
239:
240: if ((m=(binomial)malloc(sizeof(struct bin_tag)))==0){
241: fprintf(stderr,"binomial_new(): first malloc failed\n");
242: return 0;
243: }
244: if ((m->E=(int *)malloc(2*ring_N*sizeof(int)))==0){
245: fprintf(stderr,"binomial_new(): seccond malloc failed\n");
246: free((void*)m);
247: return 0;
248: }
249: m->exps1=m->E;
250: m->exps2=m->E+ring_N;
251: m->ff=UNKNOWN;
252: m->bf=BINOMIAL;
253: m->next=0;
254: for(i=0;i<ring_N;i++) { m->exps1[i]=(m->exps2[i]=0);}
255: return m;
256: }
257:
258:
259: /*
260: ** void binomial_free(binomial m):
261: ** reclaim space allocated by binomial_new()
262: **
263: */
264: void binomial_free(binomial m){
265: if (m->E!=0)free((void *)(m->E));
266: free((void *)m);
267: }
268:
269: /*
270: ** void binomial_read(FILE *is, binomial b):
271: ** read data into binomial
272: **
273: */
274: void binomial_read(FILE *is, binomial b){
275: char c, d='-';
276: eatwhite(is);
277: c=getc(is);
278: /* read in status of facet flag if present*/
279: if (c=='#') b->ff=FACET;
280: else if (c=='!') b->ff=NONFACET;
281: else if (c=='?') b->ff=UNKNOWN;
282: else ungetc(c,is);
283:
284: /* read in leading and trailing monomials if present*/
285: eatwhite(is);
286: if ((c=getc(is))=='-') d='+';
287: else ungetc(c,is);
288: get_monomial(is,b->exps1);
289: eatwhite(is);
290: if ((c=getc(is))==d){
291: get_monomial(is,b->exps2);
292: b->bf=BINOMIAL;
293: }else{
294: b->bf=MONOMIAL;
295: ungetc(c,is);
296: }
297: if (d=='+'){
298: if (b->bf==MONOMIAL){
299: fprintf(stderr,"WARNING: binomial_read() monomial with negative coef\n");
300: }
301: binomial_flip(b);
302: }
303: }
304:
305:
306: /*
307: ** void binomial_print(FILE *,binomial):
308: ** print binomial
309: **
310: */
311: void binomial_print(FILE *of, binomial b){
312: if (b==0) fprintf(of,"(NULL)");
313: else {
314: if (b->ff==FACET) fprintf(of,"# ");
315: else if (b->ff==NONFACET) fprintf(of,"! ");
316: /*else if (b->ff==UNKNOWN) fprintf(of,"? ");*/
317: print_monomial(of,b->exps1);
318: if (b->bf==BINOMIAL){
319: fprintf(of,"-");
320: print_monomial(of,b->exps2);
321: }
322: }
323: }
324:
325: /*
326: ** void binomial_copy(binomial src,binomial dest):
327: ** copy source binomials exponents and flags to those of dest binomial
328: **
329: */
330: void binomial_copy(binomial src,binomial dest){
331: int i;
332:
333: for(i=0;i<ring_N;i++){
334: binomial_lead(dest)[i]=binomial_lead(src)[i];
335: binomial_trail(dest)[i]=binomial_trail(src)[i];
336: }
337: binomial_next(dest)=0;
338: dest->bf=src->bf;
339: }
340:
341: /*
342: ** void binomial_flip(binomial b):
343: ** interchange leading and trailing terms of binomial
344: **
345: */
346: void binomial_flip(binomial b){
347: int *tmp;
348: tmp=binomial_lead(b);
349: binomial_lead(b)=binomial_trail(b);
350: binomial_trail(b)=tmp;
351: }
352:
353: /*
354: ** binomial monomial_spair(binomial b,monomial m):
355: ** form spair between marked binomial and marked monomial
356: ** and reduce result by binomial as many times as possible
357: **
358: */
359: binomial monomial_spair(binomial b,monomial m){
360: binomial S = 0;
361: int i,tmp;
362:
363: S=binomial_new();
364:
365: /* form monomial spair i.e. trailing term times factor required to */
366: /* bump lead term up to monomial -- also set trail terms to zero */
367: for(i=0;i<ring_N;i++){
368: (binomial_lead(S))[i]=(binomial_trail(b))[i];
369: (binomial_trail(S))[i]=0;
370: if ((tmp=(m[i]-(binomial_lead(b))[i]))>0)((binomial_lead(S))[i])+=tmp;
371: }
372:
373: /*reduce S by b as many times as possible*/
374: while(monomial_divides(binomial_lead(b),binomial_lead(S))==TRUE){
375: for(i=0;i<ring_N;i++){
376: tmp=(binomial_lead(S))[i]-(binomial_lead(b))[i];
377: (binomial_lead(S))[i]=(binomial_trail(b))[i]+tmp;
378: }
379:
380: }
381: return S;
382: }
383:
384: /*
385: ** int binomial_spair(binomial b1, binomial b2):
386: ** replace b1, by spair of b1 and b2 -- the spair is taken with respect
387: ** to the markings of b1 and b2, but it's marking is not meaningfull.
388: ** returns TRUE if result is equivalent to zero monomial.
389: **
390: */
391: int binomial_spair(binomial b1, binomial b2){
392: int i,jk,jk2=TRUE;
393: for(i=0;i<ring_N;i++){
394: jk=binomial_lead(b1)[i]-binomial_lead(b2)[i];
395: if (binomial_lead(b1)[i]!=0 && binomial_lead(b2)!=0) jk2=FALSE;
396: binomial_lead(b1)[i]=binomial_trail(b2)[i];
397: if (jk>0) binomial_lead(b1)[i]+=jk;
398: else binomial_trail(b1)[i]-=jk;
399: }
400: jk=monomial_equal(binomial_lead(b1),binomial_trail(b1));
401: if (jk==TRUE || jk2==TRUE) return TRUE;
402: else return FALSE;
403: }
404:
405:
406: /*
407: **void binomial_bumpto(binomial b1, binomial b2):
408: ** assuming that the leading term of b1 divides that of b2
409: ** multiply both sides of b1 so that the leading terms of b1
410: ** and b2 are equal
411: **
412: */
413: void binomial_bumpto(binomial b1, binomial b2){
414: int i,tmp;
415: for(i=0;i<ring_N;i++){
416: tmp=(binomial_lead(b2))[i]-(binomial_lead(b1))[i];
417: (binomial_lead(b1))[i]+=tmp;
418: (binomial_trail(b1))[i]+=tmp;
419: }
420: }
421:
422: /*
423: ** void reducetrail(binomial b1, binomial b2):
424: ** assuming that lead(b2) divides trail(b1) change b1 into binomial
425: ** resulting from bumping b2 so that lead(b2)=trail(b1) and adding
426: ** result to b1.
427: **
428: */
429: void reducetrail(binomial b1, binomial b2){
430: int i;
431: for( i=0;i<ring_N;i++){
432: binomial_trail(b1)[i]+=binomial_trail(b2)[i]-binomial_lead(b2)[i];
433: }
434: }
435:
436: /*
437: ** int binomial_compair(binomial b1,binomial b2):
438: ** return TRUE if the lead term of b1 is lexicographically greater than
439: ** that of b2 or if they tie return true if the trail term b1 is greater
440: ** or equal to that of b2 otherwise return false
441: */
442: int binomial_compair(binomial b1,binomial b2){
443: int tmp;
444:
445: tmp=monomial_lexcomp(binomial_lead(b1),binomial_lead(b2));
446: if (tmp>0) return TRUE;
447: else if (tmp<0) return FALSE;
448: else {
449: tmp=monomial_lexcomp(binomial_trail(b1),binomial_trail(b2));
450: if (tmp>=0) return TRUE;
451: else return FALSE;
452: }
453: }
454:
455:
456:
457:
458:
459:
460:
461:
462:
463:
464:
465:
466:
467:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to binomial.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / TiGERS_0.9