Annotation of OpenXM/src/Ti/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: extern int NT;
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: if (NT == 1) {
102: fprintf(of,"[");
103: for (i=0; i<ring_N; i++) {
104: fprintf(of,"%d",exps[i]);
105: if (i != ring_N-1) fprintf(of," , ");
106: }
107: fprintf(of,"]");
108: }else{
109: for(i=0;i<ring_N;i++){
110: if (exps[i]!=0){
111: if (jk!=0) fprintf(of,"*");
112: else jk=1;
113: ring_putvar(of,i);
114: if (exps[i]>1) fprintf(of,"^%d",exps[i]);
115: }
116: }
117: }
118: }
119:
120: /*
121: ** void get_monomial(FILE *,int *):
122: ** Read ascii representation of monomial as power product -- convert
123: ** to exponent vector stored in *exps.
124: **
125: */
126: void get_monomial(FILE *is,int *exps){
127: char c='*';
128: int v,e;
129:
130: for(v=0;v<ring_N;v++)exps[v]=0;
131: while(c=='*'){
132: eatwhite(is);
133: v=ring_getvar(is);
134: eatwhite(is);
135: c=fgetc(is);
136: if (c=='^'){
137: fscanf(is,"%d",&e);
138: c=fgetc(is);
139: }
140: else {
141: e=1;
142: }
143: exps[v]=e;
144: }
145: ungetc(c,is);
146: }
147:
148: /*
149: ** int monomial_divides(monomial m1, monomial m2):
150: ** return TRUE if m1 divides m2
151: **
152: */
153: int monomial_divides(monomial m1, monomial m2){
154: int i;
155: for(i=0;i<ring_N;i++)if (m1[i]>m2[i]) return FALSE;
156: return TRUE;
157: }
158:
159: /*
160: ** int monomial_rel_prime(monomial m1, monomial m2):
161: ** return TRUE if m1 is reletively prime to m2.
162: **
163: */
164: int monomial_rel_prime(monomial m1, monomial m2){
165: int i;
166: for(i=0;i<ring_N;i++) if (m1[i]!=0 && m2[i]!=0) return FALSE;
167: return TRUE;
168: }
169:
170: /*
171: ** int monomial_equal(monomial m1, monomial m2):
172: ** return TRUE if monomials m1 and m2 have same exponents
173: **
174: */
175: int monomial_equal(monomial m1, monomial m2){
176: int i;
177: for(i=0;i<ring_N;i++) if (m1[i]-m2[i]!=0) return FALSE;
178: return TRUE;
179: }
180:
181: /*
182: **int monomial_lcm(monomial m1, monomial m2, monomial m3):
183: ** copy least common multiple of m1 and m2 into m3.
184: **
185: */
186: int monomial_lcm(monomial m1, monomial m2, monomial m3){
187: int i;
188: for(i=0;i<ring_N;i++){
189: if (m1[i]>=m2[i]) m3[i]=m1[i];
190: else m3[i]=m2[i];
191: }
192: }
193:
194: /*
195: **
196: **
197: */
198: int monomial_stddegree(monomial m){
199: int i,d=0;
200: for(i=0;i<ring_N;i++) d+=m[i];
201: return d;
202: }
203:
204:
205: /*
206: ** int monomial_lexcomp(monomial m1, monomial m2):
207: ** return >0 if m1 is lexicographically greater then m2
208: **
209: */
210: int monomial_lexcomp(monomial m1, monomial m2){
211: int i,jk;
212: for(i=0;i<ring_N;i++) if ((jk=m1[i]-m2[i])!=0) break;
213: return jk;
214: }
215:
216: /*
217: ** int monomial_rlexcomp(monomial m1, monomial m2):
218: ** return >0 if m1 is degree reverse lexicographically greater than m2
219: ** with variables taken in order a>b...., except that variable ring_lv is
220: ** taken last (if it is >=0).
221: **
222: */
223: int monomial_rlexcomp(monomial m1, monomial m2){
224: int i,jk,d=0;
225: /* compute diference of total degrees */
226: for(i=0;i<ring_N;i++) d+=ring_weight[i]*(m1[i]-m2[i]);
227: if (d!=0) return d;
228: /* test ring_lv entree first */
229: if (ring_lv>=0 && ((jk=m1[ring_lv]-m2[ring_lv])!=0)) return -1*jk;
230: /* now test entrees in reverse order */
231: for(i=ring_N-1;i>=0;i--) if ((jk=m1[i]-m2[i])!=0) return -1*jk;
232: return 0;
233: }
234:
235: /*
236: **----------binomial defs------------------------------------------------
237: **
238: **
239: */
240:
241: /*
242: ** binomial binomial_new(): allocate storage for a new binomial--
243: ** initialize flags and set exponents to zero;
244: */
245: binomial binomial_new(){
246: binomial m;
247: int i,*E=0,*E2=0;
248:
249: if ((m=(binomial)malloc(sizeof(struct bin_tag)))==0){
250: fprintf(stderr,"binomial_new(): first malloc failed\n");
251: return 0;
252: }
253: if ((m->E=(int *)malloc(2*ring_N*sizeof(int)))==0){
254: fprintf(stderr,"binomial_new(): seccond malloc failed\n");
255: free((void*)m);
256: return 0;
257: }
258: m->exps1=m->E;
259: m->exps2=m->E+ring_N;
260: m->ff=UNKNOWN;
261: m->bf=BINOMIAL;
262: m->next=0;
263: for(i=0;i<ring_N;i++) { m->exps1[i]=(m->exps2[i]=0);}
264: return m;
265: }
266:
267:
268: /*
269: ** void binomial_free(binomial m):
270: ** reclaim space allocated by binomial_new()
271: **
272: */
273: void binomial_free(binomial m){
274: if (m->E!=0)free((void *)(m->E));
275: free((void *)m);
276: }
277:
278: /*
279: ** void binomial_read(FILE *is, binomial b):
280: ** read data into binomial
281: **
282: */
283: void binomial_read(FILE *is, binomial b){
284: char c, d='-';
285: eatwhite(is);
286: c=getc(is);
287: /* read in status of facet flag if present*/
288: if (c=='#') b->ff=FACET;
289: else if (c=='!') b->ff=NONFACET;
290: else if (c=='?') b->ff=UNKNOWN;
291: else ungetc(c,is);
292:
293: /* read in leading and trailing monomials if present*/
294: eatwhite(is);
295: if ((c=getc(is))=='-') d='+';
296: else ungetc(c,is);
297: get_monomial(is,b->exps1);
298: eatwhite(is);
299: if ((c=getc(is))==d){
300: get_monomial(is,b->exps2);
301: b->bf=BINOMIAL;
302: }else{
303: b->bf=MONOMIAL;
304: ungetc(c,is);
305: }
306: if (d=='+'){
307: if (b->bf==MONOMIAL){
308: fprintf(stderr,"WARNING: binomial_read() monomial with negative coef\n");
309: }
310: binomial_flip(b);
311: }
312: }
313:
314:
315: /*
316: ** void binomial_print(FILE *,binomial):
317: ** print binomial
318: **
319: */
320: void binomial_print(FILE *of, binomial b){
321: if (NT == 1) {
322: if (b==0) fprintf(of,"[ ]");
323: else {
324: /* if (b->ff==FACET) fprintf(of,"# ");
325: else if (b->ff==NONFACET) fprintf(of,"! "); */
326: /*else if (b->ff==UNKNOWN) fprintf(of,"? ");*/
327: fprintf(of,"[");
328: print_monomial(of,b->exps1);
329: if (b->bf==BINOMIAL){
330: fprintf(of," , ");
331: print_monomial(of,b->exps2);
332: fprintf(of,"]");
333: }
334: }
335: }else{
336: if (b==0) fprintf(of,"(NULL)");
337: else {
338: if (b->ff==FACET) fprintf(of,"# ");
339: else if (b->ff==NONFACET) fprintf(of,"! ");
340: /*else if (b->ff==UNKNOWN) fprintf(of,"? ");*/
341: print_monomial(of,b->exps1);
342: if (b->bf==BINOMIAL){
343: fprintf(of,"-");
344: print_monomial(of,b->exps2);
345: }
346: }
347: }
348: }
349:
350: /*
351: ** void binomial_copy(binomial src,binomial dest):
352: ** copy source binomials exponents and flags to those of dest binomial
353: **
354: */
355: void binomial_copy(binomial src,binomial dest){
356: int i;
357:
358: for(i=0;i<ring_N;i++){
359: binomial_lead(dest)[i]=binomial_lead(src)[i];
360: binomial_trail(dest)[i]=binomial_trail(src)[i];
361: }
362: binomial_next(dest)=0;
363: dest->bf=src->bf;
364: }
365:
366: /*
367: ** void binomial_flip(binomial b):
368: ** interchange leading and trailing terms of binomial
369: **
370: */
371: void binomial_flip(binomial b){
372: int *tmp;
373: tmp=binomial_lead(b);
374: binomial_lead(b)=binomial_trail(b);
375: binomial_trail(b)=tmp;
376: }
377:
378: /*
379: ** binomial monomial_spair(binomial b,monomial m):
380: ** form spair between marked binomial and marked monomial
381: ** and reduce result by binomial as many times as possible
382: **
383: */
384: binomial monomial_spair(binomial b,monomial m){
385: binomial S = 0;
386: int i,tmp;
387:
388: S=binomial_new();
389:
390: /* form monomial spair i.e. trailing term times factor required to */
391: /* bump lead term up to monomial -- also set trail terms to zero */
392: for(i=0;i<ring_N;i++){
393: (binomial_lead(S))[i]=(binomial_trail(b))[i];
394: (binomial_trail(S))[i]=0;
395: if ((tmp=(m[i]-(binomial_lead(b))[i]))>0)((binomial_lead(S))[i])+=tmp;
396: }
397:
398: /*reduce S by b as many times as possible*/
399: while(monomial_divides(binomial_lead(b),binomial_lead(S))==TRUE){
400: for(i=0;i<ring_N;i++){
401: tmp=(binomial_lead(S))[i]-(binomial_lead(b))[i];
402: (binomial_lead(S))[i]=(binomial_trail(b))[i]+tmp;
403: }
404:
405: }
406: return S;
407: }
408:
409: /*
410: ** int binomial_spair(binomial b1, binomial b2):
411: ** replace b1, by spair of b1 and b2 -- the spair is taken with respect
412: ** to the markings of b1 and b2, but it's marking is not meaningfull.
413: ** returns TRUE if result is equivalent to zero monomial.
414: **
415: */
416: int binomial_spair(binomial b1, binomial b2){
417: int i,jk,jk2=TRUE;
418: for(i=0;i<ring_N;i++){
419: jk=binomial_lead(b1)[i]-binomial_lead(b2)[i];
420: if (binomial_lead(b1)[i]!=0 && binomial_lead(b2)!=0) jk2=FALSE;
421: binomial_lead(b1)[i]=binomial_trail(b2)[i];
422: if (jk>0) binomial_lead(b1)[i]+=jk;
423: else binomial_trail(b1)[i]-=jk;
424: }
425: jk=monomial_equal(binomial_lead(b1),binomial_trail(b1));
426: if (jk==TRUE || jk2==TRUE) return TRUE;
427: else return FALSE;
428: }
429:
430:
431: /*
432: **void binomial_bumpto(binomial b1, binomial b2):
433: ** assuming that the leading term of b1 divides that of b2
434: ** multiply both sides of b1 so that the leading terms of b1
435: ** and b2 are equal
436: **
437: */
438: void binomial_bumpto(binomial b1, binomial b2){
439: int i,tmp;
440: for(i=0;i<ring_N;i++){
441: tmp=(binomial_lead(b2))[i]-(binomial_lead(b1))[i];
442: (binomial_lead(b1))[i]+=tmp;
443: (binomial_trail(b1))[i]+=tmp;
444: }
445: }
446:
447: /*
448: ** void reducetrail(binomial b1, binomial b2):
449: ** assuming that lead(b2) divides trail(b1) change b1 into binomial
450: ** resulting from bumping b2 so that lead(b2)=trail(b1) and adding
451: ** result to b1.
452: **
453: */
454: void reducetrail(binomial b1, binomial b2){
455: int i;
456: for( i=0;i<ring_N;i++){
457: binomial_trail(b1)[i]+=binomial_trail(b2)[i]-binomial_lead(b2)[i];
458: }
459: }
460:
461: /*
462: ** int binomial_compair(binomial b1,binomial b2):
463: ** return TRUE if the lead term of b1 is lexicographically greater than
464: ** that of b2 or if they tie return true if the trail term b1 is greater
465: ** or equal to that of b2 otherwise return false
466: */
467: int binomial_compair(binomial b1,binomial b2){
468: int tmp;
469:
470: tmp=monomial_lexcomp(binomial_lead(b1),binomial_lead(b2));
471: if (tmp>0) return TRUE;
472: else if (tmp<0) return FALSE;
473: else {
474: tmp=monomial_lexcomp(binomial_trail(b1),binomial_trail(b2));
475: if (tmp>=0) return TRUE;
476: else return FALSE;
477: }
478: }
479:
480:
481:
482:
483:
484:
485:
486:
487:
488:
489:
490:
491:
492:
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 / src / Ti