1 #include "Global.h"
2 #include "Subpoly.h"
3 #include "Rat.h"
4 /* NB mod 2^32, works for #poly<2^32 */
5
6 /* #include <types.h> -> defines _ILP32 (32-bit programs)
7 [ -> isa_defs.h ] or _LP64 (long and pointer 64 bits)
8 #include <stdint.h> (on Linux systems?) ... uintptr_t */
9
10 /* #include <limits.h> LONG_MAX = 2147483647 vs. 2^63-1 */
11
12 /* L->dbname points at constant string => don't change, use:
13 * dbname, which has "File_Ext_NCmax" extra characters allocated !!! */
14
15 /* uchar=unsigned char uint=unsigned int
16 * NP=#Polys NB=#Bytes #files=#nv's with NP>0 #lists=#(nv,nuc) with NP>0
17 *
18 * uchar rd k_1 ... k_rd // ... not in data base !!
19 *
20 * uchar dim #files nv_max nuc_max
21 * uint #lists hNF hSM hNM hNB slNF slSM slNM slNB
22 * uchar v1 #nuc's with v1
23 * uchar nuc1 uint #NF(v1,nuc1) uchar nuc2 uint #NV(v1,nuc2) ...
24 * uchar v2 #nuc's with v2
25 * uchar nuc1 uint #NF(v2,nuc1) uchar nuc2 uint #NV(v2,nuc2) ...
26 * uchar "all hNF honest nf's"
27 * uchar "all slNF sublattice {nv nuc nf[]}'s"
28 */
29
Small_Make_Dual(PolyPointList * _P,VertexNumList * _V,EqList * _E)30 void Small_Make_Dual(PolyPointList *_P, VertexNumList *_V, EqList *_E){
31 int i;
32 EqList AE=*_E;
33 VNL_to_DEL(_P, _V, _E);
34 assert(EL_to_PPL(&AE, _P, &_P->n));
35 _V->nv=_P->np;
36 for (i=0;i<_P->np;i++) _V->v[i]=i;
37 }
38
Polyi_2_DBo(char * polyi,char * dbo)39 void Polyi_2_DBo(char *polyi,char *dbo)
40 { char *dbnames = (char *) malloc(1+strlen(dbo)+File_Ext_NCmax), *fx;
41 FILE *F=fopen(polyi,"rb"), *Finfo, *Fv, *Fsl;
42 time_t Tstart=time(NULL); FInfoList L; UPint tNF=0;
43 int d, v, nu, i, j, list_num, sl_nNF, sl_SM, sl_NM, sl_NB; Along tNB=0;
44
45 if (!*polyi) {puts("With -do you require -pi or -di and -pa"); exit(0);}
46
47 if(F==NULL) {printf("Input file %s not found\n",polyi); exit(0);}
48 strcpy(dbnames,dbo); fx=&dbnames[strlen(dbo)+1];
49 strcat(dbnames,".info"); Finfo=fopen(dbnames,"w"); assert(Finfo!=NULL);
50 printf("Read %s (",polyi); fflush(stdout);
51
52 Init_FInfoList(&L); /* start reading the file */
53 d=fgetc(F); assert(d==0); /* for(i=0;i<d;i++) fgetc(F); */
54 if(d){printf("Recursion depth %d forbidden in DB !!\n",d);Finfo=stdout;}
55
56 d=fgetc(F); L.nV=fgetc(F); L.nVmax=fgetc(F); L.NUCmax=fgetc(F);
57 list_num=fgetUI(F);
58 L.nNF=fgetUI(F); L.nSM=fgetUI(F); L.nNM=fgetUI(F); L.NB=fgetUI(F);
59 sl_nNF=fgetUI(F); sl_SM=fgetUI(F); sl_NM=fgetUI(F); sl_NB=fgetUI(F);
60
61 for(i=0;i<L.nV;i++)
62 { v=fgetc(F); L.nNUC[v]=fgetc(F); /* read #nuc's per #Vert */
63 for(j=0;j<L.nNUC[v];j++)
64 { L.NFnum[v][nu=fgetc(F)]=fgetUI(F); /* read nuc and #NF(v,nu)*/
65 tNF+=L.NFnum[v][nu]; tNB+=nu*L.NFnum[v][nu];
66 }
67 } assert(tNF==L.nNF); assert(0== (unsigned int)(tNB-L.NB)); L.NB=tNB;
68
69 printf("%lldpoly +%dsl %lldb) write %s.* (%d files) ",
70 2*L.nNF-L.nSM-L.nNM, 2*sl_nNF-sl_SM-sl_NM, L.NB+sl_NB,
71 dbo,L.nV+1+(sl_nNF>0)); fflush(stdout);
72
73 fprintf(Finfo, /* write Finfo */
74 "%d %d %d %d %d %lld %d %lld %lld %d %d %d %d\n\n",
75 d, L.nV, L.nVmax, L.NUCmax,list_num,
76 L.nNF,L.nSM,L.nNM,L.NB, sl_nNF, sl_SM, sl_NM, sl_NB);
77
78 for(v=d+1;v<=L.nVmax;v++) if(L.nNUC[v]) /* honest info */
79 { i=0; fprintf(Finfo,"%d %d\n",v,L.nNUC[v]); /* v #nuc's */
80 for(nu=1;nu<=L.NUCmax;nu++) if(L.NFnum[v][nu])
81 { fprintf(Finfo,"%d %d%s",nu,L.NFnum[v][nu], /* nuc #NF(v,nuc) */
82 (++i<L.nNUC[v]) ? " " : "\n");
83 }
84 } if(Finfo==stdout) exit(0);
85 if(ferror(Finfo)) {printf("File error in %s\n",dbnames);exit(0);}
86 fclose(Finfo); fflush(stdout);
87
88 for(v=d+1;v<=L.nVmax;v++) if(L.nNUC[v]) /* write honest polys */
89 { char ext[4]={'v',0,0,0};
90 ext[1]='0' + v / 10; ext[2]='0' + v % 10;
91 strcpy(fx,ext); Fv=fopen(dbnames,"wb"); assert(Fv!=NULL);
92
93 for(nu=1;nu<=L.NUCmax;nu++) if(L.NFnum[v][nu])
94 { int vnuNB=nu*L.NFnum[v][nu];
95 for(i=0;i<vnuNB;i++) fputc(fgetc(F),Fv);
96 }
97 if(ferror(Fv)) {printf("File error in %s\n",dbnames);exit(0);}
98 fclose(Fv);
99 }
100
101 if(sl_nNF) /* write sublattice polys */
102 { strcpy(fx,"sl"); Fsl=fopen(dbnames,"wb");
103 assert(Fsl!=NULL); for(i=0;i<sl_NB;i++) fputc(fgetc(F),Fsl);
104 if(ferror(Fsl)) {printf("File error in %s\n",dbnames);exit(0);}
105 fclose(Fsl);
106 }
107
108 printf("done (%ds)\n",(int) difftime(time(NULL),Tstart));
109
110 if(ferror(F)) {printf("File error in %s\n",polyi);exit(0);} fclose(F);
111 }
112
Init_DB(NF_List * _NFL)113 void Init_DB(NF_List *_NFL){
114 /* Read the database, create RAM_poly;
115 for given, nv, nuc the matching is as follows:
116 DB: |0|1|...|B-1|B|...|2B|...|((n-1)/B)*B|...|n-2|n-1|
117 RAM: |0 | 1|...| (n-1)/B)-1|
118 (n...nNF[nv][nuc], B...BLOCK_LENGTH, the offsets are DB->Fv_pos[v][nuc]
119 and DB->RAM_pos[v][nuc], respectively;
120 each entry |x| corresponds to nuc unsigned characters) */
121
122 time_t Tstart=time(NULL);
123 char *dbname = (char *) malloc(1+strlen(_NFL->dbname)+File_Ext_NCmax), *fx;
124 DataBase *DB=&_NFL->DB;
125 int d, v, nu, i, j, list_num, sl_nNF, sl_SM, sl_NM, sl_NB,
126 RAM_pos=0; Along RAM_size=0;
127
128 strcpy(dbname,_NFL->dbname);
129 printf("Reading data-base %s: ",dbname);
130 strcat(dbname,".info");
131 fx=&dbname[strlen(_NFL->dbname)+1];
132
133 /* read the info-file: */
134 DB->Finfo=fopen(dbname,"r");
135 assert(DB->Finfo!=NULL);
136 fscanf(DB->Finfo, "%d %d %d %d %d %lld %d %lld %lld %d %d %d %d",
137 &d, &DB->nV, &DB->nVmax, &DB->NUCmax, &list_num,
138 &DB->nNF, &DB->nSM, &DB->nNM, &DB->NB, &sl_nNF, &sl_SM, &sl_NM, &sl_NB);
139 printf("%lld+%dsl %lldnf %lldb",
140 2*(DB->nNF)-DB->nSM-DB->nNM,
141 2*sl_nNF-sl_SM-sl_NM, DB->nNF+sl_nNF, DB->NB+sl_NB);
142 /* if( _FILE_OFFSET_BITS < 64 ) assert(DB->NB <= LONG_MAX); Along */
143 if(_NFL->d) assert(d==_NFL->d); else _NFL->d=d;
144
145 for (v=1;v<VERT_Nmax;v++){
146 DB->nNUC[v]=0;
147 for(nu=0; nu<NUC_Nmax; nu++) DB->NFnum[v][nu]=0;}
148
149 for(i=0; i<DB->nV; i++){
150 fscanf(DB->Finfo,"%d",&v);
151 fscanf(DB->Finfo,"%d",&(DB->nNUC[v]));
152 for(j=0;j<DB->nNUC[v];j++){
153 fscanf(DB->Finfo,"%d", &nu);
154 fscanf(DB->Finfo,"%d", &(DB->NFnum[v][nu]));
155 RAM_size+=nu*((DB->NFnum[v][nu]-1)/BLOCK_LENGTH); } }
156
157 if(ferror(DB->Finfo)) {printf("File error in %s\n",dbname); exit(0);}
158 fclose(DB->Finfo);
159 fflush(stdout); assert(RAM_size<=INT_MAX);
160
161 DB->RAM_NF=(unsigned char *) malloc(RAM_size);
162 assert(DB->RAM_NF!=NULL);
163
164 /* read the DB-files and create RAM_NF: */
165 for (v=2;v<=DB->nVmax;v++) if(DB->nNUC[v]){
166 char ext[4]={'v',0,0,0};
167 ext[1]='0' + v / 10; ext[2]='0' + v % 10;
168 strcpy(fx,ext);
169 DB->Fv[v]=fopen(dbname,"rb");
170 assert(DB->Fv[v]!=NULL);
171 FSEEK(DB->Fv[v],0,SEEK_END);
172
173 DB->Fv_pos[v][0]=0;
174 for (nu=0;nu<=DB->NUCmax;nu++){
175 DB->RAM_pos[v][nu]=RAM_pos;
176 for (i=0; i<(DB->NFnum[v][nu]-1)/BLOCK_LENGTH; i++){
177 FSEEK(DB->Fv[v], DB->Fv_pos[v][nu]+(i+1)*BLOCK_LENGTH*nu, SEEK_SET);
178 for (j=0; j<nu; j++) DB->RAM_NF[RAM_pos++]=fgetc(DB->Fv[v]); }
179 DB->Fv_pos[v][nu+1]=DB->Fv_pos[v][nu]+DB->NFnum[v][nu]*nu; } }
180
181 printf(" done (%ds)\n",(int) difftime(time(NULL),Tstart));
182 fflush(stdout);
183 }
184
Compare_Poly(int * nuc,unsigned char * uc1,unsigned char * uc2)185 char Compare_Poly(int *nuc, unsigned char *uc1, unsigned char *uc2){
186 /* uc1 always encodes a single poly,
187 uc2 might encode a single poly or a mirror pair;
188 returns: 'm': if uc1 isn't in uc2, but its mirror is;
189 'i': if uc1 is in uc2;
190 'l': if uc1 is less than uc2;
191 'g': if uc1 is greater than uc2; */
192 switch (RIGHTminusLEFT(uc1,uc2,nuc)){
193 case -1: return 'g';
194 case 1: return 'l';
195 case 0: {if ((*uc1 + *uc2) % 4 == 3) return 'm'; else return 'i';}
196 default: puts("Sth. wrong in Compare_Poly!!!"); exit(0);} return 0;
197 }
198
Is_in_DB(int * nv,int * nuc,unsigned char * uc,NF_List * _NFL)199 int Is_in_DB(int *nv, int *nuc, unsigned char *uc, NF_List *_NFL){
200 /* Uses the following basic strategy for identifying the position of some
201 object x w.r.t. entries of a list l of length n:
202 min_pos=-1;
203 max_pos=n;
204 while (max_pos-min_pos>1){
205 int pos=(max_pos+min_pos)/2;
206 switch (Compare(x,l[pos]){
207 case 'equal': return sth.;
208 case 'less': {max_pos=pos; continue:}
209 case 'greater': min_pos=pos;}}
210 results in a return or l[min_pos] < x < l[max_pos=min_pos+1]
211 Applied here to locate x=uc=encoded polyhedron first w.r.t. RAM_poly and
212 then w.r.t. the location in the database */
213
214 int pos, min_pos=-1, max_RAM_pos, i, max_Fv_piece; Along Fv_pos;
215 unsigned char Aux_poly[BLOCK_LENGTH*NUC_Nmax];
216
217 DataBase *DB=&_NFL->DB; if (!DB->NFnum[*nv][*nuc]) return 0;
218
219 /* Analyse the position of uc w.r.t. DB->RAM_NF */
220 max_RAM_pos=(DB->NFnum[*nv][*nuc]-1)/BLOCK_LENGTH;
221 while (max_RAM_pos-min_pos>1){
222 pos=(max_RAM_pos+min_pos)/2;
223 switch
224 (Compare_Poly(nuc, uc,
225 &(DB->RAM_NF[DB->RAM_pos[*nv][*nuc]+(*nuc)*pos]))){
226 case 'm': return 0;
227 case 'i': return 1;
228 case 'l': {max_RAM_pos=pos; continue;}
229 case 'g': min_pos=pos;}}
230
231 /* Look for uc in the database: */ Fv_pos=max_RAM_pos;
232 if (Fv_pos==(DB->NFnum[*nv][*nuc]-1)/BLOCK_LENGTH)
233 max_Fv_piece=DB->NFnum[*nv][*nuc]-Fv_pos*BLOCK_LENGTH;
234 else max_Fv_piece=BLOCK_LENGTH;
235 min_pos=-1;
236 if (FSEEK(DB->Fv[*nv],
237 DB->Fv_pos[*nv][*nuc]+Fv_pos*
238 (Along) ( (*nuc)*BLOCK_LENGTH) , SEEK_SET)){
239 printf("Error in fseek in Is_in_DB!"); exit(0);}
240 for (i=0;i<(*nuc)*(max_Fv_piece);i++) Aux_poly[i]=fgetc(DB->Fv[*nv]);
241 while (max_Fv_piece-min_pos>1){
242 pos=(max_Fv_piece+min_pos)/2;
243 switch(Compare_Poly(nuc, uc, &(Aux_poly[(*nuc)*pos]))){
244 case 'm': {return 0;}
245 case 'i': {return 1;}
246 case 'l': {max_Fv_piece=pos; continue;}
247 case 'g': {min_pos=pos;} } }
248 return 0;
249 }
250
Add_Polya_2_DBi(char * dbi,char * polya,char * dbo)251 void Add_Polya_2_DBi(char *dbi,char *polya,char *dbo)
252 { FInfoList FIi, FIa, FIo; Along Apos, HIpos, HApos, Inp, tnb=0, tNF=0;
253 unsigned char ucI[NUC_Nmax],ucA[NUC_Nmax],*ucSL=NULL,*uc;int SLp[SL_Nmax];
254 int d, vI, nuI, IslNF, IslSM, IslNM, nu; unsigned Ili, u, Oli=0, IslNB;
255 int v, vA, nuA, AslNF, AslSM, AslNM, i; unsigned Ali, a; Along AslNB;
256 int s, slNF=0, slSM=0, slNM=0, slNB=0, slNP=0; UPint Anp;
257 int AmI=00,ms, newout=strcmp(dbi,dbo) && (*dbo),j=1+strlen(SAVE_FILE_EXT);
258 char *Ifx, *Ifn = (char *) malloc(j+strlen(dbi)+File_Ext_NCmax), *Ofx,
259 *Ofn = (char *) malloc(j+strlen(newout ? dbo : dbi)+File_Ext_NCmax);
260 FILE *FI, *FA, *FO; if(*polya==0) {puts("-pa file required"); exit(0);}
261 Init_FInfoList(&FIi); Init_FInfoList(&FIa);
262 strcpy(Ifn,dbi); Ifx=&Ifn[strlen(dbi)]; strcpy(Ifx,".info");
263 if(*dbo==0)dbo=dbi;
264 strcpy(Ofn,dbo); Ofx=&Ofn[strlen(dbo)];
265
266 if(NULL==(FI=fopen(Ifn,"r"))) {printf("Cannot open %s",Ifn);exit(0);}
267 if(NULL==(FA=fopen(polya,"rb"))){printf("Cannot open %s",polya);exit(0);}
268 fscanf(FI,"%d%d%d%d%d%lld%d%lld %lld %d%d%d%d",&d,&i,&j,&nu,&Ili,
269 &FIi.nNF,&FIi.nSM,&FIi.nNM,&FIi.NB,&IslNF,&IslSM,&IslNM,&IslNB);
270 FIi.nV=i; FIi.nVmax=j; FIi.NUCmax=nu;
271 /* printf("%d %d %d %d %d %d %d %d %lld %d %d %d %d\n",
272 d,FIi.nV,FIi.nVmax,FIi.NUCmax,Ili,FIi.nNF,FIi.nSM,
273 FIi.nNM,FIi.NB,IslNF,IslSM,IslNM,IslNB); */
274 for(i=0;i<FIi.nV;i++)
275 { fscanf(FI,"%d",&v); fscanf(FI,"%d",&j); FIi.nNUC[v]=j;
276 for(j=0;j<FIi.nNUC[v];j++)
277 { fscanf(FI,"%d",&nu); fscanf(FI,"%d",&FIi.NFnum[v][nu]);
278 tNF+=FIi.NFnum[v][nu];
279 }
280 } assert(tNF==FIi.nNF); tNF=0;
281 assert(!fgetc(FA)); /* rd==0 (recursion depth::no aux-file) */
282 Read_Bin_Info(FA,&j,&Ali,&AslNF,&AslSM,&AslNM,&AslNB,&FIa);assert(d==j);
283 strcpy(Ifx,".sl"); if(IslNF) assert(NULL != (FI=fopen(Ifn,"rb")));
284 if((IslNB+AslNB))
285 assert(NULL != (ucSL = (unsigned char *) malloc( (IslNB+AslNB) )));
286 HApos=FTELL(FA); FSEEK(FA,0,SEEK_END); Apos=FTELL(FA);
287 Inp=2*FIi.nNF-FIi.nSM-FIi.nNM; Anp=2*FIa.nNF-FIa.nSM-FIa.nNM;
288 printf("Data on %s: %lld+%dsl %lldb (%dd)\n", dbi,Inp,
289 /* Islp= */ 2*IslNF-IslNM-IslSM,FIi.NB+slNB,d);
290 printf("Data on %s: %u+%dsl %lldb (%dd)\n", polya,Anp,
291 /* Aslp= */ 2*AslNF-AslNM-AslSM, Apos,d);
292 assert(HApos+FIa.NB+AslNB==Apos); FSEEK(FA,-AslNB,SEEK_CUR);
293 s=0; if(s<AslNF) AuxGet_vn_uc(FA,&vA,&nuA,ucA);
294 for(i=0;i<IslNF;i++)
295 { AuxGet_vn_uc(FI,&vI,&nuI,ucI); uc = &ucSL[ 2 + (SLp[slNF++]=slNB) ];
296 while(s<AslNF)
297 { if(!(AmI=vA-vI)) if(!(AmI=nuA-nuI))
298 AmI=RIGHTminusLEFT(ucI,ucA,&nuI);
299 if(AmI<0) /* put A */
300 { uc[-2]=vA; uc[-1]=nuA; for(j=0;j<nuA;j++)uc[j]=ucA[j];
301 slNB+=2+nuA; if( (ms=(*uc%4)) ) {if(ms<3) slNM++;} else slSM++;
302 if((++s)<AslNF) AuxGet_vn_uc(FA,&vA,&nuA,ucA);
303 uc = &ucSL[ 2 + (SLp[slNF++]=slNB) ];
304 }
305 else break;
306 }
307 if((s<AslNF)&&(AmI==0)) /* put I==A */
308 { uc[-2]=vI; uc[-1]=nuI; for(j=0;j<nuI;j++) uc[j]=ucI[j];
309 if((*ucI%4)!=(*ucA%4)) *uc = 3+4*(*uc/4);
310 if( (ms=(*uc%4)) ) {if(ms<3) slNM++;} else slSM++;
311 if((++s)<AslNF) AuxGet_vn_uc(FA,&vA,&nuA,ucA);
312 tnb+=2+nuI;
313 }
314 else /* put I */
315 { uc[-2]=vI; uc[-1]=nuI; for(j=0;j<nuI;j++) uc[j]=ucI[j];
316 if( (ms=(*uc%4)) ) {if(ms<3) slNM++;} else slSM++;
317 } slNB+=2+nuI;
318 }
319 while(s<AslNF) /* put A */
320 { uc = &ucSL[ 2 + (SLp[slNF++]=slNB) ]; slNB+=2+nuA;
321 uc[-2]=vA; uc[-1]=nuA; for(j=0;j<nuA;j++)uc[j]=ucA[j];
322 if( (ms=(*uc%4)) ) {if(ms<3) slNM++;} else slSM++;
323 if((++s)<AslNF) AuxGet_vn_uc(FA,&vA,&nuA,ucA);
324 } assert(tnb+slNB==IslNB+AslNB); /* SL done */
325
326 printf("SL: %dnf %dsm %dnm %db -> ",slNF,slSM,slNM,slNB);
327 if(IslNF)
328 { HIpos=FTELL(FI); assert(HIpos==IslNB); assert(!ferror(FI));
329 fclose(FI); if(!newout) remove(Ifn);
330 } /* SL file done */
331
332 FSEEK(FA,HApos,SEEK_SET); Init_FInfoList(&FIo);
333 FIo.nVmax=max(FIi.nVmax,FIa.nVmax);
334 FIo.NUCmax=max(FIi.NUCmax,FIa.NUCmax); /* Tnb=0; */
335 for(v=d+1;v<=FIo.nVmax;v++) for(nu=1;nu<=FIo.NUCmax;nu++)
336 if( (FIo.NFnum[v][nu]=FIi.NFnum[v][nu]+FIa.NFnum[v][nu]) )
337 { FIo.nNUC[v]++; Oli++; /* Tnb+=FIo.NFnum[v][nu]; */
338 } FIo.nV=0; for(v=d+1;v<=FIo.nVmax;v++) if(FIo.nNUC[v]) FIo.nV++;
339
340 for(v=d+1;v<=FIo.nVmax;v++) if(FIo.nNUC[v])
341 { char vxt[5]; strcpy(vxt,".v"); vxt[2]=v/10+'0'; vxt[3]=v%10+'0';
342 vxt[4]=0; strcpy(Ofx,vxt); strcpy(Ifx,vxt);
343 if(FIi.nNUC[v])
344 { if(!newout) {strcat(Ifx,SAVE_FILE_EXT); assert(!rename(Ofn,Ifn));}
345 if(NULL==(FI=fopen(Ifn,"rb"))){printf("Ifn %s failed",Ifn);exit(0);}
346 } if(NULL==(FO=fopen(Ofn,"wb"))){printf("Ofn %s failed",Ofn);exit(0);}
347
348 for(nu=1;nu<=FIo.NUCmax;nu++) if(FIo.NFnum[v][nu])
349 { unsigned int I_NF=FIi.NFnum[v][nu], A_NF=FIa.NFnum[v][nu], O_NF=0;
350 UPint neq=0, peq=0, pa=0; Along pi=0, po=0;
351
352 a=0; if(0<A_NF) {AuxGet_uc(FA,&nu,ucA); pa += 1+(((*ucA)%4)==3);}
353 for(u=0;u<I_NF;u++)
354 { AuxGet_uc(FI,&nu,ucI); pi += 1+(((*ucI)%4)==3);
355 while(a<A_NF)
356 { AmI=RIGHTminusLEFT(ucI,ucA,&nu);
357 if(AmI<0) /* put A */
358 { po += 1+(((*ucA)%4)==3);
359 AuxPut_hNF(FO,&v,&nu,ucA,&FIo,
360 &slNF,&slSM,&slNM,&slNB,ucSL,SLp); O_NF++; a++;
361 if(a<A_NF){AuxGet_uc(FA,&nu,ucA);pa += 1+(((*ucA)%4)==3);}
362 }
363 else break;
364 }
365 if((a<A_NF)&&(AmI==0)) /* put I==A */
366 { int mm=10*(*ucI%4) + (*ucA%4); switch(mm) {
367 case 00: peq++; break; case 11: peq++; break;
368 case 12: *ucI+=2; break; case 13: *ucI+=2;peq++; break;
369 case 21: *ucI+=1; break; case 22: peq++; break;
370 case 23: *ucI+=1;peq++; break; case 31: peq++; break;
371 case 32: peq++; break; case 33: peq+=2; break;
372 default: puts("inconsistens mirror flags");exit(0);}
373 AuxPut_hNF(FO,&v,&nu,ucI,&FIo,&slNF,&slSM,&slNM,&slNB,ucSL,
374 SLp); po += 1+(((*ucI)%4)==3); neq++; a++;
375 if(a<A_NF) {AuxGet_uc(FA,&nu,ucA); pa += 1+(((*ucA)%4)==3);}
376 }
377 else
378 { AuxPut_hNF(FO,&v,&nu,ucI,&FIo,&slNF,
379 &slSM,&slNM,&slNB,ucSL,SLp); po += 1+(((*ucI)%4)==3);
380 } O_NF++;
381 }
382 while(a<A_NF) /* put A */
383 { O_NF++; po += 1+(((*ucA)%4)==3);
384 AuxPut_hNF(FO,&v,&nu,ucA,&FIo,&slNF,&slSM,&slNM,&slNB,ucSL,SLp);
385 ++a; if(a<A_NF) {AuxGet_uc(FA,&nu,ucA); pa += 1+(((*ucA)%4)==3);}
386 } assert(pi+pa==peq+po); /* checksum(v,nu) */
387 FIo.NFnum[v][nu]=O_NF; FIo.nNF+=O_NF; FIo.NB+=O_NF*nu;
388 assert(O_NF+neq==I_NF+A_NF);
389 /*
390 {static int list;printf("#%d v=%d nu=%d Inf=%d Anf=%d ",++list,v,nu,
391 I_NF,A_NF);printf("Onf=%d pi=%d pa=%d po=%d\n",O_NF,pi,pa,po);
392 }*/
393 }
394 if(FIi.nNUC[v])
395 { assert(!ferror(FI)); fclose(FI); if(!newout) remove(Ifn);
396 } assert(!ferror(FO)); fclose(FO);
397 } tnb=0;
398
399 if(slNF)
400 { strcpy(Ofx,".sl"); assert(NULL != (FO=fopen(Ofn,"wb")));
401 for(i=0;i<slNF;i++) /* write SL */
402 { uc=&ucSL[SLp[i]+2]; v=uc[-2]; nu=uc[-1]; tnb+=nu+2;
403 assert(uc[-2]<VERT_Nmax); fputc(uc[-2],FO); slNP+=1+(((*uc)%4)==3);
404 fputc(nu,FO); for(s=0;s<nu;s++) fputc(uc[s],FO);
405 } assert(tnb==slNB);
406 assert(slNP==2*slNF-slNM-slSM);
407 assert(!ferror(FO)); fclose(FO);
408 }
409 printf("\nd=%d v%d v<=%d n<=%d vn%d %lld %d %lld %lld %d %d %d %d\n",
410 d, FIo.nV, FIo.nVmax, FIo.NUCmax,Oli,
411 FIo.nNF,FIo.nSM,FIo.nNM,FIo.NB, slNF, slSM, slNM, slNB);
412
413 strcpy(Ofx,".info"); assert(NULL != (FO=fopen(Ofn,"w")));
414 fprintf(FO, /* write FO.info */
415 "%d %d %d %d %d %lld %d %lld %lld %d %d %d %d\n\n",
416 d, FIo.nV, FIo.nVmax, FIo.NUCmax,Oli,
417 FIo.nNF,FIo.nSM,FIo.nNM,FIo.NB, slNF, slSM, slNM, slNB);
418
419 for(v=d+1;v<=FIo.nVmax;v++) if(FIo.nNUC[v]) /* honest info */
420 { i=0; fprintf(FO,"%d %d\n",v,FIo.nNUC[v]); /* v #nuc's */
421 for(nu=1;nu<=FIo.NUCmax;nu++) if(FIo.NFnum[v][nu])
422 { fprintf(FO,"%d %d%s",nu,FIo.NFnum[v][nu], /* nuc #NF(v,nuc) */
423 (++i<FIo.nNUC[v]) ? " " : "\n");
424 }
425 }
426 printf("Writing %s: %lld+%dsl %lldm+%ds %lldb",dbo,
427 2*FIo.nNF-FIo.nNM-FIo.nSM,
428 slNP,/*Tnb=*/ FIo.nNF-FIo.nNM-FIo.nSM,FIo.nSM,FIo.NB+slNB);
429 /* if(tnb>99)
430 { long long tnp=(2*FIo.nNF-FIo.nNM-FIo.nSM)/10; tnp*=tnp;
431 tnp/=20; tnp/=Tnb; printf(" [p^2/2m=%ldk]",tnp);
432 }
433 */ Print_Expect(&FIo);
434 puts("");
435 assert(ferror(FA)==0); fclose(FA);
436 assert(ferror(FO)==0); fclose(FO);
437 }
Check_sl_order(int * v,int * nu,unsigned char * uc)438 int Check_sl_order(int *v, int *nu, unsigned char *uc)
439 { static int n, V, NU; static unsigned char UC[NUC_Nmax]; if(n)
440 { if(V>(*v)) return 0; if((V==(*v))&&(NU>(*nu))) return 0;
441 if((V==(*v))&&(NU==(*nu))&&(RIGHTminusLEFT(UC,uc,&NU)<=0)) return 0;
442 } V=*v; NU=*nu; for(n=0;n<NU;n++) UC[n]=uc[n]; return 1;
443 }
Check_hnf_order(int * v,int * nu,unsigned char * uc)444 int Check_hnf_order(int *v, int *nu, unsigned char *uc)
445 { static int n, V, NU; static unsigned char UC[NUC_Nmax]; if(n)
446 { if(V>(*v)) return 0; if((V==(*v))&&(NU>(*nu))) return 0;
447 if((V==(*v))&&(NU==(*nu))&&(RIGHTminusLEFT(UC,uc,&NU)<=0)) return 0;
448 } V=*v; NU=*nu; for(n=0;n<NU;n++) UC[n]=uc[n]; return 1;
449 }
Print_NF(FILE * F,int * d,int * v,Long NF[POLY_Dmax][VERT_Nmax])450 void Print_NF(FILE *F,int *d, int *v,Long NF[POLY_Dmax][VERT_Nmax])
451 { int i,j; fprintf(F,"%d %d\n",*d,*v); for(i=0;i<*d;i++)
452 for(j=0;j<*v;j++)fprintf(F,"%d%s",(int)NF[i][j],(*v==j+1) ? "\n" : " ");
453 }
Print_Missing_Mirror(int * d,int * v,int * nu,unsigned char * uc,PolyPointList * _P)454 void Print_Missing_Mirror(int *d, int *v, int *nu, unsigned char *uc,
455 PolyPointList *_P)
456 { int I,J,MS; Long NF[POLY_Dmax][VERT_Nmax];
457 VertexNumList V; EqList E;
458 UCnf2vNF(d,v,nu,uc,NF,&MS); MS %= 4; _P->n=*d; _P->np=*v;
459 for(I=0;I<*v;I++) for(J=0;J<*d;J++) _P->x[I][J]=NF[J][I];
460 if(MS==2) Print_NF(outFILE,d,v,NF);
461 else if(MS==1)
462 { IP_Check(_P,&V,&E); Small_Make_Dual(_P,&V,&E);
463 Make_Poly_NF(_P,&V,&E,NF); Print_NF(outFILE,d,&(V.nv),NF);
464 }
465 else {puts("Only use Print_Missing_Mirror for MM!");exit(0);}
466 }
467
468 /* cF->{1::c 2::C (extended output)} cF->{-1::M (missing mirrors)} */
Check_NF_Order(char * polyi,char * dbi,int cF,PolyPointList * _P)469 void Check_NF_Order(char *polyi,char *dbi, int cF, PolyPointList *_P)/* 1=MM */
470 { FILE *F=NULL; FInfoList L; /* time_t Tstart=time(NULL); */
471 unsigned int rd, i, j, list_num, tln=0, tSM=0; int d, nu, v, si;
472 int sl_nNF, sl_SM, sl_NM, sl_NB; Along tNF=0, tNB=0, tNM=0, SLpos, Hpos=0;
473 char *Ifx=NULL, *Ifn = (char *) malloc(1+strlen(dbi)+File_Ext_NCmax);
474 if((*polyi)&&(*dbi)) puts("only give one of -pi FILE or -di FILE");
475 if((*polyi==0)&&(*dbi==0)) puts("I need one of: -pi FILE or -di FILE");
476 if((*polyi==0)+(*dbi==0)!=1) exit(0);
477
478 if(*polyi) /* read INFO PART */
479 { printf("Checking consistency of Aux/InFile %s\n",polyi);
480 F=fopen(polyi,"rb"); if(F==NULL) {puts("File not found");exit(0);}
481 Init_FInfoList(&L); /* start reading the file */
482
483 rd=fgetc(F); if(rd>127) rd=128*(rd-128)+fgetc(F);
484 assert(rd <= MAX_REC_DEPTH); printf("rd=%d: ",rd);
485 for(i=0;i<rd;i++) printf(" %d",fgetc(F));
486 if(rd) puts("");
487
488 d=fgetc(F); L.nV=fgetc(F); L.nVmax=fgetc(F); L.NUCmax=fgetc(F);
489 list_num=fgetUI(F);
490 printf("d=%d nV=%d nVmax=%d NUCmax=%d #lists=%d ",
491 d,L.nV,L.nVmax,L.NUCmax,list_num); fflush(stdout);
492 L.nNF=fgetUI(F); L.nSM=fgetUI(F); L.nNM=fgetUI(F); L.NB=fgetUI(F);
493 sl_nNF=fgetUI(F); sl_SM=fgetUI(F); sl_NM=fgetUI(F); sl_NB=fgetUI(F);
494
495 for(i=0;i<L.nV;i++)
496 { v=fgetc(F); L.nNUC[v]=fgetc(F);
497 if(cF==2) printf("v%dn%d: ",v,L.nNUC[v]);
498 fflush(stdout);
499 for(j=0;j<L.nNUC[v];j++)
500 { L.NFnum[v][nu=fgetc(F)]=fgetUI(F); /* read nuc and #NF(v,nu)*/
501 tNF+=L.NFnum[v][nu];
502 if(cF==2) printf("%d:%d ",nu,L.NFnum[v][nu]);
503 tln++; tNB+=L.NFnum[v][nu]*(Along)nu;
504 } if(cF==2) puts("");
505 } assert( 0 == (unsigned int) (tNB-L.NB) ); L.NB=tNB;
506 if(cF==2) printf("np=%lld+%dsl %lldb %d files\n",2*L.nNF-L.nSM-L.nNM,
507 2*sl_nNF-sl_SM-sl_NM,L.NB+sl_NB,L.nV+1+(sl_nNF>0));
508 if(cF==2) printf("%lldnf %dsm %lldnm %lldb sl: %d %d %d %d\n",
509 L.nNF,L.nSM,L.nNM,L.NB,sl_nNF,sl_SM,sl_NM,sl_NB); fflush(stdout);
510 }
511 if(*dbi)
512 { printf("Checking consistency of DataBase %s:\n",dbi);
513 strcpy(Ifn,dbi); Ifx=&Ifn[strlen(dbi)]; strcpy(Ifx,".info");
514 F=fopen(Ifn,"r"); if(F==NULL) {puts("Info File not found");exit(0);}
515 Init_FInfoList(&L); /* start reading the file */
516 fscanf(F,"%d%d%d%d%d%lld%d%lld %lld %d%d%d%d",
517 &d,&i,&j,&nu,&list_num,&L.nNF,&L.nSM,&L.nNM,&L.NB,
518 &sl_nNF,&sl_SM,&sl_NM,&sl_NB); L.nV=i; L.nVmax=j; L.NUCmax=nu;
519 printf("d=%d nV=%d nVmax=%d NUCmax=%d #lists=%d ",
520 d,L.nV,L.nVmax,L.NUCmax,list_num);
521 printf("np=%lld+%dsl %lldb %d files\n",2*L.nNF-L.nSM-L.nNM,
522 2*sl_nNF-sl_SM-sl_NM,L.NB+sl_NB,L.nV+1+(sl_nNF>0));
523 printf("%lldnf %dsm %lldnm %lldb sl: %d %d %d %d\n",
524 L.nNF,L.nSM,L.nNM,L.NB,sl_nNF,sl_SM,sl_NM,sl_NB); fflush(stdout);
525
526 for(i=0;i<L.nV;i++)
527 { fscanf(F,"%d",&v); fscanf(F,"%d",&j); L.nNUC[v]=j;
528 if(cF==2) printf("v%dn%d: ",v,L.nNUC[v]);
529 fflush(stdout);
530 for(j=0;j<L.nNUC[v];j++)
531 { fscanf(F,"%d",&nu); fscanf(F,"%d",&L.NFnum[v][nu]); tln++;
532 tNF+=L.NFnum[v][nu];
533 if(cF==2) printf("%d:%d ",nu,L.NFnum[v][nu]);
534 } if(cF==2) puts("");
535 } assert(tln==list_num);
536 } printf("#hNF=%lld sum=%lld %s\n",L.nNF,tNF,
537 (tNF==L.nNF) ? "o.k." : "Error");
538 if(tNF!=L.nNF) exit(0); assert(!ferror(F)); tNF=0;
539 if(tln!=list_num){printf("ERROR: #li=%d != %d\n",list_num,tln);exit(0);}
540
541 { /* long long np=2*L.nNF-L.nSM-L.nNM, pp2m=L.nNF-L.nSM-L.nNM;
542 j=2*sl_nNF-sl_SM-sl_NM; */
543 printf("np=%lld+%dsl %lldnf %dsm %lldnm %lldb sl: %d %d %d %d\n",
544 2*L.nNF-L.nSM-L.nNM,2*sl_nNF-sl_SM-sl_NM,L.nNF,L.nSM,L.nNM,L.NB,
545 sl_nNF,sl_SM,sl_NM,sl_NB);
546 fflush(stdout); /* if(pp2m) tln=(np/10)*(np/10)/20/pp2m;else tln=0;*/
547 }
548 printf("sl: "); /* check SL PART */
549 if(*polyi)
550 { Hpos=FTELL(F); FSEEK(F,-sl_NB,SEEK_END); SLpos=FTELL(F);
551 if(SLpos-Hpos==L.NB) printf("NB o.k. ");
552 }
553 else if(sl_NB)
554 { strcpy(Ifx,".sl"); fclose(F);
555 if(NULL==(F=fopen(Ifn,"rb"))){printf("Open %s failed",Ifn);exit(0);}
556 } else puts("no .sl file");
557 for(si=0;si<sl_nNF;si++)
558 { unsigned char uc[NUC_Nmax]; v=fgetc(F); assert(v<=VERT_Nmax);
559 nu=fgetc(F); /* assert(nu<=L.NUCmax); */
560 AuxGet_uc(F,&nu,uc); assert(!ferror(F));
561 if((*uc%4)==0)tSM++;
562 else if((*uc%4)<3)tNM++;
563 Test_ucNF(&d,&v,&nu,uc,_P);
564 assert(Check_sl_order(&v,&nu,uc));
565 } if(sl_NB&&(*dbi))
566 { SLpos=FTELL(F); assert(sl_NB==SLpos); printf("NB o.k. ");
567 } printf("sm=%d=%d nm=%d=%lld",sl_SM,tSM,sl_NM,tNM); fflush(stdout);
568 if((sl_SM!=(int)tSM)||(sl_NM!=(int)tNM)){puts("ERROR!!");exit(0);}
569 tSM=tNM=0; /* if(tln>1)printf(" p^2/2m=%ldkCY",tln); */
570 Print_Expect(&L);
571
572 printf("\nv:"); if(*polyi) FSEEK(F,Hpos,SEEK_SET); /* check h-order */
573 if(cF<0) puts("");
574 for(v=d+1;v<=L.nVmax;v++) /* if(cF) */ if(L.nNUC[v])
575 { Along nbsum=0; if(cF>0){printf(" %d",v);fflush(stdout);} if(*dbi)
576 { char vxt[5]; strcpy(vxt,".v"); vxt[2]=v/10+'0'; vxt[3]=v%10+'0';
577 vxt[4]=0; strcpy(Ifx,vxt); fclose(F);
578 if(NULL==(F=fopen(Ifn,"rb"))){printf("Ifn %s failed",Ifn);exit(0);}
579 }
580 for(nu=1;nu<=L.NUCmax;nu++) if(L.NFnum[v][nu])
581 { unsigned char uc[NUC_Nmax]; for(i=0;i<L.NFnum[v][nu];i++)
582 { AuxGet_uc(F,&nu,uc); Check_hnf_order(&v,&nu,uc);
583 switch(*uc%4)
584 { case 0: tSM++; break; case 1:; case 2:tNM++; case 3:;}
585
586 if(cF<0)if((*uc%4)%3)Print_Missing_Mirror(&d,&v,&nu,uc,_P);
587
588 } nbsum+=nu*L.NFnum[v][nu];
589 }
590 if(*dbi){Hpos=FTELL(F); assert(Hpos==nbsum);} assert(!ferror(F));
591 }
592 printf(" sm=%d nm=%lld",tSM,tNM); printf(" order o.k.");
593 if((L.nSM!=(int)tSM)||(L.nNM!=tNM)) puts(" CheckSum ERROR!");
594 if(cF<0) return;
595
596 printf("\nv:"); if(*polyi) FSEEK(F,Hpos,SEEK_SET); /* check h-NF */
597 for(v=d+1;v<=L.nVmax;v++) if(L.nNUC[v])
598 { int nbsum=0; printf(" %d",v);fflush(stdout); if(*dbi)
599 { char vxt[5]; strcpy(vxt,".v"); vxt[2]=v/10+'0'; vxt[3]=v%10+'0';
600 vxt[4]=0; strcpy(Ifx,vxt); fclose(F);
601 if(NULL==(F=fopen(Ifn,"rb"))){printf("Ifn %s failed",Ifn);exit(0);}
602 }
603 for(nu=1;nu<=L.NUCmax;nu++) if(L.NFnum[v][nu])
604 { unsigned char uc[NUC_Nmax]; for(i=0;i<L.NFnum[v][nu];i++)
605 { AuxGet_uc(F,&nu,uc); Test_ucNF(&d,&v,&nu,uc,_P);
606 switch(*uc%4)
607 { case 0: tSM++; break; case 1:; case 2:tNM++; case 3:;}
608 } nbsum+=nu*L.NFnum[v][nu];
609 }
610 if(*dbi){Hpos=FTELL(F); assert(Hpos==nbsum);} assert(!ferror(F));
611 } printf(" NF o.k.\n"); assert(!ferror(F)); fclose(F);
612 }
613
614 #define SUBTRACT_H_FROM_SL
615
616 /* A=h1&h2 B=h1&s2 C=s1&h2 D=s1&s2 1n=1-2 2n=2-1 1i=1-1n 2i=2-2n
617 * 1n=(h1-A, s1-C-D) 2n=(h2-A, s2-B-D) 1i=(A,C+D) 2i=(A,B+D)
618
619 ln -s ../zzu.47 zzu.1 ; ln -s ../zzu.58 zzu.2
620 make && class.x -pi zzu.1 -ps zzu.2 -po zzu.1n
621 make && class.x -pi zzu.2 -ps zzu.1 -po zzu.2n
622 class.x -pi zzu.1 -ps zzu.1n -po zzu.1i
623 class.x -pi zzu.2 -ps zzu.2n -po zzu.2i
624
625 class.x -pi zzu.1i -ps zzu.2n -po zzu.AD1
626 class.x -pi zzu.2i -ps zzu.1n -po zzu.AD2
627 class.x -pi zzu.1i -ps zzu.2i -po zzu.0C
628 class.x -pi zzu.2i -ps zzu.1i -po zzu.0B
629 * 1 + 2 = 1n ^ 2n ^ AD */
630
Reduce_Aux_File(char * polyi,char * polys,char * dbsub,char * polyo)631 void Reduce_Aux_File(char *polyi,char *polys,char *dbsub,char *polyo)
632 { FILE *FI=fopen(polyi,"rb"), *FS, *FO=fopen(polyo,"wb");
633 FInfoList FIi, FIs, FIo;
634 Along Ipos,Spos=00, HIpos,HSpos=00, Opos,HOpos, IslNB, SslNB, tnb=0;
635 unsigned char ucI[NUC_Nmax], ucS[NUC_Nmax], *ucSL, *uc; int SLp[SL_Nmax];
636 int IslNF, IslSM, IslNM, db=0, vI, nuI, j, i; unsigned Ili, u; UPint Inp;
637 int SslNF, SslSM, SslNM, dv=0, vS, nuS, s, d; unsigned Sli; Along Snp;
638 int slNF=0, slSM=0, slNM=0, slNB=0, slNP=0, SmI=00, nu,ms, v; UPint Oli=0;
639 char *Sfx=NULL, *Sfn = (*polys) ? (char *) NULL :
640 (char *) malloc(1+strlen(dbsub)+File_Ext_NCmax);
641 Along HIPli[VERT_Nmax][NUC_Nmax], HSPli[VERT_Nmax][NUC_Nmax];
642
643 if((*polys==0) != (*dbsub==0)) db=(*dbsub!=0); else
644 { printf("Need ONE of ps=%s and ds=%s\n",polys,dbsub); exit(0); }
645 if(!*polyi||!*polyo) {
646 puts("With -ps or -ds you have to specify I/O files via -pi and -po");
647 exit(0);}
648
649 if(NULL==FI) {printf("Cannot open %s",polyi);exit(0);}
650 if(NULL==FO) {printf("Cannot open %s",polyo);exit(0);}
651 ucSL = (unsigned char *) malloc( SL_Nmax * CperR_MAX * sizeof(char) );
652 assert(ucSL!=NULL); Init_FInfoList(&FIs);
653
654 if(db)
655 { unsigned tln=0; Along tNF=0; strcpy(Sfn,dbsub);
656 Sfx=&Sfn[strlen(dbsub)]; strcpy(Sfx,".info"); FS=fopen(Sfn,"r");
657 if(FS==NULL) {puts("Info File not found");exit(0);}
658 {Along sNF,sNM; /* start reading the file */
659 fscanf(FS,"%d%d%d%d%d%lld%d%lld %lld %d%d%d%lld",&d,&i,&j,&nu,&Sli,
660 &sNF,&FIs.nSM,&sNM,&FIs.NB,&SslNF,&SslSM,&SslNM,&SslNB);
661 FIs.nNF=sNF; FIs.nNM=sNM; } FIs.nV=i; FIs.nVmax=j; FIs.NUCmax=nu;
662 for(i=0;i<FIs.nV;i++)
663 { fscanf(FS,"%d",&v); fscanf(FS,"%d",&j); FIs.nNUC[v]=j;
664 for(j=0;j<FIs.nNUC[v];j++)
665 { fscanf(FS,"%d",&nu);fscanf(FS,"%d",&FIs.NFnum[v][nu]);
666 tln++; tNF+=FIs.NFnum[v][nu];
667 }
668 } assert(tln==Sli); s=d;
669 assert(tNF==FIs.nNF); assert(!ferror(FS)); tNF=0;
670 if(tln!=Sli){printf("ERROR: #li=%d != %d\n",Sli,tln);exit(0);}
671 }
672 else
673 { if(NULL==(FS=fopen(polys,"rb")))
674 {printf("Cannot open %s",polys);exit(0);}
675 if(fgetc(FS)) {puts("don't subtract aux files!");exit(0);}
676 Read_Bin_Info(FS,&s,&Sli,&SslNF,&SslSM,&SslNM,&SslNB,&FIs);
677 HSpos=FTELL(FS); FSEEK(FS,0,SEEK_END); Spos=FTELL(FS);
678 assert(HSpos+FIs.NB+SslNB==Spos); FSEEK(FS,-SslNB,SEEK_CUR);
679 }
680 d=fgetc(FI); if(d>127) d=128*(d-128)+fgetc(FI); Init_FInfoList(&FIi);
681 if(d) printf("RecursionDepth=%d on %s\n",d,polyi);
682 if(d < 128) fputc(d,FO); else {fputc(128+(d)/128,FO); fputc(d%128,FO);}
683 for(i=0;i<d;i++) {j=fgetc(FI); fputc(j,FO);}
684 Read_Bin_Info(FI,&d,&Ili,&IslNF,&IslSM,&IslNM,&IslNB,&FIi); assert(d==s);
685 HIpos=FTELL(FI); FSEEK(FI,0,SEEK_END); Ipos=FTELL(FI);
686 Inp=2*FIi.nNF-FIi.nSM-FIi.nNM; Snp=2*FIs.nNF-FIs.nSM-FIs.nNM;
687
688 printf("Data on %s: %d+%dsl %lldb (%dd)\n", polyi,Inp,
689 2*IslNF-IslNM-IslSM,Ipos,s);
690 printf("Data on %s: %lld+%dsl %lldb (%dd)\n", db ? dbsub : polys,Snp,
691 2*SslNF-SslNM-SslSM,db ? FIs.NB+SslNB : Spos,d);
692 assert(HIpos+FIi.NB+IslNB==Ipos); FSEEK(FI,-IslNB,SEEK_CUR);
693
694 if(db && SslNF)
695 { strcpy(Sfx,".sl"); fclose(FS);
696 if(NULL==(FS=fopen(Sfn,"rb"))){printf("Open %s failed",Sfn);exit(0);}
697 }
698 s=0; if(0<SslNF) AuxGet_vn_uc(FS,&vS,&nuS,ucS);
699 for(v=0;v<IslNF;v++)
700 { AuxGet_vn_uc(FI,&vI,&nuI,ucI); uc = &ucSL[ 2 + (SLp[slNF++]=slNB) ];
701 while(s<SslNF)
702 { if(!(SmI=vS-vI)) if(!(SmI=nuS-nuI))
703 SmI=RIGHTminusLEFT(ucI,ucS,&nuI);
704 if(SmI<0) /* next S */
705 { if((++s)<SslNF) AuxGet_vn_uc(FS,&vS,&nuS,ucS);
706 }
707 else break;
708 }
709 if((s<SslNF)&&(SmI==0)) /* I==S */
710 { ms=(*ucS%4); if( (ms%3) && ((*ucI%4)!=ms) ) /* put I==S */
711 { uc[-2]=vI; uc[-1]=nuI; for(nu=0;nu<nuI;nu++) uc[nu]=ucI[nu];
712 *uc=(3-ms)+4*(*uc/4); slNM++; slNB+=2+nuI; tnb+=2+nuI;
713 }
714 else slNF--;
715 if((++s)<SslNF) AuxGet_vn_uc(FS,&vS,&nuS,ucS);
716 }
717 else /* put I */
718 { uc[-2]=vI; uc[-1]=nuI; for(i=0;i<nuI;i++) uc[i]=ucI[i];
719 if( (ms=(*uc%4)) ) {if(ms<3) slNM++;} else slSM++; slNB+=2+nuI;
720 }
721 } /* assert(tnb+slNB==IslNB+AslNB); */ /* SL done */
722
723 printf("SL: %dnf %dsm %dnm %db -> ",slNF,slSM,slNM,slNB); tnb=0;
724 fflush(stdout);
725
726 if(!db) FSEEK(FS,HSpos,SEEK_SET); FSEEK(FI,HIpos,SEEK_SET); Opos=HIpos;
727
728 for(v=d+1;v<=FIi.nVmax;v++) if(FIi.nNUC[v]) /* init HIPli */
729 for(nu=1;nu<=FIi.NUCmax;nu++) if(FIi.NFnum[v][nu])
730 { HIPli[v][nu]=Opos; Opos+=nu*FIi.NFnum[v][nu];
731 } assert(Opos==HIpos+FIi.NB); Opos=HSpos;
732 for(v=d+1;v<=FIs.nVmax;v++) if(FIs.nNUC[v]) /* init HSPli */
733 for(nu=1;nu<=FIs.NUCmax;nu++) if(FIs.NFnum[v][nu])
734 { if(db) if(v>dv) {Opos=0; dv=v;}
735 HSPli[v][nu]=Opos; Opos+=nu*FIs.NFnum[v][nu];
736 } if(!db) assert(Opos==HSpos+FIs.NB); dv=0; Init_FInfoList(&FIo);
737
738 for(v=d+1;v<=FIi.nVmax;v++) if(FIi.nNUC[v]) /* FIo.NFnum[v][nu]?=0 */
739 for(nu=1;nu<=FIi.NUCmax;nu++) if(FIi.NFnum[v][nu])
740 { if(FIs.NFnum[v][nu]) /* check for add polys */
741 { volatile unsigned int *_nf=&FIo.NFnum[v][nu];
742 /* (*_nf)=#NFnum(v,n) ?= 0 */
743 unsigned int n, I_NF=FIi.NFnum[v][nu], S_NF=FIs.NFnum[v][nu];
744 if(db) if(v>dv)
745 { char vxt[5]; strcpy(vxt,".v"); vxt[2]=v/10+'0';
746 vxt[3]=v%10+'0'; vxt[4]=0; strcpy(Sfx,vxt);
747 assert(!ferror(FS));fclose(FS); if(NULL==(FS=fopen(Sfn,"rb")))
748 {printf("%s open failed",Sfn);exit(0);} dv=v;
749 }
750 FSEEK(FI,HIPli[v][nu],SEEK_SET); FSEEK(FS,HSPli[v][nu],SEEK_SET);
751 u=0; if(0<S_NF) AuxGet_uc(FS,&nu,ucS);
752 for(n=0;n<I_NF;n++)
753 { AuxGet_uc(FI,&nu,ucI);
754 while(u<S_NF)
755 { SmI=RIGHTminusLEFT(ucI,ucS,&nu);
756 if(SmI<0)
757 { if((++u)<S_NF)AuxGet_uc(FS,&nu,ucS); /* get S */
758 }
759 else break;
760 }
761 if((u<S_NF)&&(SmI==0)) /* found I==S */
762 { ms=(*ucS)%4; if((ms%3) && (((*ucI)%4)!=ms)) (*_nf)++;
763 if((++u)<S_NF)AuxGet_uc(FS,&nu,ucS);
764 }
765 else (*_nf)++;
766 if(*_nf) break;
767 }
768 }
769 else FIo.NFnum[v][nu]=FIi.NFnum[v][nu];
770 if(FIo.NFnum[v][nu])
771 { FIo.nVmax=max(FIo.nVmax,v); FIo.NUCmax=max(FIo.NUCmax,nu);
772 FIo.nNUC[v]++; Oli++;
773 }
774 } dv=0;
775 for(v=d+1;v<=FIo.nVmax;v++) if(FIo.nNUC[v]) FIo.nV++;
776
777 fputc(d,FO); fputc(FIo.nV,FO); fputc(FIo.nVmax,FO);
778 fputc(FIo.NUCmax,FO); fputUI(Oli,FO); Opos=FTELL(FO);
779 j=InfoSize(0,Oli,&FIo)-5-sizeof(int); for(i=0;i<j;i++) fputc(0,FO);
780
781 HOpos=FTELL(FO); /* printf("Opos=%d HOpos=%d\n\n\n",Opos,HOpos); */
782
783 for(v=d+1;v<=FIo.nVmax;v++) if(FIo.nNUC[v])
784 for(nu=1;nu<=FIo.NUCmax;nu++) if(FIo.NFnum[v][nu])
785 { unsigned I_NF=FIi.NFnum[v][nu], S_NF=FIs.NFnum[v][nu], O_NF=0, n;
786 UPint neq=0, peq=0, pi=0, po=0;
787 if(v>tnb)
788 { printf(" %d",v); tnb=v; fflush(stdout); if(FIs.nNUC[v]) if(db)
789 { char vxt[5]; strcpy(vxt,".v"); vxt[2]=v/10+'0';
790 vxt[3]=v%10+'0'; vxt[4]=0; strcpy(Sfx,vxt);
791 assert(!ferror(FS));fclose(FS); if(NULL==(FS=fopen(Sfn,"rb")))
792 {printf("%s open failed",Sfn);exit(0);}
793 }
794 }
795 FSEEK(FI,HIPli[v][nu],SEEK_SET); FSEEK(FS,HSPli[v][nu],SEEK_SET);
796 u=0; if(u<S_NF) AuxGet_uc(FS,&nu,ucS);
797 for(n=0;n<I_NF;n++)
798 { AuxGet_uc(FI,&nu,ucI); pi += 1+(((*ucI)%4)==3);
799 while(u<S_NF)
800 { SmI=RIGHTminusLEFT(ucI,ucS,&nu);
801 if(SmI<0) /* get S */
802 { u++; if(u<S_NF) AuxGet_uc(FS,&nu,ucS);
803 }
804 else break;
805 }
806 if((u<S_NF)&&(SmI==0)) /* put I-S */
807 { int k, mm=10*(*ucI%4) + (*ucS%4); switch(mm) {
808 case 33: peq++; case 22:; case 11:; case 00: peq++; break;
809 case 31:; case 32: *ucI-=(*ucS%4); peq++;
810 case 12:; case 21: neq--; po++; O_NF++;
811 FIo.nNM++; for(k=0;k<nu;k++) fputc(ucI[k],FO); break;
812 case 13:; case 23: peq++; break;
813 default: puts("inconsistens mirror flags");exit(0);}
814 neq++; u++; if(u<S_NF) AuxGet_uc(FS,&nu,ucS);
815 }
816 else
817 { int k;for(k=0;k<nu;k++)fputc(ucI[k],FO);O_NF++;k=(*ucI)%4;
818 po += 1+(k==3); if(k) {if(k<3)FIo.nNM++;} else FIo.nSM++;
819 }
820 } assert(pi-peq==po); assert(O_NF+neq==I_NF); /* checksum(v,nu) */
821 FIo.NFnum[v][nu]=O_NF; FIo.nNF+=O_NF; FIo.NB+=O_NF*nu;
822 } tnb=0;
823
824 #ifdef SUBTRACT_H_FROM_SL
825 uc = & (ucSL[SLp[i=0]]); dv=0;
826 for(v=d+1;v<=FIs.nVmax;v++) if(FIs.nNUC[v]) /* subtract S from SL */
827 for(nu=1;nu<=FIs.NUCmax;nu++) if(FIs.NFnum[v][nu])
828 { while( (uc[0]<v) || ((uc[0]==v)&&(uc[1]<nu)))
829 { uc = & (ucSL[SLp[++i]]); if(i>=slNF) goto END_SL;
830 } /* go up to (v,nu) in SL-list */
831 if(db) if(v>dv)
832 { char vxt[5]; strcpy(vxt,".v"); vxt[2]=v/10+'0';
833 vxt[3]=v%10+'0'; vxt[4]=0; strcpy(Sfx,vxt);
834 assert(!ferror(FS));fclose(FS); if(NULL==(FS=fopen(Sfn,"rb")))
835 {printf("%s open failed",Sfn);exit(0);} dv=v;
836 }
837
838 if((uc[0]==v)&&(uc[1]==nu)) /* read from file and sort out */
839 { FSEEK(FS,HSPli[v][nu],SEEK_SET);
840 for(u=0;u<FIs.NFnum[v][nu];u++)
841 { int HmSL; AuxGet_uc(FS,&nu,ucS);/* Test_ucNF(&d,&v,&nu,ucS);*/
842 while(0 < (HmSL=RIGHTminusLEFT(&uc[2],ucS,&nu))) /* next SL */
843 { uc = & (ucSL[SLp[++i]]); if(i>=slNF) goto END_SL;
844 if((uc[0]!=v)||(uc[1]!=nu)) goto END_VN;
845 }
846 if(HmSL==0) /* remove SL */
847 { int k, sms=uc[2]%4, hms=(*ucS)%4;
848 switch(10*hms + sms){
849 case 31: case 32: case 00: case 11: case 22: case 33:
850 for(k=i+1;k<slNF;k++) SLp[k-1]=SLp[k];
851 slNB -= nu+2;
852 if(sms==0) --slSM; else if(sms<3) --slNM;
853 --slNF;
854 if(i>=slNF) goto END_SL;
855 uc = & (ucSL[SLp[i]]);
856 if((uc[0]!=v)||(uc[1]!=nu)) goto END_VN;
857 break;
858 case 12: case 21: break;
859 case 13: case 23: ++slNM; uc[2]-=hms; break;
860 default: puts("inconsistent MS flags in SL-H"); exit(0);}
861 } /* else (SL>H): hence next H */
862 } assert(!ferror(FS));
863 }
864 else END_VN: assert( (uc[0]>v) || ((uc[0]==v)&&(uc[1]>nu)) );
865 } END_SL: ;
866 #endif
867
868 for(i=0;i<slNF;i++) /* write SL */
869 { uc=&ucSL[SLp[i]+2]; v=uc[-2]; nu=uc[-1]; tnb+=nu+2;
870 /* printf("#%d: SLp=%d v=%d nu=%d\n",i,SLp[i],v,nu); */
871 assert(uc[-2]<VERT_Nmax); fputc(uc[-2],FO); slNP+=1+(((*uc)%4)==3);
872 fputc(nu,FO); for(s=0;s<nu;s++) fputc(uc[s],FO);
873 } assert(tnb==slNB);
874 assert(slNP==2*slNF-slNM-slSM);
875 printf("\nd=%d v%d v<=%d n<=%d vn%d %lld %d %lld %lld %d %d %d %d\n",
876 d, FIo.nV, FIo.nVmax, FIo.NUCmax,Oli,
877 FIo.nNF,FIo.nSM,FIo.nNM,FIo.NB, slNF, slSM, slNM, slNB);
878
879 FSEEK(FO,Opos,SEEK_SET);
880 fputUI(FIo.nNF,FO); /* write info */
881 fputUI(FIo.nSM,FO); fputUI(FIo.nNM,FO); fputUI(FIo.NB,FO);
882 fputUI(slNF,FO); fputUI(slSM,FO); fputUI(slNM,FO); fputUI(slNB,FO);
883 for(v=d+1;v<=FIo.nVmax;v++) if(FIo.nNUC[v])
884 { i=0; for(nu=1;nu<=FIo.NUCmax;nu++) if(FIo.NFnum[v][nu]) i++;
885 fputc(v,FO); fputc(i,FO); /* v #nuc(v) */
886 for(nu=1;nu<=FIo.NUCmax;nu++) if(FIo.NFnum[v][nu])
887 { fputc(nu,FO); fputUI(FIo.NFnum[v][nu],FO); /* nuc #NF */
888 }
889 }
890
891 printf("Writing %s: %lld+%dsl %lldm+%ds %lldb",polyo,2*FIo.nNF-FIo.nNM
892 -FIo.nSM,slNP,tnb=FIo.nNF-FIo.nNM-FIo.nSM,FIo.nSM,FIo.NB+slNB);
893 /* if(tnb>99)
894 { long long tnp=(2*FIo.nNF-FIo.nNM-FIo.nSM)/1000; tnp*=tnp;
895 tnp/=(2*tnb); printf(" [p^2/2m=%ldM]",tnp);
896 } */
897 Print_Expect(&FIo);
898 puts("");
899 assert(ferror(FI)==0); fclose(FI);
900 assert(ferror(FS)==0); fclose(FS);
901 assert(HOpos==FTELL(FO)); assert(ferror(FO)==0); fclose(FO);
902 }
Bin2a(char * polyi,int max,PolyPointList * _P)903 void Bin2a(char *polyi, int max, PolyPointList *_P)
904 { FILE *F=fopen(polyi,"rb"); FInfoList L; UPint list_num,tNF=0; Along tNB=0;
905 int d, v, s, sl_nNF, sl_SM, sl_NM, sl_NB, mc=0, MS, nu; unsigned i, j;
906 unsigned char uc[POLY_Dmax*VERT_Nmax]; VertexNumList V; EqList E;
907 Long NF[POLY_Dmax][VERT_Nmax]; Init_FInfoList(&L);
908
909 if(F==NULL) {printf("Input file %s not found\n",polyi); exit(0);}
910 d=fgetc(F); assert(d==0); /* for(i=0;i<d;i++) fgetc(F); */
911 d=fgetc(F); L.nV=fgetc(F); L.nVmax=fgetc(F); L.NUCmax=fgetc(F);
912 list_num=fgetUI(F);
913 L.nNF=fgetUI(F); L.nSM=fgetUI(F); L.nNM=fgetUI(F); L.NB=fgetUI(F);
914 sl_nNF=fgetUI(F); sl_SM=fgetUI(F); sl_NM=fgetUI(F); sl_NB=fgetUI(F);
915
916 for(i=0;i<L.nV;i++)
917 { v=fgetc(F); L.nNUC[v]=fgetc(F); /* read #nuc's per #Vert */
918 for(j=0;j<L.nNUC[v];j++)
919 { L.NFnum[v][nu=fgetc(F)]=fgetUI(F); /* read nuc and #NF(v,nu)*/
920 tNF+=L.NFnum[v][nu]; tNB+=L.NFnum[v][nu]*(Along)nu;
921 }
922 } assert( 0 == (unsigned int)(tNB-L.NB) ); L.NB=tNB; assert(tNF==L.nNF);
923
924 for(v=d+1;v<=L.nVmax;v++) if(L.nNUC[v]) /* write honest polys */
925 { int I,J; for(nu=1;nu<=L.NUCmax;nu++) for(j=0;j<L.NFnum[v][nu];j++)
926 { for(s=0;s<nu;s++) uc[s]=fgetc(F);
927 UCnf2vNF(&d,&v,&nu,uc,NF,&MS); MS %= 4; _P->n=d; _P->np=v;
928 for(I=0;I<v;I++) for(J=0;J<d;J++) _P->x[I][J]=NF[J][I];
929 assert(Ref_Check(_P,&V,&E));
930
931 if(MS!=2) /* if(MS!=2) print NF */
932 if(!max || Poly_Max_check(_P,&V,&E))
933 { mc++; Print_NF(outFILE,&d,&v,NF);}
934 if(MS>1) /* if(MS>1); print Mirror */
935 if(!max || Poly_Min_check(_P,&V,&E))
936 { mc++; Small_Make_Dual(_P,&V,&E);
937 Make_Poly_NF(_P,&V,&E,NF); Print_NF(outFILE,&d,&(V.nv),NF);
938 }
939 }
940 }
941 printf("np=%lld+%dsl ",2*L.nNF-L.nSM-L.nNM,2*sl_nNF-sl_SM-sl_NM);
942 printf( /* write Finfo */
943 "%dd %dv<=%d n<=%d %dnv %lld %d %lld %lld %d %d %d %d",
944 d, L.nV, L.nVmax, L.NUCmax,list_num,
945 L.nNF,L.nSM,L.nNM,L.NB, sl_nNF, sl_SM, sl_NM, sl_NB);
946 if(max) printf(" r-max=%d",mc);
947 puts("");
948 }
DB_fromVF_toVT(DataBase * DB,int vf,int vt)949 void DB_fromVF_toVT(DataBase *DB,int vf,int vt)
950 { int v,n; DB->nNF=0;
951 for(v=vf;v<=vt;v++) for(n=0;n<NUC_Nmax;n++) DB->nNF+=DB->NFnum[v][n];
952 if(!DB->nNF){fprintf(stderr,"No NF with %d<=v<=%d\n",vf,vt);exit(0);}
953 while(0==DB->nNUC[vt]) {vt--;assert(vf<=vt);} DB->v=vf; DB->nVmax=vt;
954 }
Bin2aDBsl(char * dbi,int max,int vf,int vt,PolyPointList * _P)955 void Bin2aDBsl(char *dbi, int max, int vf, int vt, PolyPointList *_P)
956 { FILE *F; FInfoList L;
957 int d, v, nu, i, j, list_num, mc=0, MS,
958 sl_nNF, sl_SM, sl_NM, sl_NB, tSM=0, tNM=0;
959 char *Ifx, *Ifn = (char *) malloc(1+strlen(dbi)+File_Ext_NCmax);
960 Long NF[POLY_Dmax][VERT_Nmax];
961 VertexNumList V;EqList E;
962 strcpy(Ifn,dbi); Ifx=&Ifn[strlen(dbi)]; strcpy(Ifx,".info");
963 F=fopen(Ifn,"r"); if(F==NULL) {puts("Info File not found");exit(0);}
964 Init_FInfoList(&L); /* start reading the file */
965 fscanf(F,"%d%d%d%d%d%lld%d%lld %lld %d%d%d%d",
966 &d,&i,&j,&nu,&list_num,&L.nNF,&L.nSM,&L.nNM,&L.NB,
967 &sl_nNF,&sl_SM,&sl_NM,&sl_NB); L.nV=i; L.nVmax=j; L.NUCmax=nu;
968 if(sl_NB && (vf==2)&&(vt==VERT_Nmax-1))
969 { strcpy(Ifx,".sl"); fclose(F);
970 if(NULL==(F=fopen(Ifn,"rb"))){printf("Open %s failed",Ifn);exit(0);}
971 }
972 else /* puts("no .sl file"); */
973 { DataBase *DB; Open_DB(dbi, &DB, 0); if((DB->v<vf)||(DB->nVmax>vt))
974 DB_fromVF_toVT(DB,vf,vt); /* read only vf <= v <= vt :: */
975 for(i=0; Read_H_poly_from_DB(DB,_P); i++) Print_PPL(_P,"");
976 printf("#poly=%d\n",i); Close_DB(DB);
977 }
978 for(i=0;i<sl_nNF;i++)
979 { int I,J; unsigned char uc[NUC_Nmax]; v=fgetc(F); assert(v<=VERT_Nmax);
980 nu=fgetc(F); AuxGet_uc(F,&nu,uc); assert(!ferror(F));
981 if((*uc%4)==0)tSM++; else if((*uc%4)<3)tNM++;
982 Test_ucNF(&d,&v,&nu,uc,_P);
983
984 UCnf2vNF(&d,&v,&nu,uc,NF,&MS); MS %= 4; _P->n=d; _P->np=v;
985 for(I=0;I<v;I++) for(J=0;J<d;J++) _P->x[I][J]=NF[J][I];
986 assert(Ref_Check(_P,&V,&E));
987 if(MS!=2) /* if(MS!=2) print NF */
988 if(!max || Poly_Max_check(_P,&V,&E))
989 { mc++; Print_NF(outFILE,&d,&v,NF);}
990 if(MS>1) /* if(MS>1); print Mirror */
991 if(!max || Poly_Min_check(_P,&V,&E))
992 { mc++; Small_Make_Dual(_P,&V,&E);
993 Make_Poly_NF(_P,&V,&E,NF); Print_NF(outFILE,&d,&(V.nv),NF);
994 }
995 }
996 printf("np=%lld+%dsl ",2*L.nNF-L.nSM-L.nNM,2*sl_nNF-sl_SM-sl_NM);
997 printf( /* write Finfo */
998 "%dd %dv<=%d n<=%d %dnv ... %d %d %d %d", /* %d %d %d %lld */
999 d, L.nV, L.nVmax, L.NUCmax,list_num,
1000 /* L.nNF,L.nSM,L.nNM,L.NB, */ sl_nNF, sl_SM, sl_NM, sl_NB);
1001 if(max) printf(" r-max=%d",mc);
1002 puts("");
1003 }
Bin_2_ascii(char * polyi,char * dbin,int max,int vf,int vt,PolyPointList * P)1004 void Bin_2_ascii(char *polyi,char *dbin,int max,int vf,int vt,PolyPointList *P)
1005 { if(*polyi) Bin2a(polyi, max, P);
1006 else if(*dbin) Bin2aDBsl(dbin, max, vf, vt, P);
1007 else puts("With -b[2a] you have to specify input via -pi or -di");
1008 }
1009
1010 /* ====================================================================== */
1011 /* ========== ========== */
1012 /* ========== Hodge-database routines ========== */
1013 /* ========== ========== */
1014 /* ====================================================================== */
1015
1016 #define Hod_Dif_max 480
1017 #define Hod_Min_max 251
1018
1019 #if (POLY_Dmax < 6)
DB_to_Hodge(char * dbin,char * dbout,int vfrom,int vto,PolyPointList * _P)1020 void DB_to_Hodge(char *dbin, char *dbout, int vfrom, int vto,
1021 PolyPointList *_P){
1022
1023 /* Read the database, write the Hodge numbers */
1024
1025 DataBase DB;
1026 VertexNumList V;
1027 Long VPM[EQUA_Nmax][VERT_Nmax];
1028 static EqList E;
1029 time_t Tstart;
1030 char *dbname = (char *) malloc(1+strlen(dbin)+File_Ext_NCmax), *fx;
1031 char *dbhname = (char *) malloc(6+strlen(dbout)+File_Ext_NCmax), *fhx;
1032 unsigned char uc_poly[NUC_Nmax];
1033 int d, v, nu, i, j, list_num, sl_nNF, sl_SM, sl_NM, sl_NB, dh,
1034 nnf_vd[VERT_Nmax][Hod_Dif_max+1],
1035 nnf_v[VERT_Nmax];
1036 BaHo BH;
1037 FILE *Faux[Hod_Dif_max+1];
1038 FILE *Fvinfo;
1039 PolyPointList *_PD = (PolyPointList *) malloc(sizeof(PolyPointList));
1040 assert(_PD!=NULL);
1041
1042 if (!*dbin||!*dbout){
1043 puts("You have to specify I/O database names via -di and -do");
1044 exit(0);}
1045
1046 for(i=0;i<=Hod_Dif_max;i++) for(j=0;j<VERT_Nmax;j++) nnf_vd[j][i]=0;
1047 strcpy(dbname,dbin);
1048 strcpy(dbhname,dbout);
1049 strcat(dbname,".info");
1050 strcat(dbhname,".vinfo");
1051 fx=&dbname[strlen(dbin)+1];
1052 fhx=&dbhname[strlen(dbout)+1];
1053 Fvinfo=fopen(dbhname,"a");
1054
1055 printf("Reading %s: ",dbname); fflush(0);
1056
1057 /* read the info-file: */
1058 DB.Finfo=fopen(dbname,"r");
1059 assert(DB.Finfo!=NULL);
1060 fscanf(DB.Finfo, "%d %d %d %d %d %lld %d %lld %lld %d %d %d %d",
1061 &d, &DB.nV, &DB.nVmax, &DB.NUCmax, &list_num,
1062 &DB.nNF, &DB.nSM, &DB.nNM, &DB.NB, &sl_nNF, &sl_SM, &sl_NM, &sl_NB);
1063 printf("%lldp (%dsl) %lldnf %lldb\n",
1064 2*(DB.nNF)-DB.nSM-DB.nNM+2*sl_nNF-sl_SM-sl_NM,
1065 2*sl_nNF-sl_SM-sl_NM, DB.nNF+sl_nNF, DB.NB+sl_NB);
1066
1067 for (v=1;v<VERT_Nmax;v++){
1068 DB.nNUC[v]=0;
1069 for(nu=0; nu<NUC_Nmax; nu++) DB.NFnum[v][nu]=0;}
1070
1071 for(i=0; i<DB.nV; i++){
1072 fscanf(DB.Finfo,"%d",&v);
1073 fscanf(DB.Finfo,"%d",&(DB.nNUC[v]));
1074 for(j=0;j<DB.nNUC[v];j++){
1075 fscanf(DB.Finfo,"%d", &nu);
1076 fscanf(DB.Finfo,"%d", &(DB.NFnum[v][nu])); } }
1077
1078 if(ferror(DB.Finfo)) {printf("File error in %s\n",dbname); exit(0);}
1079 fclose(DB.Finfo);
1080 fflush(stdout);
1081
1082 printf("Reading DB-files, calculating Hodge numbers, writing aux-files:\n");
1083
1084 /* read the DB-files and calculate Hodge numbers */
1085 for (v=vfrom;v<=vto;v++) if(DB.nNUC[v]){
1086 int nd=0;
1087 char ext[4], aext[8];
1088 ext[0]='v'; ext[1]='0' + v / 10; ext[2]='0' + v % 10; ext[3]=0;
1089 aext[0]='v'; aext[1]='0' + v / 10; aext[2]='0' + v % 10; aext[3]='d';
1090 aext[4]=aext[5]=aext[6]=aext[7]=0;
1091 Tstart=time(NULL);
1092 printf("v=%d: ", v); fflush(0);
1093 strcpy(fx,ext);
1094 DB.Fv[v]=fopen(dbname,"rb");
1095 assert(DB.Fv[v]!=NULL);
1096 for (nu=0;nu<=DB.NUCmax;nu++) for (i=0; i<DB.NFnum[v][nu]; i++){
1097 int mirror=0;
1098 int MS;
1099 unsigned char c;
1100 for (j=0; j<nu; j++) uc_poly[j]=fgetc(DB.Fv[v]);
1101 uc_nf_to_P(_P, &MS, &d, &v, &nu, uc_poly);
1102 assert(Ref_Check(_P, &V, &E));
1103 assert(V.nv==v);
1104 Make_VEPM(_P,&V,&E,VPM);
1105 Complete_Poly(VPM,&E,V.nv,_P);
1106 Make_Dual_Poly(_P,&V,&E,_PD);
1107 RC_Calc_BaHo(_P,&V,&E,_PD,&BH);
1108 if (BH.h1[1]<BH.h1[2]) mirror=1;
1109 if (mirror) dh=BH.h1[2]-BH.h1[1];
1110 else dh=BH.h1[1]-BH.h1[2];
1111 if (!nnf_vd[v][dh]||(dh>250)){
1112 aext[4]='0'+dh/100; aext[5]='0'+(dh/10)%10; aext[6]='0'+dh%10;
1113 strcpy(fhx,aext);
1114 Faux[dh]=fopen(dbhname,"ab"); }
1115 if (!nnf_vd[v][dh]) nd++;
1116 c=BH.h1[2-mirror]; fputc(c,Faux[dh]);
1117 c=BH.mp/256+4*BH.mv; fputc(c,Faux[dh]);
1118 c=BH.mp%256; fputc(c,Faux[dh]);
1119 c=BH.np/256+4*BH.nv; fputc(c,Faux[dh]);
1120 c=BH.np%256; fputc(c,Faux[dh]);
1121 c=nu+64*mirror; fputc(c,Faux[dh]);
1122 for (j=0;j<nu;j++) fputc(uc_poly[j],Faux[dh]);
1123 if (dh>250) fclose(Faux[dh]);
1124 (nnf_vd[v][dh])++;
1125 (nnf_v[v])++; }
1126
1127 if(ferror(DB.Fv[v])) {printf("File error in %s\n",dbname); exit(0);}
1128 fclose(DB.Fv[v]);
1129 for (dh=0;dh<=250;dh++) if (nnf_vd[v][dh]) {
1130 if(ferror(Faux[dh])) {printf("File error at dh=%d\n",dh); exit(0);}
1131 fclose(Faux[dh]);}
1132
1133 fprintf(Fvinfo, "%d %d %d\n", v, nd, nnf_v[v]);
1134 for (dh=0;dh<=Hod_Dif_max;dh++) if (nnf_vd[v][dh])
1135 fprintf(Fvinfo, "%d %d ", dh, nnf_vd[v][dh]);
1136 fprintf(Fvinfo, "\n");
1137 printf(" %d NF (%ds)\n", nnf_v[v], (int) difftime(time(NULL),Tstart));
1138 fflush(0); }
1139 free(_PD);
1140 fclose(Fvinfo);
1141 }
1142
Sort_Hodge(char * dbaux,char * dbout)1143 void Sort_Hodge(char *dbaux, char *dbout){
1144 /* Sort from v-chi-format to chi-h12-format */
1145
1146 time_t Tstart;
1147 char *dbaname = (char *) malloc(6+strlen(dbaux)+File_Ext_NCmax), *fax;
1148 char *dbhname =
1149 (char *) malloc(6+strlen(*dbout ? dbout : dbaux)+File_Ext_NCmax), *fhx;
1150 int v, i, j, dh, nd,
1151 nnf_d[Hod_Dif_max+1], nnf_vd[VERT_Nmax][Hod_Dif_max+1],
1152 nnf_h[Hod_Min_max+1], nnf_v[VERT_Nmax];
1153 FILE *Fh[Hod_Min_max+1];
1154 FILE *Fchia, *Fvinfo, *Fhinfo;
1155
1156 for(i=0;i<=Hod_Dif_max;i++){
1157 nnf_d[i]=0;
1158 for(j=0;j<VERT_Nmax;j++) nnf_vd[j][i]=0;}
1159
1160 if(!*dbout) dbout=dbaux;
1161 strcpy(dbaname,dbaux);
1162 strcat(dbaname,".vinfo");
1163 strcpy(dbhname,dbout);
1164 strcat(dbhname,".hinfo");
1165 fhx=&dbhname[strlen(dbout)+1];
1166 fax=&dbaname[strlen(dbaux)+1];
1167
1168 printf("Reading %s\n",dbaname); fflush(0);
1169
1170 /* read the info-file: */
1171 Fvinfo=fopen(dbaname,"r");
1172 assert(Fvinfo!=NULL);
1173 while ((fscanf(Fvinfo, "%d",&v))!=EOF){
1174 fscanf(Fvinfo, "%d %d", &nd, &(nnf_v[v]));
1175 /* printf("%d %d %d ", v, nd, nnf_v[v] ); fflush(0); */
1176 for (i=0;i<nd;i++) {
1177 fscanf(Fvinfo, "%d", &dh);
1178 fscanf(Fvinfo, "%d", &(nnf_vd[v][dh]));
1179 nnf_d[dh]+=nnf_vd[v][dh];}}
1180 if(ferror(Fvinfo)) {printf("File error in %s\n",dbaname); exit(0);}
1181 fclose(Fvinfo);
1182
1183 printf("Sorting and writing Hodge-files:\n"); fflush(0);
1184 Tstart=time(NULL);
1185
1186 Fhinfo=fopen(dbhname,"w");
1187 assert(Fhinfo!=NULL);
1188
1189 /* Sort the Hodge&Poly-Data */
1190 for (dh=0;dh<=Hod_Dif_max;dh++) if (nnf_d[dh]){
1191 char aext[8], hext[9];
1192 int h12, nh=0;
1193 unsigned char c, nuc;
1194 aext[0]='v'; aext[3]='d'; aext[4]='0'+dh/100;
1195 aext[5]='0'+(dh/10)%10; aext[6]='0'+dh%10; aext[7]=0;
1196 hext[0]='d'; hext[1]='0'+dh/100; hext[2]='0'+(dh/10)%10;
1197 hext[3]='0'+dh%10; hext[4]='h'; hext[5]=hext[6]=hext[7]=hext[8]=0;
1198
1199 printf("dh=%d: %dNF...", dh, nnf_d[dh]); fflush(0);
1200 for(j=0;j<=Hod_Min_max;j++) nnf_h[j]=0;
1201
1202 for (v=2;v<VERT_Nmax;v++) if(nnf_vd[v][dh]){
1203 aext[1]='0'+v/10; aext[2]='0'+v%10;
1204 strcpy(fax,aext);
1205 Fchia=fopen(dbaname,"rb");
1206 for(i=0;i<nnf_vd[v][dh];i++){
1207 h12=fgetc(Fchia);
1208 if (!nnf_h[h12]){
1209 hext[5]='0'+h12/100; hext[6]='0'+(h12/10)%10; hext[7]='0'+h12%10;
1210 strcpy(fhx,hext);
1211 Fh[h12]=fopen(dbhname,"wb");
1212 nh++;}
1213 nnf_h[h12]++;
1214 c=fgetc(Fchia); fputc(c,Fh[h12]);
1215 if (v!=c/4) {printf("v=%d, hp.mv=%d", v, (int) (c/4)); exit(0);}
1216 /* for (j=0;j<4;j++) {c=fgetc(Fchia); fputc(c,Fh[h12]);}
1217 for (j=0;j<c%64;j++) fputc(fgetc(Fchia),Fh[h12]);}*/
1218 for (j=0;j<3;j++) fputc(fgetc(Fchia),Fh[h12]);
1219 nuc=fgetc(Fchia); fputc(nuc,Fh[h12]);
1220 c=fgetc(Fchia);
1221 { int ms=c%4; if(ms) c+=3-ms; }
1222 fputc(c,Fh[h12]);
1223 for (j=0;j<nuc%64-1;j++) fputc(fgetc(Fchia),Fh[h12]);}
1224 if(ferror(Fchia)) {
1225 printf("File error in Fchia at dh=%d v=%d\n",dh,v); exit(0);}
1226 fclose(Fchia); }
1227
1228 fprintf(Fhinfo, "%d %d %d\n", dh, nh, nnf_d[dh]);
1229 for (h12=0;h12<=Hod_Min_max;h12++) if (nnf_h[h12]) {
1230 fprintf(Fhinfo, "%d %d ", h12, nnf_h[h12]);
1231 nnf_d[dh]-=nnf_h[h12];
1232 if(ferror(Fh[h12])) {
1233 printf("File error at dh=%d h12=%d\n",dh, h12); exit(0);}
1234 fclose(Fh[h12]);}
1235 fprintf(Fhinfo, "\n");
1236 if (nnf_d[dh])
1237 {printf("nnf_d[%d]!=sum nnf_dh[%d][h12]!", dh, dh); exit(0);}
1238 printf(" sorted\n"); }
1239
1240 if(ferror(Fhinfo)) {printf("File error in Fhinfo\n"); exit(0);}
1241 fclose(Fhinfo);
1242 printf(" done (%ds)\n",(int) difftime(time(NULL),Tstart));
1243 fflush(stdout);
1244 }
1245
Test_Hodge_db(char * dbname)1246 void Test_Hodge_db(char *dbname){
1247
1248 time_t Tstart;
1249 char *filename = (char *) malloc(6+strlen(dbname)+File_Ext_NCmax), *fhx;
1250 int i, j, dh, h12, nh, nnf_sum,
1251 nnf_d[Hod_Dif_max+1], nnf_dh[Hod_Dif_max+1][Hod_Min_max+1];
1252 FILE *Fh;
1253 FILE *Fhinfo;
1254
1255 for(i=0;i<=Hod_Dif_max;i++){
1256 nnf_d[i]=0;
1257 for(j=0;j<=Hod_Min_max;j++) nnf_dh[i][j]=0;}
1258
1259 strcpy(filename,dbname);
1260 strcat(filename,".hinfo");
1261 fhx=&filename[strlen(dbname)+1];
1262
1263 printf("Reading %s\n",filename); fflush(0);
1264
1265 /* read the info-file: */
1266 Fhinfo=fopen(filename,"r");
1267 assert(Fhinfo!=NULL);
1268 while ((fscanf(Fhinfo, "%d",&dh))!=EOF){
1269 fscanf(Fhinfo, "%d %d", &nh, &(nnf_d[dh]));
1270 nnf_sum=0;
1271 for (i=0;i<nh;i++) {
1272 fscanf(Fhinfo, "%d", &h12);
1273 fscanf(Fhinfo, "%d", &(nnf_dh[dh][h12]));
1274 nnf_sum+=nnf_dh[dh][h12];}
1275 if (nnf_sum!=nnf_d[dh]) {
1276 printf("nnf_d[%d]!=sum nnf_dh[%d][h12]!", dh, dh); exit(0);}}
1277 if(ferror(Fhinfo)) {printf("File error in %s\n",filename); exit(0);}
1278 fclose(Fhinfo);
1279
1280 printf("Analysing Hodge-files:\n"); fflush(0);
1281 Tstart=time(NULL);
1282
1283 /* Analyse Hodge&Poly-Data */
1284 for (dh=0;dh<=Hod_Dif_max;dh++) if (nnf_d[dh]){
1285 char hext[9];
1286 hext[0]='d'; hext[1]='0'+dh/100; hext[2]='0'+(dh/10)%10;
1287 hext[3]='0'+dh%10; hext[4]='h'; hext[5]=hext[6]=hext[7]=hext[8]=0;
1288
1289 printf("dh=%d: %dNF...", dh, nnf_d[dh]); fflush(0);
1290 nnf_sum=0;
1291
1292 for (h12=0;h12<=Hod_Min_max;h12++) if(nnf_dh[dh][h12]){
1293 /* unsigned char uc_poly[NUC_Nmax]; */
1294 int c1, /* c2, */ nuc;
1295 hext[5]='0'+h12/100; hext[6]='0'+(h12/10)%10; hext[7]='0'+h12%10;
1296 strcpy(fhx,hext);
1297 Fh=fopen(filename,"rb");
1298 assert (Fh!=0);
1299 while((c1=fgetc(Fh))!=EOF){
1300 nnf_sum++;
1301 /* c2= */ fgetc(Fh);
1302 /* mv=c1/4; mp=(c1%4)*256+c2; */
1303 c1=fgetc(Fh); /* c2= */ fgetc(Fh);
1304 /* nv=c1/4; np=(c1%4)*256+c2; */
1305 c1=fgetc(Fh);
1306 /* mirror=c1/64; */
1307 nuc=c1%64;
1308 for (j=0;j<nuc;j++) /* uc_poly[j]= */ fgetc(Fh);
1309 /* uc_nf_to_P(_P, &MS, &d, &mv, &nuc, uc_poly); */}
1310 if(ferror(Fh)) {
1311 printf("File error in Fh at dh=%d h12=%d\n",dh,h12); exit(0);}
1312 fclose(Fh); }
1313
1314 printf("%d\n",nnf_sum); }
1315
1316 printf(" done (%ds)\n",(int) difftime(time(NULL),Tstart));
1317 fflush(stdout);
1318 }
1319
Extract_from_Hodge_db(char * dbname,char * x_string,PolyPointList * _P)1320 void Extract_from_Hodge_db(char *dbname, char *x_string, PolyPointList *_P){
1321
1322 time_t Tstart;
1323 char c=*x_string, hext[9], com[35],
1324 *filename = (char *) malloc(6+strlen(dbname)+File_Ext_NCmax), *fhx;
1325 unsigned char uc_poly[NUC_Nmax];
1326 int E=998, H1=0, H2=0, M=0, V=0, N=0, F=0, L=1000, i=0, j, dh, h12, nh,
1327 nnf_sum, nnf_d[Hod_Dif_max+1], nnf_dh[Hod_Dif_max+1][Hod_Min_max+1],
1328 HD_from=0, HD_to=Hod_Dif_max, HM_from=0, HM_to=Hod_Min_max, c1, c2, nuc,
1329 mirror, nv, np, mv, mp, d=4, MS, is_poly, is_dual, true_H1, true_H2,
1330 max_mv=34;
1331 FILE *Fh;
1332 FILE *Fhinfo;
1333
1334 hext[0]='d'; hext[4]='h'; hext[8]=0;
1335
1336 while(c){
1337 if(c=='E'){
1338 int neg=0;
1339 c=x_string[++i];
1340 if(c=='-') {neg=1; c=x_string[++i];}
1341 if((c-'0'>=0)&&(c-'0'<=9)) E=0;
1342 while((c-'0'>=0)&&(c-'0'<=9)) {E=10*E+c-'0'; c=x_string[++i];}
1343 if(E%2) {puts("The Euler number E number must be even!"); return;}
1344 if(neg) E*=-1;}
1345 else if(c=='H'){
1346 c=x_string[++i];
1347 while((c-'0'>=0)&&(c-'0'<=9)) {H1=10*H1+c-'0'; c=x_string[++i];}}
1348 else if(c==':'){
1349 c=x_string[++i];
1350 while((c-'0'>=0)&&(c-'0'<=9)) {H2=10*H2+c-'0'; c=x_string[++i];}}
1351 else if(c=='M'){
1352 c=x_string[++i];
1353 while((c-'0'>=0)&&(c-'0'<=9)) {M=10*M+c-'0'; c=x_string[++i];}}
1354 else if(c=='V'){
1355 c=x_string[++i];
1356 while((c-'0'>=0)&&(c-'0'<=9)) {V=10*V+c-'0'; c=x_string[++i];}}
1357 else if(c=='N'){
1358 c=x_string[++i];
1359 while((c-'0'>=0)&&(c-'0'<=9)) {N=10*N+c-'0'; c=x_string[++i];}}
1360 else if(c=='F'){
1361 c=x_string[++i];
1362 while((c-'0'>=0)&&(c-'0'<=9)) {F=10*F+c-'0'; c=x_string[++i];}}
1363 else if(c=='L'){
1364 c=x_string[++i];
1365 L=0;
1366 while((c-'0'>=0)&&(c-'0'<=9)) {L=10*L+c-'0'; c=x_string[++i];}}
1367 else {printf("`%c' is not valid input", c); return;}}
1368
1369 if(H1&&H2){E=2*(H1-H2);}
1370 else if(H1&&(E!=998)) H2=H1-E/2;
1371 else if(H2&&(E!=998)) H1=H2+E/2;
1372
1373 if( H1>491 || H2>491 || H1+H2>502 || (abs(E)!=998 && abs(E)>960) ){
1374 puts("#NF: 0 Note the range for Hodge numbers:");
1375 puts(" h11,h12<=491, h11+h12<=502, |E|<=960.");
1376 return;}
1377
1378 if( (V&&(V<5||V>33)) || (F&&(F<5||F>33))){
1379 puts("#NF: 0 Note the range [5,33] for vertex/facet numbers!");
1380 return;}
1381
1382 if( (M&&(M<6||M>680)) || (N&&(N<6||N>680)) ){
1383 puts("#NF: 0 Note the range [6,680] for point numbers!");
1384 return;}
1385
1386 if(V&&(V<max_mv)) max_mv=V;
1387 if(F&&(F<max_mv)) max_mv=F;
1388 if(M&&(M-1<max_mv)) max_mv=M-1;
1389 if(N&&(N-1<max_mv)) max_mv=N-1;
1390
1391 for(i=0;i<=Hod_Dif_max;i++){
1392 nnf_d[i]=0;
1393 for(j=0;j<=Hod_Min_max;j++) nnf_dh[i][j]=0;}
1394
1395 strcpy(filename,dbname);
1396 strcat(filename,".hinfo");
1397 fhx=&filename[strlen(dbname)+1];
1398
1399 /* printf("Reading %s\n",filename); fflush(0); */
1400
1401 /* read the info-file: */
1402 Fhinfo=fopen(filename,"r");
1403 assert(Fhinfo!=NULL);
1404 while ((fscanf(Fhinfo, "%d",&dh))!=EOF){
1405 fscanf(Fhinfo, "%d %d", &nh, &(nnf_d[dh]));
1406 nnf_sum=0;
1407 for (i=0;i<nh;i++) {
1408 fscanf(Fhinfo, "%d", &h12);
1409 fscanf(Fhinfo, "%d", &(nnf_dh[dh][h12]));
1410 nnf_sum+=nnf_dh[dh][h12];}
1411 if (nnf_sum!=nnf_d[dh]) {
1412 printf("nnf_d[%d]!=sum nnf_dh[%d][h12]!", dh, dh); exit(0);}}
1413 if(ferror(Fhinfo)) {printf("File error in %s\n",filename); exit(0);}
1414 fclose(Fhinfo);
1415
1416 /* printf("Analysing Hodge-files:\n"); fflush(0); */
1417 Tstart=time(NULL);
1418
1419 /* Analyse Hodge&Poly-Data */
1420
1421 nnf_sum=0;
1422 if (H1&&H2) {HM_from=(H1<H2 ? H1 : H2); HM_to=HM_from;}
1423 else if (H1||H2) {HM_to=max(H1,H2); HM_from=max(0,HM_to-Hod_Dif_max);}
1424
1425 for (h12=HM_from; h12<=HM_to; h12++) {
1426
1427 if(abs(E)!=998) {HD_from=abs(E/2); HD_to=abs(E/2);}
1428 else if ((H1||H2)&&(h12<HM_to)) {HD_from=HM_to-h12; HD_to=HM_to-h12;}
1429 else {HD_from=0; HD_to=Hod_Dif_max;}
1430
1431 for (dh=HD_from; dh<=HD_to; dh++) if (nnf_dh[dh][h12]){
1432
1433 hext[1]='0'+dh/100; hext[2]='0'+(dh/10)%10; hext[3]='0'+dh%10;
1434 hext[5]='0'+h12/100; hext[6]='0'+(h12/10)%10; hext[7]='0'+h12%10;
1435 strcpy(fhx,hext);
1436 Fh=fopen(filename,"rb");
1437 assert (Fh!=0);
1438 while((c1=fgetc(Fh))!=EOF){
1439 mv=c1/4;
1440 if(mv>max_mv) break;
1441 c2=fgetc(Fh);
1442 mp=(c1%4)*256+c2;
1443 c1=fgetc(Fh); c2=fgetc(Fh);
1444 nv=c1/4; np=(c1%4)*256+c2;
1445 c1=fgetc(Fh);
1446 mirror=c1/64;
1447 nuc=c1%64;
1448 for (j=0;j<nuc;j++) uc_poly[j]=fgetc(Fh);
1449
1450 if (mirror) {true_H1=h12; true_H2=dh+h12;}
1451 else {true_H2=h12; true_H1=dh+h12;}
1452 is_poly=0;
1453 if ((abs(E)==998)||(E==2*(true_H1-true_H2))) if (!H1||(true_H1==H1))
1454 if (!H2||(true_H2==H2)) if(!M||(M==mp)) if (!V||(V==mv))
1455 if (!N||(N==np)) if (!F||(F==nv)) is_poly=1;
1456 is_dual=0;
1457 if ((abs(E)==998)||(E==2*(true_H2-true_H1))) if (!H1||(true_H2==H1))
1458 if (!H2||(true_H1==H2)) if(!M||(M==np)) if (!V||(V==nv))
1459 if (!N||(N==mp)) if (!F||(F==mv)) is_dual=1;
1460 if (!is_poly&&!is_dual) continue;
1461
1462 uc_nf_to_P(_P, &MS, &d, &mv, &nuc, uc_poly);
1463
1464 if (is_poly){
1465 /* if(!MS) puts("!MS");
1466 if(!mirror) puts("!mirror"); */
1467 if (++nnf_sum>L) {printf("Exceeded limit of %d\n",L); return;}
1468 sprintf(com,"M:%d %d N:%d %d H:%d,%d [%d]",
1469 mp,mv,np,nv,true_H1,true_H2,2*(true_H1-true_H2));
1470 Print_PPL(_P,com); }
1471 if (is_dual&&MS){
1472 VertexNumList VNL; EqList EL;
1473 Long NF[POLY_Dmax][VERT_Nmax];
1474 if (++nnf_sum>L) {printf("Exceeded limit of %d\n",L); return;}
1475 sprintf(com,"M:%d %d N:%d %d H:%d,%d [%d]",
1476 np,nv,mp,mv,true_H2,true_H1,2*(true_H2-true_H1));
1477 Find_Equations(_P,&VNL,&EL);
1478 Small_Make_Dual(_P, &VNL, &EL);
1479 Make_Poly_NF(_P, &VNL, &EL, NF);
1480 Print_Matrix(NF, _P->n, VNL.nv, com); }
1481 }
1482 if(ferror(Fh)) {
1483 printf("File error in Fh at dh=%d h12=%d\n",dh,h12); exit(0);}
1484 fclose(Fh); }
1485 }
1486 printf("#NF: %d\n",nnf_sum);
1487 printf(" done (%ds)\n",(int) difftime(time(NULL),Tstart));
1488 fflush(stdout);
1489 }
1490
Test_Hodge_file(char * filename,PolyPointList * _P)1491 void Test_Hodge_file(char *filename, PolyPointList *_P){
1492 FILE *Ft=fopen(filename,"rb");
1493 unsigned char uc_poly[NUC_Nmax];
1494
1495 int nuc, mirror, nv, np, mv, mp, j, c1, c2, d=4, MS;
1496 assert (Ft!=0);
1497 while((c1=fgetc(Ft))!=EOF){
1498 c2=fgetc(Ft);
1499 mv=c1/4; mp=(c1%4)*256+c2;
1500 c1=fgetc(Ft); c2=fgetc(Ft);
1501 nv=c1/4; np=(c1%4)*256+c2;
1502 c1=fgetc(Ft);
1503 mirror=c1/64; nuc=c1%64;
1504 for (j=0;j<nuc;j++) uc_poly[j]=fgetc(Ft);
1505 uc_nf_to_P(_P, &MS, &d, &mv, &nuc, uc_poly);
1506 Print_PPL(_P,"");
1507 printf("nuc=%d mirror=%d mv=%d mp=%d nv=%d np=%d MS=%d\n",
1508 nuc,mirror,mv,mp,nv,np,MS); }
1509 fclose(Ft);
1510 }
1511 #endif
1512
1513 /* ===================================================================== */
Open_DB(char * dbin,DataBase ** _DB,int info)1514 void Open_DB(char *dbin, DataBase **_DB, int info)
1515 { int i,j,v,nu; DataBase *DB; char *dbname,*fx,ext[4];
1516 if(*dbin==0){*_DB=NULL;return;}
1517 dbname = (char *) malloc(1+strlen(dbin)+File_Ext_NCmax);
1518 DB = (DataBase *) malloc(sizeof(DataBase)); assert(DB != NULL); *_DB=DB;
1519 strcpy(dbname,dbin); strcat(dbname,".info"); fx=&dbname[strlen(dbin)+1];
1520 if(info) {printf("Reading %s: ",dbname); fflush(0);}
1521 DB->Finfo=fopen(dbname,"r"); assert(DB->Finfo!=NULL);
1522 fscanf(DB->Finfo, "%d %d %d %d %d %lld %d %lld %lld %d %d %d %d",
1523 &DB->d,&DB->nV,&DB->nVmax,&DB->NUCmax,&DB->list_num, &DB->nNF,&DB->nSM,
1524 &DB->nNM, &DB->NB, &DB->sl_nNF, &DB->sl_SM, &DB->sl_NM, &DB->sl_NB);
1525 if(info) printf("%lldp (%dsl) %lldnf %lldb\n",
1526 2*(DB->nNF)-DB->nSM-DB->nNM+2*DB->sl_nNF-DB->sl_SM-DB->sl_NM, 2*
1527 DB->sl_nNF-DB->sl_SM-DB->sl_NM, DB->nNF+DB->sl_nNF,DB->NB+DB->sl_NB);
1528 for (v=1;v<VERT_Nmax;v++){ DB->nNUC[v]=0;
1529 for(nu=0; nu<NUC_Nmax; nu++) DB->NFnum[v][nu]=0;}
1530 for(i=0; i<DB->nV; i++){
1531 fscanf(DB->Finfo,"%d",&v); fscanf(DB->Finfo,"%d",&(DB->nNUC[v]));
1532 for(j=0;j<DB->nNUC[v];j++){
1533 fscanf(DB->Finfo,"%d", &nu);
1534 fscanf(DB->Finfo,"%d", &(DB->NFnum[v][nu])); } }
1535 if(ferror(DB->Finfo)) {printf("File error in %s\n",dbname); exit(0);}
1536 fclose(DB->Finfo); ext[0]='v'; ext[3]=0; DB->v=DB->p=DB->nu=0;
1537 for(v=DB->d+1; v<=DB->nVmax; v++) if(DB->nNUC[v])
1538 { ext[1]='0' + v / 10; ext[2]='0' + v % 10; strcpy(fx,ext);
1539 DB->Fv[v]=fopen(dbname,"rb"); assert(DB->Fv[v]!=NULL);
1540 if(0==DB->v) {DB->v=v; while(0==DB->NFnum[v][DB->nu]) (DB->nu)++;}
1541 }
1542 }
Close_DB(DataBase * DB)1543 void Close_DB(DataBase *DB)
1544 { int v; if(DB==NULL) return;
1545 for(v=1+DB->d; v<=DB->nVmax; v++) if(DB->nNUC[v]) {
1546 if(ferror(DB->Fv[v])) {printf("File error at v=%d\n",v); exit(0);}
1547 fclose(DB->Fv[v]);} free(DB);
1548 }
1549 #undef TEST
Read_H_ucNF_from_DB(DataBase * DB,unsigned char * uc)1550 int Read_H_ucNF_from_DB(DataBase *DB, unsigned char *uc)/* p=next read pos */
1551 { static Along totNF; int rest; assert(DB!=NULL);
1552 rest = DB->NFnum[DB->v][DB->nu] - DB->p;
1553 #ifdef TEST
1554 {static int vNF, xvNF, nuNF=DB->NFnum[DB->v][DB->nu];
1555 if(totNF==0){int i,j;for(i=1;i<=DB->NUCmax;i++)vNF+=DB->NFnum[DB->v][i];
1556 for(i=DB->v;i<=DB->nVmax;i++) for(j=1;j<=DB->NUCmax;j++)
1557 xvNF+=DB->NFnum[i][j];
1558 if(xvNF!=DB->nNF) {printf("xvNF=%d nNF=%d",xvNF,DB->nNF);exit(0);}
1559 xvNF=vNF;
1560 }}
1561 #endif
1562 assert(0<=rest); if(rest==0) /* search next v/nu */
1563 { int v=DB->v, nu=DB->nu; /* search nu(v): */
1564 while(DB->nu < DB->NUCmax) if(DB->NFnum[DB->v][++(DB->nu)]) break;
1565 if((nu==DB->nu)||(DB->NFnum[DB->v][DB->nu]==0)) /* search v: */
1566 { if(DB->v < DB->nVmax)
1567 { while(0==DB->nNUC[++DB->v]) assert(DB->v < DB->nVmax); DB->nu=0;
1568 while(DB->nu < DB->NUCmax) if(DB->NFnum[DB->v][++(DB->nu)]) break;
1569 } else { assert(totNF==DB->nNF); return 0; } }
1570 if((v<DB->v)||(nu<DB->nu)) {rest=DB->NFnum[DB->v][DB->nu]; DB->p=0;}
1571 #ifdef TEST
1572 if((v<DB->v)||(nu<DB->nu)){assert(totNF==nuNF);nuNF+=DB->NFnum[DB->v][DB->nu];}
1573 if(v<DB->v) {
1574 if(vNF){printf("vNF[%d]=%d was=%d v_new=%d\n",v,xvNF,vNF,DB->v);exit(0);}
1575 for(i=1;i<=DB->NUCmax;i++) vNF+=DB->NFnum[DB->v][i];}
1576 #endif
1577
1578 } assert(DB->p<=DB->NFnum[DB->v][DB->nu]); assert(rest);
1579 assert(totNF < DB->nNF);
1580 #ifdef TEST
1581 printf("return 1 at totNF=%d v=%d nu=%d p=%d\n",
1582 totNF,DB->v,DB->nu,DB->p); assert(0<vNF--);
1583 #endif
1584 AuxGet_uc(DB->Fv[DB->v],&DB->nu,uc); ++DB->p; totNF++; return 1;
1585 }
1586
Read_SLucNF_from_DB()1587 int Read_SLucNF_from_DB( /* DataBase *DB,unsigned char *uc */ )
1588 { puts("Read_SLucNF_from_DB: to be implemented"); exit(0); return 0;
1589 }
1590
Read_SLpoly_from_DB()1591 int Read_SLpoly_from_DB( /* DataBase *DB,unsigned char *uc */ )
1592 { puts("Read_SLpoly_from_DB: to be implemented"); exit(0); return 0;
1593 }
1594
Read_H_poly_from_DB(DataBase * DB,PolyPointList * P)1595 int Read_H_poly_from_DB(DataBase *DB,PolyPointList *P)
1596 { static unsigned char uc[NUC_Nmax]; static int ms3;
1597 int i,j, MS; Long NF[POLY_Dmax][VERT_Nmax]; VertexNumList V; EqList E;
1598 if(0==ms3) { if(0==Read_H_ucNF_from_DB(DB,uc)) return 0;} else --uc[0];
1599 UCnf2vNF(&DB->d,&DB->v,&DB->nu,uc,NF,&MS); P->n=DB->d; P->np=DB->v;
1600 for(i=0;i<P->np;i++) for(j=0;j<P->n;j++) P->x[i][j]=NF[j][i];
1601 MS%=4;
1602 ms3=(MS==3); if(MS==2)
1603 { assert(Ref_Check(P,&V,&E)); P->np=E.ne;
1604 for(i=0;i<P->np;i++)for(j=0;j<P->n;j++) P->x[i][j]=E.e[i].a[j];
1605 } return 1;
1606 }
1607
Read_H_poly_from_DB_or_inFILE(DataBase * DB,PolyPointList * P)1608 int Read_H_poly_from_DB_or_inFILE(DataBase *DB,PolyPointList *P)
1609 { if((DB==NULL)||(inFILE==NULL)) {CWS cws; return Read_CWS_PP(&cws,P);}
1610 else return Read_H_poly_from_DB(DB,P);
1611 }
1612 /* ===================================================================== */
1613
1614
1615 /* ===================================================================== */
1616 /* ============== from Subpoly.c: ============= */
1617 /* void Make_All_Sublat(NF_List *_L, int n, int v, subl_int diag[POLY_Dmax],
1618 subl_int u[][VERT_Nmax], char *mFlag) */
1619 /* create all inequivalent decompositions diag=s*t into upper
1620 triangular matrices s,t;
1621 t*(first lines of u) becomes the poly P;
1622 choose the elements of t (columns rising, line # falling),
1623 calculate corresponding element of s at the same time;
1624 finally calculate P and add to list */
1625 /* void MakePolyOnSublat(NF_List *_L, subl_int x[VERT_Nmax][VERT_Nmax],
1626 int v, int f, int *max_order, char *mFlag) */
1627 /* Decompose the VPM x as x=w*diag*u, where w and u are SL(Z);
1628 the first lines of u are the coordinates on the coarsest lattices */
1629 /* ============== End of "from Subpoly.c" ============= */
1630
1631 /* =================== Sublattice: phv =================== */
AuxGxP(Long * Gi,Long * V,int * d)1632 Long AuxGxP(Long *Gi,Long *V,int *d)
1633 { Long x=0; int j; for(j=0;j<*d;j++) x+=Gi[j]*V[j]; return x;
1634 }
1635
1636 /* Glz x (lincomb(Points)) -> Diag: if(index>1) print vertices of G*P & D */
1637 #define TEST
1638 #define TEST_OUT
Aux_Print_CoverPoly(int * I,int * d,int * N,Long * X[POLY_Dmax],Long G[][POLY_Dmax],Long * D,int * x,Long Z[][VERT_Nmax],Long * M,int r)1639 void Aux_Print_CoverPoly(int *I,int *d, int *N,Long *X[POLY_Dmax],
1640 Long G[][POLY_Dmax],Long *D,int *x,Long Z[][VERT_Nmax],Long *M,int r)
1641 { int i,j,dia=1, err=0; fprintf(outFILE,"%d %d index=%d D=%ld",*d,*N,
1642 *I,D[0]); for(i=1;i<*d;i++) printf(" %ld",D[i]);
1643 for(i=0;i<r;i++)
1644 { fprintf(outFILE," /Z%ld:",M[i]);
1645 for(j=0;j<*N;j++) fprintf(outFILE," %ld",Z[i][j]);
1646 } printf(" #%d\n",*x);
1647 for(i=0;i<*d;i++){
1648 for(j=0;j<*N;j++){
1649 Long Xij=AuxGxP(G[i],X[j],d);
1650 if(0!=(Xij % D[i])) err=1;
1651 printf("%ld ",Xij);}
1652 puts("");
1653 }
1654 #ifdef TEST
1655 if((dia==0)||err){
1656 int i,j;printf("D=");for(i=0;i<*d;i++)printf("%ld ",D[i]);
1657 printf(" index[%d]=%d\n",*x,*I);
1658 for(i=0;i<*d;i++) {
1659 printf("G=");
1660 for(j=0;j<*d;j++) printf("%2ld ",G[i][j]);
1661 printf(" G.P=");
1662 for(j=0;j<*N;j++) printf("%2ld ",AuxGxP(G[i],X[j],d));
1663 puts("");}
1664 exit(0);}
1665 #endif
1666 }
Aux_Print_SLpoly(int * I,int * d,int * N,Long * X[POLY_Dmax],Long G[][POLY_Dmax],Long * D,int * x)1667 void Aux_Print_SLpoly(int *I,int *d, int *N,Long *X[POLY_Dmax],
1668 Long G[][POLY_Dmax],Long *D,int *x)
1669 { int i,j,dia=1, err=0; printf("%d %d index=%d D=%ld",*d,*N,*I,D[0]);
1670 for(i=1;i<*d;i++) printf(" %ld",D[i]);
1671 printf(" #%d\n",*x);
1672 for(i=0;i<*d;i++)
1673 { /* int z=1;*/ for(j=0;j<*N;j++)
1674 { Long Xij=AuxGxP(G[i],X[j],d);
1675 #ifdef TEST_STUFF
1676 if(Xij){if(z&&dia)dia=(labs(Xij)==D[i]);z=0;} /* SUFFICIENT only !!!*/
1677 if(dia==0)if(i==(*d-1)){Long g=Xij,T; int J; for(J=j+1;J<*N;J++){T=
1678 labs(AuxGxP(G[i],X[J],d)); if(T) g=Fgcd(g,T);} if(g==D[i]) dia=1;}
1679 #endif
1680 if(0!=(Xij % D[i])) err=1;
1681 printf("%ld ",Xij/D[i]);}
1682 puts("");
1683 }
1684 #ifdef TEST
1685 if((dia==0)||err){
1686 int i,j;printf("D=");for(i=0;i<*d;i++)printf("%ld ",D[i]);
1687 printf(" index[%d]=%d\n",*x,*I);for(i=0;i<*d;i++)
1688 {printf("G=");for(j=0;j<*d;j++)printf("%2ld ",G[i][j]);
1689 printf(" G.P=");
1690 for(j=0;j<*N;j++) printf("%2ld ",AuxGxP(G[i],X[j],d));
1691 puts("");}
1692 exit(0);}
1693 #endif
1694 }
Aux_Make_Dual(PolyPointList * P,VertexNumList * V,EqList * E)1695 void Aux_Make_Dual(PolyPointList *P, VertexNumList *V, EqList *E)
1696 { Long VM[VERT_Nmax][POLY_Dmax]; int i,j, d=P->n, e=E->ne, v=V->nv;
1697 assert(e<=VERT_Nmax); P->np=V->nv=e; E->ne=v;
1698 for(i=0;i<v;i++)for(j=0;j<d;j++) VM[i][j]=P->x[V->v[i]][j];
1699 for(i=0;i<e;i++){for(j=0;j<d;j++)P->x[i][j]=E->e[i].a[j]; V->v[i]=i;}
1700 for(i=0;i<v;i++){for(j=0;j<d;j++)E->e[i].a[j]=VM[i][j]; E->e[i].c=1;}
1701 assert(Ref_Check(P,V,E));
1702 }
1703 void PrintVPHMusage();
1704 int Make_Lattice_Basis(int d, int p, Long *P[POLY_Dmax], /* index=det(D) */
1705 Long G[][POLY_Dmax], Long *D);/* G x P generates diagonal lattice D */
PH_Sublat_Polys(char * dbin,int omitFIP,PolyPointList * _P,char sF)1706 void PH_Sublat_Polys(char *dbin, int omitFIP, PolyPointList *_P, char sF)
1707 { static EqList E; VertexNumList V; int x=0, K, B, I=1; /* index>I only */
1708 Long *RelPts[POINT_Nmax];DataBase *DB=NULL; if(*dbin)Open_DB(dbin,&DB,0);
1709 K=((sF!='P')&&(sF!='H')&&(sF!='Q')&&(sF!='B')); /* K=='CoverPoly' */
1710 if(('1'<sF)&&(sF<='9')) I=sF-'0';
1711 B=(omitFIP==2); /* 'q' for index >I */
1712 while(Read_H_poly_from_DB_or_inFILE(DB,_P))
1713 { Long D[POLY_Dmax],G[POLY_Dmax][POLY_Dmax],PM[VERT_Nmax][VERT_Nmax];
1714 int index, N=0; /* if(!Ref_Check(_P,&V,&E)) Print_PPL(_P,""); */
1715 if(B)assert(_P->n==4);/* b:Brower group only for CY hypersurface d=4 */
1716 assert(Ref_Check(_P,&V,&E));
1717 /* Aux_Make_Dual(_P,&V,&E); */ /* don't dualize: take M-lattice poly */
1718 Make_VEPM(_P,&V,&E,PM); _P->np=V.nv; Complete_Poly(PM,&E,V.nv,_P); ++x;
1719 if(omitFIP==0) for(N=0;N<_P->np;N++) RelPts[N]=_P->x[N];
1720 else if(omitFIP<3 ) /* Omit facet-IPs */
1721 { int p; for(p=0;p<_P->np;p++)
1722 { int e, z=0;
1723 for(e=0;e<E.ne;e++)
1724 if(0==Eval_Eq_on_V(&E.e[e],_P->x[p],_P->n)) z++;
1725 if(z>omitFIP) RelPts[N++]=_P->x[p];
1726 } assert(V.nv<=N); /* count <E,.>=0; if(n>1) add_to_RelPts; */
1727 }
1728 else if(omitFIP==3) /* Omit all non-vertices */
1729 { for(N=0;N<V.nv;N++) RelPts[N]=_P->x[N];
1730 /* { int tnv=V.nv; Find_Equations(_P,&V,&E);assert(V.nv==tnv);
1731 for(tnv=0;tnv<V.nv;tnv++)assert(V.v[tnv]<V.nv);} */
1732 }
1733 else {puts("something wrong in PH_Sublat_Polys"); exit(0);}
1734 if(K)
1735 { Long Z[POLY_Dmax][VERT_Nmax], M[POLY_Dmax]; int r;
1736 index=Sublattice_Basis(_P->n,N,RelPts,Z,M,&r,G,D);
1737 if(B)if(D[2]==1)continue;
1738 if(omitFIP==3)if(D[1]==1)continue;
1739 assert(index>0); if(index<=I)continue; assert(Ref_Check(_P,&V,&E));
1740 Aux_Print_CoverPoly(&index,&_P->n,&N,RelPts,G,D,&x,Z,M,r);
1741 }
1742 else
1743 { index=Make_Lattice_Basis(_P->n,N,RelPts,G,D); if(1==index)continue;
1744 if(B)if(D[2]==1)continue;
1745 if(omitFIP==3)if(D[1]==1)continue;
1746 assert(index>0); assert(Ref_Check(_P,&V,&E)); Print_VL(_P,&V,"");
1747 Aux_Print_SLpoly(&index,&_P->n,&N,RelPts,G,D,&x);
1748 }
1749 } if(*dbin) Close_DB(DB);
1750 }
1751 #undef TEST
1752 /* Lattice generated by vertices; UT-decomp of diag */
V_Sublat_Polys(char mr,char * dbin,char * polyi,char * polyo,PolyPointList * _P)1753 void V_Sublat_Polys(char mr,char *dbin,char *polyi,char *polyo,
1754 PolyPointList *_P)
1755 { NF_List *_L=(NF_List *) malloc(sizeof(NF_List));
1756 int max_order=1;
1757 static EqList E; VertexNumList V; int x=0;
1758 Long *RelPts[VERT_Nmax];DataBase *DB=NULL; if(*dbin)Open_DB(dbin,&DB,0);
1759 assert(_L!=NULL);
1760 if(!(*polyo)) {
1761 puts("You have to specify an output file via -po in -sv-mode.");
1762 printf("For more help use option `-h'\n");
1763 exit(0);}
1764 _L->of=0; _L->rf=0; _L->iname=polyi; _L->oname=polyo; _L->dbname=dbin;
1765 Init_NF_List(_L); _L->SL=0;
1766 while(Read_H_poly_from_DB_or_inFILE(DB,_P))
1767 { Long D[POLY_Dmax],G[POLY_Dmax][POLY_Dmax];
1768 int index, N; assert(Ref_Check(_P,&V,&E));
1769 for(N=0;N<V.nv;N++) RelPts[N]=_P->x[V.v[N]];
1770 ++x;
1771 index=Make_Lattice_Basis(_P->n,N,RelPts,G,D); if(1==index) continue;
1772 assert(index>0); if(index>max_order) max_order=index;
1773 #ifdef TEST
1774 Aux_Print_SLpoly(&index,&_P->n, &N,RelPts,G,D,&x);
1775 #endif
1776 { int i,j; subl_int diag[POLY_Dmax], U[POLY_Dmax][VERT_Nmax];
1777 for(i=0;i<_P->n;i++){diag[i]=D[i]; for(j=0;j<V.nv;j++) {int k;
1778 U[i][j]=0; for(k=0;k<_P->n;k++) U[i][j]+=G[i][k]*RelPts[j][k];
1779 assert(0==(U[i][j]%D[i])); U[i][j]/=D[i];}}
1780 Make_All_Sublat(_L, _P->n, V.nv, diag, U, &mr, _P);
1781 }
1782 } if(*dbin) Close_DB(DB);
1783 printf("max_order=%d\n", max_order); Write_List_2_File(polyo,_L);
1784 _L->TIME=time(NULL); fputs(ctime(&_L->TIME),stdout);
1785 free(_L);
1786 }
VPHM_Sublat_Polys(char sFlag,char mr,char * dbin,char * polyi,char * polyo,PolyPointList * _P)1787 void VPHM_Sublat_Polys(char sFlag,char mr,char *dbin,char *polyi,char *polyo,
1788 PolyPointList *_P)
1789 { switch(sFlag){ /* if(dbin=0) read from inFILE; */
1790 case 'p': case 'P': PH_Sublat_Polys(dbin,0,_P,sFlag); break ;
1791 case 'h': case 'H': PH_Sublat_Polys(dbin,1,_P,sFlag); break ;
1792 case 'b': case 'B': PH_Sublat_Polys(dbin,2,_P,sFlag); break ;
1793 case 'q': case 'Q': PH_Sublat_Polys(dbin,3,_P,sFlag); break ;
1794 case 'v': case 'V': V_Sublat_Polys(mr,dbin,polyi,polyo,_P); break ;
1795 case 'm': case 'M': Find_Sublat_Polys(mr,dbin,polyi,polyo,_P); break ;
1796 default:if(('1'<sFlag)&&(sFlag<='9'))PH_Sublat_Polys(dbin,3,_P,sFlag);
1797 else{puts("-s# requires that # is in {v,p,h,b,m,q}");PrintVPHMusage();
1798 } }
1799 }
PrintVPHMusage()1800 void PrintVPHMusage()
1801 {puts(" -sh ... gen by codim>1 points (omit IPs of facets)");
1802 puts(" -sp ... gen by all points");
1803 puts(" -sb ... generated by dim<=1 (edges), print if rank=2 ");
1804 puts(" -sq ... generated by vertices, print if rank=3 ");
1805 puts(" q,b currently assume that dim=4");
1806 exit(0);
1807 }
1808
1809
1810
Bin_2_ANF(char * polyi,int max,PolyPointList * _P)1811 void Bin_2_ANF(char *polyi, int max, PolyPointList *_P)
1812 { FILE *F=fopen(polyi,"rb"); FInfoList L; UPint list_num,tNF=0; Along tNB=0;
1813 int d, v, s, sl_nNF, sl_SM, sl_NM, sl_NB, mc=0, MS, nu; unsigned i, j;
1814 unsigned char uc[POLY_Dmax*VERT_Nmax]; VertexNumList V; EqList E;
1815 Long NF[POLY_Dmax][VERT_Nmax]; Init_FInfoList(&L);
1816
1817 if(F==NULL) {printf("Input file %s not found\n",polyi); exit(0);}
1818 d=fgetc(F); assert(d==0); /* for(i=0;i<d;i++) fgetc(F); */
1819 d=fgetc(F); L.nV=fgetc(F); L.nVmax=fgetc(F); L.NUCmax=fgetc(F);
1820 list_num=fgetUI(F);
1821 L.nNF=fgetUI(F); L.nSM=fgetUI(F); L.nNM=fgetUI(F); L.NB=fgetUI(F);
1822 sl_nNF=fgetUI(F); sl_SM=fgetUI(F); sl_NM=fgetUI(F); sl_NB=fgetUI(F);
1823
1824 for(i=0;i<L.nV;i++)
1825 { v=fgetc(F); L.nNUC[v]=fgetc(F); /* read #nuc's per #Vert */
1826 for(j=0;j<L.nNUC[v];j++)
1827 { L.NFnum[v][nu=fgetc(F)]=fgetUI(F); /* read nuc and #NF(v,nu)*/
1828 tNF+=L.NFnum[v][nu]; tNB+=L.NFnum[v][nu]*(Along)nu;
1829 }
1830 } assert( 0 == (unsigned int)(tNB-L.NB) ); L.NB=tNB; assert(tNF==L.nNF);
1831
1832 for(v=d+1;v<=L.nVmax;v++) if(L.nNUC[v]) /* write honest polys */
1833 { int I,J; for(nu=1;nu<=L.NUCmax;nu++) for(j=0;j<L.NFnum[v][nu];j++)
1834 { for(s=0;s<nu;s++) uc[s]=fgetc(F);
1835 UCnf_2_ANF(&d,&v,&nu,uc,NF,&MS); MS %= 4; _P->n=d; _P->np=v;
1836 for(I=0;I<v;I++) for(J=0;J<d;J++) _P->x[I][J]=NF[J][I];
1837 /* assert(Ref_Check(_P,&V,&E)); */ assert(MS==1);
1838
1839 if(MS!=2) /* if(MS!=2) print NF */
1840 if(!max || Poly_Max_check(_P,&V,&E))
1841 { mc++; Print_NF(outFILE,&d,&v,NF);}
1842 if(MS>1) /* if(MS>1); print Mirror */
1843 if(!max || Poly_Min_check(_P,&V,&E))
1844 { mc++; Small_Make_Dual(_P,&V,&E);
1845 Make_Poly_NF(_P,&V,&E,NF); Print_NF(outFILE,&d,&(V.nv),NF);
1846 }
1847 }
1848 }
1849 printf("np=%lld+%dsl ",2*L.nNF-L.nSM-L.nNM,2*sl_nNF-sl_SM-sl_NM);
1850 printf( /* write Finfo */
1851 "%dd %dv<=%d n<=%d %dnv %lld %d %lld %lld %d %d %d %d",
1852 d, L.nV, L.nVmax, L.NUCmax,list_num,
1853 L.nNF,L.nSM,L.nNM,L.NB, sl_nNF, sl_SM, sl_NM, sl_NB);
1854 if(max) printf(" r-max=%d",mc);
1855 puts("");
1856 }
1857
Bin_2_ANF_DBsl(char * dbi,int max,int vf,int vt,PolyPointList * _P)1858 void Bin_2_ANF_DBsl(char *dbi, int max, int vf, int vt, PolyPointList *_P)
1859 { FILE *F; FInfoList L;
1860 int d, v, nu, i, j, list_num, mc=0, MS,
1861 sl_nNF, sl_SM, sl_NM, sl_NB, tSM=0, tNM=0;
1862 char *Ifx, *Ifn = (char *) malloc(1+strlen(dbi)+File_Ext_NCmax);
1863 Long NF[POLY_Dmax][VERT_Nmax];
1864 VertexNumList V;EqList E;
1865 strcpy(Ifn,dbi); Ifx=&Ifn[strlen(dbi)]; strcpy(Ifx,".info");
1866 F=fopen(Ifn,"r"); if(F==NULL) {puts("Info File not found");exit(0);}
1867 Init_FInfoList(&L); /* start reading the file */
1868 fscanf(F,"%d%d%d%d%d%lld%d%lld %lld %d%d%d%d",
1869 &d,&i,&j,&nu,&list_num,&L.nNF,&L.nSM,&L.nNM,&L.NB,
1870 &sl_nNF,&sl_SM,&sl_NM,&sl_NB); L.nV=i; L.nVmax=j; L.NUCmax=nu;
1871 if(sl_NB && (vf==2)&&(vt==VERT_Nmax-1))
1872 { strcpy(Ifx,".sl"); fclose(F);
1873 if(NULL==(F=fopen(Ifn,"rb"))){printf("Open %s failed",Ifn);exit(0);}
1874 }
1875 else /* puts("no .sl file"); */
1876 { DataBase *DB; Open_DB(dbi, &DB, 0); if((DB->v<vf)||(DB->nVmax>vt))
1877 DB_fromVF_toVT(DB,vf,vt); /* read only vf <= v <= vt :: */
1878 for(i=0; Read_H_poly_from_DB(DB,_P); i++)
1879 { int c, l, p=_P->np-1, off=_P->x[p][0];
1880 if(off) for(l=0;l<d;l++) for(c=l;c<=p;c++) _P->x[c][l]-=off;
1881 /* if(off) Print_PPL(_P,"off"); else */
1882 Print_PPL(_P,"");
1883 } printf("#poly=%d\n",i); Close_DB(DB);
1884 }
1885 for(i=0;i<sl_nNF;i++)
1886 { int I,J; unsigned char uc[NUC_Nmax]; v=fgetc(F); assert(v<=VERT_Nmax);
1887 nu=fgetc(F); AuxGet_uc(F,&nu,uc); assert(!ferror(F));
1888 if((*uc%4)==0)tSM++; else if((*uc%4)<3)tNM++;
1889
1890 UCnf_2_ANF(&d,&v,&nu,uc,NF,&MS); MS %= 4; _P->n=d; _P->np=v;
1891 for(I=0;I<v;I++) for(J=0;J<d;J++) _P->x[I][J]=NF[J][I];
1892 if(MS!=2) /* if(MS!=2) print NF */
1893 if(!max || Poly_Max_check(_P,&V,&E))
1894 { mc++; Print_NF(outFILE,&d,&v,NF);}
1895 if(MS>1) /* if(MS>1); print Mirror */
1896 if(!max || Poly_Min_check(_P,&V,&E))
1897 { mc++; Small_Make_Dual(_P,&V,&E);
1898 Make_Poly_NF(_P,&V,&E,NF); Print_NF(outFILE,&d,&(V.nv),NF);
1899 }
1900 }
1901 printf("np=%lld+%dsl ",2*L.nNF-L.nSM-L.nNM,2*sl_nNF-sl_SM-sl_NM);
1902 printf( /* write Finfo */
1903 "%dd %dv<=%d n<=%d %dnv ... %d %d %d %d", /* %d %d %d %lld */
1904 d, L.nV, L.nVmax, L.NUCmax,list_num,
1905 /* L.nNF,L.nSM,L.nNM,L.NB, */ sl_nNF, sl_SM, sl_NM, sl_NB);
1906 if(max) printf(" r-max=%d",mc);
1907 puts("");
1908 }
1909
Gen_Bin_2_ascii(char * pi,char * dbi,int max,int vf,int vt,PolyPointList * P)1910 void Gen_Bin_2_ascii(char *pi,char *dbi,int max,int vf,int vt,PolyPointList *P)
1911 { if(*pi) Bin_2_ANF(pi, max, P);
1912 else if(*dbi) Bin_2_ANF_DBsl(dbi, max, vf, vt, P);
1913 else puts("With -B[2A] you have to specify input via -pi or -di");
1914 }
1915
1916