1 #include "Global.h"
2 #include "Nef.h"
3 #include "Rat.h"
4
5 /* =============== Typedefs and Headers =================== */
6 #define min(a,b) (((a)<(b)) ? (a):(b))
7 #define max(a,b) (((a)>(b)) ? (a):(b))
8
9 void Sort_PPL(PolyPointList *_P, VertexNumList *_V);
10 void part_nef(PolyPointList *, VertexNumList *, EqList *, PartList *,
11 int *, NEF_Flags *);
12
13 void IP_Fiber_Data(PolyPointList *, PolyPointList *,int nv,
14 Long G[VERT_Nmax][POLY_Dmax][POLY_Dmax],int fd[VERT_Nmax],
15 int *nf,int CD);
16
17 /* =============== End of Typedefs and Headers =================== */
18
19 /* =============== Begin of DYNamical Complete =================== */
20
Print_DYN_PPL(DYN_PPL * _MP,const char * comment)21 void Print_DYN_PPL( DYN_PPL *_MP, const char *comment){
22 int i,j;
23 if(_MP->np > 20){
24 fprintf(outFILE,"%d %d %s\n", (int) _MP->np,_MP->n, comment);
25 for(i = 0; i < _MP->np; i++) {
26 for(j = 0; j < _MP->n; j++)
27 fprintf(outFILE,"%d ",(int) _MP->L[i].x[j]);
28 fprintf(outFILE,"\n");
29 }
30 }
31 else {
32 fprintf(outFILE,"%d %d %s\n",_MP->n, (int)_MP->np, comment);
33 for(i = 0;i < _MP->n; i++) {
34 for(j = 0; j < _MP->np; j++)
35 fprintf(outFILE," %4d",(int) _MP->L[j].x[i]);
36 fprintf(outFILE,"\n");
37 }
38 }
39 }
40
DYNadd_for_completion(Long * yDen,Long Den,EqList * _E,DYN_PPL * _CP)41 void DYNadd_for_completion(Long *yDen, Long Den, EqList *_E, DYN_PPL *_CP){
42 int i,n=_CP->n;
43 Long yold[POLY_Dmax];
44
45 if(Den>1) for(i=0;i<n;i++) {
46 if(yDen[i]%Den) return;
47 yold[i]=yDen[i]/Den;}
48 else for(i=0;i<n;i++) yold[i]=yDen[i];
49 for (i=0;i<_E->ne;i++)
50 if (Eval_Eq_on_V(&(_E->e[i]), yold, n) < 0) return;
51 if(!(_CP->np < _CP->NP_max)){
52 _CP->NP_max += 1000000;
53 if((_CP->L=(Vector *)realloc(_CP->L,_CP->NP_max * sizeof(Vector))) == NULL)
54 Die("Unable to realloc space for _CP->L");
55 }
56 for(i=0;i<n;i++) _CP->L[_CP->np].x[i]=yold[i];
57 _CP->np++;
58 }
59
DYNComplete_Poly(Long VPM[][VERT_Nmax],EqList * _E,int nv,DYN_PPL * _CP)60 void DYNComplete_Poly(Long VPM[][VERT_Nmax], EqList *_E, int nv, DYN_PPL *_CP){
61 int i,j,k,l,InsPoint,rank=0,n=_CP->n;
62 Long MaxDist[EQUA_Nmax], InvMat[POLY_Dmax][POLY_Dmax], Den=1;
63 Long yDen[POLY_Dmax];
64 int OrdFac[VERT_Nmax],
65 BasFac[POLY_Dmax], one[POLY_Dmax], position[POLY_Dmax];
66 LRat ind[POLY_Dmax][POLY_Dmax], x[POLY_Dmax], y[POLY_Dmax], f,
67 PInvMat[POLY_Dmax][POLY_Dmax];
68
69 _CP->np=0;
70
71 /* Calculate maximal distances from facets of Delta^* (Vertices of Delta) */
72
73 for (i=0;i<_E->ne;i++) {
74 MaxDist[i]=0;
75 for (j=0;j<nv;j++)
76 if (MaxDist[i]<VPM[i][j]) MaxDist[i]=VPM[i][j];}
77
78 /* Order facets of Delta^* (Vertices of Delta) w.r.t. MaxDist */
79
80 OrdFac[0]=0;
81 for (i=1;i<_E->ne;i++){
82 InsPoint=i;
83 while (InsPoint&&(MaxDist[i]<MaxDist[OrdFac[InsPoint-1]])) InsPoint--;
84 for (j=i;j>InsPoint;j--) OrdFac[j]=OrdFac[j-1];
85 OrdFac[InsPoint]=i; }
86
87 /* Find first POLY_Dmax linearly independent facets + Inverse Matrix */
88
89 for (i=0;i<n;i++) for (j=0;j<n;j++) PInvMat[i][j]=LrI(0);
90 for (i=0;i<n;i++) PInvMat[i][i]=LrI(1);
91 i=0;
92 while (rank<n){
93 for (j=0;j<n;j++) x[j]=LrI(_E->e[OrdFac[i]].a[j]);
94 for (j=0;j<n;j++) y[j]=LrI(0);
95 y[rank]=LrI(1);
96 for (j=0;j<rank;j++) {
97 f=x[one[j]];
98 for (k=0;k<n;k++) {
99 x[k]=LrD(x[k],LrP(f,ind[j][k]));
100 y[k]=LrD(y[k],LrP(f,PInvMat[j][k])); } }
101 one[rank]=-1;
102 for (l=0;(l<n)&&(one[rank]==-1);l++) if (x[l].N) one[rank]=l;
103 if(one[rank]>-1){
104 for (k=0;k<n;k++) {
105 ind[rank][k]=LrQ(x[k],x[one[rank]]);
106 PInvMat[rank][k]=LrQ(y[k],x[one[rank]]); }
107 for (j=0;j<rank;j++) {
108 f=ind[j][one[rank]];
109 for (k=0;k<n;k++) {
110 ind[j][k]=LrD(ind[j][k],LrP(ind[rank][k],f));
111 PInvMat[j][k]=LrD(PInvMat[j][k],LrP(PInvMat[rank][k],f)); } }
112 BasFac[rank]=OrdFac[i];
113 rank++; }
114 i++; }
115 for (i=0;i<n;i++) for (j=0;j<n;j++)
116 Den=(Den/LFgcd(Den,PInvMat[i][j].D))*PInvMat[i][j].D;
117 for (i=0;i<n;i++) for (j=0;j<n;j++)
118 InvMat[one[i]][j]=(Den/PInvMat[i][j].D)*PInvMat[i][j].N;
119
120 for (i=0;i<n;i++){
121 for (j=0;j<n;j++) {
122 long long s=0;
123 for(k=0;k<n;k++) s+=((long long) (InvMat[k][i]))*
124 ((long long) (_E->e[BasFac[j]].a[k]));
125 if (s!=Den*(i==j)) {
126 puts("something wrong in Make_Dual_Poly");
127 exit(0);}}}
128
129 /* Examine all integer points of parallelogram: */
130 /* The basic structure of the algorithm is:
131 for (k=0;k<n-1;k++) position[k]=-1; / * sets k=n-1; important! *
132 position[n-1]=-2; / * starting point just outside the parallelogram *
133 while(k>=0){
134 position[k]++;
135 DO AT position;
136 for(k=n-1;((position[k]==MaxDist[BasFac[k]]-1)&&(k>=0));k--)
137 position[k]=-1; }
138 / * sets k to the highest value where pos.[k] wasn't the max value;
139 resets the following max values to min values */
140 /* Quantities linear in position can be changed with every change of
141 position (here: yDen) */
142
143 for(i=0;i<n;i++) yDen[i]=0;
144 for (k=0;k<n-1;k++) { /* sets k=n-1; important! */
145 position[k]=-_E->e[BasFac[k]].c;
146 for(i=0;i<n;i++) yDen[i]-=_E->e[BasFac[k]].c*InvMat[i][k]; }
147 position[n-1]=-_E->e[BasFac[n-1]].c-1;
148 for(i=0;i<n;i++) yDen[i]-=(_E->e[BasFac[k]].c+1)*InvMat[i][n-1];
149 while(k>=0){
150 position[k]++;
151 for(i=0;i<n;i++) yDen[i]+=InvMat[i][k];
152 DYNadd_for_completion(yDen, Den, _E, _CP);
153 for(k=n-1;(k>=0);k--){
154 if (position[k]!=MaxDist[BasFac[k]]-_E->e[BasFac[k]].c) break;
155 position[k]=-_E->e[BasFac[k]].c;
156 for (i=0;i<n;i++) yDen[i]-=MaxDist[BasFac[k]]*InvMat[i][k]; }}
157 }
158
DYNMake_VEPM(DYN_PPL * _P,VertexNumList * _V,EqList * _E,Long PM[][VERT_Nmax])159 void DYNMake_VEPM(DYN_PPL *_P, VertexNumList *_V, EqList *_E,
160 Long PM[][VERT_Nmax]){
161 int i, j;
162 for (i=0;i<_E->ne;i++) for (j=0;j<_V->nv;j++)
163 PM[i][j]=Eval_Eq_on_V(&_E->e[i],_P->L[_V->v[j]].x,_P->n);
164 }
165
166 /* =============== End of DYNamical Complete =================== */
167 /* =============== Begin of FIBRATIONS =================== */
PRINT_APL(AmbiPointList * _AP,const char * comment)168 void PRINT_APL(AmbiPointList *_AP, const char *comment){
169 int i,j;
170
171 fprintf(outFILE,"%d %d %s\n", _AP->N, _AP->np, comment);
172 for(i = 0; i < _AP->N; i++){
173 for(j = 0; j < _AP->np; j++)
174 fprintf(outFILE,(_AP->np>20) ? " %2d" : " %4d", (int) _AP->x[j][i]);
175 fprintf(outFILE,"\n");
176 }
177 }
178
PRINT_MATRIX(Long * M,int l,int c,int C)179 void PRINT_MATRIX(Long *M, int l, int c, int C){
180 int i,j;
181 for(j = 0; j < l; j++){
182 for(i = 0; i < c; i++)
183 fprintf(outFILE,(c>20) ? " %3d" : " %4d", (int) *(M+i+C*j));
184 fprintf(outFILE,"\n");
185 }
186 }
187
PRINT_TMATRIX(Long * M,int l,int c,int C)188 void PRINT_TMATRIX(Long *M, int l, int c, int C){
189 int i,j;
190
191 for(i = 0; i < c; i++){
192 for(j = 0; j < l; j++)
193 fprintf(outFILE,(l>20) ? " %3d" : " %4d", (int) *(M+i+C*j));
194 fprintf(outFILE,"\n");
195 }
196 }
197
PRINT_PL(PolyPointList * _P,const char * comment)198 void PRINT_PL(PolyPointList *_P, const char *comment){
199
200 fprintf(outFILE,"%d %d %s\n", _P->n, _P->np, comment);
201 PRINT_TMATRIX(&_P->x[0][0], _P->np, _P->n, POLY_Dmax);
202 }
203
PRINT_GORE(PolyPointList * _P,int codim,int n,const char * comment)204 void PRINT_GORE(PolyPointList *_P, int codim, int n, const char *comment){
205
206 PolyPointList *_P_AUX;
207 VertexNumList *_V_AUX;
208 EqList *_E_AUX;
209 int i, j, Z=0;
210
211 _P_AUX = (PolyPointList *) malloc(sizeof(PolyPointList));
212 if (_P_AUX == NULL) Die("Unable to alloc space for PolyPointLis _P_AUX");
213 _V_AUX = (VertexNumList *) malloc(sizeof(VertexNumList));
214 if (_V_AUX == NULL) Die("Unable to alloc space for VertexNumList _V_AUX");
215 _E_AUX = (EqList *) malloc(sizeof(EqList));
216 if (_E_AUX == NULL) Die("Unable to alloc space for EqList _E_AUX");
217
218 _P_AUX->n = _P->n - codim + 1;
219
220 /* troubles if IP==0 not in last position:
221 _P_AUX->np = _P->np - codim + 1;
222 for(i = 0; i < _P_AUX->np; i++)
223 for(j = 0; j < _P_AUX->n; j++)
224 _P_AUX->x[i][j] = _P->x[i][j + codim - 1]; ** replaced by: */
225 _P_AUX->np = _P->np - codim;
226 for(i = 0; i < _P->np; i++){ int z=0; j=0;
227 while((j < _P_AUX->n)&&(!z)) {if(_P->x[i][j + codim - 1]) z=1; j++;};
228 if(z){ for(j = 0; j < _P_AUX->n; j++)
229 _P_AUX->x[i-Z][j] = _P->x[i][j + codim - 1]; }
230 else Z++;}
231 for(j = 0; j < _P_AUX->n; j++) _P_AUX->x[_P_AUX->np][j] = 0;
232 _P_AUX->np++;
233 /* Print_PPL(_P_AUX,"Ref_Check input"); */
234 assert(Ref_Check(_P_AUX, _V_AUX, _E_AUX));
235 /* Find_Equations(_P_AUX, _V_AUX, _E_AUX); ... redundant */
236 if(n == 0){
237 Sort_PPL(_P_AUX, _V_AUX);
238 fprintf(outFILE,"%d %d %s (nv=%d)\n",_P_AUX->n,_P_AUX->np,comment,_V_AUX->nv);
239 for(i = 0; i < _P_AUX->n; i++){
240 for(j = 0; j < _P_AUX->np; j++)
241 fprintf(outFILE,(_P_AUX->np>20) ? " %3d" : " %4d", (int) _P_AUX->x[j][i]);
242 fprintf(outFILE,"\n");
243 } /*PRINT_PL(_P_AUX, comment);*/
244 }
245 else{
246 Sort_VL(_V_AUX);
247 Sort_PPL(_P, _V_AUX);
248 if(n == 1)
249 fprintf(outFILE,"%d %d %s (nv=%d)\n",_P->n,_P->np,comment,_V_AUX->nv );
250 if(n == 2){
251 int o;
252 fprintf(outFILE,"%d %d %s (nv=%d)\n",_P->n+1,_P->np,comment,_V_AUX->nv );
253 for(j = 0; j < _P->np; j++){
254 o = 1; i = 0;
255 while(o && (i < codim - 1)){
256 if(_P->x[j][i] == 1)
257 o = 0;
258 i++;
259 }
260 fprintf(outFILE,(_P->np>20) ? " %3d" : " %4d", o);
261 }
262 fprintf(outFILE,"\n");
263 }
264 for(i = 0; i < _P->n; i++){
265 for(j = 0; j < _P->np; j++)
266 fprintf(outFILE,(_P->np>20) ? " %3d" : " %4d", (int) _P->x[j][i]);
267 fprintf(outFILE,"\n");
268 }
269 }
270 free(_P_AUX);free(_V_AUX);free(_E_AUX);
271 }
272
G_x_P(Long * Gi,Long * V,int * d)273 Long G_x_P(Long *Gi, Long *V, int *d){
274
275 Long x=0; int j;
276
277 for(j=0;j<*d;j++)
278 x+=Gi[j]*V[j];
279 return x;
280 }
281
PRINT_Fibrations(VertexNumList * _V,PolyPointList * _P,Flags * _F)282 void PRINT_Fibrations(VertexNumList *_V, PolyPointList *_P, Flags *_F
283 /* ,PartList *_PTL */ ){
284
285 Long G[VERT_Nmax][POLY_Dmax][POLY_Dmax];
286 int s[VERT_Nmax], CD=_F->f, n, c, i, j, fib, nf, nv, np, dim[VERT_Nmax];
287 PolyPointList *_P_AUX;
288 VertexNumList *_V_AUX;
289 EqList *_E_AUX;
290 char C[VERT_Nmax];
291
292 _P_AUX = (PolyPointList *) malloc(sizeof(PolyPointList));
293 if (_P_AUX == NULL) Die("Unable to alloc space for PolyPointLis _P_AUX");
294 _V_AUX = (VertexNumList *) malloc(sizeof(VertexNumList));
295 if (_V_AUX == NULL) Die("Unable to alloc space for VertexNumList _V_AUX");
296 _E_AUX = (EqList *) malloc(sizeof(EqList));
297 if (_E_AUX == NULL) Die("Unable to alloc space for EqList _E_AUX");
298
299 if (_P->np >= VERT_Nmax) Die("Need _P->np < VERT_Nmax in PRINT_Fibrations");
300 IP_Fiber_Data(_P, _P_AUX, _V->nv, G, dim, &nf, CD);
301 if (nf >= VERT_Nmax) Die("Need nf < VERT_Nmax in PRINT_Fibrations");
302
303
304
305 if(nf){
306 for(j = 0; j < _P->np - 1; j++)
307 fprintf(outFILE,(_P->np > 20) ? "----" :"-----");
308 fprintf(outFILE," #fibrations=%d\n",nf);
309 }
310 for(n = 0; n < nf; n++){
311 c = 0;
312 for(i = 0; i < (_P->np - 1) ; i++){
313 j = dim[n]; fib = 1;
314 while((j < _P->n) && fib){
315 if (G_x_P(G[n][j], _P->x[i], &_P->n))
316 fib = 0;
317 j++;
318 }
319 if(fib){
320 for(j = 0; j < dim[n]; j++)
321 _P_AUX->x[c][j] = G_x_P(G[n][j], _P->x[i], &_P->n);
322 s[c++] = i;
323 }
324 }
325 _P_AUX->np = c; _P_AUX->n = dim[n];
326 assert(Ref_Check(_P_AUX, _V_AUX, _E_AUX));
327 for(i = 0; i < (_P->np-1); i++) C[i]='_';
328 for(i = 0; i < c; i++)
329 C[s[i]] = 'p';
330 for(i = 0; i < _V_AUX->nv; i++)
331 C[s[_V_AUX->v[i]]] = 'v';
332 for(i = 0; i < _P->np-1; i++)
333 fprintf(outFILE,"%s%c", (_P->np > 20) ? " " : " ",C[i]);
334 nv = _V_AUX->nv; np = (_P_AUX->np + 1);
335 EL_to_PPL(_E_AUX, _P_AUX, &dim[n]);
336 assert(Ref_Check(_P_AUX, _V_AUX, _E_AUX));
337 {
338 Long X[VERT_Nmax][VERT_Nmax];
339 Make_VEPM(_P_AUX, _V_AUX, _E_AUX, X);
340 Complete_Poly(X, _E_AUX, _V_AUX->nv, _P_AUX);
341 /*Dim_Fib_CI(dim, n, _PTL, C);*/
342 fprintf(outFILE," cd=%d m:%3d %2d n:%2d %d\n",(_P->n - dim[n]),
343 _P_AUX->np, _V_AUX->nv, np, nv);
344 }
345
346 }
347 free(_P_AUX);free(_V_AUX);free(_E_AUX);
348 }
349
350 /* =============== End of FIBRATIONS =================== */
351
Die(char * comment)352 void Die(char *comment){
353 printf("\n%s\n",comment); exit(0);
354 }
355
Time_Info(time_t * _Tstart,clock_t * _Cstart,const char * comment)356 void Time_Info(time_t *_Tstart, clock_t *_Cstart, const char *comment){
357
358 fprintf(outFILE, "%s %dsec %dcpu\n", comment, (int)
359 /* CLOCKS_PER_SEC::10^6 */
360 difftime(time(NULL), *_Tstart), (int) ((((long long)
361 clock() - *_Cstart) / (long long) CLOCKS_PER_SEC))); fflush(0);
362 }
363
Print_Nefinfo(PartList * _PTL,time_t * _Tstart,clock_t * _Cstart)364 void Print_Nefinfo(PartList *_PTL, /* Flags *_F,*/ time_t *_Tstart,
365 clock_t *_Cstart){
366
367 int i, d = 0, p = 0;
368
369 for(i = 0; i < _PTL->n; i++){
370 if(_PTL->DirProduct[i]) d++;
371 }
372 for(i = 0; i < _PTL->n; i++){
373 if(_PTL->Proj[i]) p++;
374 }
375 fprintf(outFILE, "np=%d d:%d p:%d %4dsec %4dcpu\n", _PTL->n-d-p, d, p,
376 /* CLOCKS_PER_SEC::10^6 */ (int) difftime(time(NULL), *_Tstart), (int)
377 ((((long long) clock() - *_Cstart) / (long long) CLOCKS_PER_SEC)));
378 fflush(0);
379 }
380
N_Part(PartList * _PTL)381 int N_Part(PartList *_PTL){
382
383 int i;
384 for(i = 0; i < _PTL->n; i++)
385 if((!_PTL->DirProduct[i]) && (!_PTL->Proj[i]))
386 return 1;
387 return 0;
388 }
389
Print_L(LInfo * _L,int p,int v)390 void Print_L(LInfo *_L, int p, int v)
391 {
392 int codim, D, i, j, N;
393
394 if(v && p) Die("only -Lp OR -Lv !");
395
396 N = (v ? _L->nv : (_L->nv + 1));
397 if(v){
398 assert(FIB_POINT_Nmax >= N);
399 fprintf(outFILE,"%d %d Vertices in N-lattice:\n",_L->d, N);
400 for(j = 0; j < _L->d; j++){
401 for(i = 0; i < N; i++)
402 fprintf(outFILE,(_L->nv > 20) ? " %3d" : " %4d", (int) _L->VM[i][j]);
403 fprintf(outFILE,"\n");
404 }
405 }
406 for(j = 0; j < N; j++)
407 fprintf(outFILE,(_L->nv > 20) ? "----" :"-----");
408 fprintf(outFILE,"\n");
409 assert(FIB_Nmax >= _L->nw); assert(FIB_POINT_Nmax >= _L->nv);
410 for(i = 0; i < _L->nw; i++){
411 D = 0; codim = 0;
412 for(j = 0; j < _L->nv; j++){
413 D += _L->W[i][j];
414 if(_L->W[i][j] == 0)
415 codim++;
416 }
417 for(j = 0; j < _L->nv; j++)
418 fprintf(outFILE,(N > 20) ? " %3d" : " %4d",(int) _L->W[i][j]);
419 fprintf(outFILE," d=%d ",D);
420 fprintf(outFILE,"codim=%d\n",codim+_L->d-_L->nv+1);
421 }
422 }
423
Make_L(PolyPointList * _P,VertexNumList * _V,LInfo * _L,int p,int v)424 void Make_L(PolyPointList *_P, VertexNumList *_V, LInfo *_L , int p, int v){
425 int i,j;
426
427 if(v){
428 for(i = 0; i < _V->nv; i++)
429 for(j = 0; j < _P->n; j++)
430 _L->VM[i][j] = _P->x[_V->v[i]][j];
431 _L->nv = _V->nv;
432 }
433 if(p){
434 assert(_P->np <= VERT_Nmax);
435 for(i = 0; i < _P->np; i++)
436 for(j = 0; j < _P->n; j++)
437 _L->VM[i][j] = _P->x[i][j];
438 _L->nv = (_P->np-1);
439 }
440 _L->Wmax = FIB_Nmax; _L->nw = 0; _L->d = _P->n;
441 IP_Simplex_Decomp(_L->VM, _L->nv, _P->n, &_L->nw, _L->W, _L->Wmax,0);
442 }
443
IntSqrt(int q)444 int IntSqrt(int q)
445 { /* sqrt(q) => r=1; r'=(q+r*r)/(2r); */
446 if (q == 0)
447 return 0;
448 if (q < 4)
449 return 1;
450 else { /* troubles: e.g. 9408 */
451 long long r = (q + 1) / 2, n;
452 while (r > (n = (q + r * r) / (2 * r)))
453 r = n;
454 if (q < r * r)
455 r--;
456 if ((r * r <= q) && (q < (r + 1) * (r + 1)))
457 return (int) r;
458 else{
459 printf("Error in sqrt(%d)=%d\n", q, (int) n);
460 exit(0);
461 }
462 }
463 return 0;
464
465 }void INCI_TO(int I[], INCI * _X, int *_n)
466 {
467 /* INCI X -> (0,0,1,0,.....,0,1,0);*/
468 int i;
469
470 INCI Y = *_X;
471 for (i = 0; i < *_n; i++) {
472 if (INCI_M2(Y))
473 I[i] = 1;
474 else
475 I[i] = 0;
476 Y = INCI_D2(Y);
477 }
478 }
479
Num_Pos(Cone * _C)480 int Num_Pos(Cone *_C){
481
482 int d, n=2;
483
484 for (d = 1; d < _C->dim; d++)
485 n += _C->nface[d];
486 return n;
487 }
488
Make_PosetList(Cone * _C,Poset_Element_List * _PEL)489 void Make_PosetList(Cone *_C, Poset_Element_List *_PEL){
490
491 int n = 0, d, i;
492
493 for (d = 0; d <= _C->dim; d++)
494 for (i = 0; i < _C->nface[d]; i++) {
495 _PEL->L[n].dim = d;
496 _PEL->L[n].num = i;
497 n++;
498 }
499 }
500
Interval_Check(int rank,Poset_Element * _x,Poset_Element * _y,Cone * _C)501 int Interval_Check(int rank, Poset_Element * _x, Poset_Element * _y,
502 Cone * _C)
503 {
504 int flag = 0;
505
506 if (_x->dim == (_y->dim - rank))
507 if (INCI_LE(_C->edge[_x->dim][_x->num], _C->edge[_y->dim][_y->num]))
508 flag = 1;
509 return flag;
510 }
511
Make_Intervallist(Interval_List * _IL,Poset_Element_List * _PEL,Cone * _C)512 void Make_Intervallist(Interval_List *_IL, Poset_Element_List *_PEL, Cone *_C){
513
514 int d, i, j;
515
516 _IL->n = 0;
517 for (d = 0; d <= _C->dim; d++)
518 for (i = 0; i < _PEL->n; i++)
519 for (j = 0; j <= i; j++)
520 if (Interval_Check(d, &_PEL->L[j], &_PEL->L[i], _C) == 1) {
521 _IL->L[_IL->n].min = j;
522 _IL->L[_IL->n].max = i;
523 _IL->n++;
524 }
525 }
526
Make_Mirror(EPoly * _EP,int h[][POLY_Dmax],int D,int dim)527 int Make_Mirror(EPoly *_EP, int h[][POLY_Dmax], int D, int dim)
528 {
529 int i, j, k, u, v, chi = 0, H;
530
531 for (u = 0; u < 4*(Pos_Max); u++)
532 for (v = 0; v < 3*(Pos_Max); v++){
533 if(((-2 * D + u) > dim) || ((-2 * D + u) < 0) ||
534 ((-2 * D + v) > dim) || ((-2 * D + v) < 0)){
535 if (_EP->E[u][v] != 0)
536 Die("Something wrong with E poly");
537 }
538 else {
539 h[dim - (u - 2*D)][v - 2*D] = _EP->E[u][v];
540 chi += _EP->E[u][v];
541 if(((u - 2*D + v - 2*D) %2 ) != 0)
542 h[dim - (u - 2*D)][v - 2*D] *= -1;
543 }
544 }
545 if((dim %2) != 0) chi = chi*(-1);
546 H = - dim*chi;
547 for (i = 0; i <= dim; i++)
548 for (j = 0; j <= dim; j++){
549 if(((i+j) %2) == 0)
550 k = 1;
551 else
552 k = -1;
553 H += 3*k*(2*i - dim)*(2*i - dim)*h[i][j];
554 }
555 if (H != 0) Die("paper: CY 4-folds and toric fibrations; equation (8) don't hold");
556 if (dim == 4){
557 if((-h[2][2] + 44*h[0][0] + 4 * h[1][1] - 2 * h[1][2] + 4 * h[1][3] +
558 20 * h[0][2] - 52 * h[0][1]) != 0)
559 Die("paper: CY 4-folds and toric fibrations; equation (9) don't hold");
560 }
561 return chi;
562 }
563
Max_S(PartList * _PTL,int * _n)564 int Max_S(PartList * _PTL, int *_n)
565 {
566 int i, k, m = 0, nv, Nv = 0;
567 for (i = 0; i < _PTL->codim; i++) {
568 nv = 0;
569 for (k = 0; k < _PTL->nv; k++)
570 if (_PTL->S[*_n][k] == i)
571 nv++;
572 if (Nv < nv) {
573 Nv = nv;
574 m = i;
575 }
576 }
577 return m;
578 }
579
PrintDegrees(LInfo * _L,PartList * _PTL,int m,int S[VERT_Nmax])580 void PrintDegrees( /*Flags *_F,*/ LInfo *_L, PartList *_PTL, int m,
581 /*int n, */ int S[VERT_Nmax])
582 {
583 int d[POLY_Dmax], i, j;
584
585 for (i = 0; i < _L->nw; i++){
586 fprintf(outFILE," (");
587 for (j = 0; j < _PTL->codim; j++)
588 d[j] = 0;
589 for (j = 0; j < _L->nv; j++)
590 d[S[j]] += _L->W[i][j];
591 for (j = 0; j < _PTL->codim; j++)
592 if(j != m)
593 fprintf(outFILE,"%d ",d[j]);
594 fprintf(outFILE,"%d",d[m]);
595 fprintf(outFILE,")");
596 }
597 }
598
PrintWeights(CWS * _W)599 void PrintWeights(CWS * _W)
600 {
601 int i,j;
602 for (i = 0; i < _W->nw; i++) {
603 fprintf(outFILE, "%ld", (long) _W->d[i]);
604 for (j = 0; j < _W->N; j++)
605 fprintf(outFILE, " %ld", (long) _W->W[i][j]);
606 if (i != (_W->nw - 1))
607 fprintf(outFILE, " ");
608 }
609 }
610
PrintDiamond(int h[][POLY_Dmax],int dim)611 void PrintDiamond(int h[][POLY_Dmax], int dim)
612 {
613 int i,j;
614
615 fprintf(outFILE, "\n\n"); fflush(0);
616 for (i = 0; i <= dim; i++) {
617 fprintf(outFILE, " ");
618 for (j = 0; j <= (dim - i); j++)
619 fprintf(outFILE, " ");
620 for (j = 0; j <= i; j++)
621 fprintf(outFILE, " h%2d%2d ", i - j, j);
622 fprintf(outFILE, "\n\n");
623 fflush(0);
624 }
625 for (i = 1; i <= dim; i++) {
626 fprintf(outFILE, " ");
627 for (j = 0; j <= i; j++)
628 fprintf(outFILE, " ");
629 for (j = i; j <= dim; j++)
630 fprintf(outFILE, " h%2d%2d ", dim - j + i, j);
631 fprintf(outFILE, "\n\n");
632 fflush(0);
633 }
634 fprintf(outFILE, "\n\n");fflush(0);
635 for (i = 0; i <= dim; i++) {
636 fprintf(outFILE, " ");
637 for (j = 0; j <= (dim - i); j++)
638 fprintf(outFILE, " ");
639 for (j = 0; j <= i; j++)
640 fprintf(outFILE, "%10d", h[i - j][j]);
641 fprintf(outFILE, "\n\n");
642 }
643 for (i = 1; i <= dim; i++) {
644 fprintf(outFILE, " ");
645 for (j = 0; j <= i; j++)
646 fprintf(outFILE, " ");
647 for (j = i; j <= dim; j++)
648 fprintf(outFILE, "%10d", h[dim - j + i][j]);
649 fprintf(outFILE, "\n\n");
650 }
651 fflush(0);
652 }
653
Print_Points(PolyPointList * _P,int c,int nv,int S[VERT_Nmax])654 void Print_Points(PolyPointList *_P, int c, int nv, int S[VERT_Nmax]){
655
656 int i, P=0;
657
658 for(i = nv; i < (_P->np - 1); i++)
659 if(S[i] == c){
660 fprintf(outFILE,"%d ",i);
661 P=1;
662 }
663 if(P)
664 fprintf(outFILE," ");
665 }
666
Output(PolyPointList * _P,PolyPointList * _P_D,CWS * _W,EPoly * _EP,VertexNumList * _V,int * _n,PartList * _PTL,int * _codim,FILE * outFILE,Flags * _F,int * _D,LInfo * _L)667 void Output(PolyPointList * _P, /* PolyPointList * _DP,*/
668 PolyPointList * _P_D, CWS * _W, EPoly *_EP, /* EqList * _E,*/
669 VertexNumList * _V, int *_n, PartList *_PTL,
670 int *_codim, FILE *outFILE, Flags * _F, int *_D, LInfo *_L)
671 {
672 int i, j, k, m, chi = 0, D = (_P_D->n + 1), dim = (_P->n - *_codim);
673 int h[POLY_Dmax][POLY_Dmax] = {{0}, {0}}, S[VERT_Nmax];
674
675
676 m = Max_S(_PTL, _n);
677 if (((!_PTL->DirProduct[*_n]) || _F->Dir) && ((!_PTL->Proj[*_n]) || _F->Proj)) {
678 if (_F->H == 0){
679 if (_F->w){
680 PrintWeights(_W); fprintf(outFILE, " ");
681 }
682 #ifdef WRITE_CWS
683 if (_F->Msum == 1)
684 fprintf(outFILE, " d=%d %d", (int) _D[0], (int) _D[1]);
685 #endif
686 if (!_F->p){
687 chi = Make_Mirror(_EP, h, D, dim);
688 fprintf(outFILE, "H:");
689 for (i = 1; i < dim; i++)
690 fprintf(outFILE, "%d ", h[1][i]);
691 fprintf(outFILE, "[%d]", chi);
692 for (i = 1; i <= dim/2; i++)
693 if (h[0][i] != 0)
694 fprintf(outFILE, " h%d=%d", i, h[0][i]);
695 if(h[0][0] != 1){
696 if (_PTL->DirProduct[*_n])
697 fprintf(outFILE, " h%d=%d", 0, h[0][0]);
698 else
699 Die("\nh00 not 1 !!!\n");
700 }
701 }
702 if((_P_D->np - *_codim) <= VERT_Nmax)
703 for(i = 0; i < (_P_D->np - *_codim); i++){
704 S[i] = 0;
705 for(j = 0; j < (*_codim - 1); j++)
706 if(_P_D->x[i][j])
707 S[i] = (j+1);
708 }
709 if (*_codim == 2) {
710 fprintf(outFILE, " P:%d V:", *_n);
711 i = 0;
712 if (m == 0)
713 i = 1;
714 for (j = 0; j < _PTL->nv; j++)
715 if (_PTL->S[*_n][j] == i)
716 fprintf(outFILE, "%d ", j);
717 fprintf(outFILE," ");
718 if((_P_D->np - *_codim) <= VERT_Nmax)
719 Print_Points(_P, i, _V->nv, S);
720 else
721 fprintf(outFILE, " _P->np > VERT_Nmax! ");
722 }
723 else {
724 fprintf(outFILE, " P:%d ", *_n);
725 j = 0;
726 for (i = 0; i < *_codim; i++)
727 if (i != m) {
728 fprintf(outFILE, "V%d:", j);
729 for (k = 0; k < _PTL->nv; k++)
730 if (_PTL->S[*_n][k] == i)
731 fprintf(outFILE, "%d ", k);
732 fprintf(outFILE," ");
733 j++;
734 if((_P_D->np - *_codim) <= VERT_Nmax)
735 Print_Points(_P, i, _V->nv, S);
736 else
737 fprintf(outFILE, " _P->np > VERT_Nmax! ");
738 }
739 }
740 if(_PTL->DProj[*_n])
741 fprintf(outFILE, "DP ");
742 if(_PTL->DirProduct[*_n])
743 fprintf(outFILE, " D");
744 if(_F->Lv || (_F->Lp && (_P_D->np - *_codim) <= VERT_Nmax))
745 PrintDegrees(/*_F,*/ _L, _PTL, m, /* *_n,*/ S);
746 fflush(0);
747 }
748 else{
749 chi = Make_Mirror(_EP, h, D, dim);
750 PrintDiamond(h, dim);
751 }
752 }
753 }
754
755
Min_Dim(int n,int T_flag)756 int Min_Dim(int n, int T_flag)
757 {
758 if (!T_flag){
759 int i;
760 i = n / 2;
761 if((n % 2) != 0)
762 i += 1;
763 return i;
764 } else
765 return n;
766 }
767
Init_ST(SPoly * _S,SPoly * _T,Poset_Element_List * _PEL)768 void Init_ST(SPoly *_S, SPoly *_T, Poset_Element_List *_PEL)
769 {
770 int i, d;
771
772 for(i = 0; i < _PEL->n; i++){
773 for(d = 1; d <= _PEL->L[i].dim; d++){
774 _S[i].S[d] = 0;
775 _T[i].S[d] = 0;
776 }
777 _S[i].S[0] = 1; _T[i].S[0] = 0;
778 }
779 }
780
Eval_Eq_on_x(Long * _x,Equation * _E,int dim)781 Long Eval_Eq_on_x(Long *_x, Equation *_E, int dim){
782
783 Long c = _E->c;
784 int d;
785
786 for(d = 0; d < dim; d++)
787 c += _x[d] * _E->a[d];
788 return c;
789 }
790
INCI_to_x(int n,DYN_PPL * _P,EqList * _E)791 INCI INCI_to_x(int n, DYN_PPL *_P, EqList *_E){
792
793 INCI X = INCI_0();
794 int i;
795
796 for(i = 0; i < _E->ne; i++)
797 X = INCI_PN(X,Eval_Eq_on_x(_P->L[n].x, &_E->e[i], _P->n));
798 return X;
799 }
800
Poly_To_ST(DYN_PPL * _P,EqList * _E,Cone * _C,SPoly * _S,SPoly * _T,Poset_Element_List * _PEL,int l,int T_flag)801 void Poly_To_ST(DYN_PPL *_P, EqList *_E, Cone *_C, SPoly *_S, SPoly *_T,
802 Poset_Element_List *_PEL, int l, int T_flag)
803
804 {
805 int min, i, dim_flag, in_flag, n, dim, num;
806 INCI X;
807
808 min = Min_Dim(l, T_flag);
809 for(i = 0; i < _P->np; i++){
810 dim_flag = 1; in_flag = 0;
811 n = (_PEL->n - 1);
812 X = INCI_to_x(i, _P, _E);
813 while(dim_flag && !in_flag){
814 dim = (_C->dim - _PEL->L[n].dim);
815 num = (_C->nface[dim] - _PEL->L[n].num -1);
816 if (INCI_LE(_C->edge[dim][num], X)){
817 _S[n].S[l] += 1;
818 if (INCI_LE(X, _C->edge[dim][num])){
819 _T[n].S[l] += 1;
820 in_flag = 1;
821 }
822 }
823 n--;
824 if (_PEL->L[n].dim < min)
825 dim_flag = 0;
826 }
827 }
828 }
829
New_CPVE(PolyPointList * _P,DYN_PPL * _CP,VertexNumList * _V,VertexNumList * _CV,EqList * _E,EqList * _CE,int n)830 void New_CPVE(PolyPointList *_P, DYN_PPL *_CP, VertexNumList *_V,
831 VertexNumList *_CV, EqList *_E, EqList *_CE, int n)
832 {
833 int j, d, nv = 0;
834
835 for(j = 0; j < _E->ne; j++){
836 _CE->e[j].c = n*(_E->e[j].c);
837 for(d = 0; d < _P->n; d++)
838 _CE->e[j].a[d] = _E->e[j].a[d];
839 }
840 for(j = 0; j < _V->nv; j++){
841 for(d = 0; d < _P->n; d++)
842 _CP->L[nv].x[d] = n*(_P->x[_V->v[j]][d]);
843 _CV->v[nv] = nv;
844 nv++;
845 }
846 _CV->nv = _V->nv; _CP->np = _V->nv; _CP->n = _P->n; _CE->ne = _E->ne;
847 }
848
Poly_To_DYNPoly(DYN_PPL * _CP,PolyPointList * _P)849 void Poly_To_DYNPoly(DYN_PPL *_CP, PolyPointList *_P){
850 int i, d;
851
852 assert(_P->np <= _CP->NP_max);
853 for(i = 0; i < _P->np; i++)
854 for(d = 0; d < _P->n; d++)
855 _CP->L[i].x[d] = _P->x[i][d];
856 _CP->n = _P->n; _CP->np = _P->np;
857 }
858
Make_S_Poly(Cone * _C,VertexNumList * _V,EqList * _E,PolyPointList * _P,Poset_Element_List * _PEL,SPoly * _S,int SINFO,int CHECK_SERRE)859 void Make_S_Poly(Cone *_C, VertexNumList *_V, EqList *_E, PolyPointList *_P,
860 Poset_Element_List *_PEL, SPoly *_S, int SINFO, int CHECK_SERRE)
861 {
862 Long PM[EQUA_Nmax][VERT_Nmax];
863 DYN_PPL CP;
864 VertexNumList *_CV;
865 EqList *_CE;
866 SPoly *_T;
867 int i=1, j, d, min;
868
869 CP.NP_max = NP_Max;
870 min = Min_Dim(_C->dim,CHECK_SERRE);
871
872 CP.L = (Vector *) calloc(CP.NP_max, sizeof(Vector));
873 if (CP.L == NULL) Die("Unable to alloc space for PolyPointList _CP.L");
874 _CV = (VertexNumList *) malloc(sizeof(VertexNumList));
875 if (_CV == NULL) Die("Unable to alloc space for VertexNumList _CV");
876 _CE = (EqList *) malloc(sizeof(EqList));
877 if (_CE == NULL) Die("Unable to alloc space for EqList _CE");
878 _T = (SPoly *) calloc(_PEL->n, sizeof(SPoly ));
879 if (_T == NULL) Die("Unable to alloc space for SPoly _T");
880
881 Poly_To_DYNPoly(&CP, _P);
882 Init_ST(_S, _T, _PEL);
883 Poly_To_ST(&CP, _E, _C, _S, _T, _PEL, i, CHECK_SERRE);
884
885 if(SINFO){
886 printf("\n\n#points in largest cone:\n");
887 printf("layer: %2d #p: %8d #ip: %8d\n", i=1, (int) _S[_PEL->n-1].S[1],
888 (int)_T[_PEL->n-1].S[1]);
889 }
890 for (i = 2; i <= min; i++) {
891 New_CPVE(_P, &CP, _V, _CV, _E, _CE, i);
892 DYNMake_VEPM(&CP,_CV,_CE, PM);
893 DYNComplete_Poly(PM, _CE, _V->nv, &CP);
894 Poly_To_ST(&CP, _CE, _C, _S, _T, _PEL, i, CHECK_SERRE);
895 if(SINFO)
896 printf("layer: %2d #p: %8d #ip: %8d\n",i, (int)_S[_PEL->n-1].S[i],
897 (int) _T[_PEL->n-1].S[i]);
898 }
899 for(i = 0; i < _PEL->n; i++){
900 for (j = 0; j < _PEL->L[i].dim; j++){
901 min = Min_Dim(_PEL->L[i].dim, CHECK_SERRE);
902 for (d = min; d > 0; d--) {
903 _S[i].S[d] += -_S[i].S[d - 1];
904 _T[i].S[d] += -_T[i].S[d - 1];
905 }
906 }
907 if(CHECK_SERRE){
908 for (d = 1; d < _PEL->L[i].dim; d++)
909 if(_S[i].S[d] != _T[i].S[_PEL->L[i].dim - d]){
910 puts("No Serre Duality!");
911 for(d = 1; d < _PEL->L[i].dim; d++)
912 printf("S[%d]: %3d T[%d]: %3d\n",d, (int) _S[i].S[d],
913 (_PEL->L[i].dim - d), (int) _T[i].S[_PEL->L[i].dim - d]);
914 exit(0);
915 }
916 }
917 if (!CHECK_SERRE){
918 for (d = 1; d < _PEL->L[i].dim / 2; d++)
919 _S[i].S[_PEL->L[i].dim - d] = _T[i].S[d];
920 if (_PEL->L[i].dim != 0)
921 _S[i].S[_PEL->L[i].dim] = 0;
922 }
923 }
924 free(_T);free(_CV);free(_CE);free(CP.L);
925 }
926
SB_To_E(EPoly * _EP,Cone * _C,Poset_Element_List * _PEL,BPoly * _BL,Interval_List * _IL,SPoly * _S_D,SPoly * _S_N,int * _codim)927 void SB_To_E(EPoly *_EP, Cone *_C, Poset_Element_List *_PEL, BPoly * _BL,
928 Interval_List * _IL, SPoly *_S_D, SPoly *_S_N, int *_codim){
929 /* E[0][0] = E[2*_C->dim][2*_C->dim] */
930
931 int i, j, x, y, u, v, s, dim, dmax, dmin;
932
933 for (u = 0; u < 4*(Pos_Max); u++)
934 for (v = 0; v < 4*(Pos_Max); v++)
935 _EP->E[u][v] = 0;
936
937 for (i = 0; i < _IL->n; i++) {
938 dmin = _PEL->L[_IL->L[i].min].dim;
939 dmax = _PEL->L[_IL->L[i].max].dim;
940 dim = dmax - dmin;
941 s = 1;
942 for (j = 0; j < dmin; j++)
943 s *= -1;
944 for (y = 0; y < (_C->dim - dmax + 1); y++)
945 for (x = 0; x < (dmin + 1); x++)
946 for (u = 0; u < (dim + 1); u++)
947 for (v = 0;v <= (dim / 2); v++)
948 _EP->E[2*_C->dim - x + y + dmax - u - *_codim]
949 [2*_C->dim + y + x + v - *_codim] +=
950 _S_D[_IL->L[i].min].S[x] * _S_N[_PEL->n - _IL->L[i].max - 1].S[y]
951 * _BL[i].B[u][v] * s;
952 }
953 }
954
Make_Cone(Cone * _C_D,Cone * _C_N,FaceInfo * _I,PolyPointList * _P)955 void Make_Cone(Cone * _C_D, Cone * _C_N, FaceInfo * _I, PolyPointList * _P)
956 {
957 int d, i;
958
959 _C_D->dim = (_P->n + 1); _C_N->dim = (_P->n + 1);
960 _C_D->nface[0] = 1; _C_N->nface[0] = 1;
961 _C_D->nface[_P->n + 1] = 1; _C_N->nface[_P->n + 1] = 1;
962 _C_D->edge[0][0] = INCI_0(); _C_N->edge[0][0] = INCI_0();
963 _C_D->edge[_C_D->dim][0] = INCI_0();
964 for (i = 0; i < _I->nf[0]; i++)
965 _C_D->edge[_C_D->dim][0] =
966 INCI_OR(_C_D->edge[_C_D->dim][0], _I->v[0][i]);
967 _C_N->edge[_C_N->dim][0] = INCI_0();
968 for (i = 0; i < _I->nf[_P->n - 1]; i++)
969 _C_N->edge[_C_N->dim][0] =
970 INCI_OR(_C_N->edge[_C_N->dim][0], _I->f[_P->n - 1][i]);
971 for (d = 0; d < _P->n; d++) {
972 _C_D->nface[d + 1] = _I->nf[d];
973 _C_N->nface[_C_D->dim - d -1] = _I->nf[d];
974 for (i = 0; i < _I->nf[d]; i++) {
975 _C_D->edge[d + 1][i] = _I->v[d][i];
976 _C_N->edge[_C_D->dim - d -1][_I->nf[d] - i -1] = _I->f[d][i];
977 }
978 }
979 }
980
Make_EN(PolyPointList * _P,VertexNumList * _V,EqList * _EN,int * _codim)981 void Make_EN(PolyPointList * _P, VertexNumList * _V, EqList * _EN, int *_codim)
982 {
983 int x, i, j;
984
985 _EN->ne = _V->nv;
986 for (i = 0; i < _V->nv; i++) {
987 x = 1; j = 0;
988 while((x == 1) && (j < (*_codim - 1))){
989 if(_P->x[_V->v[i]][j] == 1)
990 x = 0;
991 j++;
992 }
993 _EN->e[i].c = x;
994 for (j = 0; j < _P->n; j++){
995 if(j < *_codim - 1)
996 _EN->e[i].a[j] = _P->x[_V->v[i]][j] - x;
997 else
998 _EN->e[i].a[j] = _P->x[_V->v[i]][j];
999 }
1000 }
1001 }
1002
Remove_Proj(PolyPointList * _P,int * _codim)1003 int Remove_Proj(PolyPointList * _P, /* int *_n,*/ int *_codim)
1004 {
1005 int nv, Nv=0, j, i = 0, proj_flag = 0;
1006
1007 while (!proj_flag && (i < *_codim - 1)) {
1008 nv=0;
1009 for(j = 0; j < _P->np; j++)
1010 if(_P->x[j][i] == 1)
1011 nv++;
1012 Nv += nv;
1013 if(nv == 2)
1014 proj_flag = 1;
1015 i++;
1016 }
1017 if((_P->np - Nv) == 2)
1018 proj_flag = 1;
1019 return proj_flag;
1020 }
1021
Make_Gore_Poly(PolyPointList * _P,PolyPointList * _DP,PolyPointList * _P_D,PolyPointList * _P_N,VertexNumList * _V,PartList * _PTL,int * _codim,int * _n)1022 void Make_Gore_Poly(PolyPointList * _P, PolyPointList * _DP,
1023 PolyPointList * _P_D, PolyPointList * _P_N,
1024 VertexNumList * _V, PartList * _PTL, int *_codim,
1025 int *_n)
1026 {
1027 /*_P_N from _DP (in M-lattice), _P_D from _P (in N-lattice) */
1028
1029 int i, l, k, d, sum, ip, c;
1030
1031 if(*_codim <= 0)
1032 Die("Need Codim > 0");
1033 _P_N->n = _DP->n + *_codim - 1;
1034 _P_N->np = 0;
1035 for (l = 0; l < _DP->np; l++) {
1036 for (i = 0; i < *_codim; i++) {
1037 ip = 1;
1038 k = 0;
1039 while (ip && (k < _PTL->nv)) {
1040 sum = 0;
1041 for (d = 0; d < _DP->n; d++)
1042 sum += _DP->x[l][d] * _P->x[_V->v[k]][d];
1043 if (((_PTL->S[*_n][k] == i) && (sum < -1)) ||
1044 ((_PTL->S[*_n][k] != i) && (sum < 0)))
1045 ip = 0;
1046 k++;
1047 }
1048 if (ip == 1) {
1049 assert(_P_N->np < POINT_Nmax);
1050 for (d = 0; d < _DP->n; d++)
1051 _P_N->x[_P_N->np][d + *_codim - 1] = _DP->x[l][d];
1052 for (d = 1; d < *_codim; d++) {
1053 if (d == i)
1054 _P_N->x[_P_N->np][d - 1] = 1;
1055 else
1056 _P_N->x[_P_N->np][d - 1] = 0;
1057 }
1058 _P_N->np++;
1059 }
1060 }
1061 }
1062 _P_D->n = _P->n + *_codim - 1;
1063 _P_D->np = 0;
1064 for (l = 0; l < _P->np; l++) {
1065 for (i = 0; i < *_codim; i++) {
1066 ip = 1;
1067 k = 0;
1068 while (ip && (k < _P_N->np)) {
1069 sum = 0;
1070 for (d = 0; d < _P->n; d++)
1071 sum += _P->x[l][d] * _P_N->x[k][d + *_codim - 1];
1072 if (i > 0) {
1073 if (((_P_N->x[k][i - 1] == 1) && (sum < -1)) ||
1074 ((_P_N->x[k][i - 1] == 0) && (sum < 0)))
1075 ip = 0;
1076 } else {
1077 c = -1;
1078 for (d = 0; d < *_codim - 1; d++)
1079 if (_P_N->x[k][d] == 1)
1080 c = 0;
1081 if (sum < c)
1082 ip = 0;
1083 }
1084 k++;
1085 }
1086 if (ip == 1) {
1087 assert(_P_D->np < POINT_Nmax);
1088 for (d = 0; d < _P->n; d++)
1089 _P_D->x[_P_D->np][d + *_codim - 1] = _P->x[l][d];
1090 for (d = 1; d < *_codim; d++) {
1091 if (d == i)
1092 _P_D->x[_P_D->np][d - 1] = 1;
1093 else
1094 _P_D->x[_P_D->np][d - 1] = 0;
1095 }
1096 _P_D->np++;
1097 }
1098 }
1099 }
1100 }
1101
Int_Pos(int i,Interval_List * _IL,Poset_Element_List * _PEL)1102 Poset Int_Pos(int i, Interval_List * _IL, Poset_Element_List * _PEL)
1103 {
1104 Poset P;
1105
1106 P.x = _PEL->L[_IL->L[i].min];
1107 P.y = _PEL->L[_IL->L[i].max];
1108 return P;
1109 }
1110
M_To_B(BPoly * _BP,BPoly * _MP,int d,int rho)1111 void M_To_B(BPoly *_BP, BPoly *_MP, int d, int rho){
1112
1113 int M[Pos_Max][Pos_Max] = {{0}, {0}}, i, u, v;
1114
1115 if(rho == 0)
1116 M[0][0] = 1;
1117 else
1118 for(u = 0; u <= rho; u++)
1119 for(v = 0; v < rho/2 + (rho % 2); v++)
1120 M[rho - u][rho - v] = _MP->B[u][v];
1121 for(i = 1; i <= d - rho; i++){
1122 for(u = rho + i; u > 0; u--) /* (u degree != 0) && (v degree != 0) */
1123 for(v = rho + i; v > 0; v--)
1124 M[u][v] = - M[u-1][v] + M[u][v-1];
1125 for(v = rho + i; v > 0; v--) /* (u degree == 0) && (v degree != 0) */
1126 M[0][v] = M[0][v-1];
1127 for(u = rho + i; u > 0; u--) /* (u degree != 0) && (v degree == 0) */
1128 M[u][0] = - M[u-1][0];
1129 M[0][0] = 0; /* (u degree == 0) && (v degree == 0) */
1130 }
1131 for(u = 0; u <= d; u++)
1132 for(v = 0; v < d/2 + (d % 2); v++)
1133 _BP->B[u][v] += M[u][v];
1134 }
1135
N_To_B(BPoly * _BP,BPoly * _NP,int d,int rho)1136 void N_To_B(BPoly *_BP, BPoly *_NP, int d, int rho){
1137
1138 int N[Pos_Max][Pos_Max] = {{0}, {0}}, i, u, v;
1139
1140 if(rho == 0)
1141 N[0][0] = 1;
1142 else
1143 for(u = 0; u <= rho; u++)
1144 for(v = 0; v < rho/2 + (rho % 2); v++)
1145 N[u][v] = _NP->B[u][v];
1146 for(i = 1; i <= d - rho; i++){
1147 for(u = rho + i; u > 0; u--) /* (u degree != 0) && (v degree != 0) */
1148 for(v = rho / 2 + i; v > 0; v--)
1149 N[u][v] = N[u-1][v-1] - N[u][v];
1150 for(v = rho/2 + i; v > 0; v--) /* (u degree == 0) && (v degree != 0) */
1151 N[0][v] *= -1;
1152 for(u = rho + i; u >= 0; u--) /* (u degree >= 0) && (v degree == 0) */
1153 N[u][0] *= -1;
1154 }
1155 for(u = 0; u <= d; u++)
1156 for(v = 0; v < d/2 + (d % 2); v++)
1157 _BP->B[u][v] -= N[u][v];
1158 }
1159
Make_B_Poly(Cone * _C,Poset_Element_List * _PEL,Interval_List * _IL,BPoly * _BL)1160 void Make_B_Poly(Cone * _C, Poset_Element_List * _PEL, Interval_List * _IL,
1161 BPoly *_BL)
1162 {
1163 int D, i, j, k, d, n, n_;
1164 Poset Pos, pos;
1165
1166 for (i = 0; i < _IL->n; i++) {
1167 Pos = Int_Pos(i, _IL, _PEL);
1168 d = (Pos.y.dim - Pos.x.dim);
1169 for (j = 0; j <= d; j++)
1170 for (k = 0; k <= d / 2; k++)
1171 _BL[i].B[j][k] = 0;
1172 if (d == 0)
1173 _BL[i].B[0][0] = 1;
1174 }
1175 for (D = 1; D <= _C->dim; D++)
1176 for (n = 0; n < _IL->n; n++) {
1177 Pos = Int_Pos(n, _IL, _PEL);
1178 if ((Pos.y.dim - Pos.x.dim) == D) {
1179 for (d = 0; d < D; d++)
1180 for (n_ = 0; n_ < n; n_++) {
1181 pos = Int_Pos(n_, _IL, _PEL);
1182 if (Interval_Check(0, &pos.x, &Pos.x, _C))
1183 if (Interval_Check((D - d), &pos.y, &Pos.y, _C))
1184 M_To_B(&_BL[n], &_BL[n_], D, d);
1185 if (Interval_Check(0, &pos.y, &Pos.y, _C))
1186 if (Interval_Check((D - d), &Pos.x, &pos.x, _C))
1187 N_To_B(&_BL[n], &_BL[n_], D, d);
1188 }
1189 }
1190 }
1191 }
1192
Compute_E_Poly(EPoly * _EP,PolyPointList * _P_D,VertexNumList * _V_D,EqList * _E_D,PolyPointList * _P_N,VertexNumList * _V_N,EqList * _E_N,int * _codim,Flags * _F,time_t * _Tstart,clock_t * _Cstart)1193 void Compute_E_Poly(EPoly *_EP,
1194 PolyPointList * _P_D, VertexNumList * _V_D, EqList * _E_D,
1195 PolyPointList * _P_N, VertexNumList * _V_N, EqList * _E_N,
1196 int *_codim, Flags * _F, time_t *_Tstart, clock_t *_Cstart){
1197
1198 /* Should compute _EP from the rest. Probably requires alignment of
1199 _E_D with _V_N and of _V_D with _E_N */
1200
1201 Interval_List IL;
1202 SPoly *_S_D = NULL, *_S_N = NULL;
1203 BPoly *_BL = NULL;
1204 Poset_Element_List PEL_D, PEL_N;
1205 FaceInfo *_I_D;
1206 Cone *_C_D, *_C_N;
1207
1208 _I_D = (FaceInfo *) malloc(sizeof(FaceInfo));
1209 if (_I_D == NULL) Die("Unable to alloc space for FaceInfo _I_D");
1210 _C_D = (Cone *) malloc(sizeof(Cone));
1211 if (_C_D == NULL) Die("Unable to alloc space for Cone _C_D");
1212 _C_N = (Cone *) malloc(sizeof(Cone));
1213 if (_C_N == NULL) Die("Unable to alloc space for Cone _C_N");
1214
1215 Make_Incidence(_P_D, _V_D, _E_D, _I_D);
1216
1217 Make_Cone(_C_D, _C_N, _I_D, _P_D);
1218 PEL_D.n = Num_Pos(_C_D); PEL_N.n = Num_Pos(_C_N);
1219
1220 PEL_D.L = (Poset_Element *) calloc(PEL_D.n, sizeof(Poset_Element));
1221 if (PEL_D.L == NULL) Die("Unable to alloc space for PEL_D.L");
1222 PEL_N.L = (Poset_Element *) calloc(PEL_D.n, sizeof(Poset_Element));
1223 if (PEL_N.L == NULL) Die("Unable to alloc space for PEL_N.L");
1224
1225 Make_PosetList(_C_D, &PEL_D); Make_PosetList(_C_N, &PEL_N);
1226 if(_F->t) Time_Info(_Tstart, _Cstart, " BEGIN S-Poly");
1227
1228 _S_D = (SPoly *) calloc(PEL_D.n, sizeof(SPoly ));
1229 if (_S_D == NULL) Die("Unable to alloc space for SPoly _S_D");
1230 _S_N = (SPoly *) calloc(PEL_D.n, sizeof(SPoly ));
1231 if (_S_N == NULL) Die("Unable to alloc space for SPoly _S_N");
1232
1233 Make_S_Poly(_C_N, _V_D, _E_D, _P_D, &PEL_D, _S_D, _F->S, _F->T);
1234 Make_S_Poly(_C_D, _V_N, _E_N, _P_N, &PEL_N, _S_N, _F->S, _F->T);
1235
1236 if(_F->t) Time_Info(_Tstart, _Cstart, " BEGIN B-Poly");
1237
1238 IL.L = (Interval *) calloc(((1 + PEL_D.n)/2 + 1)*PEL_D.n,
1239 sizeof(Interval));
1240 if (IL.L == NULL) Die("Unable to alloc space for IL.L");
1241
1242 Make_Intervallist(&IL, &PEL_D, _C_D);
1243
1244 _BL = (BPoly *) calloc(IL.n, sizeof(BPoly));
1245 if (_BL == NULL) Die("Unable to alloc space for _BL");
1246
1247 Make_B_Poly(_C_D, &PEL_D, &IL, _BL);
1248
1249 if(_F->t) Time_Info(_Tstart, _Cstart, " BEGIN E-Poly");
1250 SB_To_E(_EP, _C_D, &PEL_D, _BL, &IL, _S_D, _S_N, _codim);
1251
1252 free(PEL_D.L); free(PEL_N.L); free(_S_D); free(_S_N); free(IL.L);
1253 free(_BL); free(_I_D); free(_C_D); free(_C_N);
1254 }
1255
Make_E_Poly(FILE * outFILE,CWS * _W,PolyPointList * _CP,VertexNumList * _CV,EqList * _CE,int * _codim,Flags * _F,int * _D)1256 void Make_E_Poly(FILE * outFILE, CWS * _W, PolyPointList * _CP,
1257 VertexNumList * _CV, EqList * _CE, int *_codim,
1258 Flags * _F, int *_D)
1259 {
1260 time_t Tstart;
1261 clock_t Cstart;
1262 int n;
1263 /* Interval_List IL;
1264 SPoly *_S_D = NULL, *_S_N = NULL;
1265 BPoly *_BL = NULL; */
1266 EPoly EP;
1267 /* Poset_Element_List PEL_D, PEL_N;*/
1268 PartList *_PTL;
1269 PolyPointList *_P = NULL, *_DP = NULL, *_P_D, *_P_N;
1270 VertexNumList *_V = NULL, *_DV = NULL, *_V_D, *_V_N;
1271 EqList *_E = NULL, *_DE = NULL, *_E_D, *_E_N;
1272 /*FaceInfo *_I_D;
1273 Cone *_C_D, *_C_N;*/
1274 LInfo *_L = NULL;
1275
1276 /* =============== Begin of Static Allocation =================== */
1277 _PTL = (PartList *) malloc(sizeof(PartList));
1278 if (_PTL == NULL) Die("Unable to alloc space for PartList _PTL");
1279 _P_D = (PolyPointList *) malloc(sizeof(PolyPointList));
1280 if (_P_D == NULL) Die("Unable to alloc space for PolyPointLis _P_D");
1281 _P_N = (PolyPointList *) malloc(sizeof(PolyPointList));
1282 if (_P_N == NULL) Die("Unable to alloc space for PolyPointLis _P_N");
1283 _V_D = (VertexNumList *) malloc(sizeof(VertexNumList));
1284 if (_V_D == NULL) Die("Unable to alloc space for VertexNumList _V_D");
1285 _V_N = (VertexNumList *) malloc(sizeof(VertexNumList));
1286 if (_V_N == NULL) Die("Unable to alloc space for VertexNumList _V_N");
1287 _E_D = (EqList *) malloc(sizeof(EqList));
1288 if (_E_D == NULL) Die("Unable to alloc space for EqList _E_D");
1289 _E_N = (EqList *) malloc(sizeof(EqList));
1290 if (_E_N == NULL) Die("Unable to alloc space for EqList _E_N");
1291 if (_F->Lv || _F->Lp){
1292 _L = (LInfo *) malloc(sizeof(LInfo));
1293 if (_L == NULL) Die("Unable to alloc space for LInfo _L");
1294 }
1295 if(_F->N){
1296 _DP = (PolyPointList *) malloc(sizeof(PolyPointList));
1297 if (_DP == NULL) Die("Unable to alloc space for PolyPointLis _DP");
1298 _DV = (VertexNumList *) malloc(sizeof(VertexNumList));
1299 if (_DV == NULL) Die("Unable to alloc space for VertexNumList _DV");
1300 _DE = (EqList *) malloc(sizeof(EqList));
1301 if (_DE == NULL) Die("Unable to alloc space for EqList _DE");
1302 }
1303 else{
1304 _P = (PolyPointList *) malloc(sizeof(PolyPointList));
1305 if (_P == NULL) Die("Unable to alloc space for PolyPointLis _P");
1306 _V = (VertexNumList *) malloc(sizeof(VertexNumList));
1307 if (_V == NULL) Die("Unable to alloc space for VertexNumList _V");
1308 _E = (EqList *) malloc(sizeof(EqList));
1309 if (_E == NULL) Die("Unable to alloc space for EqList _E");
1310 }
1311 /* =============== End of Static Allocation =================== */
1312
1313 if(_F->N){
1314 _P=_CP; _E=_CE; _V=_CV;
1315 Make_Dual_Poly(_P, _V, _E, _DP);
1316 Find_Equations(_DP, _DV, _DE);
1317 Sort_PPL(_DP, _DV);
1318 }
1319 else{
1320 _DP=_CP; _DE=_CE; _DV=_CV;
1321 Make_Dual_Poly(_DP, _DV, _DE, _P);
1322 Find_Equations(_P, _V, _E);
1323 Sort_PPL(_P, _V);
1324 }
1325 Tstart = time(NULL); Cstart = clock(); _F->Test = 0;fflush(0);
1326 { NEF_Flags NF; NF.Sym=_F->Sym;
1327 NF.noconvex=_F->noconvex; NF.Test=0; NF.Sort=_F->Sort;
1328 part_nef(_P, _V, _E, _PTL, _codim, &NF);
1329 if(_F->y && (_PTL->n != 0) && (_W->nw == 0))
1330 fprintf(outFILE, "%d %d Vertices of Poly in M-lattice: ", _DP->n, _DV->nv);
1331 if(!(_F->y && (_PTL->n==0))){
1332 fprintf(outFILE, "M:%d %d N:%d %d ", _DP->np, _E->ne, _P->np, _V->nv);
1333 fprintf(outFILE, " codim=%d", *_codim);
1334 fprintf(outFILE, " #part="); fflush(0);
1335 fprintf(outFILE, "%d\n",_PTL->n);fflush(0);
1336 }
1337 }
1338 if(_F->y && (_W->nw == 0) && (_PTL->n != 0))
1339 PRINT_TMATRIX(&_DP->x[0][0], _DV->nv, _DP->n, POLY_Dmax);
1340 if(_F->V)
1341 Print_VL(_P, _V, "Vertices of P:");
1342 if((!_F->Lv) && (!_F->n) && _F->Lp)
1343 PRINT_PL(_P, "Points of Poly in N-Lattice:");
1344 if(_F->Lv || _F->Lp){
1345 Make_L(_P, _V, _L ,_F->Lp,_F->Lv);
1346 Print_L(_L,_F->Lp,_F->Lv);
1347 }
1348 if(_F->f){
1349 if(_F->Lv)
1350 PRINT_PL(_P, "Points of Poly in N-Lattice:");
1351 PRINT_Fibrations(_V, _P, _F /*, _PTL*/ );
1352 }
1353 for (n = 0; n < _PTL->n; n++) {
1354 /*if ((!_PTL->DirProduct[n]) || _F->Dir){*/
1355 if(_F->Dir==2) if(!_PTL->DirProduct[n]) continue ;
1356 Tstart = time(NULL);
1357 Cstart = clock();
1358 Make_Gore_Poly(_P, _DP, _P_D, _P_N, _V, _PTL, _codim, &n);
1359 _PTL->Proj[n] = Remove_Proj(_P_D, /*&n,*/ _codim);
1360 _PTL->DProj[n] = Remove_Proj(_P_N, /*&n,*/ _codim);
1361 Find_Equations(_P_D, _V_D, _E_D);
1362 Find_Equations(_P_N, _V_N, _E_N);
1363 Sort_PPL(_P_N, _V_N);
1364
1365 if (((!_PTL->Proj[n]) || _F->Proj) && ((!_PTL->DirProduct[n]) || _F->Dir)
1366 && !_F->p && !_F->y) {
1367 Make_EN(_P_D, _V_D, _E_N, _codim);
1368 Compute_E_Poly(&EP, _P_D, _V_D, _E_D, _P_N, _V_N, _E_N, _codim, _F,
1369 &Tstart, &Cstart);}
1370 /*for (i=0;i<=10;i++){
1371 for (j=0;j<=10;j++) fprintf(outFILE,"%d ",EP.E[i][j]);
1372 fprintf(outFILE,"\n");}*/
1373
1374 if(!_F->n && !_F->g && !_F->d && !_F->y){
1375 Output( _P, /*_DP,*/ _P_D, _W, &EP, /*_E,*/
1376 _V, &n, _PTL, _codim, outFILE, _F, _D, _L);
1377 if (((!_PTL->Proj[n]) || _F->Proj) && ((!_PTL->DirProduct[n]) || _F->Dir))
1378 Time_Info(&Tstart, &Cstart, "");
1379 }
1380 if (((!_PTL->Proj[n]) || _F->Proj) && ((!_PTL->DirProduct[n]) || _F->Dir)){
1381 if(_F->g) PRINT_GORE(_P_D, *_codim, _F->gd, "Points of PG:");
1382 if(_F->d) PRINT_GORE(_P_N, *_codim, _F->dd, "Points of dual PG:");
1383 }
1384 }
1385 if(!_F->n && !_F->g && !_F->d && !_F->y) Print_Nefinfo(_PTL, /* _F,*/
1386 &Tstart, &Cstart);
1387 if(_F->n && N_Part(_PTL)) PRINT_PL(_P, "Points of Poly in N-Lattice:");
1388 /* =============== Begin of FREE Static Allocation =================== */
1389 free(_PTL);free(_P_D); free(_P_N); free(_V_D); free(_V_N); free(_E_D);
1390 free(_E_N); free(_L);
1391 if(_F->N){free(_DP); free(_DE); free(_DV); }
1392 else{free(_P); free(_E); free(_V); }
1393 /* =============== End of FREE Static Allocation =================== */
1394 }
1395
1396 void SL2Z_Make_Poly_UTriang(PolyPointList *P);
1397
AnalyseGorensteinCone(CWS * _CW,PolyPointList * _P,VertexNumList * _V,EqList * _E,int * _codim,Flags * _F)1398 void AnalyseGorensteinCone(CWS *_CW, PolyPointList *_P, VertexNumList *_V,
1399 EqList *_E, int *_codim, Flags * _F){
1400 /* _P should be the Gorenstein-polytope in M - what is called _P_N in
1401 certain other parts of the program */
1402
1403 time_t Tstart;
1404 clock_t Cstart;
1405 EPoly EP;
1406 int i, j, k, dim = _P->n - (*_codim * 2) + 1, chi, r=1;
1407 int h[POLY_Dmax][POLY_Dmax] = {{0}, {0}};
1408 PolyPointList *_P_D = (PolyPointList *) malloc(sizeof(PolyPointList));
1409 VertexNumList *_V_D = (VertexNumList *) malloc(sizeof(VertexNumList));
1410 EqList *_E_D = (EqList *) malloc(sizeof(EqList));
1411 /* the Gorenstein-polytope in N-lattice and its vertices and equations */
1412 EqList *_new_E_D = (EqList *) malloc(sizeof(EqList));
1413 /* the equations of _P_D in an order corresponding to the vertices of _P */
1414 PairMat VPM, VPM_D;
1415
1416 if ((_P_D == NULL)||(_V_D == NULL)||(_E_D == NULL)||(_new_E_D == NULL))
1417 Die("Unable to allocate space for _P_D in AnalyseGorensteinCone");
1418
1419 if (POLY_Dmax < _P->n + 1){/* only relevant if index == 1 */
1420 printf("Please increase POLY_Dmax to at least %d = %d + 1\n",
1421 (_P->n + 1), _P->n);
1422 printf("(POLY_Dmax >= dim(cone) = dim(support) + 1 required)\n");
1423 assert(*_codim == 1);
1424 exit(0);}
1425 /* Print_PPL(_P, "_P before sorting");fflush(0); */
1426 Find_Equations(_P,_V,_E);
1427 Make_VEPM(_P, _V, _E, VPM);
1428 Complete_Poly(VPM, _E, _V->nv, _P);
1429 Sort_PPL(_P,_V);
1430 /* Print_PPL(_P, "_P after sorting"); */
1431 /* Print_VL(_P,_V, "_V"); */
1432 /* Print_EL(_E, &_P->n, 0, "_E"); */
1433
1434 /* Create _P_D from the equations of _P:
1435 _P_D is the hull of the generators of the cone dual to the one over _P,
1436 Gorenstein <-> _P_D lies in a plane at distance 1 from the origin,
1437 index = modulus of the last component of the equation of this plane */
1438 _P_D->n = _P->n + 1;
1439 _P_D->np = _E->ne;
1440 for (i=0; i<_P_D->np; i++){
1441 for (j=0; j<_P_D->n; j++)
1442 _P_D->x[i][j] = _E->e[i].a[j];
1443 _P_D->x[i][_P->n] = _E->e[i].c;}
1444 /* Print_PPL(_P_D, "_P_D "); */
1445 Find_Equations(_P_D,_V_D,_E_D);
1446 /* Print_EL(_E_D, &_P_D->n, 0, "_E_D"); */
1447 if (_E_D->ne == 1){
1448 if ((_E_D->e[0].c != 1)&&(_E_D->e[0].c != -1)) r = 0;
1449 else if (abs((int) _E_D->e[0].a[_P->n]) != *_codim){
1450 printf("Warning: Input has index %d, should be %d! ",
1451 abs(_E_D->e[0].a[_P->n]), *_codim);
1452 r = 0;}}
1453 else {assert(_E_D->ne > _P_D->n); r = 0;}
1454
1455 Tstart = time(NULL); Cstart = clock(); _F->Test = 0;
1456 for (i = 0; i < _CW->nw; i++) {
1457 fprintf(outFILE, "%d ", (int) _CW->d[i]);
1458 for (j = 0; j < _CW->N; j++)
1459 fprintf(outFILE, "%d ", (int) _CW->W[i][j]);
1460 if (i + 1 < _CW->nw)
1461 fprintf(outFILE, " "); }
1462 fflush(0);
1463
1464 if (r) {
1465 /* shift _P_D into a plane through the origin: */
1466 for (i=0; i<_P_D->np; i++){
1467 for (j=0; j<_P_D->n; j++)
1468 _P_D->x[i][j] -= _P_D->x[_P_D->np-1][j];}
1469 /* Print_PPL(_P_D, "_P_D before Make_Poly_UTriang"); */
1470 SL2Z_Make_Poly_UTriang(_P_D);
1471 /* Print_PPL(_P_D, "_P_D after Make_Poly_UTriang"); */
1472 _P_D->n--;
1473 for (i=0; i<_P_D->np; i++) assert(_P_D->x[i][_P->n] == 0);
1474 /* Print_PPL(_P_D, "_P_D after reduction"); */
1475 Find_Equations(_P_D,_V_D,_E_D);
1476 /* Print_EL(_E_D, &_P_D->n, 0, "_E_D"); */
1477 Sort_VL(_V_D);
1478 /* Make_VEPM(_P, _V, _E, VPM); */
1479 /* Print_Matrix(VPM, _E->ne, _V->nv, "VPM"); */
1480 Make_VEPM(_P_D, _V_D, _E_D, VPM_D);
1481 /* Print_Matrix(VPM_D, _E_D->ne, _V_D->nv, "VPM_D"); */
1482 Complete_Poly(VPM_D, _E_D, _V_D->nv, _P_D);
1483 /* Print_PPL(_P_D, "_P_D after Complete_Poly"); */
1484 assert(_E_D->ne = _V->nv);
1485 assert(_E->ne = _V_D->nv);
1486 /* Compute _new_E_D->ne by comparing VPM and VPM_D:
1487 if (VPM[j][i] == VPM_D[k][j]) for all j then _new_E_D->e[i] = _E_D[k] */
1488 _new_E_D->ne = _E_D->ne;
1489 for (i=0;i<_V->nv;i++){
1490 for(k=0;k<_V->nv;k++){
1491 for (j=0;j<_V_D->nv;j++) if (VPM[j][i] != VPM_D[k][j]) break;
1492 if (j==_V_D->nv) {_new_E_D->e[i] = _E_D->e[k]; break;} }
1493 if (k >= _V->nv){
1494 printf("k = %d, _V->nv = %d\n", k, _V->nv);
1495 Print_Matrix(VPM, _E->ne, _V->nv, "VPM");
1496 Print_Matrix(VPM_D, _E_D->ne, _V_D->nv, "VPM_D");}
1497 assert(k<_V->nv);}
1498 /* Print_EL(_new_E_D, &_P_D->n, 0, "_new_E_D");
1499 Make_VEPM(_P_D, _V_D, _new_E_D, VPM_D);
1500 Print_Matrix(VPM_D, _E_D->ne, _V_D->nv, "VPM_D"); */
1501 if (_F->N){ /* swap M and N */
1502 PolyPointList *_auxP = _P_D;
1503 VertexNumList *_auxV = _V_D;
1504 EqList *_auxE = _new_E_D;
1505 _P_D = _P; _V_D = _V; _new_E_D = _E;
1506 _P = _auxP; _V = _auxV; _E = _auxE;}
1507 fprintf(outFILE,"M:%d %d ",_P->np,_V->nv);
1508 fprintf(outFILE,"N:%d %d ",_P_D->np,_V_D->nv);
1509 if (!_F->g && !_F->d){
1510 Compute_E_Poly(&EP, _P_D, _V_D, _new_E_D, _P, _V, _E, _codim, _F,
1511 &Tstart, &Cstart);
1512 chi = Make_Mirror(&EP, h, _P_D->n + 1, dim);
1513 fprintf(outFILE, "H:");
1514 for (i = 1; i < dim; i++)
1515 fprintf(outFILE, "%d ", h[1][i]);
1516 fprintf(outFILE, "[%d]", chi);
1517 for (i = 1; i <= dim/2; i++)
1518 if (h[0][i] != 0)
1519 fprintf(outFILE, " h%d=%d", i, h[0][i]);
1520 if(h[0][0] != 1) fprintf(outFILE, " h%d=%d", 0, h[0][0]);
1521 puts("");
1522 if (_F->H) PrintDiamond(h, dim);}
1523 else puts("");
1524 if(_F->V)
1525 Print_VL(_P_D, _V_D, "Vertices of support in N:");
1526 if (_F->g) Print_PPL(_P_D, "Points of support in N:");
1527 if (_F->d) Print_PPL(_P, "Points of support in M:");
1528 if (_F->t) Time_Info(&Tstart, &Cstart, "");
1529 if (_F->N){ /* revert M-N-swap (necessary because of alloc/free) */
1530 PolyPointList *_auxP = _P_D;
1531 VertexNumList *_auxV = _V_D;
1532 EqList *_auxE = _new_E_D;
1533 _P_D = _P; _V_D = _V; _new_E_D = _E;
1534 _P = _auxP; _V = _auxV; _E = _auxE;}}
1535 else{
1536 if (_F->N) fprintf(outFILE,"N:%d %d ",_P->np,_V->nv);
1537 else fprintf(outFILE,"M:%d %d ",_P->np,_V->nv);
1538 fprintf(outFILE,"F:%d ",_E->ne);
1539 puts("");
1540 if ((_F->Rv) || ((_F->V)&&(_F->N)))
1541 Print_VL(_P, _V, "Vertices of input polytope:");}
1542 free(_P_D); free(_E_D); free(_V_D); free(_new_E_D);
1543 }
1544
1545