Annotation of OpenXM/src/Ti/gset.c, Revision 1.1
1.1 ! maekawa 1: /*
! 2: ** gset.c Birk Huber, 4/99
! 3: ** -- definitions and basic operations on gsets
! 4: **
! 5: **
! 6: ** TiGERS, Toric Groebner Basis Enumeration by Reverse Search
! 7: ** copyright (c) 1999 Birk Huber
! 8: **
! 9: */
! 10: #include<stdio.h>
! 11: #include<stdlib.h>
! 12: #include "utils.h"
! 13: #include "gset.h"
! 14: #include "matrices.h"
! 15: #include "Rsimp.h"
! 16: extern int NT;
! 17:
! 18: /* uncomment the following to cause use sainity checks (not fully tested)*/
! 19: /*
! 20: #define GTEST 0
! 21: #include "gtesting.h"
! 22: */
! 23:
! 24: /*
! 25: ** create an empty gset, initialize all its pointers to null.
! 26: */
! 27: gset gset_new(){
! 28: gset g=0;
! 29: if ((g=(gset)malloc(sizeof(struct gset_tag)))==0){
! 30: fprintf(stderr,"gset_new(): malloc failed\n");
! 31: return 0;
! 32: }
! 33: gset_first(g)=0;
! 34: gset_cache_vtx(g)=0;
! 35: gset_cache_edge(g)=0;
! 36: return g;
! 37: }
! 38:
! 39: /*
! 40: ** delete a gset and all binomials in it.
! 41: */
! 42: void gset_free(gset g){
! 43: binomial b;
! 44: while((b=gset_first(g))!=0){
! 45: gset_first(g)=binomial_next(b);
! 46: binomial_free(b);
! 47: }
! 48: free((void *)g);
! 49: }
! 50:
! 51:
! 52: /*
! 53: ** read a gset from an input stream
! 54: ** (return TRUE if gset read in )
! 55: */
! 56: int gset_read(FILE *is,gset g){
! 57: char c;
! 58: binomial b;
! 59:
! 60: eatwhite(is);
! 61: c=getc(is);
! 62: if (c!='{') return FALSE;
! 63: c=',';
! 64: while(c==','){
! 65: b=binomial_new();
! 66: binomial_read(is,b);
! 67: gset_insert(g,b);
! 68: eatwhite(is);
! 69: c=getc(is);
! 70: }
! 71: if (c!='}'){
! 72: fprintf(stderr,"gset_read(): gset should end with a '}'\n");
! 73: return FALSE;
! 74: }
! 75: ungetc(c,is);
! 76: return TRUE;
! 77: }
! 78:
! 79: /*
! 80: ** Display a gset
! 81: */
! 82: void gset_print(FILE *of,gset g){
! 83: binomial b;
! 84:
! 85: if (NT == 1) {
! 86: fprintf(of,"[");
! 87: for(b=g->bottom;b!=0;b=b->next){
! 88: binomial_print(of,b);
! 89: if (b->next!=0) fprintf(of," , ");
! 90: }
! 91: fprintf(of,"]");
! 92: }else{
! 93: fprintf(of,"{");
! 94: for(b=g->bottom;b!=0;b=b->next){
! 95: binomial_print(of,b);
! 96: if (b->next!=0) fprintf(of,", ");
! 97: }
! 98: fprintf(of,"}");
! 99: }
! 100: }
! 101:
! 102: /*
! 103: ** Display a gsets initial ideal
! 104: */
! 105: void gset_init_print(FILE *of,gset g){
! 106: binomial b;
! 107: fprintf(of,"{");
! 108: for(b=g->bottom;b!=0;b=b->next){
! 109: print_monomial(of,binomial_lead(b));
! 110: if (b->next!=0) fprintf(of,", ");
! 111: }
! 112: fprintf(of,"}");
! 113: }
! 114:
! 115: /*
! 116: ** Display a facet binomials
! 117: */
! 118: void gset_facet_print(FILE *of,gset g){
! 119: binomial b;
! 120: int tog=0;
! 121: fprintf(of,"{");
! 122: for(b=g->bottom;b!=0;b=b->next){
! 123: if (b->ff==FACET){
! 124: if (tog==1) fprintf(of,",");
! 125: else tog=1;
! 126: binomial_print(of,b);
! 127: }
! 128: }
! 129: fprintf(of,"}");
! 130: }
! 131:
! 132: /*
! 133: ** place binomial in gset in decreasing lexicographic order.
! 134: ** (using binomial_compair(,) to determine reletive order or
! 135: ** binomials)
! 136: */
! 137: void gset_insert(gset g, binomial b){
! 138: binomial t1,t2;
! 139: if (gset_first(g)==0 || binomial_compair(gset_first(g),b)==FALSE){
! 140: binomial_next(b)=gset_first(g);
! 141: gset_first(g)=b;
! 142: }
! 143: else{
! 144: t1=gset_first(g);
! 145: while(t1!=0){
! 146: t2=binomial_next(t1);
! 147: if (t2==0 || binomial_compair(t2,b)==FALSE){
! 148: binomial_next(t1)=b;
! 149: binomial_next(b)=t2;
! 150: break;
! 151: }
! 152: t1=t2;
! 153: }
! 154: }
! 155: }
! 156:
! 157: /*
! 158: ** return first facet binomial which is mis-marked
! 159: */
! 160: binomial gset_downedge(gset g){
! 161: binomial ptr=0;
! 162: for(ptr=gset_first(g);ptr!=0;ptr=binomial_next(ptr)){
! 163: if (binomial_ordered(ptr)==FALSE && gset_isfacet(g,ptr)==TRUE)break;
! 164: }
! 165: return ptr;
! 166: }
! 167:
! 168:
! 169: /*
! 170: ++void bmprint(FILE *of,binomial B, binomial M):
! 171: ++ simple output function for a generating set consisting of a binomial and
! 172: ++ a list of monomials --- usefull mainly for debugging. -- gset output
! 173: ++ routing may be modified in such a way as to use code which could double
! 174: ++ for this purpose.
! 175: ++
! 176: */
! 177: void bmprint(FILE *of,binomial B, binomial M){
! 178: fprintf(of,"{ ");
! 179: if (B!=0) binomial_print(of,B);
! 180: while(M!=0){
! 181: fprintf(of,",");
! 182: binomial_print(of,M);
! 183: M=binomial_next(M);
! 184: }
! 185: fprintf(of,"}\n");
! 186: }
! 187:
! 188: /*
! 189: **binomial bmleadreduce(binomial m1, binomial B, binomial mlist):
! 190: ** "lift" monomial m1 to a binomial using bmlist B,mlist.
! 191: */
! 192: binomial bmleadreduce(binomial m1, binomial B, binomial mlist){
! 193: binomial res,ptr;
! 194: /* case 1: the binomial is used in the decomp:*/
! 195: if (monomial_divides(binomial_lead(B),binomial_lead(m1))==TRUE){
! 196: res=binomial_new();
! 197: binomial_copy(B,res);
! 198: binomial_bumpto(res,m1);
! 199:
! 200: /* now reduce by lt(b) as many times as possible */
! 201: while(monomial_divides(binomial_lead(B),binomial_trail(res))==TRUE){
! 202: reducetrail(res,B);
! 203: }
! 204:
! 205: /* now find a monomial deviding whats left and reduce by it*/
! 206: for(ptr=mlist;ptr!=0;ptr=ptr->next){
! 207: if (monomial_divides(binomial_lead(ptr),binomial_trail(res))==TRUE){
! 208: reducetrail(res,ptr);
! 209: return res;
! 210: }
! 211: }
! 212: }
! 213: else { /* case 2 -- our monomial is a multiple of some monomial*/
! 214: for(ptr=mlist;ptr!=0;ptr=ptr->next){
! 215: if (monomial_divides(binomial_lead(ptr),binomial_lead(m1))==TRUE){
! 216: res=binomial_new();
! 217: binomial_copy(ptr,res);
! 218: binomial_set(res);
! 219: binomial_bumpto(res,m1);
! 220: return res;
! 221: }
! 222: }
! 223: }
! 224: fprintf(stderr, "Warning BMreduce does not get to zero\n");
! 225: }
! 226:
! 227:
! 228:
! 229: /*
! 230: ++ binomial find_divisor(binomial S,binomial L):
! 231: ++ search list starting at L for a monomial dividing S
! 232: ++ return this binomial or zero if no such binomial exists on list.
! 233: ++
! 234: */
! 235: binomial find_divisor(binomial S,binomial L){
! 236: binomial ptr;
! 237: for(ptr=L;ptr!=0;ptr=ptr->next){
! 238: if(ptr!=S&&monomial_divides(binomial_lead(ptr),binomial_lead(S))==TRUE){
! 239: break;
! 240: }
! 241: }
! 242: return ptr;
! 243: }
! 244:
! 245:
! 246: /*
! 247: ++ void remove_multiples(binomial S, binomial *L):
! 248: ++ Pass through list starting at *L, removing all multiples of *S
! 249: ++ encountered (except S if it is on the list).
! 250: ++
! 251: */
! 252: void remove_multiples(binomial S, binomial *L){
! 253: binomial tlist=*L,ptr;
! 254: *L=0;
! 255: while(tlist!=0){
! 256: ptr=tlist; tlist=tlist->next; ptr->next=0; /* pop tlist*/
! 257: if (S==ptr || monomial_divides(binomial_lead(S),binomial_lead(ptr))==FALSE) {
! 258: ptr->next=*L; /* push L */
! 259: *L=ptr;
! 260: }
! 261: else binomial_free(ptr);
! 262: }
! 263: }
! 264:
! 265: /*
! 266: ** void bmrgb(binomial B, binomial *ms):
! 267: ** transform a marked binomial and a list of monomials into the
! 268: ** reduced grobner basis for the ideal the generate.
! 269: **
! 270: */
! 271: void bmrgb(binomial B, binomial *ms){
! 272: int i,jk;
! 273: binomial DL,SL,tlist,S,tmp,ptr;
! 274: int szero;
! 275:
! 276: DL=*ms;
! 277: SL=0;
! 278:
! 279: while(DL!=0){
! 280: tmp=DL; DL=DL->next; /* pop DL */
! 281: tmp->next=SL; SL=tmp; /* push SL */
! 282:
! 283: /* calculate normal form of S pair with respect to {B}*/
! 284: S=monomial_spair(B,binomial_lead(tmp));
! 285:
! 286: /* check if anything on DL or on SL divides S*/
! 287: if (find_divisor(S,SL)!=0||find_divisor(S,DL)!=0)binomial_free(S);
! 288: else {
! 289: remove_multiples(S,&SL); /* remove any multiples of S from SL*/
! 290: remove_multiples(S,&DL); /* remove any multiples of S from DL*/
! 291: S->next=DL; DL=S; /* finally push S onto DL */
! 292: }
! 293: }
! 294:
! 295: remove_multiples(B,&SL); /* now that we have a grobner basis
! 296: we can remove multiples of lead(B)
! 297: from the monomial list */
! 298:
! 299: // replace list bottom pointer
! 300: *ms=SL;
! 301:
! 302: }
! 303:
! 304: /*
! 305: ** gset gset_flip(gset g1, binomial b):
! 306: ** Take a gset and one of its facet binomials and produce the new gset
! 307: ** corresponding to the grobner cone which shares the given facet binomial.
! 308: **
! 309: */
! 310: gset gset_flip(gset g1, binomial b){
! 311: binomial old=0,new=0,tmp,ptr,B,Bold;
! 312: gset gres=0, gtest=0;
! 313:
! 314: /* sainity check -- remove eventually*/
! 315: #ifdef GTEST
! 316: if (lp_isfacet(g1,b)!=FACET){
! 317: fprintf(stdout,"WARNING flip asked to flip non-binomial\n");
! 318: }
! 319: #endif
! 320:
! 321: /* copy all binomials of g's list ecept b to old and new lists*/
! 322: for(ptr=gset_first(g1);ptr!=0;ptr=ptr->next){
! 323: if (ptr!=b){
! 324: tmp=binomial_new();
! 325: binomial_copy(ptr,tmp);
! 326: monomial_set(tmp);
! 327: binomial_next(tmp)=new;
! 328: new=tmp;
! 329: tmp=binomial_new();
! 330: binomial_copy(ptr,tmp);
! 331: monomial_set(tmp);
! 332: binomial_next(tmp)=old;
! 333: old=tmp;
! 334: }
! 335: }
! 336: Bold=binomial_new();
! 337: binomial_copy(b,Bold);
! 338:
! 339: /* make fliped copy of b*/
! 340: B=binomial_new();
! 341: binomial_copy(b,B);
! 342: binomial_flip(B);
! 343:
! 344: /* calculate the reduced GB of the binomial b and the monomials in M */
! 345: bmrgb(B,&new);
! 346:
! 347: /* create new gset and insert flipped copy of b */
! 348: gres=gset_new();
! 349: gset_insert(gres,B);
! 350:
! 351: /* do lifting of elements in mlist and add them to new gset*/
! 352: while((tmp=new)!=0){
! 353: new=new->next;
! 354: ptr=bmleadreduce(tmp,Bold,old);
! 355: gset_insert(gres,ptr);
! 356: binomial_free(tmp);
! 357: }
! 358:
! 359: /* sainity check -- remove eventually*/
! 360: #ifdef GTEST
! 361: if (gset_is_gb(gres)==FALSE){
! 362: fprintf(stderr,"flipped failed to find grobner basis \n");
! 363: fprintf(stderr,"g1=");gset_print(stderr,g1);fprintf(stderr,"\n");
! 364: fprintf(stderr,"b=");binomial_print(stderr,b);fprintf(stderr,"\n");
! 365: exit(1);
! 366: }
! 367: #endif
! 368:
! 369: binomial_free(Bold);
! 370: while(old!=0){
! 371: ptr=old;
! 372: old=old->next;
! 373: binomial_free(ptr);
! 374: }
! 375:
! 376: gset_autoreduce(gres);
! 377:
! 378: /* sainity check -- remove eventually */
! 379: #ifdef GTEST
! 380: if ( gset_is_reduced(gres)==FALSE){
! 381: fprintf(stderr,"flipped failed to find reduced grobner basis \n");
! 382: fprintf(stderr,"g1=");gset_print(stderr,g1);fprintf(stderr,"\n");
! 383: fprintf(stderr,"b=");binomial_print(stderr,b);fprintf(stderr,"\n");
! 384: exit(1);
! 385: }
! 386: #endif
! 387:
! 388: return gres;
! 389: }
! 390:
! 391:
! 392: /*
! 393: ** void gset_autoreduce(gset g):
! 394: ** autoreduce g according to its marking
! 395: */
! 396: void gset_autoreduce(gset g){
! 397: binomial b1,b2;
! 398: int changed;
! 399:
! 400: /* First remove redundant elements from the grobner basis*/
! 401: while(find_divisor(gset_first(g),gset_first(g))!=0){
! 402: gset_first(g)=(b1=gset_first(g))->next;
! 403: binomial_free(b1);
! 404: }
! 405: b1=gset_first(g);
! 406: while(b1!=0 && binomial_next(b1)!=0){
! 407: if (find_divisor(binomial_next(b1),gset_first(g))!=0) {
! 408: b2=binomial_next(b1);
! 409: binomial_next(b1)=binomial_next(b2);
! 410: binomial_free(b2);
! 411: }
! 412: else b1=binomial_next(b1);
! 413: }
! 414:
! 415: /* now reduce trailing terms */
! 416: b1=gset_first(g);
! 417: while(b1!=0){
! 418: changed=1;
! 419: while(changed==1){
! 420: changed=0;
! 421: b2=gset_first(g);
! 422: while(b2!=0){
! 423: if (b2!=b1 &&
! 424: monomial_divides(binomial_lead(b2),binomial_trail(b1))==TRUE){
! 425: reducetrail(b1,b2);
! 426: changed=1;
! 427: }
! 428: b2=binomial_next(b2);
! 429: }
! 430: }
! 431: b1=binomial_next(b1);
! 432: }
! 433:
! 434: /* now reinsert terms to ensure list is in order */
! 435: b1=gset_first(g);
! 436: gset_first(g)=0;
! 437: while((b2=b1)!=0) {
! 438: b1=binomial_next(b1);
! 439: gset_insert(g,b2);
! 440: }
! 441:
! 442: }
! 443:
! 444: /*
! 445: ** void gset_rgb(gset g, int (*comp)(monomial,monomial))
! 446: ** Use buchberger's algorithm to transform g into groebner basis
! 447: ** wrt the term order implemented by comp
! 448: */
! 449: #define ordered(b) ((comp(binomial_lead(b),binomial_trail(b))>=0)?TRUE:FALSE)
! 450: int changed;
! 451: void gset_rgb(gset g, int (*comp)(monomial,monomial)){
! 452: void process_pair(binomial,binomial,binomial,binomial,binomial *,int (*comp)());
! 453: binomial ptr,b,DL=0,SL=0,NL=0;
! 454: int passes=0;
! 455:
! 456: do{
! 457: DL=0; SL=0; NL=0;
! 458: changed=FALSE;
! 459:
! 460: /* reorient g and put binomial list on DL*/
! 461: ptr=gset_first(g);
! 462: while((b=ptr)!=0){
! 463: ptr=binomial_next(ptr);
! 464: if (ordered(b)==FALSE) binomial_flip(b);
! 465: binomial_next(b)=DL;
! 466: DL=b;
! 467: }
! 468:
! 469: /* */
! 470: while(DL!=0){
! 471: NL=0;
! 472: b=DL; DL=DL->next; /* pop DL*/
! 473: b->next=SL; SL=b; /* push ptr onto SL*/
! 474: for(ptr=SL;ptr!=0;ptr=binomial_next(ptr)) process_pair(b,ptr,DL,SL,&NL,comp);
! 475: for(ptr=DL;ptr!=0;ptr=binomial_next(ptr)){
! 476: process_pair(b,ptr,DL,SL,&NL,comp);
! 477: if (binomial_next(ptr)==0){ binomial_next(ptr)=NL; NL=0;}
! 478: }
! 479: }
! 480:
! 481: /* reinsert elements in SL into gset */
! 482: ptr=SL;
! 483: gset_first(g)=0;
! 484: while(ptr!=0){
! 485: ptr=(b=ptr)->next;
! 486: gset_insert(g,b);
! 487: }
! 488:
! 489:
! 490: gset_autoreduce(g);
! 491:
! 492: if (passes>1){
! 493: if (changed==TRUE ){fprintf(stderr,"looping pass (%d)\n",passes);}
! 494: else {fprintf(stderr,"stopping after %d passes\n",passes);}
! 495: }
! 496: passes++;
! 497: }while(changed==TRUE);
! 498:
! 499: }
! 500:
! 501: /*
! 502: -- process_pair(b1,b2,DL,SL,NL,comp):
! 503: -- process spair defined by binomials b1 and b2.
! 504: -- taking spair wrt union of DL SL and *NL lists
! 505: -- if non-zero add result to NL list.
! 506: -- comp is a monomial compairison function used throuout.
! 507: */
! 508: void process_pair(binomial b1,binomial b2,binomial DL,binomial SL,binomial *NL,
! 509: int (*comp)(monomial,monomial)){
! 510: binomial S,ptr,find_divisor(binomial,binomial);
! 511: int i,tmp;
! 512:
! 513: /* test is lead terms are rel prime [buchberger's first criterion]*/
! 514: if (monomial_rel_prime(binomial_lead(b1),binomial_lead(b2))==TRUE)return;
! 515:
! 516: /* buchberger's second criterion ?*/
! 517:
! 518:
! 519: S=binomial_new();
! 520: binomial_copy(b1,S);
! 521: if (binomial_spair(S,b2)==TRUE){
! 522: binomial_free(S);
! 523: return;
! 524: }
! 525: if (ordered(S)==FALSE) binomial_flip(S);
! 526: while(((ptr=find_divisor(S,DL))!=0)||
! 527: ((ptr=find_divisor(S,SL))!=0)||
! 528: ((ptr=find_divisor(S,*NL))!=0)){
! 529: if (binomial_spair(S,ptr)==TRUE){
! 530: binomial_free(S);
! 531: return;
! 532: }
! 533: if (ordered(S)==FALSE) binomial_flip(S);
! 534: }
! 535: changed=TRUE;
! 536: binomial_next(S)=*NL;
! 537: *NL=S;
! 538: }
! 539:
! 540: /*
! 541: ** void gset_setfacets(gset g):
! 542: ** Use repeated calls to gset_isfacet to ensure that all facet flags
! 543: ** are set.
! 544: ** Also sets nfacets, nelts and deg counters (these contain garbage except
! 545: ** right after a call to set facets).
! 546: **
! 547: ** NOTE 4/20 the deg counter was just added and the code to load it is
! 548: ** cludged in here -- logically it should be part of the
! 549: ** binomial.c file -- for now it is here.
! 550: */
! 551: #define max(a,b) (((a)>(b))?(a):(b))
! 552: void gset_setfacets(gset g){
! 553: binomial b;
! 554: monomial ml,mt;
! 555: int i,d=0,dl,dt;
! 556: gset_nfacets(g)=0;
! 557: gset_nelts(g)=0;
! 558: for(b=gset_first(g);b!=0;b=binomial_next(b)){
! 559: if (gset_isfacet(g,b)==FACET) gset_nfacets(g)++;
! 560: gset_nelts(g)++;
! 561: ml=binomial_lead(b); dl=0;
! 562: mt=binomial_trail(b); dt=0;
! 563: for(i=0;i<ring_N;i++){ dl+=ml[i]; dt+=mt[i];}
! 564: d=max(d,max(dl,dt));
! 565: }
! 566: gset_deg(g)=d;
! 567: }
! 568:
! 569: /*
! 570: ** int gset_isfacet(gset g,binomial b):
! 571: ** check whether b is known to be a facet or not. if not call lp_isfacet()
! 572: ** to use a linear program to set b's facet flag.
! 573: **
! 574: */
! 575: int lptest=1; /* 1 test only LP, 2 test only flipability, 3 test both */
! 576: int gset_isflippable(gset g,binomial b);
! 577: int gset_isfacet(gset g,binomial b){
! 578: int rval=UNKNOWN;
! 579: rval=b->ff;
! 580: if (rval!=FACET && rval!=NONFACET){
! 581: switch(lptest){
! 582: case 1: rval=(b->ff=lp_isfacet(g,b));
! 583: break;
! 584: case 2: if (gset_isflippable(g,b)==FALSE) rval=(b->ff=NONFACET);
! 585: else rval=(b->ff=FACET);
! 586: break;
! 587: case 3: if (gset_isflippable(g,b)==FALSE) rval=(b->ff=NONFACET);
! 588: else rval=(b->ff=lp_isfacet(g,b));
! 589: break;
! 590: default: fprintf(stderr,"Warning lptest has invalid value in gset_isflippable()\n");
! 591: exit(4);
! 592: }
! 593: }
! 594: return rval;
! 595: }
! 596:
! 597: /*
! 598: **
! 599: */
! 600: int gset_isflippable(gset g,binomial b){
! 601: int ctg=0,ctn=1,val=TRUE;
! 602: binomial new=0,tmp,ptrg,ptrn,btmp;
! 603: /* copy all binomials of g's list ecept b to old and new lists*/
! 604: for(ptrg=gset_first(g);ptrg!=0;ptrg=ptrg->next){
! 605: ctg++;
! 606: if (ptrg!=b){
! 607: tmp=binomial_new();
! 608: binomial_copy(ptrg,tmp);
! 609: monomial_set(tmp);
! 610: binomial_next(tmp)=new;
! 611: new=tmp;
! 612: }
! 613: }
! 614: btmp=binomial_new();
! 615: binomial_copy(b,btmp);
! 616:
! 617:
! 618: /* calculate the reduced GB of the binomial b and the monomials in M */
! 619: bmrgb(btmp,&new);
! 620:
! 621: /* test if new and old have same leading terms */
! 622: for(ptrn=new ;ptrn!=0;ptrn=ptrn->next) ctn++;
! 623: if (ptrn!=ptrg) val=FALSE; /* initial ideals have diferent numbers of generators */
! 624: else {
! 625: for(ptrn=new;ptrn!=0;ptrn=ptrn->next){
! 626: for(ptrg=gset_first(g);ptrg!=0;ptrg=ptrg->next){
! 627: if (monomial_equal(binomial_lead(ptrn),binomial_lead(ptrg))==TRUE)break;
! 628: }
! 629: if (ptrg==0)break; /* no match found */
! 630: }
! 631: if (ptrn!=0)val=FALSE; /* not all terms had matches found */
! 632: }
! 633:
! 634: /* clean up */
! 635: binomial_free(btmp);
! 636: while(new!=0){
! 637: ptrn=new;
! 638: new=new->next;
! 639: binomial_free(ptrn);
! 640: }
! 641:
! 642: return val;
! 643: }
! 644:
! 645: /*
! 646: ** int lp_isfacet(gset g,binomial b):
! 647: ** set up and run linear program to set b's facet flag.
! 648: **
! 649: */
! 650: int lp_isfacet(gset g,binomial b){
! 651: binomial ptr;
! 652: int D=ring_N,M=0;
! 653: int idx=0,i;
! 654: double dtmp0,dtmp;
! 655:
! 656: /* need cheeper way to determine the size of g */
! 657: for(ptr=gset_first(g);ptr!=0;ptr=binomial_next(ptr)){
! 658: if (ptr->ff!=NONFACET)M++;
! 659: }
! 660:
! 661: /* make sure the lp structures have enough space for the program*/
! 662: LP_get_space(M,D+2*M);
! 663:
! 664: /* Set up linear program. */
! 665: ptr=gset_first(g);
! 666: idx=0;
! 667: while(ptr!=0){
! 668: if (ptr->ff!=NONFACET){
! 669: for(i=0;i<D;i++){ /* make sure supports of trail and lead disjoint? */
! 670: LP_A(idx,i)=binomial_lead(ptr)[i]-binomial_trail(ptr)[i];
! 671: }
! 672: for(i=0;i<M;i++){
! 673: LP_A(idx,D+i)=0;
! 674: LP_A(idx,D+M+i)=0;
! 675: }
! 676: if (ptr==b) {for(i=0;i<D;i++) LP_A(idx,i)*=-1;}
! 677: LP_A(idx,D+M+idx)=1;
! 678: LP_A(idx,D+idx)=-1;
! 679: LP_B[idx]=1;
! 680: idx++;
! 681: }
! 682: ptr=ptr->next;
! 683: }
! 684: for(i=0;i<D+M;i++){
! 685: LP_C[i]=0;
! 686: LP_X[i]=0;
! 687: LP_NonBasis[i]=i;
! 688: }
! 689: for(i=0;i<M;i++){
! 690: LP_Basis[i]=D+M+i;
! 691: LP_X[D+M+i]=1;
! 692: LP_C[D+M+i]=1;
! 693: }
! 694:
! 695: /* Run simplex algorithm */
! 696: Rsimp(LP_M,LP_N,LP_A,LP_B,LP_C,LP_X,LP_Basis,LP_NonBasis,
! 697: LP_R,LP_Q,LP_t1,LP_t2);
! 698:
! 699: /* use results to determine if the b is a facet*/
! 700: b->ff=FACET;
! 701: for(i=0;i<M;i++) if (LP_Basis[i]>=M+D) {b->ff=NONFACET;break;}
! 702:
! 703: return b->ff;
! 704: }
! 705:
! 706:
! 707:
! 708:
! 709:
! 710:
! 711:
! 712:
! 713:
! 714:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to gset.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / Ti