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>
![[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_contrib / TiGERS_0.9