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