Annotation of OpenXM/src/Ti/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: 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