1 #include "Global.h"
2 #include "Rat.h"
3
4 #define MAX_BAD_EQ (POLY_Dmax>5) /* previously 6; needed for nef !? */
5 #define SHOW_NEW_CEq 0 /* (POLY_Dmax>12)
6 tracks polytope analysis */
7 #ifndef CEQ_Nmax
8 #define CEQ_Nmax EQUA_Nmax
9 #endif
10
11 typedef struct {int ne; Equation e[CEQ_Nmax];} CEqList;
12
13 /* ====================================================================== */
14 /* ========== ========== */
15 /* ========== I N C I D E N C E S (as bit patterns) ========== */
16 /* ========== ========== */
17 /* ====================================================================== */
18
19 #if (VERT_Nmax > LONG_LONG_Nbits)
INCI_AND(INCI x,INCI y)20 INCI INCI_AND(INCI x, INCI y){
21 INCI z; int i; for (i=0;i<I_NUI;i++) z.ui[i]=(x.ui[i])&(y.ui[i]); return z;}
INCI_OR(INCI x,INCI y)22 INCI INCI_OR(INCI x, INCI y){
23 INCI z; int i; for (i=0;i<I_NUI;i++) z.ui[i]=(x.ui[i])|(y.ui[i]); return z;}
INCI_XOR(INCI x,INCI y)24 INCI INCI_XOR(INCI x, INCI y){
25 INCI z; int i; for (i=0;i<I_NUI;i++) z.ui[i]=(x.ui[i])^(y.ui[i]); return z;}
INCI_EQ(INCI x,INCI y)26 int INCI_EQ(INCI x, INCI y){
27 int i; for (i=0;i<I_NUI;i++) if ((x.ui[i])!=(y.ui[i])) return 0; return 1;}
INCI_EQ_0(INCI x)28 int INCI_EQ_0(INCI x){
29 int i; for (i=0;i<I_NUI;i++) if (x.ui[i]) return 0; return 1;}
INCI_LE(INCI x,INCI y)30 int INCI_LE(INCI x, INCI y){
31 int i;
32 /* for (i=0;i<I_NUI;i++) if ((x.ui[i]|y.ui[i])!=(y.ui[i])) return 0; */
33 unsigned int *X=x.ui, *Y=y.ui;
34 for (i=0;i<I_NUI;i++) if ((X[i]|Y[i])!=(Y[i])) return 0;
35 return 1;}
INCI_0()36 INCI INCI_0(){
37 INCI z; int i; for (i=0;i<I_NUI;i++) z.ui[i]=0; return z;}
INCI_1()38 INCI INCI_1(){
39 INCI z; int i; z.ui[0]=1; for (i=1;i<I_NUI;i++) z.ui[i]=0; return z;}
INCI_PN(INCI x,Long y)40 INCI INCI_PN(INCI x, Long y){
41 INCI z; int i;
42 z.ui[0]=(x.ui[0]<<1)|(!y);
43 for (i=1;i<I_NUI;i++) z.ui[i]=(x.ui[i]<<1)|(x.ui[i-1]>>(INT_Nbits-1));
44 return z;}
INCI_D2(INCI x)45 INCI INCI_D2(INCI x){
46 INCI z; int i;
47 for (i=0;i<I_NUI-1;i++) z.ui[i]=(x.ui[i]>>1)|(x.ui[i+1]<<(INT_Nbits-1));
48 z.ui[I_NUI-1]=x.ui[I_NUI-1]>>1;
49 return z;}
INCI_lex_GT(INCI * x,INCI * y)50 int INCI_lex_GT(INCI *x, INCI *y){
51 int i=I_NUI; while(i--) if(x->ui[i]>y->ui[i]) return 1;
52 else if(x->ui[i]<y->ui[i]) return 0; return 0; }
INCI_LmR(INCI * x,INCI * y)53 int INCI_LmR(INCI *x, INCI *y){ puts("Implement INCI_LmR"); exit(0); }
54 #else
INCI_lex_GT(INCI * x,INCI * y)55 int INCI_lex_GT(INCI *x, INCI *y){ return (*x > *y) ? 1 : 0 ; }
INCI_LmR(INCI * x,INCI * y)56 int INCI_LmR(INCI *x, INCI *y){ return (*x>*y) ? 1 : (*x<*y) ? -1 : 0; }
57 /* int Lead_Vert(INCI x){int i=0; while(!(x%2)) {i++; x/=2;} return i;} */
58 #endif
59
INCI_abs(INCI X)60 int INCI_abs(INCI X){
61 int abs=0; while(!INCI_EQ_0(X)) {abs+=INCI_M2(X); X=INCI_D2(X);} return abs;
62 }
63
Eq_To_INCI(Equation * _Eq,PolyPointList * _P,VertexNumList * _V)64 INCI Eq_To_INCI(Equation *_Eq, PolyPointList *_P, VertexNumList *_V){
65 int j; INCI X=INCI_0();
66 for (j=0;j<_V->nv;j++) X=INCI_PN(X,Eval_Eq_on_V(_Eq,_P->x[_V->v[j]],_P->n));
67 return X;
68 }
69
Print_INCI(INCI X)70 int Print_INCI(INCI X) {
71 int i=0;
72 while(!INCI_EQ_0(X)) {printf("%d", (int) INCI_M2(X)); X=INCI_D2(X); i++;}
73 return i; /* number of printed bits */
74 }
75
Print_FaceInfo(int M,FaceInfo * _I)76 void Print_FaceInfo(int M,FaceInfo *_I){
77 int i,j, k, l;
78 M--;
79 printf("Incidences as binary numbers [F-vector=(%d",_I->nf[0]);
80 for(i=1;i<=M;i++) printf(" %d",_I->nf[i]);
81 puts(")]:");
82 puts("v[d][i]: sum_j Incidence(i'th dim-d-face, j-th vertex) x 2^j");
83 for(i=0;i<=M;i++) {
84 printf("v[%d]: ",i);
85 for(j=0;j<_I->nf[i];j++){
86 k=Print_INCI(_I->v[i][j]);
87 for (l=k;l<_I->nf[0];l++) printf("0");
88 printf(" ");}
89 puts(""); }
90 puts("f[d][i]: sum_j Incidence(i'th dim-d-face, j-th facet) x 2^j");
91 for(i=0;i<=M;i++) {
92 printf("f[%d]: ",i);
93 for(j=0;j<_I->nf[i];j++){
94 k=Print_INCI(_I->f[i][j]);
95 for (l=k;l<_I->nf[M];l++) printf("0");
96 printf(" ");}
97 puts(""); }
98 }
99
100 void Make_CD2Faces(PolyPointList *_P, VertexNumList *_V, EqList *_E,
101 FaceInfo *_I);
102
Make_Incidence(PolyPointList * _P,VertexNumList * _V,EqList * _E,FaceInfo * _I)103 void Make_Incidence(PolyPointList *_P, VertexNumList *_V, EqList *_E,
104 FaceInfo *_I)
105 /* The incidence relations for faces are stored on the structure FaceInfo:
106 *
107 * int nf[d] == #faces(dim.=d) == #dual faces[dim.=n-d-1]
108 * INCI v[d][i] :: vertices on i-th dim=d face
109 * INCI f[d][i] :: dual vertices on face dual to i-th dim=n-d-1 face
110 * Long nip[d][i] :: #(IPs of i-th dim=d face)
111 * Long dip[d][i] :: #(IPs of i-th dim=n-d-1 face on dual)
112 *
113 * .v: compute vertices on facets; make intersections `&' of faces with facets
114 * while keeping only maximal intersections;
115 * .f: take analogous union; if same intersection again make union `|' of f's
116 */
117 { int i, j, M=_P->n-1, d=M, D;
118 assert(_E->ne<=VERT_Nmax);
119 _I->nf[M]=_E->ne;
120 for(i=0;i<_E->ne;i++)_I->v[M][i]=Eq_To_INCI(&_E->e[i],_P,_V); /* init .v*/
121 assert(i>0);
122 _I->f[M][--i]=INCI_1();
123 while(i--) _I->f[M][i]=INCI_PN(_I->f[M][i+1],1); /* init .f */
124 while((D=d--))
125 { int *n=&_I->nf[d]; *n=0; /* #(d-faces) */
126 for(i=0; i<_I->nf[D]; i++) for(j=0; j<_I->nf[M]; j++)
127 { int k; INCI x=INCI_AND(_I->v[D][i],_I->v[M][j]); /* x=candidate */
128 INCI u=INCI_OR(_I->f[D][i],_I->f[M][j]);
129 if( (!INCI_EQ(x,_I->v[D][i]))&&(INCI_abs(x)>d) )/* x!=vD & |x|>d */
130 { for(k=0;k<*n;k++)
131 { INCI *y=&_I->v[d][k],*v=&_I->f[d][k]; /* (*y)==v[d][k] */
132 if(INCI_LE(*y,x))
133 { if(INCI_EQ(x,*y))
134 { *v=INCI_OR(*v,u); break; /* x=y :: .f&=... */
135 }
136 else {int l=k; *y=x; _I->f[d][k]=u;/* x>y :: y=x;f= */
137 while((++l) < (*n))
138 if(INCI_LE(_I->v[d][l],x)) {(*n)--;
139 if(l<(*n)){_I->v[d][l]=_I->v[d][*n]; _I->f[d][l]=_I->f[d][*n];}}
140 else assert(!INCI_LE(x,_I->v[d][l]));
141
142 /* for(k++;k<*n;k++) if(!INCI_LE(x,_I->v[d][k])&&!INCI_LE(_I->v[d][k],x))
143 {Print_PPL(_P,"FACE_Info trouble");assert(0);} */
144 break;}
145 }
146 else if(INCI_LE(x,*y)) break; /* x<y :: break */
147 }
148 if(*n==k) /* non-comparable => new face */
149 { assert(k<FACE_Nmax); _I->v[d][k]=x; _I->f[d][k]=u; (*n)++;
150 }
151 }
152 }
153 } d=_P->n; M=0; for(i=0;i<d;i++) M += _I->nf[i] * (1-2*(i%2));
154 if(M!=2*(d%2)){for(i=0;i<d;i++)printf("%3d",_I->nf[i]);puts("=F-vector");
155 Print_PPL(_P,"PPL for incidence error");Print_FaceInfo(_P->n,_I);
156 printf("d=%d euler=%d\n",d,M); assert(M==2*(d%2));}
157 }
158
159 /* ====================================================================== */
160 /* ========== ========== */
161 /* ========== G E N E R A L P U R P O S E R O U T I N E S ========== */
162 /* ========== ========== */
163 /* ====================================================================== */
164
swap(int * i,int * j)165 void swap(int* i,int* j) {register int k; k=*i; *i=*j; *j=k;}
166
diff(const void * a,const void * b)167 int diff(const void *a, const void *b){return *((int *) a) - *((int *) b);}
168
Sort_VL(VertexNumList * _V)169 void Sort_VL(VertexNumList *_V){qsort(_V->v, _V->nv, sizeof(int), &diff);}
170
Make_VEPM(PolyPointList * _P,VertexNumList * _V,EqList * _E,PairMat PM)171 void Make_VEPM(PolyPointList *_P, VertexNumList *_V, EqList *_E,
172 PairMat PM){
173 int i, j;
174 for (i=0;i<_E->ne;i++) for (j=0;j<_V->nv;j++)
175 PM[i][j]=Eval_Eq_on_V(&_E->e[i],_P->x[_V->v[j]],_P->n);
176 }
177
Transpose_PM(PairMat PM,PairMat DPM,int nv,int ne)178 int Transpose_PM(PairMat PM, PairMat DPM, int nv, int ne){
179 int i, j;
180 if((nv>EQUA_Nmax)||(ne>VERT_Nmax)) return 0;
181 for (i=0;i<ne;i++) for (j=0;j<nv;j++) DPM[j][i]=PM[i][j];
182 return 1;
183 }
184
VNL_to_DEL(PolyPointList * _P,VertexNumList * _V,EqList * _DE)185 int VNL_to_DEL(PolyPointList *_P, VertexNumList *_V, EqList *_DE){
186 int i, j;
187 if (_V->nv>EQUA_Nmax) return 0;
188 _DE->ne=_V->nv;
189 for (i=0;i<_V->nv;i++){
190 for (j=0;j<_P->n;j++) _DE->e[i].a[j]=_P->x[_V->v[i]][j];
191 _DE->e[i].c=1; }
192 return 1;
193 }
194
195 #define LLong_EEV (1) /* 1 @ [4662 4 20 333 422 1554 2329] */
196 #define TEST_EEV (0) /* compare Long to LLong EEV */
197
EEV_To_Equation(Equation * _E1,Equation * _E2,Long * _V,int n)198 Equation EEV_To_Equation(Equation *_E1, Equation *_E2, Long *_V, int n){
199 /* Calculate the equation spanned by _V and the intersection of _E1, _E2 */
200 int i; Long l, m, g; Equation Eq;
201 l=Eval_Eq_on_V(_E2,_V,n);
202 m=Eval_Eq_on_V(_E1,_V,n);
203 g=NNgcd(l,m); assert(g); l/=g; m/=g;
204 #if ((!(LLong_EEV))||(TEST_EEV)) /* Long version */
205 for(i=0;i<n;i++) Eq.a[i]=l*_E1->a[i]-m*_E2->a[i];
206 { int gcd=Eq.c=l*_E1->c-m*_E2->c;
207 for(i=0;i<n;i++) gcd=NNgcd(gcd,Eq.a[i]); assert(gcd);
208 if (gcd!=1) { for(i=0;i<n;i++) Eq.a[i]/=gcd; Eq.c/=gcd;}}
209 #endif
210 #if ((LLong_EEV)||(TEST_EEV)) /* LLong version */
211 { LLong A[POLY_Dmax], C, G; for(i=0;i<n;i++)
212 A[i]=((LLong) l)*((LLong)_E1->a[i])-((LLong)m)*((LLong)_E2->a[i]);
213 G=C=((LLong) l)*((LLong)_E1->c)-((LLong)m)*((LLong)_E2->c);
214 for(i=0;i<n;i++) G=LNNgcd(G,A[i]); assert(G);
215 if(G!=1) {C/=G; for(i=0;i<n;i++) A[i]/=G;}
216 #if (TEST_EEV) /* Compare */
217 { int e=(Eq.c!=C); for(i=0;i<n;i++) if(Eq.a[i]!=A[i]) e=1;
218 if(e) { printf("Error in EEV: l=%d m=%d g=%d\n",l,m,g);
219 for(i=0;i<n;i++)printf("%d ",_E1->a[i]);printf(" %d = E1\n",_E1->c);
220 for(i=0;i<n;i++)printf("%d ",_E2->a[i]);printf(" %d = E2\n",_E2->c);
221 for(i=0;i<n;i++)printf("%d ",Eq.a[i]);printf(" %d = Eq\n",Eq.c);
222 for(i=0;i<n;i++)printf("%d ",A[i]);printf(" %d = LL_Eq\n",C);
223 exit(0); }
224 }
225 #else
226 Eq.c=C; for(i=0;i<n;i++) Eq.a[i]=A[i];
227 #endif
228 }
229 #endif
230 return Eq;
231 }
232
Eval_Eq_on_V(Equation * E,Long * V,int i)233 Long Eval_Eq_on_V(Equation *E, Long *V, int i){
234 Long p=E->c; while(i--) p+=V[i]*E->a[i];
235 return p;
236 }
237
DualBraP1(Long * X,Long * Y,int n)238 Long DualBraP1(Long *X, Long *Y, int n){
239 Long p=1; while(n--) p+=X[n] * Y[n];
240 return (Long) p;
241 }
242
CompareEq(Equation * X,Equation * Y,int n)243 Long CompareEq(Equation *X, Equation *Y, int n) { /* return "X-Y"; */
244 Long d; while(n--) if((d=X->a[n]-Y->a[n])) return d;
245 return X->c-Y->c;
246 }
247
Vec_Greater_Than(Long * X,Long * Y,int i)248 int Vec_Greater_Than(Long *X, Long *Y, int i){ /* return 1 iff `X > Y' */
249 while(i--) {if(X[i]>Y[i]) return 1; if(X[i]<Y[i]) return 0;}
250 puts("Identical points in Vec_Greater_Than !!"); exit(0); return 0;
251 }
252
Vec_is_zero(Long * X,int i)253 int Vec_is_zero(Long *X, int i){ /* return 1 iff `X == 0' */
254 while(i--) if(X[i]) return 0;
255 return 1;
256 }
257
Vec_Equal(Long * X,Long * Y,int i)258 int Vec_Equal(Long *X, Long *Y, int i){ /* return 1 iff `X == Y' */
259 while(i--) if(X[i]!=Y[i]) return 0;
260 return 1;
261 }
262
Swap_Vecs(Long * X,Long * Y,int i)263 void Swap_Vecs(Long *X, Long *Y, int i){ /* `X <--> Y' */
264 Long z; while(i--) {z=X[i]; X[i]=Y[i]; Y[i]=z;}
265 }
266
Ref_Equations(EqList * _F)267 int Ref_Equations(EqList *_F){
268 int i; for (i=0;i<_F->ne;i++) if (_F->e[i].c!=1) return 0; return 1;
269 }
270
Span_Check(EqList * _F,EqList * H,int * n)271 int Span_Check(EqList *_F, EqList *H, int *n){
272 int i, j;
273 for(i=0;i<H->ne;i++){
274 for(j=0;j<_F->ne;j++) if (!CompareEq(&(H->e[i]),&(_F->e[j]),*n)) break;
275 if(j==_F->ne) return 0; } /* didn't find H[i] among the F's */
276 return 1;
277 }
278
IsGoodCEq(Equation * _E,PolyPointList * _P,VertexNumList * _V)279 int IsGoodCEq(Equation *_E, PolyPointList *_P, VertexNumList *_V){
280 int i=_V->nv;
281 Long s;
282 while(!(s=Eval_Eq_on_V(_E, _P->x[_V->v[--i]], _P->n)));
283 if(s < 0) { int j=_P->n; while(j--) _E->a[j]=-_E->a[j]; _E->c=-_E->c; }
284 while(i) if(Eval_Eq_on_V(_E, _P->x[_V->v[--i]], _P->n) < 0) return 0;
285 return 1;
286 }
287
Search_New_Vertex(Equation * _E,PolyPointList * _P)288 int Search_New_Vertex(Equation *_E, PolyPointList *_P){
289 int i, v=0;
290 Long *X=_P->x[0], x=Eval_Eq_on_V(_E,X,(_P->n));
291 for(i=1;i<_P->np;i++) {
292 Long *Y=_P->x[i], y=Eval_Eq_on_V(_E,Y,(_P->n));
293 if(y>x) continue;
294 if(y==x) if(Vec_Greater_Than(X,Y,_P->n)) continue;
295 v=i; X=Y; x=y; }
296 return v;
297 }
298
Sort_PPL(PolyPointList * _P,VertexNumList * _V)299 void Sort_PPL(PolyPointList *_P, VertexNumList *_V){
300 /* Vertices first, IP last */
301 int i,j;
302 for (i=0;i<_V->nv;i++){
303 Swap_Vecs(_P->x[i], _P->x[_V->v[i]], _P->n);
304 for (j=i+1;j<_V->nv;j++) if (_V->v[j]==i) _V->v[j]=_V->v[i];
305 _V->v[i]=i;}
306 for (i=_V->nv;i<_P->np-1;i++) if(Vec_is_zero(_P->x[i], _P->n)){
307 Swap_Vecs(_P->x[i], _P->x[_P->np-1], _P->n);
308 return;}
309 }
310
EL_to_PPL(EqList * _E,PolyPointList * _P,int * n)311 int EL_to_PPL(EqList *_E, PolyPointList *_P, int *n){
312 int i,j;
313 for (i=0;i<_E->ne;i++){
314 if (_E->e[i].c!=1) {_P->np=0; return 0;}
315 for (j=0;j<*n;j++) _P->x[i][j]=_E->e[i].a[j];}
316 _P->np=_E->ne;
317 _P->n=*n;
318 return 1;
319 }
320
321 /* ====================================================================== */
322 /* ========== ========== */
323 /* ========== S T A R T -- S I M P L E X ========== */
324 /* ========== ========== */
325 /* ====================================================================== */
326
327 /* return 0 <=> max.dim., E.ne==P.n+1, made Simplex of Vertices of P; *
328 * return (P.n-E.ne) == codim. > 0 <=> E.ne defining equations on E; */
329
330 /* #define NEW_START_SIMPLEX (1) 1 @ [16644 1 38 439 2315 5548 8303] */
331
332 #define VERT_WITH_MAX_DISTANCE (0) /* 0 @ [1845 2 15 97 247 610 874] */
333 #define LONG_EQ_FIRST (0) /* 0 @ [3425 2 7 137 429 1141 1709] */
334 #define TEST_GLZ_EQ (0) /* trace StartSimplex EQs */
335
VZ_to_Base(Long * V,int * d,Long M[POLY_Dmax][POLY_Dmax])336 Long VZ_to_Base(Long *V,int *d,Long M[POLY_Dmax][POLY_Dmax]) /* 0 iff V=0 */
337 { int p[POLY_Dmax], i, j, J=0; Long g=0, W[POLY_Dmax], *G[POLY_Dmax];
338 for(i=0;i<*d;i++) if(V[i]) {W[J]=V[i]; G[J]=M[i]; p[J++]=i;}
339 else for(j=0;j<*d;j++) M[i][j]=(i==j);
340 if(J) if(p[0]) { G[0]=M[0]; for(j=0;j<*d;j++) M[p[0]][j]=(j==0);}
341 if(J>1) g=W_to_GLZ(W,&J,G); else if(J){g=*W; M[0][0]=0; M[0][p[0]]=1;}
342 if(J>1)
343 { for(i=0;i<J;i++) { int I=J;
344 for(j=*d-1;j>=0;j--) G[i][j] = (V[j]) ? G[i][--I] : 0; assert(I==0);}
345 } return g;
346 }
347
OrthBase_red_by_V(Long * V,int * d,Long A[][POLY_Dmax],int * r,Long B[][POLY_Dmax])348 int OrthBase_red_by_V(Long *V, int *d, Long A[][POLY_Dmax], int *r,
349 Long B[][POLY_Dmax])
350 { int i, j, k; Long W[POLY_Dmax], G[POLY_Dmax][POLY_Dmax];
351 for(i=0;i<*r;i++) {int j; W[i]=0; for(j=0;j<*d;j++) W[i]+=A[i][j]*V[j];}
352 assert( VZ_to_Base(W,r,G) );
353 for(i=0;i<*r-1;i++) for(k=0;k<*d;k++)
354 { B[i][k]=0; for(j=0;j<*r;j++) B[i][k]+=G[i+1][j]*A[j][k];
355 }
356 #if (TEST_GLZ_EQ)
357 printf("A -> B ... V = "); for(k=0;k<*d;k++) printf(" %5d",V[k]);
358 printf(" W=");for(k=0;k<*r;k++)printf(" %5d",W[k]);puts(""); {int
359 a,b; for(a=0;a<*r-1;a++){for(b=0;b<*d;b++)printf(" %5d",A[a][b]);
360 printf(" =A B= ");for(b=0;b<*d;b++)printf(" %5d",B[a][b]);
361 puts("");}for(b=0;b<*d;b++)printf(" %5d",A[a][b]);printf(" =A\n");}
362 #endif
363 return (*r)--;
364 }
365
New_Start_Vertex(Long * V0,Long * Ea,PolyPointList * P,int * v)366 int New_Start_Vertex(Long *V0,Long *Ea, PolyPointList *P,int *v) /* P.x[v] */
367 { Equation E; int i, n=0, p=0; Long d, dn=0, dp=0, *Xn=P->x[0], *Xp=Xn;
368 for(i=0;i<P->n;i++) E.a[i]=Ea[i];
369 E.c=0; E.c=-Eval_Eq_on_V(&E,V0,P->n);
370 d=Eval_Eq_on_V(&E,P->x[0],P->n); if(d>0) dp=d; if(d<0) dn=d;
371 for(i=1;i<P->np;i++)
372 { d=Eval_Eq_on_V(&E,P->x[i],P->n); if(d==0) continue;
373 if(d==dp) if(Vec_Greater_Than(P->x[i],Xp,P->n)) Xp=P->x[p=i];
374 if(d>dp) {dp=d; Xp=P->x[p=i];}
375 if(d==dn) if(Vec_Greater_Than(P->x[i],Xn,P->n)) Xn=P->x[n=i];
376 if(d<dn) {dn=d; Xn=P->x[n=i];}
377 }
378 if(dp) if(dn) /* points on both sides */
379 #if (VERT_WITH_MAX_DISTANCE)
380 {if(dp+dn>0) *v=p; else *v=n;}
381 #else
382 {if(dp+dn>0) *v=n; else *v=p;}
383 #endif
384 else *v=p; /* d >=0 */
385 else if(dn) *v=n; /* d <=0 */
386 else return 0;
387 /* for(i=0;i<P->n;i++) printf(" %d ",Xp[i]); printf(" = Xp Xn =");
388 for(i=0;i<P->n;i++) printf(" %d ",Xn[i]); printf(" \n");
389 */
390 return 1;
391 }
392
393 /* ========= GLZ_Start_Simplex() => return codimension ======= */
GLZ_Start_Simplex(PolyPointList * _P,VertexNumList * _V,CEqList * _C)394 int GLZ_Start_Simplex(PolyPointList *_P, VertexNumList *_V, CEqList *_C)
395 { int i, x=0, y=0, *VN=_V->v, *d=&_P->n, r=*d, b[POLY_Dmax];
396 Long *X=_P->x[x], *Y=_P->x[y], XX=0, YY=0,
397 B[(POLY_Dmax*(POLY_Dmax+1))/2][POLY_Dmax], W[POLY_Dmax]; if(_P->np<2)
398 { for(x=0;x<_P->n;x++) for(y=0;y<_P->n;y++) _C->e[x].a[y]=(x==y);
399 assert(_P->np>0); for(x=0;x<_P->n;x++) _C->e[x].c=-_P->x[0][x];
400 return _C->ne=_P->n;
401 }
402 for(i=1; i<_P->np; i++)
403 { Long *Z=_P->x[i];
404 if(Vec_Greater_Than(X,Z,_P->n)) X=_P->x[x=i]; /* (x_n)-max: VN[0] */
405 if(Vec_Greater_Than(Z,Y,_P->n)) Y=_P->x[y=i]; /* (x_n)-min: VN[1] */
406 } assert(x!=y); /* at this point I need two different vertices */
407 for(i=0;i<*d;i++) { Long Xi=(X[i]>0) ? X[i]: -X[i],
408 Yi=(Y[i]>0) ? Y[i]: -Y[i]; if(Xi>XX)XX=Xi; if(Yi>YY)YY=Yi;}
409 if(YY<XX) {VN[0]=y;VN[1]=x;} else {VN[0]=x;VN[1]=y;} _V->nv=2; y=VN[1];
410 X=_P->x[VN[0]];Y=_P->x[VN[1]]; for(i=0;i<*d;i++) b[i]=(i*(2*(*d)-i+1))/2;
411 for(x=0;x<*d;x++) for(i=0;i<*d;i++) B[x][i]=(x==i); /* b[i+1]-b[i]=d-i */
412 for(x=1;x<*d;x++)
413 { for(i=0;i<*d;i++) W[i]=_P->x[y][i]-X[i];
414 OrthBase_red_by_V(W,d,&B[b[x-1]],&r,&B[b[x]]); for(i=0;i<r;i++)
415 #if (LONG_EQ_FIRST)
416 if(New_Start_Vertex(X,B[b[x]+r-i-1],_P,&y)) break;
417 #else
418 if(New_Start_Vertex(X,B[b[x]+i],_P,&y)) break;
419 #endif
420 if(i==r) break;
421 _V->v[_V->nv++]=y; /* x = dim(span) < d */
422 }
423 if(x<*d)
424 { for(y=0;y<r;y++)
425 { Equation *E=&_C->e[y]; Long *Z=B[b[x]+y];
426 E->c=0; for(i=0;i<*d;i++) E->a[i]=Z[i];
427 E->c=-Eval_Eq_on_V(E,X,_P->n);
428 } return _C->ne=r;
429 }
430 else
431 { Equation *E=_C->e; Long *Z=B[b[*d-1]]; E->c=0; _C->ne=2;
432 for(i=0;i<*d;i++)E->a[i]=Z[i];
433 E->c=-Eval_Eq_on_V(E,X,_P->n);
434 if(Eval_Eq_on_V(E,_P->x[_V->v[*d]],_P->n)<0) {
435 for(i=0;i<*d;i++)E->a[i]=-Z[i];
436 E->c*=-1;}
437 X=_P->x[_V->v[r=*d]];
438 for(x=1;x<*d;x++) /* now the 2nd equation */
439 { Y=_P->x[_V->v[x-1]]; for(i=0;i<*d;i++) W[i]=X[i]-Y[i];
440 OrthBase_red_by_V(W,d,&B[b[x-1]],&r,&B[b[x]]);
441 }
442 E=&_C->e[1]; E->c=0;
443 for(i=0;i<*d;i++)E->a[i]=Z[i];
444 E->c=-Eval_Eq_on_V(E,X,_P->n);
445 assert(XX=Eval_Eq_on_V(E,_P->x[_V->v[*d-1]],_P->n));
446 if(XX<0) {for(i=0;i<*d;i++)E->a[i]=-Z[i]; E->c*=-1;}
447 for(x=*d-2;x>=0;x--) /* omit vertex #x */
448 { r=*d-x; for(y=x+1;y<*d;y++)
449 { Y=_P->x[_V->v[y]]; for(i=0;i<*d;i++) W[i]=X[i]-Y[i];
450 OrthBase_red_by_V(W,d,&B[b[y-1]],&r,&B[b[y]]);
451 }
452 E=&_C->e[(_C->ne)++]; E->c=0;
453 for(i=0;i<*d;i++)E->a[i]=Z[i];
454 E->c=-Eval_Eq_on_V(E,X,_P->n);
455 assert(XX=Eval_Eq_on_V(E,_P->x[_V->v[x]],_P->n));
456 if(XX<0) {for(i=0;i<*d;i++)E->a[i]=-Z[i]; E->c*=-1;}
457 }
458 }
459 assert(*d+1==_C->ne); for(x=0;x<_C->ne;x++) for(i=0;i<=*d;i++)
460 assert((x==i)==(0!=Eval_Eq_on_V(&_C->e[x],_P->x[_V->v[*d-i]],_P->n)));
461 return 0;
462 }
463
464
465 /* ====================================================================== */
466 /* ========== ========== */
467 /* ========== P O L Y H E D R O N A N A L Y S I S ========== */
468 /* ========== ========== */
469 /* ====================================================================== */
470
Make_New_CEqs(PolyPointList * _P,VertexNumList * _V,CEqList * _C,EqList * _F,INCI * CEq_I,INCI * F_I)471 void Make_New_CEqs(PolyPointList *_P, VertexNumList *_V, CEqList *_C,
472 EqList *_F, INCI *CEq_I, INCI *F_I){
473 int i,j, Old_C_ne=_C->ne;
474 static CEqList Bad_C;
475 static INCI Bad_C_I[CEQ_Nmax];
476
477 #if (SHOW_NEW_CEq)
478 static int init; static clock_t CLOCK1; static time_t DATE1;
479 if(!init){init=1; CLOCK1=clock(); DATE1=time(NULL);}
480 printf("V=%d F=%d: Ceq=%d",_V->nv,_F->ne,_C->ne);fflush(stdout);
481 #endif
482
483 Bad_C.ne=_C->ne=0;
484 for (i=0;i<Old_C_ne;i++){
485 Long dist = Eval_Eq_on_V(&_C->e[i],_P->x[_V->v[_V->nv-1]],_P->n);
486 CEq_I[i]=INCI_PN(CEq_I[i],dist);
487 if (dist<0) {Bad_C.e[Bad_C.ne]=_C->e[i]; Bad_C_I[Bad_C.ne++]=CEq_I[i];}
488 else {_C->e[_C->ne]=_C->e[i]; CEq_I[_C->ne++]=CEq_I[i];}}
489
490 #if (SHOW_NEW_CEq)
491 printf("=%dg+%db",_C->ne,Bad_C.ne);fflush(stdout);
492 #endif
493
494 Old_C_ne=_C->ne;
495 for (i=0;i<_F->ne;i++) F_I[i]=
496 INCI_PN(F_I[i],Eval_Eq_on_V(&_F->e[i],_P->x[_V->v[_V->nv-1]],_P->n));
497 for (j=0;j<_F->ne;j++) if (!INCI_M2(F_I[j]))
498 for (i=0;i<Bad_C.ne;i++){
499 INCI New_Face=INCI_AND(Bad_C_I[i],F_I[j]);
500 int k;
501 if (INCI_abs(New_Face)<_P->n-1) continue;
502 for (k=0;k<Bad_C.ne;k++) if (INCI_LE(New_Face,Bad_C_I[k]))
503 if (k!=i) break;
504 if (k!=Bad_C.ne) continue;
505 for (k=0;k<Old_C_ne;k++) if (INCI_LE(New_Face,CEq_I[k])) break;
506 if (k!=Old_C_ne) continue;
507 for (k=0;k<_F->ne;k++) if (INCI_LE(New_Face,F_I[k])) if (k!=j) break;
508 if (k!=_F->ne) continue;
509 assert(_C->ne<CEQ_Nmax);
510 CEq_I[_C->ne]=INCI_PN(INCI_D2(New_Face),0);
511 _C->e[_C->ne]=EEV_To_Equation(&(Bad_C.e[i]),&(_F->e[j]),
512 _P->x[_V->v[_V->nv-1]],_P->n);
513 assert(IsGoodCEq(&(_C->e[_C->ne++]),_P,_V));}
514 for (j=0;j<Old_C_ne;j++) if (!INCI_M2(CEq_I[j]))
515 for (i=Bad_C.ne-1;i>=0;i--){
516 INCI New_Face=INCI_AND(Bad_C_I[i],CEq_I[j]);
517 int k;
518 if (INCI_abs(New_Face)<_P->n-1) continue;
519 for (k=0;k<Bad_C.ne;k++) if (INCI_LE(New_Face,Bad_C_I[k]))
520 if (k!=i) break;
521 if (k!=Bad_C.ne) continue;
522 for (k=0;k<Old_C_ne;k++) if (INCI_LE(New_Face,CEq_I[k]))
523 if (k!=j) break;
524 if (k!=Old_C_ne) continue;
525 for (k=0;k<_F->ne;k++) if (INCI_LE(New_Face,F_I[k])) break;
526 if (k!=_F->ne) continue;
527 assert(_C->ne<CEQ_Nmax);
528 CEq_I[_C->ne]=INCI_PN(INCI_D2(New_Face),0);
529 _C->e[_C->ne]=EEV_To_Equation(&(Bad_C.e[i]),&(_C->e[j]),
530 _P->x[_V->v[_V->nv-1]],_P->n);
531 assert(IsGoodCEq(&(_C->e[_C->ne++]),_P,_V));}
532
533 #if (SHOW_NEW_CEq)
534 {time_t DATE2=time(NULL); char sm[2]={'s',0};
535 int Rs= (int)difftime(DATE2,DATE1);if(Rs>999){Rs/=60; *sm='m';}
536 printf(" done: C.ne=%d %d%s\n",_C->ne,Rs,sm);fflush(0);}
537 #endif
538 }
539
540
541 #if MAX_BAD_EQ
542
IP_Search_Bad_Eq(CEqList * _C,EqList * _F,INCI * CEq_I,INCI * F_I,PolyPointList * _P,int * _IP)543 int IP_Search_Bad_Eq(CEqList *_C, EqList *_F, INCI *CEq_I, INCI *F_I,
544 PolyPointList *_P, int *_IP){ /* return 0 :: no bad eq. */
545 while(_C->ne--) {
546 int j, M=_C->ne; /* INCI_LmR INCI_lex_GT */
547 for(j=0;j<_C->ne;j++) if(INCI_lex_GT(&CEq_I[j],&CEq_I[M])) M=j;
548 for(j=0;j<_P->np;j++)
549 if(Eval_Eq_on_V(&(_C->e[M]),_P->x[j],_P->n) < 0) {
550 INCI AI=CEq_I[M]; Equation AE=_C->e[M];
551 CEq_I[M]=CEq_I[_C->ne]; _C->e[M]=_C->e[_C->ne];
552 CEq_I[_C->ne]=AI; _C->e[_C->ne]=AE; return ++_C->ne;}
553 if(_C->e[M].c < 1) {*_IP=0; return 1;}
554 assert(_F->ne<EQUA_Nmax);
555 /* printf("#Feq=%d #Ceq=%d\n",_F->ne,_C->ne); fflush(stdout); */
556 _F->e[_F->ne]=_C->e[M]; F_I[_F->ne++]=CEq_I[M];
557 if(M<_C->ne) {_C->e[M]=_C->e[_C->ne]; CEq_I[M]=CEq_I[_C->ne];} }
558 return 0;
559 }
560
FE_Search_Bad_Eq(CEqList * _C,EqList * _F,INCI * CEq_I,INCI * F_I,PolyPointList * _P,int * _IP)561 int FE_Search_Bad_Eq(CEqList *_C, EqList *_F, INCI *CEq_I, INCI *F_I,
562 PolyPointList *_P, int *_IP){ /* return 0 :: no bad eq. */
563 while(_C->ne--) {
564 int j, M=_C->ne; /* INCI_LmR INCI_lex_GT */
565 for(j=0;j<_C->ne;j++) if(INCI_lex_GT(&CEq_I[j],&CEq_I[M])) M=j;
566 for(j=0;j<_P->np;j++)
567 if(Eval_Eq_on_V(&(_C->e[M]),_P->x[j],_P->n) < 0) {
568 INCI AI=CEq_I[M]; Equation AE=_C->e[M];
569 CEq_I[M]=CEq_I[_C->ne]; _C->e[M]=_C->e[_C->ne];
570 CEq_I[_C->ne]=AI; _C->e[_C->ne]=AE; return ++_C->ne;}
571 if(_C->e[M].c < 1) *_IP=0;
572 assert(_F->ne<EQUA_Nmax);
573 /* printf("#Feq=%d #Ceq=%d\n",_F->ne,_C->ne); fflush(stdout); */
574 _F->e[_F->ne]=_C->e[M]; F_I[_F->ne++]=CEq_I[M];
575 if(M<_C->ne) {_C->e[M]=_C->e[_C->ne]; CEq_I[M]=CEq_I[_C->ne];}
576 }
577 return 0;
578 }
579
580 #else
581
IP_Search_Bad_Eq(CEqList * _C,EqList * _F,INCI * CEq_I,INCI * F_I,PolyPointList * _P,int * _IP)582 int IP_Search_Bad_Eq(CEqList *_C, EqList *_F, INCI *CEq_I, INCI *F_I,
583 PolyPointList *_P, int *_IP){ /* return 0 :: no bad eq. */
584 while(_C->ne--) {
585 int j;
586 for(j=0;j<_P->np;j++)
587 if(Eval_Eq_on_V(&(_C->e[_C->ne]),_P->x[j],_P->n) < 0) return ++_C->ne;
588 if(_C->e[_C->ne].c < 1) { *_IP=0; return 1;}
589 assert(_F->ne<EQUA_Nmax);
590 _F->e[_F->ne]=_C->e[_C->ne];
591 F_I[_F->ne++]=CEq_I[_C->ne];}
592 return 0;
593 }
594
FE_Search_Bad_Eq(CEqList * _C,EqList * _F,INCI * CEq_I,INCI * F_I,PolyPointList * _P,int * _IP)595 int FE_Search_Bad_Eq(CEqList *_C, EqList *_F, INCI *CEq_I, INCI *F_I,
596 PolyPointList *_P, int *_IP){ /* return 0 :: no bad eq. */
597 while(_C->ne--) {
598 int j;
599 for(j=0;j<_P->np;j++)
600 if(Eval_Eq_on_V(&(_C->e[_C->ne]),_P->x[j],_P->n) < 0) return ++_C->ne;
601 if(_C->e[_C->ne].c < 1) *_IP=0;
602 assert(_F->ne<EQUA_Nmax);
603 _F->e[_F->ne]=_C->e[_C->ne];
604 F_I[_F->ne++]=CEq_I[_C->ne];}
605 return 0;
606 }
607
608 #endif
609
REF_Search_Bad_Eq(CEqList * _C,EqList * _F,INCI * CEq_I,INCI * F_I,PolyPointList * _P,int * _REF)610 int REF_Search_Bad_Eq(CEqList *_C, EqList *_F, INCI *CEq_I, INCI *F_I,
611 PolyPointList *_P, int *_REF){ /* return 0 :: no bad eq. */
612 while(_C->ne--) {
613 int j;
614 for(j=0;j<_P->np;j++)
615 if(Eval_Eq_on_V(&(_C->e[_C->ne]),_P->x[j],_P->n) < 0) return ++_C->ne;
616 if(_C->e[_C->ne].c != 1) { *_REF=0; return 1;}
617 assert(_F->ne<EQUA_Nmax);
618 _F->e[_F->ne]=_C->e[_C->ne];
619 F_I[_F->ne++]=CEq_I[_C->ne];}
620 return 0;
621 }
622
Finish_Find_Equations(PolyPointList * _P,VertexNumList * _V,EqList * _F,CEqList * _CEq,INCI * F_I,INCI * CEq_I)623 int Finish_Find_Equations(PolyPointList *_P, VertexNumList *_V,
624 EqList *_F, CEqList *_CEq, INCI *F_I, INCI *CEq_I){
625 int IP=1;
626 while(0<=_CEq->ne) if (FE_Search_Bad_Eq(_CEq,_F,CEq_I,F_I,_P,&IP)){
627 assert(_V->nv<VERT_Nmax);
628 _V->v[_V->nv++]=Search_New_Vertex(&(_CEq->e[_CEq->ne-1]),_P);
629 Make_New_CEqs(_P,_V,_CEq,_F,CEq_I,F_I); }
630 return IP;
631 }
632
Find_Equations(PolyPointList * _P,VertexNumList * _V,EqList * _F)633 int Find_Equations(PolyPointList *_P, VertexNumList *_V, EqList *_F){
634 /* return: IP, finds Vertices and Equations for _P even if not IP */
635 int i;
636 CEqList *CEq = (CEqList *) malloc(sizeof(CEqList));
637 INCI *CEq_I = (INCI *) malloc(sizeof(INCI)*CEQ_Nmax);
638 INCI *F_I = (INCI *) malloc(sizeof(INCI)*EQUA_Nmax);
639 CEq->ne=0;
640 if((CEq==NULL)||(CEq_I==NULL)||(F_I==NULL)) {
641 printf("Allocation failure in Find_Equations\n"); exit(0);}
642 if (GLZ_Start_Simplex(_P, _V, CEq)) {
643 _F->ne=CEq->ne;
644 for(i=0;i<_F->ne;i++) _F->e[i]=CEq->e[i];
645 free(CEq); free(CEq_I); free(F_I);
646 return 0;}
647 _F->ne=0;
648 for (i=0;i<CEq->ne;i++)
649 if(INCI_abs(CEq_I[i]=Eq_To_INCI(&(CEq->e[i]),_P,_V))<_P->n)
650 {fprintf(outFILE,"Bad CEq in Find_Equations"); exit(0);}
651 i=Finish_Find_Equations(_P, _V, _F, CEq, F_I, CEq_I);
652 free(CEq); free(CEq_I); free(F_I);
653 return i;
654 }
655
Finish_IP_Check(PolyPointList * _P,VertexNumList * _V,EqList * _F,CEqList * _CEq,INCI * F_I,INCI * CEq_I)656 int Finish_IP_Check(PolyPointList *_P, VertexNumList *_V, EqList *_F,
657 CEqList *_CEq, INCI *F_I, INCI *CEq_I){
658 int IP=1;
659 while(0<=_CEq->ne) if (IP_Search_Bad_Eq(_CEq,_F,CEq_I,F_I,_P,&IP)){
660 if(!IP) return 0; /* found d<=0 */
661 assert(_V->nv<VERT_Nmax);
662 _V->v[_V->nv++]=Search_New_Vertex(&(_CEq->e[_CEq->ne-1]),_P);
663 Make_New_CEqs(_P,_V,_CEq,_F,CEq_I,F_I); }
664 return 1;
665 }
666
IP_Check(PolyPointList * _P,VertexNumList * _V,EqList * _F)667 int IP_Check(PolyPointList *_P, VertexNumList *_V, EqList *_F){
668 int i;
669 CEqList *CEq = (CEqList *) malloc(sizeof(CEqList));
670 INCI *CEq_I = (INCI *) malloc(sizeof(INCI)*CEQ_Nmax);
671 INCI *F_I = (INCI *) malloc(sizeof(INCI)*EQUA_Nmax);
672 if((CEq==NULL)||(CEq_I==NULL)||(F_I==NULL)) {
673 printf("Allocation failure in IP_Check\n"); exit(0);}
674 if (GLZ_Start_Simplex(_P, _V, CEq)) {
675 free(CEq); free(CEq_I); free(F_I); return 0;}
676 for (i=0;i<CEq->ne;i++)
677 if(INCI_abs(CEq_I[i]=Eq_To_INCI(&(CEq->e[i]),_P,_V))<_P->n)
678 {fprintf(outFILE,"Bad CEq in IP_Check"); exit(0);}
679 _F->ne=0;
680 i=Finish_IP_Check(_P, _V, _F, CEq, F_I, CEq_I);
681 free(CEq); free(CEq_I); free(F_I);
682 return i;
683 }
684
Finish_REF_Check(PolyPointList * _P,VertexNumList * _V,EqList * _F,CEqList * _CEq,INCI * F_I,INCI * CEq_I)685 int Finish_REF_Check(PolyPointList *_P, VertexNumList *_V, EqList *_F,
686 CEqList *_CEq, INCI *F_I, INCI *CEq_I){
687 int REF=1;
688 while(0<=_CEq->ne) if(REF_Search_Bad_Eq(_CEq,_F,CEq_I,F_I,_P,&REF)){
689 if(!REF) return 0; /* found d!=1 */
690 assert(_V->nv<VERT_Nmax);
691 _V->v[_V->nv++]=Search_New_Vertex(&(_CEq->e[_CEq->ne-1]),_P);
692 Make_New_CEqs(_P,_V,_CEq,_F,CEq_I,F_I); }
693 return 1;
694 }
695
Ref_Check(PolyPointList * _P,VertexNumList * _V,EqList * _F)696 int Ref_Check(PolyPointList *_P, VertexNumList *_V, EqList *_F){
697 int i;
698 CEqList *CEq = (CEqList *) malloc(sizeof(CEqList));
699 INCI *CEq_I = (INCI *) malloc(sizeof(INCI)*CEQ_Nmax);
700 INCI *F_I = (INCI *) malloc(sizeof(INCI)*EQUA_Nmax);
701 if((CEq==NULL)||(CEq_I==NULL)||(F_I==NULL)) {
702 printf("Allocation failure in Ref_Check\n"); exit(0);}
703 if (GLZ_Start_Simplex(_P, _V, CEq)) {
704 free(CEq); free(CEq_I); free(F_I); return 0;}
705 for (i=0;i<CEq->ne;i++) CEq_I[i]=Eq_To_INCI(&(CEq->e[i]),_P,_V);
706 _F->ne=0;
707 i=Finish_REF_Check(_P, _V, _F, CEq, F_I, CEq_I);
708 free(CEq); free(CEq_I); free(F_I);
709 return i;
710 }
711
712
713 /* ====================================================================== */
714 /* ========== ========== */
715 /* ========== D U A L P O L Y & C O M P L E T I O N ========== */
716 /* ========== ========== */
717 /* ====================================================================== */
718
Make_Dual_Poly(PolyPointList * _P,VertexNumList * _V,EqList * _E,PolyPointList * _DP)719 void Make_Dual_Poly(PolyPointList *_P, VertexNumList *_V, EqList *_E,
720 PolyPointList *_DP){
721 EqList DE;
722 PairMat PM, DPM;
723 Make_VEPM(_P,_V,_E,PM);
724 Transpose_PM(PM, DPM, _V->nv, _E->ne);
725 VNL_to_DEL(_P,_V,&DE);
726 _DP->n=_P->n;
727 _DP->np=0;
728 Complete_Poly(DPM, &DE, _E->ne, _DP);
729 }
730
add_for_completion(Long * yDen,Long Den,EqList * _E,PolyPointList * _CP,int * old_np)731 void add_for_completion(Long *yDen, Long Den,
732 EqList *_E, PolyPointList *_CP, int *old_np){
733 int i,n=_CP->n;
734 Long yold[POLY_Dmax];
735
736 if(Den>1) for(i=0;i<n;i++) {
737 if(yDen[i]%Den) return;
738 yold[i]=yDen[i]/Den;}
739 else for(i=0;i<n;i++) yold[i]=yDen[i];
740 for (i=0;i<_E->ne;i++) if (Eval_Eq_on_V(&(_E->e[i]), yold, n) < 0) return;
741 for (i=0;i<*old_np;i++) if (Vec_Equal(_CP->x[i],yold,n)) return;
742 assert(_CP->np<POINT_Nmax);
743 for(i=0;i<n;i++) _CP->x[_CP->np][i]=yold[i];
744 _CP->np++;
745 }
746
Compute_InvMat(int n,EqList * _E,int OrdFac[VERT_Nmax],int BasFac[POLY_Dmax],Long * Den,Long InvMat[POLY_Dmax][POLY_Dmax])747 void Compute_InvMat(int n, EqList *_E, int OrdFac[VERT_Nmax],
748 int BasFac[POLY_Dmax], Long *Den,
749 Long InvMat[POLY_Dmax][POLY_Dmax]){
750 /* Find first POLY_Dmax linearly independent facets + Inverse Matrix */
751
752 LRat ind[POLY_Dmax][POLY_Dmax], x[POLY_Dmax], y[POLY_Dmax], f,
753 PInvMat[POLY_Dmax][POLY_Dmax];
754 int i, j, k, l, rank=0, one[POLY_Dmax];
755
756 for (i=0;i<n;i++) for (j=0;j<n;j++) PInvMat[i][j]=LrI(0);
757 for (i=0;i<n;i++) PInvMat[i][i]=LrI(1);
758 i=0;
759 while (rank<n){
760 for (j=0;j<n;j++) x[j]=LrI(_E->e[OrdFac[i]].a[j]);
761 for (j=0;j<n;j++) y[j]=LrI(0);
762 y[rank]=LrI(1);
763 for (j=0;j<rank;j++) {
764 f=x[one[j]];
765 for (k=0;k<n;k++) {
766 x[k]=LrD(x[k],LrP(f,ind[j][k]));
767 y[k]=LrD(y[k],LrP(f,PInvMat[j][k])); } }
768 one[rank]=-1;
769 for (l=0;(l<n)&&(one[rank]==-1);l++) if (x[l].N) one[rank]=l;
770 if(one[rank]>-1){
771 for (k=0;k<n;k++) {
772 ind[rank][k]=LrQ(x[k],x[one[rank]]);
773 PInvMat[rank][k]=LrQ(y[k],x[one[rank]]); }
774 for (j=0;j<rank;j++) {
775 f=ind[j][one[rank]];
776 for (k=0;k<n;k++) {
777 ind[j][k]=LrD(ind[j][k],LrP(ind[rank][k],f));
778 PInvMat[j][k]=LrD(PInvMat[j][k],LrP(PInvMat[rank][k],f)); } }
779 BasFac[rank]=OrdFac[i];
780 rank++; }
781 i++; }
782 for (i=0;i<n;i++) for (j=0;j<n;j++)
783 *Den=(*Den/LFgcd(*Den,PInvMat[i][j].D))*PInvMat[i][j].D;
784 for (i=0;i<n;i++) for (j=0;j<n;j++)
785 InvMat[one[i]][j]=(*Den/PInvMat[i][j].D)*PInvMat[i][j].N;
786
787 for (i=0;i<n;i++){
788 for (j=0;j<n;j++) {
789 long long s=0;
790 for(k=0;k<n;k++) s+=((long long) (InvMat[k][i]))*
791 ((long long) (_E->e[BasFac[j]].a[k]));
792 if (s!=*Den*(i==j)) {
793 puts("something wrong in Make_Dual_Poly");
794 exit(0);}}}
795 }
796
Complete_Poly(PairMat VPM,EqList * _E,int nv,PolyPointList * _CP)797 void Complete_Poly(PairMat VPM, EqList *_E, int nv,
798 PolyPointList *_CP){
799 int i, j, k, InsPoint, n=_CP->n, old_np=_CP->np;
800 Long MaxDist[EQUA_Nmax], InvMat[POLY_Dmax][POLY_Dmax], Den=1;
801 Long yDen[POLY_Dmax];
802 int OrdFac[VERT_Nmax], BasFac[POLY_Dmax], position[POLY_Dmax];
803
804 /*_CP->np=0;*/
805
806 /* Calculate maximal distances from facets of Delta^* (Vertices of Delta) */
807
808 for (i=0;i<_E->ne;i++) {
809 MaxDist[i]=0;
810 for (j=0;j<nv;j++)
811 if (MaxDist[i]<VPM[i][j]) MaxDist[i]=VPM[i][j];}
812
813 /* Order facets of Delta^* (Vertices of Delta) w.r.t. MaxDist */
814
815 OrdFac[0]=0;
816 for (i=1;i<_E->ne;i++){
817 InsPoint=i;
818 while (InsPoint&&(MaxDist[i]<MaxDist[OrdFac[InsPoint-1]])) InsPoint--;
819 for (j=i;j>InsPoint;j--) OrdFac[j]=OrdFac[j-1];
820 OrdFac[InsPoint]=i; }
821
822 Compute_InvMat(n, _E, OrdFac, BasFac, &Den, InvMat);
823
824 /* printf("Den=%ld ", Den); mostly 1 or 2!!! */
825
826 /* Examine all integer points of parallelogram: */
827 /* The basic structure of the algorithm is:
828 for (k=0;k<n-1;k++) position[k]=-1; / * sets k=n-1; important! *
829 position[n-1]=-2; / * starting point just outside the parallelogram *
830 while(k>=0){
831 position[k]++;
832 DO AT position;
833 for(k=n-1;((position[k]==MaxDist[BasFac[k]]-1)&&(k>=0));k--)
834 position[k]=-1; }
835 / * sets k to the highest value where pos.[k] wasn't the max value;
836 resets the following max values to min values */
837 /* Quantities linear in position can be changed with every change of
838 position (here: yDen) */
839
840 for(i=0;i<n;i++) yDen[i]=0;
841 for (k=0;k<n-1;k++) { /* sets k=n-1; important! */
842 position[k]=-_E->e[BasFac[k]].c;
843 for(i=0;i<n;i++) yDen[i]-=_E->e[BasFac[k]].c*InvMat[i][k]; }
844 position[n-1]=-_E->e[BasFac[n-1]].c-1;
845 for(i=0;i<n;i++) yDen[i]-=(_E->e[BasFac[k]].c+1)*InvMat[i][n-1];
846 while(k>=0){
847 position[k]++;
848 for(i=0;i<n;i++) yDen[i]+=InvMat[i][k];
849 add_for_completion(yDen, Den, _E, _CP, &old_np);
850 for(k=n-1;(k>=0);k--){
851 if (position[k]!=MaxDist[BasFac[k]]-_E->e[BasFac[k]].c) break;
852 position[k]=-_E->e[BasFac[k]].c;
853 for (i=0;i<n;i++) yDen[i]-=MaxDist[BasFac[k]]*InvMat[i][k]; }}
854 }
855
856
857
858 /* ====================================================================== */
859 /* ========== ========== */
860 /* ========== B A T Y R E V ' S F O R M U L A S ========== */
861 /* ========== ========== */
862 /* ====================================================================== */
863
RaiseDip(INCI x,FaceInfo * _I,int n,int mult)864 void RaiseDip(INCI x, FaceInfo *_I, int n, int mult){
865 int i, j;
866 for(i=0;i<n;i++)
867 for(j=0;j<_I->nf[i];j++)
868 if(INCI_EQ(x,_I->v[i][j])) {_I->dip[i][j] += mult; return;}
869 }
870
RaiseNip(INCI x,FaceInfo * _I,int n)871 void RaiseNip(INCI x, FaceInfo *_I, int n){
872 int i, j;
873 for(i=n-1;i>=0;i--)
874 for(j=0;j<_I->nf[i];j++)
875 if(INCI_EQ(x,_I->f[i][j])) {_I->nip[i][j]++; return;}
876 }
877
Make_FaceIPs(PolyPointList * _P,VertexNumList * _V,EqList * _E,PolyPointList * _DP,FaceInfo * _I)878 void Make_FaceIPs(PolyPointList *_P, VertexNumList *_V, EqList *_E,
879 PolyPointList *_DP, FaceInfo *_I){
880 /* compute IP's of faces by computing Incidences for all points and
881 * comparing with Incidences of dual faces */
882 int i, j, k;
883 INCI x;
884 for(i=0;i<_P->n;i++) for(j=0;j<_I->nf[i];j++) {
885 _I->nip[i][j]=0;
886 _I->dip[i][j]=0; }
887 for(k=0;k<_P->np;k++){
888 x=INCI_0();
889 for(i=0;i<_E->ne;i++)
890 x=INCI_PN(x,Eval_Eq_on_V(&(_E->e[i]),_P->x[k],_P->n));
891 RaiseNip(x, _I, _P->n);}
892 for (k=0;k<_DP->np;k++){
893 x=INCI_0();
894 for(i=0;i<_V->nv;i++)
895 x=INCI_PN(x,DualBraP1(_P->x[_V->v[i]],_DP->x[k],_P->n));
896 RaiseDip(x, _I, _P->n, 1);}
897 }
898
PrintFaceIPs(PolyPointList * _P,FaceInfo * _I)899 void PrintFaceIPs(PolyPointList *_P,FaceInfo *_I){
900 int i,j, M=_P->n-1;
901 for(i=0;i<=M;i++) {
902 printf("ip[%d]:",i);
903 for(j=0;j<_I->nf[i];j++) printf(" %ld",(long) _I->nip[i][j]);
904 puts(""); }
905 for(i=0;i<=M;i++) {
906 printf("dip[%d]:",i);
907 for(j=0;j<_I->nf[i];j++)printf(" %ld",(long) _I->dip[i][j]);
908 puts(""); }
909 }
910
Eval_BaHo(FaceInfo * _I,BaHo * _BH)911 void Eval_BaHo(FaceInfo *_I, BaHo *_BH){
912 /* Calculate Hodge/Picard numbers from FaceInfo */
913 int i,j, n=_BH->n;
914 int *h1;
915 _BH->cor=0;
916 h1=_BH->h1;
917 for(i=0;i<n-1;i++) h1[i]=0;
918 h1[1]+=_BH->np-n-1;
919 for (i=0;i<_I->nf[0];i++) h1[1]-=_I->dip[0][i];
920 for (j=1;j<n-1;j++) for (i=0;i<_I->nf[j];i++) {
921 h1[j]+=_I->dip[j][i]*_I->nip[j][i];
922 _BH->cor+=_I->dip[j][i]*_I->nip[j][i];}
923 if (n==3) _BH->pic=h1[1];
924 for (i=0;i<_I->nf[n-1];i++) h1[n-2]-=_I->nip[n-1][i];
925 h1[n-2]+=_BH->mp-n-1;
926 if (n==5) _BH->h22=44+4*h1[1]+4*h1[3]-2*h1[2];
927 }
928
RC_Calc_BaHo(PolyPointList * _P,VertexNumList * _V,EqList * _E,PolyPointList * _DP,BaHo * _BH)929 void RC_Calc_BaHo(PolyPointList *_P, VertexNumList *_V, EqList *_E,
930 PolyPointList *_DP, BaHo *_BH){
931 /* Needs reflexive and complete _P and _DP */
932 FaceInfo *_FI=(FaceInfo *) malloc(sizeof(FaceInfo));
933 if(_FI==NULL) {printf("RC_Calc_BaHo: Unable to allocate _FI\n"); exit(0);}
934 _BH->mp=_P->np; _BH->mv=_V->nv; _BH->nv=_E->ne; _BH->np=_DP->np;
935 _BH->n=_P->n;
936 Make_Incidence(_P, _V, _E, _FI);
937 Make_FaceIPs(_P, _V, _E, _DP, _FI);
938 Eval_BaHo(_FI, _BH);
939 free(_FI);
940 }
941
942
943
944
Qadd_for_completion(Long * yDen,Long Den,int n,int * np,EqList * _E,FaceInfo * _I)945 void Qadd_for_completion(Long *yDen, Long Den, int n, int *np, EqList *_E,
946 FaceInfo *_I){
947 int i;
948 Long y[POLY_Dmax];
949 INCI x = INCI_0();
950 for(i=0;i<n;i++) if(yDen[i]%Den) return;
951 for(i=0;i<n;i++) y[i] = yDen[i]/Den;
952 for(i=0;i<_E->ne;i++) x=INCI_PN(x,Eval_Eq_on_V(&(_E->e[i]),y,n));
953 RaiseDip(x, _I, n, 1);
954 (*np)++;
955 }
956
lastline(Long * EyD,Long * yDen,Long Den,EqList * _E,int n,int * np,Long MD,Long IME[POLY_Dmax][EQUA_Nmax],Long InvMat[POLY_Dmax][POLY_Dmax],FaceInfo * _I)957 void lastline(Long *EyD, Long *yDen, Long Den, EqList *_E, int n, int *np,
958 Long MD /*MaxDist[BasFac[n-1]*/, Long IME[POLY_Dmax][EQUA_Nmax],
959 Long InvMat[POLY_Dmax][POLY_Dmax], FaceInfo *_I){
960 int i, j, l, lmin=0, lmax=MD;
961 INCI xmin, xmax;
962 Long y[POLY_Dmax];
963 for (j=0;j<_E->ne;j++) if ((EyD[j]<0)&&(IME[n-1][j]<=0)) return;
964 for (j=0;(j<_E->ne)&&(lmax >=lmin);j++)
965 if (IME[n-1][j]<0) {
966 l = EyD[j]/(-IME[n-1][j]);
967 if (lmax > l) lmax = l;}
968 else if ((IME[n-1][j]>0)&&(EyD[j]<0)){
969 l = 1 + (-1-EyD[j])/IME[n-1][j];
970 if (lmin < l) lmin = l;}
971 if (lmax < lmin) return;
972
973 if(Den>1) {
974 for(i=0;i<n;i++)
975 yDen[i] += InvMat[i][n-1] * lmin;
976 for (l=lmin;l<=lmax;l++){
977 Qadd_for_completion(yDen, Den, n, np, _E, _I);
978 for (i=0;i<n;i++) yDen[i]+=InvMat[i][n-1];}
979 for(i=0;i<n;i++)
980 yDen[i] -= InvMat[i][n-1] * (lmax+1);
981 return;}
982
983 for(i=0;i<n;i++)
984 y[i] = yDen[i] + InvMat[i][n-1] * lmin;
985 xmin=INCI_0();
986 for(i=0;i<_E->ne;i++)
987 xmin=INCI_PN(xmin,Eval_Eq_on_V(&(_E->e[i]),y,n));
988 RaiseDip(xmin, _I, n, 1);
989 (*np)++;
990 if (lmax == lmin) return;
991 for(i=0;i<n;i++)
992 y[i] = yDen[i] + InvMat[i][n-1] * lmax;
993 xmax=INCI_0();
994 for(i=0;i<_E->ne;i++)
995 xmax=INCI_PN(xmax,Eval_Eq_on_V(&(_E->e[i]),y,n));
996 RaiseDip(xmax, _I, n, 1);
997 (*np)++;
998 if (lmax == lmin +1) return;
999 RaiseDip(INCI_AND(xmin,xmax), _I, n, lmax - lmin - 1);
1000 (*np) += lmax - lmin - 1;
1001 }
1002
QComplete_Poly(PairMat VPM,EqList * _E,int nv,int n,int * np,FaceInfo * _I)1003 void QComplete_Poly(PairMat VPM, EqList *_E, int nv, int n, int *np,
1004 FaceInfo *_I){
1005 int i, j, k, InsPoint;
1006 Long MaxDist[EQUA_Nmax], InvMat[POLY_Dmax][POLY_Dmax], Den=1;
1007 Long yDen[POLY_Dmax], EyD[EQUA_Nmax], IME[POLY_Dmax][EQUA_Nmax];
1008 /* IME[k][j] = InvMat [k] * _E[j] */
1009 int OrdFac[VERT_Nmax], BasFac[POLY_Dmax], position[POLY_Dmax];
1010
1011 for(i=0;i<n;i++) for(j=0;j<_I->nf[i];j++) _I->dip[i][j]=0;
1012
1013 (*np) = 0;
1014
1015 /* Calculate maximal distances from facets of Delta^* (Vertices of Delta) */
1016
1017 for (i=0;i<_E->ne;i++) {
1018 MaxDist[i]=0;
1019 for (j=0;j<nv;j++)
1020 if (MaxDist[i]<VPM[i][j]) MaxDist[i]=VPM[i][j];}
1021
1022 /* Order facets of Delta^* (Vertices of Delta) w.r.t. MaxDist */
1023
1024 OrdFac[0]=0;
1025 for (i=1;i<_E->ne;i++){
1026 InsPoint=i;
1027 while (InsPoint&&(MaxDist[i]<MaxDist[OrdFac[InsPoint-1]])) InsPoint--;
1028 for (j=i;j>InsPoint;j--) OrdFac[j]=OrdFac[j-1];
1029 OrdFac[InsPoint]=i; }
1030
1031 Compute_InvMat(n, _E, OrdFac, BasFac, &Den, InvMat);
1032
1033 for (k=0;k<n;k++)
1034 for (j=0;j<_E->ne;j++){
1035 IME[k][j] = 0;
1036 for(i=0;i<n;i++) IME[k][j] += InvMat[i][k] * _E->e[j].a[i];}
1037 /* Examine all integer points of parallelogram: */
1038 /* The basic structure of the algorithm is:
1039 for (k=0;k<n-1;k++) position[k]=-1; / * sets k=n-1; important! *
1040 position[n-1]=-2; / * starting point just outside the parallelogram *
1041 while(k>=0){
1042 position[k]++;
1043 DO AT position;
1044 for(k=n-1;((position[k]==MaxDist[BasFac[k]]-1)&&(k>=0));k--)
1045 position[k]=-1; }
1046 / * sets k to the highest value where pos.[k] wasn't the max value;
1047 resets the following max values to min values */
1048 /* Quantities linear in position can be changed with every change of
1049 position (here: yDen) */
1050 for(i=0;i<n;i++) yDen[i]=0;
1051 for (k=0;k<n-2;k++) { /* sets k=n-2; important! */
1052 position[k]=-_E->e[BasFac[k]].c;
1053 for(i=0;i<n;i++) yDen[i]-=_E->e[BasFac[k]].c*InvMat[i][k];}
1054 position[n-2]=-_E->e[BasFac[n-2]].c-1;
1055 for(i=0;i<n;i++) {
1056 yDen[i]-=(_E->e[BasFac[k]].c+1)*InvMat[i][k];
1057 yDen[i]-=_E->e[BasFac[n-1]].c*InvMat[i][n-1];}
1058 for (j=0;j<_E->ne;j++) {
1059 EyD[j]=_E->e[i].c * Den;
1060 for(i=0;i<n;i++) EyD[j] += yDen[i] * _E->e[j].a[i]; }
1061 while(k>=0){
1062 position[k]++;
1063 for (i=0;i<n;i++) yDen[i]+=InvMat[i][k];
1064 for (j=0;j<_E->ne;j++) EyD[j] += IME[k][j];
1065 lastline(EyD, yDen, Den, _E, n, np, MaxDist[BasFac[n-1]], IME, InvMat, _I);
1066 for (k=n-2;(k>=0);k--){
1067 if (position[k]!=MaxDist[BasFac[k]]-_E->e[BasFac[k]].c) break;
1068 position[k]=-_E->e[BasFac[k]].c;
1069 for (i=0;i<n;i++) yDen[i]-=MaxDist[BasFac[k]]*InvMat[i][k];
1070 for (j=0;j<_E->ne;j++)
1071 EyD[j] -= MaxDist[BasFac[k]] * IME[k][j];}}
1072 }
1073
Compute_nip(PolyPointList * _P,EqList * _E,FaceInfo * _I)1074 void Compute_nip(PolyPointList *_P, EqList *_E, FaceInfo *_I){
1075 /* compute IP's of faces by computing Incidences for all points and
1076 * comparing with Incidences of dual faces */
1077 int i, j, k;
1078 INCI x;
1079 for(i=0;i<_P->n;i++) for(j=0;j<_I->nf[i];j++) _I->nip[i][j]=0;
1080 for(k=0;k<_P->np;k++){
1081 x=INCI_0();
1082 for(i=0;i<_E->ne;i++)
1083 x=INCI_PN(x,Eval_Eq_on_V(&(_E->e[i]),_P->x[k],_P->n));
1084 RaiseNip(x, _I, _P->n);}
1085 }
1086
QuickAnalysis(PolyPointList * _P,BaHo * _BH,FaceInfo * _FI)1087 int QuickAnalysis(PolyPointList *_P, BaHo *_BH, FaceInfo *_FI){
1088 VertexNumList V; EqList E; EqList DE; PairMat PM; PairMat DPM; //FaceInfo FI;
1089 int i;
1090 if (!IP_Check(_P,&V,&E)) return 0;
1091 _BH->mp=_P->np; _BH->mv=V.nv; _BH->nv=E.ne; _BH->np = 0; _BH->n=_P->n;
1092 for(i = 0; i < E.ne; i++) if (E.e[i].c != 1) return 1;
1093 Make_VEPM(_P,&V,&E,PM);
1094 if (!Transpose_PM(PM, DPM, V.nv, E.ne)) {
1095 fprintf(stderr,"Transpose_PM failed because #eq=%d > VERT_Nmax\n", E.ne);
1096 exit(0);}
1097 VNL_to_DEL(_P, &V, &DE);
1098 Make_Incidence(_P, &V, &E, _FI);
1099 Compute_nip(_P, &E, _FI);
1100 QComplete_Poly(DPM, &DE, E.ne, _P->n, &_BH->np, _FI);
1101 Eval_BaHo(_FI, _BH);
1102 return 1;
1103 }
1104