Annotation of OpenXM_contrib/TiGERS_0.9/binomial.c, Revision 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