[BACK]Return to gset.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / TiGERS_0.9

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

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>