1 #include "Global.h"
2 #include "Subpoly.h"
3
4
5 /* VF_2_ucNF / UCnf2vNF compression/decompression assumes that vNF[0][0]==1
6 * (see UCnf2vNF: "{ off=NF[0][0]-1; ..."
7 */
8
9 /* Write_List_2_File */
10
11 /* BIN File Format: ==[FileINFO|data] MAX_REC_DEPTH
12 * [FileINFO]:: [RecDepth<128] or [128+rd/256,rd%256] but: .<7 !!!
13 * L->b[] rd byte <256 !!
14 * d #Vposs #Vmax NUCmax 4 byte <256 !!
15 * #li #NF #SM #NM #NB 5*4byte <2^32!! but: NB%2^32
16 * sublattice #N? #SM #NM #NB 4*4byte
17 * #v_ #nu_ #v2 nu2 ... 2 byte <256 !!
18 * -> nu nNF 1+4byte ... li=sum_{#nu(v)}
19 */
20
21
22 #define TEST_UCnf /* ==== check that decomp(compress)=id ==== */
23
24 #if ( 2 * WATCHREF > SAVE_INC )
25 #error increase SAVE_INC / WATCHREF
26 #endif
27
28 #ifdef ADD_LIST_LENGTH /* ... square root of SAVE_INC */
29 #if ( 2 * ADD_LIST_LENGTH > SAVE_INC )
30 #error increase SAVE_INC / ADD_LIST_LENGTH
31 #endif
32 #else
IntSqrt(int q)33 int IntSqrt(int q) /* sqrt(q) => r=1; r'=(q+r*r)/(2r); */
34 { assert(q>0); if(q<4) return 1; else { /* troubles: e.g. 9408 */
35 long long r=(q+1)/2,n; while(r > (n=(q+r*r)/(2*r))) r=n; if(q<r*r)r--;
36 if((r*r<=q)&&(q<(r+1)*(r+1))) return (int)r;
37 else {printf("Error in sqrt(%d)=%d\n",q,(int)n);exit(0);}} return 0;
38 }
39 #define ADD_LIST_LENGTH (IntSqrt(SAVE_INC))
40 #endif
41
42 #undef More_File_IO_Data
43
44 /* === POLY_NF on FILE: np nv nf complete lines of upper_trian_matrix === */
45 /* === POLY_NF on LIST: np nv nf lines with v00 & v[ij] (i>j) omitted === */
46
47 #define INCREMENTAL_TIME /* these flags should stay "on" */
48 #define INCREMENTAL_WRITE /* write no polys from pi-file */
49 #define ACCEL_PEntComp /* speed optimization flags */
50 #define USE_UNIT_ENCODE
51
52 /* ========== Auxiliary routines from Polynf.c ========== */
53 void Eval_Poly_NF(int *d,int *v,int *f, Long VM[POLY_Dmax][VERT_Nmax],
54 Long VPM[VERT_Nmax][VERT_Nmax], /* in */
55 Long pNF[POLY_Dmax][VERT_Nmax], int t); /* out */
56 int Init_rVM_VPM(PolyPointList *P, VertexNumList *_V,EqList *_F, /* in */
57 int *d,int *v,int *f, Long VM[POLY_Dmax][VERT_Nmax], /* out */
58 Long VPM[VERT_Nmax][VERT_Nmax]); /* return reflexive */
59
60
61 typedef struct {int base[VERT_Nmax+1][NB_MAX], nuc[VERT_Nmax+1][NB_MAX],
62 v[VERT_Nmax+1];} Base_List;
63
64 void VF_2_ucNF(PolyPointList *P, VertexNumList *V, EqList *E, /* IN */
65 int *NV, int *nUC, unsigned char *UC); /* OUT */
66 void Print_Statistics(NF_List *);
67 void Init_BaseList(Base_List **BL,int *d); /* malloc + init.; BL=&(list) */
68 void Insert_PPent_into_Pent(NF_List *S);
69 void Read_In_File(NF_List *);
70 void Read_Aux_File(NF_List *);
71
72 /* ========= End of Headers and TypeDefs ========== */
73 /* ===================================================================== */
74
Init_FInfoList(FInfoList * FI)75 void Init_FInfoList(FInfoList *FI)
76 { int i, j;
77 FI->nNF= FI->nSM= FI->nNM= FI->NB=0; FI->nVmax= FI->NUCmax= FI->nV=0;
78 for(i=0;i<=VERT_Nmax;i++) FI->nNUC[i]=0;
79 FI->NFli=NULL;
80 for(i=0;i<=VERT_Nmax;i++) for(j=0;j<NUC_Nmax;j++) FI->NFnum[i][j]=0;
81 }
82
Init_New_List(NF_List * S)83 void Init_New_List(NF_List *S)
84 { S->ANB = (SAVE_INC+SL_Nmax) * CperR_MAX; S->NewNB = S->NP = S->nSLP = 0;
85 S->RemNB=S->PEN=S->PPEN=S->SLN=0; S->peNM=S->peSM=S->slNM=S->slSM=0;
86 S->PE = (PEnt *) malloc( (SAVE_INC+SL_Nmax) * sizeof(PEnt) );
87 S->PPE = (PPEnt *) malloc( ADD_LIST_LENGTH * sizeof(PPEnt) );
88 S->SLp = (int *) malloc( SL_Nmax * sizeof(int) );
89 S->NewNF = (unsigned char *) malloc( S->ANB * sizeof(char) ); S->NC = 0;
90 assert((S->PE!=NULL)&&(S->PPE!=NULL)&&(S->SLp!=NULL)&&(S->NewNF!=NULL));
91 #ifdef __DECC
92 /* printf("NULL=%p S->PE=%p S->NewNF=%p\n",NULL,S->PE,S->NewNF); */
93 #endif
94 }
95
96 NF_List *AuxNFLptr=NULL; /* dirty trick for Xmin Xmax Xdif */
Init_NF_List(NF_List * L)97 void Init_NF_List(NF_List *L)
98 {
99 L->TIME= L->SAVE= time(NULL); fputs(ctime(&L->TIME),stdout); L->CLOCK= clock();
100 L->IP_Time = L->NF_Time = 0; L->d = L->nNF = L->nIP = L->nSLNF = 0;
101
102 Init_FInfoList(&L->In); if(*L->iname) Read_In_File(L);
103 if(*L->dbname) Init_DB(L);
104 Init_FInfoList(&L->Aux); Init_New_List(L); L->savedNP= 0; AuxNFLptr=L;
105 if(L->rf) Read_Aux_File(L); else L->rd=0;
106 }
107
PrintNumbers(NF_List * S)108 void PrintNumbers(NF_List *S)
109 { printf(
110 "NP=%lld nSLP=%d H=2*%lld - %lldnm - %dsm SL=2*%d - %dnm - %dsm\n",
111 S->NP, S->nSLP, S->Aux.nNF+S->PEN+S->PPEN,
112 S->Aux.nNM+S->peNM, S->peSM+S->Aux.nSM, S->SLN,
113 S->slNM, S->slSM); fflush(stdout);
114 }
Test_SLnbMS(NF_List * S)115 void Test_SLnbMS(NF_List *S)
116 { int i, tNB=0, tNM=0, tSM=0, testNB=!(S->PEN||S->PPEN);
117 if(testNB)assert(S->RemNB==0);
118 for(i=0;i<S->SLN;i++)
119 { unsigned char *C=&S->NewNF[S->SLp[i]]; int ms=C[2]%4;
120 if(C[0]>VERT_Nmax) {printf("v[%d]=%d\n",i,C[0]);exit(0);}
121 tNB+=C[1]+2; if(ms) {if(ms<3)tNM++;} else tSM++;
122 } if((S->slNM!=tNM) || (S->slSM!=tSM) || (testNB&&(S->NewNB!=tNB)))
123 { printf("Test_SLnbMS NM: %d=%d SM: %d=%d NM: %lld=%d",S->slNM,tNM,
124 S->slSM,tSM,S->NewNB,tNB); exit(0);
125 } assert(S->nSLP==2*S->SLN-tSM-tNM);
126 }
127 void Test_ucNF(int *d, int *tnv, int *tnuc, unsigned char *tuc,
128 PolyPointList *_P);
Test_NF_List(NF_List * S,PolyPointList * _P)129 void Test_NF_List(NF_List *S, PolyPointList *_P)
130 { int i; if(S->PPEN) Insert_PPent_into_Pent(S);
131 printf("test PEN=%d ",S->PEN);
132 for(i=0;i<S->PEN;i++)
133 { unsigned char *C=&S->NewNF[S->PE[i].c], *uc=&C[2];
134 int nv=C[0], nuc=C[1]; Test_ucNF(&S->d,&nv,&nuc,uc,_P);
135 if(i)
136 { unsigned char *oldC=&S->NewNF[S->PE[i-1].c], *olduc=&C[2];
137 unsigned int oldn=S->PE[i-1].n;
138 if(oldC[0]>nv) { printf("NV failed at i=%d\n",i);exit(0); }
139 if(oldC[0]==nv)
140 { assert(oldC[1]<=nuc);
141 if(oldC[1]==nuc) if(RIGHTminusLEFT(olduc,uc,&nuc)<0)
142 { printf("failed at i=%d\n",i);exit(0); }
143 if(oldC[1]==nuc) if(oldn>S->PE[i].n)
144 { printf("oldn>S->PE[%d].n\n",i);exit(0); }
145 }
146 }
147 } /* printf("test SLN=%d ",S->SLN);
148 for(i=0;i<S->SLN;i++)
149 { unsigned char *C=&S->NewNF[S->SLp[i]], *uc=&C[2];
150 int nv=C[0], nuc=C[1]; Test_ucNF(&S->d,&nv,&nuc,uc);
151 } puts("Test o.k."); */
152 }
Test_PPEN(NF_List * S)153 void Test_PPEN(NF_List *S)
154 { int i; /* printf("test PPEN=%d ",S->PPEN); */
155 for(i=0;i<S->PPEN;i++)
156 { unsigned char *C=&S->NewNF[S->PPE[i].pe.c], *uc=&C[2];
157 int nv=C[0], nuc=C[1]; /** / Test_ucNF(&S->d,&nv,&nuc,uc); / **/
158 if(i)
159 { unsigned char *oldC=&S->NewNF[S->PPE[i-1].pe.c], *olduc=&oldC[2];
160 if(oldC[0]>nv) { printf("NV failed at i=%d: v%d=%d v%d=%d\n",
161 i,i,*C,i-1,*oldC);exit(0); }
162 if(oldC[0]==nv)
163 { assert(oldC[1]<=nuc);
164 if(oldC[1]==nuc) if(RIGHTminusLEFT(olduc,uc,&nuc)<0)
165 { printf("failed at i=%d\n",i);exit(0);
166 }
167 }
168 }
169 }
170 }
InfoSize(int rd,int lists,FInfoList * FI)171 int InfoSize(int rd, int lists, FInfoList *FI)
172 { return (rd/128+rd+1 + 4+ 2*FI->nV+ lists) + sizeof(int) * (9 + lists);
173 }
fgetUI(FILE * F)174 unsigned int fgetUI(FILE *F) /* read unsigned int from bin file */
175 { unsigned char A,B,C,D; /* L=D+256*(C+256*(B+256*A)); */
176 fscanf(F,"%c%c%c%c",&A,&B,&C,&D);
177 return (((unsigned int)A * 256 + (unsigned int)B) * 256 +
178 (unsigned int)C) * 256 + (unsigned int)D;
179 }
Read_Bin_Info(FILE * F,int * d,unsigned * li,int * SLN,int * slSM,int * slNM,Along * NewSLnb,FInfoList * FI)180 void Read_Bin_Info(FILE *F, int *d, unsigned *li, int *SLN, int *slSM,
181 int *slNM, Along *NewSLnb, FInfoList *FI)
182 { int i,j,v; unsigned tli=0; Along tNF=0,tNB=0; *d=fgetc(F);
183 assert(*d <= POLY_Dmax);
184 FI->nV=fgetc(F); FI->nVmax=fgetc(F); FI->NUCmax=fgetc(F);
185 *li=fgetUI(F); FI->nNF=fgetUI(F); FI->nSM=fgetUI(F);
186 FI->nNM=fgetUI(F); FI->NB=fgetUI(F);
187 *SLN=fgetUI(F); *slSM=fgetUI(F); *slNM=fgetUI(F); *NewSLnb=fgetUI(F);
188 for(i=0;i<FI->nV;i++)
189 { v=fgetc(F); tli+=(FI->nNUC[v]=fgetc(F));
190 for(j=0;j<FI->nNUC[v];j++)
191 { unsigned nu, nnb=FI->NFnum[v][nu=fgetc(F)]=fgetUI(F);
192 tNF+=nnb; tNB+=nnb*nu;
193 }
194 } assert(tNF==FI->nNF); assert(tli==*li);
195 assert( 0 == (unsigned int) (tNB-FI->NB) ); FI->NB=tNB;
196 }
Read_Honest_Poly(FILE * F,FInfoList * FI,NF_List * L)197 void Read_Honest_Poly(FILE *F,FInfoList *FI,NF_List *L)
198 { int i, v, nu; unsigned li; Along tNF=0, pos=0;
199 Init_FInfoList(FI); L->rd=fgetc(F); if(128<=(L->rd))
200 { assert((L->rd-=128)<7); L->rd = 128*L->rd + fgetc(F); /* DirtyFix rd */
201 }
202 for(i=0;i<L->rd;i++) L->b[i]=fgetc(F);
203 Read_Bin_Info(F,&L->d,&li,&L->SLN,&L->slSM,&L->slNM,&L->NewNB,FI);
204
205 L->nSLP=2*L->SLN-L->slSM-L->slNM; L->PEN=L->PPEN=L->peNM=L->peSM=0;
206 L->NP=L->savedNP=L->nSLP+2*FI->nNF-FI->nSM-FI->nNM; L->RemNB=0;
207 printf("%lld+%dsl %lldm+%ds %lldb",L->NP-L->nSLP,L->nSLP,
208 FI->nNF-FI->nSM-FI->nNM,FI->nSM, FI->NB+L->NewNB);
209 if(L->rd) printf(" rd=%d",L->rd);
210 if(FI->NFli != NULL) printf("WARNing: NFli != NULL");
211 fflush(stdout);
212 FI->NFli = (unsigned char *) malloc( FI->NB * sizeof(char) );
213 if(FI->NFli==NULL) {puts("Aux.NFli allocation failed");exit(0);}
214 for(v=L->d+1;v<=FI->nVmax;v++) if(FI->nNUC[v])
215 for(nu=1;nu<=FI->NUCmax;nu++) if(FI->NFnum[v][nu])
216 { unsigned j; FI->NF[v][nu]=&FI->NFli[pos]; tNF+=FI->NFnum[v][nu];
217 for(j=0;j<FI->NFnum[v][nu];j++) for(i=0;i<nu;i++)
218 FI->NFli[pos++]=fgetc(F);
219 } assert(pos==FI->NB);
220 /* printf("\nRead: tNF=%d pos=%d\n",tNF,pos);
221 for(v=L->d+1;v<=FI->nVmax;v++) if(FI->nNUC[v])
222 for(nu=1;nu<=FI->NUCmax;nu++) for(j=0;j<FI->NFnum[v][nu];j++)
223 Test_ucNF(&L->d,&v,&nu,&FI->NF[v][nu][nu*j]); puts("o.k.");
224 */
225 }
Read_SubLat_Poly(FILE * F,NF_List * L)226 void Read_SubLat_Poly(FILE *F,NF_List *L)
227 { int i, j, pos=0; for(i=0;i<L->SLN;i++)
228 { unsigned char *C=&L->NewNF[pos]; L->SLp[i]=pos;
229 L->NewNF[pos++]=fgetc(F); L->NewNF[pos++]=fgetc(F);
230 for(j=0;j<C[1];j++) L->NewNF[pos++]=fgetc(F);
231 } assert(pos==L->NewNB);
232 }
Read_In_File(NF_List * S)233 void Read_In_File(NF_List *S)
234 { time_t Tstart=time(NULL); FILE *F=fopen(S->iname,"rb"); /* F=fopen */
235 printf("Reading In-File %s: ",S->iname);fflush(stdout);
236 if(F==NULL) {puts("Cannot open (read)!"); exit(0);}
237 Read_Honest_Poly(F,&S->In,S); S->NP=S->nSLP=0; assert(S->rd==0);
238 fclose(F); printf(" done (%ds)\n",(int) difftime(time(NULL),Tstart));
239 fflush(stdout);
240 }
Read_File_2_List(char * fn,NF_List * L)241 void Read_File_2_List(char *fn,NF_List *L) /* ... like Read_Aux_File */
242 { time_t Tstart=time(NULL); FILE *F=fopen(fn,"rb");
243 printf("Reading %s: ",fn);
244 if(F==NULL) {puts("Cannot open (read)!"); exit(0);}
245 Read_Honest_Poly(F,&L->Aux,L);
246 Read_SubLat_Poly(F,L);
247 printf(" done (%ds)\n",(int) difftime(time(NULL),Tstart));
248 fclose(F); fflush(stdout);
249 }
Read_Aux_File(NF_List * L)250 void Read_Aux_File(NF_List *L)
251 { time_t Tstart=time(NULL); FILE *F; /* F=fopen */
252 int NCalloc=strlen(L->oname)+strlen(SAVE_FILE_EXT)+ (USE_TMP_DIR ? 6 : 1);
253 char *auxfn = (char *) malloc(NCalloc);
254 if(USE_TMP_DIR) {strcpy(auxfn,"/tmp/"); strcat(auxfn,L->oname);}
255 else strcpy(auxfn,L->oname);
256 strcat(auxfn,SAVE_FILE_EXT);
257 F=fopen(auxfn,"rb"); printf("Reading %s: ",auxfn); fflush(stdout);
258 if(F==NULL) puts("No aux-file found!");
259 else
260 { Read_Honest_Poly(F,&L->Aux,L); Read_SubLat_Poly(F,L);
261 printf(" done (%ds)\n",(int) difftime(time(NULL),Tstart));
262 fclose(F); fflush(stdout);
263 #ifdef MOVE_SAVE_FILE /* inconsistent with USE_TMP_DIR !! */
264 {char *Mfn=(char *) malloc(1+strlen(L->oname)+strlen(MOVE_SAVE_FILE));
265 strcpy(Mfn,L->oname); strcat(Mfn,MOVE_SAVE_FILE);
266 assert(!rename(auxfn,Mfn)); free(Mfn);}
267 #endif
268 }
269 fflush(stdout); free(auxfn);
270 }
fputUI(unsigned int l,FILE * F)271 void fputUI(unsigned int l,FILE *F) /* write unsigned int to bin file */
272 { unsigned char A,B,C,D;
273 D=l % 256; l/=256; C=l % 256; l/=256; B=l % 256; l/=256; A=l;
274 fprintf(F,"%c%c%c%c",A,B,C,D);
275 }
TestMSbits(NF_List * S,PolyPointList * _P)276 void TestMSbits(NF_List *S, PolyPointList *_P)
277 { int i, v, nu, tc, peNF=0,peSM=0,peNM=0,
278 slNF=0, slSM=0,slNM=0,axSM=0; UPint axNF=0,axNM=0;
279 for(i=0;i<S->PEN;i++)
280 { unsigned char *C=&S->NewNF[S->PE[i].c]; int ms=C[2]%4; peNF++;
281 { /* int nv=*C,nuc=C[1]; Test_ucNF(&S->d,&nv,&nuc,&C[2]);*/ }
282 if(ms==0) (peSM)++; else if(ms!=3) (peNM)++;
283 }
284 for(i=0;i<S->PPEN;i++)
285 { unsigned char *C=&S->NewNF[S->PPE[i].pe.c]; int ms=C[2]%4; peNF++;
286 { /* int nv=*C,nuc=C[1]; Test_ucNF(&S->d,&nv,&nuc,&C[2]);*/ }
287 if(ms==0) (peSM)++; else if(ms!=3) (peNM)++;
288 }
289 for(v=S->d+1;v<=S->Aux.nVmax;v++) if(S->Aux.nNUC[v])
290 for(nu=1;nu<=S->Aux.NUCmax;nu++) if((tc=S->Aux.NFnum[v][nu]))
291 for(i=0;i<tc;i++)
292 { unsigned char *C=&S->Aux.NF[v][nu][nu*i]; int ms=C[0]%4; axNF++;
293 if(ms==0) (axSM)++; else if(ms!=3) (axNM)++;
294 }
295 printf("\nTEST: hNP=%lld hNF=%d+%lld=%d+%d nNM=%d+%d nSM=%d+%d\n",
296 S->NP-S->nSLP,S->PEN+S->PPEN,S->Aux.nNF,peNF,axNF,peNM,axNM,peSM,
297 axSM); fflush(stdout);
298 assert(S->NP-S->nSLP == 2*peNF-peNM-peSM + 2*axNF-axNM-axSM);
299 for(i=0;i<S->SLN;i++)
300 { unsigned char *C=&S->NewNF[S->SLp[i]]; int ms=C[2]%4; slNF++;
301 {int nv=*C,nuc=C[1]; Test_ucNF(&S->d,&nv,&nuc,&C[2],_P); }
302 if(ms==0) (slSM)++; else if(ms!=3) (slNM)++;
303 }
304 printf(" sl: slNP=%d slNF=%d=%d slNM=%d slSM=%d\n\n",
305 S->nSLP,S->SLN,slNF,slNM,slSM); fflush(stdout);
306 assert(S->nSLP == 2*slNF-slNM-slSM);
307 }
AuxCalcNumbers(NF_List * S,FInfoList * AI)308 void AuxCalcNumbers(NF_List *S,FInfoList *AI)
309 { int i,j; Init_FInfoList(AI); if(S->PPEN) Insert_PPent_into_Pent(S);
310 for(i=0;i<S->PEN;i++)
311 { unsigned char *C=&S->NewNF[S->PE[i].c]; int ms=C[2]%4; AI->nNF++;
312 ++(AI->NFnum[C[0]][C[1]]);
313 if(ms==0) (AI->nSM)++; else if(ms!=3) (AI->nNM)++;
314 }
315 for(i=S->d+1;i<=VERT_Nmax;i++) for(j=0;j<NUC_Nmax;j++)if(AI->NFnum[i][j])
316 { AI->nNUC[i]++; if(j>AI->NUCmax)AI->NUCmax=j;
317 AI->NB += j * AI->NFnum[i][j];
318 } /* TestMSbits(S); */
319 for(i=S->d+1;i<=VERT_Nmax;i++) if(AI->nNUC[i]) AI->nVmax=i;
320 }
321 /* uchar=unsigned char uint=unsigned int
322 * NP=#Polys NB=#Bytes #files=#nv's with NP>0 #lists=#(nv,nuc) with NP>0
323 *
324 * uchar rd k_1 ... k_rd // ... not in data base !!
325 * // big rd: "128+rd/128" "rd%128"
326 * uchar dim #files nv_max nuc_max
327 * uint #lists hNF hSM hNM hNB slNF slSM slNM slNB
328 * uchar v1 #nuc's with v1
329 * uchar nuc1 uint #NF(v1,nuc1) uchar nuc2 uint #NV(v1,nuc2) ...
330 * uchar v2 #nuc's with v2
331 * uchar nuc1 uint #NF(v2,nuc1) uchar nuc2 uint #NV(v2,nuc2) ...
332 * uchar "all hNF honest nf's"
333 * uchar "all slNF sublattice {nv nuc nf[]}'s"
334 */
Write_Bin_File(FILE * F,NF_List * L)335 void Write_Bin_File(FILE *F,NF_List *L)
336 { unsigned int i,j,fi=0,li=0,v,nu,tc=0,pen=0, tnb;
337 UPint NFnum[VERT_Nmax][NUC_Nmax];
338 int n; unsigned char *C, nVmax,NUCmax;
339 FInfoList AI, *_AI = (FInfoList *) &AI;
340
341 printf("%lld+%dsl %lldm+%ds %lldb",L->NP-L->nSLP,L->nSLP,
342 L->Aux.nNF - L->Aux.nNM - L->Aux.nSM + L->PEN + L->PPEN - L->peNM
343 - L->peSM, L->Aux.nSM + L->peSM,
344 /* nNF :: L->Aux.nNF+L->PEN+L->PPEN+L->SLN, */
345 L->Aux.NB+L->NewNB-L->RemNB-2*(L->PEN+L->PPEN));
346 if(L->rd) printf(" rd=%d",L->rd);
347 #ifdef USE_UNIT_ENCODE
348 if(L->NC>0) printf(" u%lld",L->NC*100 / (L->Aux.nNF+L->PEN+L->PPEN) );
349 #endif
350 AuxCalcNumbers(L,_AI);
351 Print_Expect(_AI);
352 fflush(stdout);
353 nVmax=max(L->Aux.nVmax,AI.nVmax); NUCmax=max(L->Aux.NUCmax,AI.NUCmax);
354 for(i=L->d+1;i<=nVmax;i++)
355 { int v_li=0;
356 for(j=1;j<=NUCmax;j++)
357 if((NFnum[i][j]=L->Aux.NFnum[i][j]+AI.NFnum[i][j])) v_li++;
358 li+=v_li;
359 if(v_li) fi++;
360 }
361 assert(L->rd <= MAX_REC_DEPTH);
362 if(L->rd < 128) fputc(L->rd,F);
363 else {fputc(128+(L->rd)/128,F); fputc((L->rd)%128,F);}
364 for(n=0;n<L->rd;n++) fputc(L->b[n],F);
365 fputc(L->d,F); fputc(fi,F); fputc(nVmax,F); fputc(NUCmax,F);
366 fputUI(li,F); fputUI(L->Aux.nNF+AI.nNF,F);
367 fputUI(L->Aux.nSM+AI.nSM,F); fputUI(L->Aux.nNM+AI.nNM,F);
368 #ifdef TEST_Write
369 printf("\nd=%d nV=%d nVmax=%d NUCmax=%d #lists=%d ",
370 L->d,fi,nVmax,NUCmax,li); fflush(stdout);
371 printf("%dnf %dsm %dnm %lldb sl: %d %d %d %lld\n",L->Aux.nNF+AI.nNF,
372 L->Aux.nSM+AI.nSM,L->Aux.nNM+AI.nNM,L->Aux.NB+AI.NB,L->SLN,
373 L->slSM,L->slNM,L->NewNB-L->RemNB-AI.NB-2*AI.nNF); fflush(stdout);
374 #endif
375 assert(L->NP-L->nSLP==2*AI.nNF-AI.nNM-AI.nSM
376 +2*L->Aux.nNF-L->Aux.nNM-L->Aux.nSM);
377 fputUI(tnb=L->Aux.NB+AI.NB,F);
378 fputUI(L->SLN,F); fputUI(L->slSM,F); fputUI(L->slNM,F);
379 assert(L->nSLP==2*L->SLN-L->slNM-L->slSM);
380 fputUI(L->NewNB-L->RemNB-AI.NB-2*AI.nNF,F);
381 for(v=L->d+1;v<=nVmax;v++) if(AI.nNUC[v]+L->Aux.nNUC[v])
382 { i=0; for(nu=1;nu<=NUCmax;nu++) if(NFnum[v][nu]) i++;
383 #ifdef TEST_Write
384 printf("v%d:n%d ",v,i);
385 #endif
386 fputc(v,F); fputc(i,F); /* v #nuc(v) */
387 for(nu=1;nu<=NUCmax;nu++) if(NFnum[v][nu])
388 { fputc(nu,F); fputUI(NFnum[v][nu],F); /* nuc #NF */
389 }
390 }
391 j=0; for(v=L->d+1;v<=nVmax;v++) if(AI.nNUC[v] || L->Aux.nNUC[v])
392 for(nu=1;nu<=NUCmax;nu++) if((tc=NFnum[v][nu]))
393 { unsigned int k; pen+=AI.NFnum[v][nu]; /* n -> i; lipos -> j */
394 for(i=0; i < (L->Aux.NFnum[v][nu]); i++)
395 { if(j<pen) while(L->PE[j].n == i)
396 { C=&L->NewNF[L->PE[j++].c+2];for(k=0;k<nu;k++)fputc(C[k],F);
397 tc--; tnb-=nu; if(j==pen) break;
398 }
399 C=&L->Aux.NF[v][nu][nu*i]; for(k=0;k<nu;k++) fputc(C[k],F);
400 tnb-=nu; tc--;
401 }
402 while(j<pen)
403 { C=&L->NewNF[L->PE[j++].c+2]; for(k=0;k<nu;k++)fputc(C[k],F);
404 tnb-=nu; tc--;
405 } assert(tc==0);
406 } assert(tnb==0);
407 for(n=0;n<L->SLN;n++) /* =S->nSLP */
408 { C=&L->NewNF[L->SLp[n]]; nu=C[1]+2; for(j=0;j<nu;j++) fputc(C[j],F);
409 tc += C[1] + 2;
410 } assert(tc==L->NewNB-L->RemNB-AI.NB-2*AI.nNF);
411 #ifdef More_File_IO_Data
412 printf("\n nv<=%d nuc<=%d files=%d lists=%d nNF=%d NB=%lld ..",
413 nVmax,NUCmax,fi,li,AI.nNF,AI.NB);
414 #endif
415 if(NULL!=L->Aux.NFli) {free(L->Aux.NFli); L->Aux.NFli=NULL;}
416 }
417
Write_Aux_File(NF_List * S)418 void Write_Aux_File(NF_List *S)
419 { time_t Tstart=time(NULL); FILE *F;
420 int NCalloc=strlen(S->oname)+strlen(SAVE_FILE_EXT)+ (USE_TMP_DIR ? 6 : 1);
421 char *auxfn = (char *) malloc(NCalloc);
422 #ifdef TEMP_FILE_EXT
423 char *tmpfn = (char *) malloc(1+strlen(S->oname)+4);
424 strcpy(tmpfn,S->oname); strcat(tmpfn,".tmp");
425 #else
426 char *tmpfn = auxfn;
427 #endif
428 if(USE_TMP_DIR) {strcpy(auxfn,"/tmp/"); strcat(auxfn,S->oname);}
429 else strcpy(auxfn,S->oname);
430 strcat(auxfn,SAVE_FILE_EXT);
431 F=fopen(tmpfn,"wb"); printf("Writing %s: ",auxfn); fflush(stdout);
432 if(F==NULL){puts("Cannot open!");exit(0);}
433 Write_Bin_File(F,S); if(ferror(F)) {puts("File ERROR!!");exit(0);}
434 fclose(F); printf(" done: %ds\n",(int) difftime(time(NULL),Tstart));
435 fflush(stdout); free(auxfn);
436 #ifdef TEMP_FILE_EXT
437 rename(tmpfn,auxfn); free(tmpfn);
438 #endif
439 }
440
Write_List_2_File(char * fn,NF_List * S)441 void Write_List_2_File(char *fn,NF_List *S)
442 { time_t Tstart=time(NULL); FILE *F=fopen(fn,"wb");
443 printf("Writing %s: ",fn); fflush(stdout);
444 if(F==NULL){puts("Cannot open!");exit(0);}
445 Write_Bin_File(F,S); if(ferror(F)) {puts("File ERROR!!");exit(0);}
446 fclose(F); printf(" done: %ds\n",(int) difftime(time(NULL),Tstart));
447 S->SAVE=time(NULL); fflush(stdout);
448 }
ReAlloc_SortList(NF_List * L)449 void ReAlloc_SortList(NF_List *L)
450 { Write_Aux_File(L); Read_Aux_File(L); /* Test_SLnbMS(L); TestMSbits(L); */
451 L->SAVE=time(NULL);
452 }
CheckLastSaveTime(NF_List * L,int maxsec)453 void CheckLastSaveTime(NF_List *L, int maxsec)
454 { if((int)difftime(time(NULL),L->SAVE) > maxsec)
455 { Print_Statistics(L); ReAlloc_SortList(L);
456 }
457 }
PEntComp(int * FIpos,int * nv,int * nuc,unsigned char * uc,PEnt * pe,NF_List * S)458 int PEntComp(int *FIpos,int *nv,int *nuc,unsigned char *uc, /* pe/pos - uc */
459 PEnt *pe,NF_List *S)
460 { unsigned char *C=&S->NewNF[pe->c];
461 int i = (int) C[0] - *nv; if(i) return i;
462 i = (int) C[1] - *nuc; if(i) return i;
463 #ifdef ACCEL_PEntComp
464 i = pe->n - *FIpos; if(i) return i;
465 #endif
466 return RIGHTminusLEFT(uc,&(C[2]),nuc);
467 }
SearchPEntList(int * nv,int * nuc,unsigned char * uc,int * PEpos,int * FIpos,NF_List * S)468 int SearchPEntList(int *nv,int *nuc,unsigned char *uc, /* 0 :: new */
469 int *PEpos, int *FIpos, NF_List *S) /* -1 exists(SL) */
470 { int i, n0, n1; /* +1 ex. honest */
471 *PEpos=0; if(!S->PEN) return 0;
472 if((i = PEntComp(FIpos,nv,nuc,uc,&S->PE[n0=0],S))) /* i = PE[0] - uc */
473 { if(i>0) { return 0;}}
474 else return 1;
475 if((i = PEntComp(FIpos,nv,nuc,uc,&S->PE[n1=S->PEN-1],S)))
476 { if(i<0) { *PEpos=S->PEN; return 0;}}
477 else { *PEpos=n1; return 1; }
478 while(n1>n0+1)
479 { *PEpos=(n0+n1)/2; i=PEntComp(FIpos,nv,nuc,uc,&S->PE[*PEpos],S);
480 if(i) if(i>0) n1=*PEpos; else n0=*PEpos;
481 else return 1;
482 }
483 *PEpos=n1; return 0; /* exists: return 1; else PEpos=where it would go */
484 }
Insert_PPent_into_Pent(NF_List * S)485 void Insert_PPent_into_Pent(NF_List *S)
486 { int Mpos=S->PEN, Spos=S->PPEN, pos=Mpos+Spos;
487 if(S->PPEN<=0)
488 { puts("This should not happen in Insert_PPent_into_Pent");exit(0);
489 } assert(pos<=SAVE_INC+SL_Nmax);
490 while(Spos--)
491 { if(S->PPE[Spos].n == Mpos) S->PE[--pos] = S->PPE[Spos].pe;
492 else break;
493 }
494 if(++Spos) while(Spos--)
495 { while(S->PPE[Spos].n < Mpos) S->PE[--pos] = S->PE[--Mpos];
496 S->PE[--pos] = S->PPE[Spos].pe;
497 }
498 S->PEN+=S->PPEN; S->PPEN=0; /* Test_NF_List(S); */
499 }
PPEntComp(int * FIpos,int * PEpos,int * nv,int * nuc,unsigned char * uc,PPEnt * ppe,NF_List * S)500 int PPEntComp(int *FIpos,int *PEpos,int *nv,int *nuc,unsigned char *uc,
501 PPEnt *ppe,NF_List *S) /* ppe/pos - uc */
502 { unsigned char *C=&S->NewNF[ppe->pe.c];
503 int i = (int) C[0] - *nv; if(i) return i;
504 i = (int) C[1] - *nuc; if(i) return i;
505 #ifdef ACCEL_PEntComp
506 i = ppe->n - *PEpos; if(i) return i; if((i=ppe->pe.n -*FIpos)) return i;
507 #endif
508 return RIGHTminusLEFT(uc,&(C[2]),nuc);
509 }
InsertPNFintoPPEntList(int * spos,int * mpos,int * lpos,int * nv,int * nuc,unsigned char * uc,NF_List * S)510 void InsertPNFintoPPEntList(int *spos, int *mpos, int *lpos,
511 int *nv, int *nuc, unsigned char *uc, NF_List *S)
512 { int l=S->PPEN++; unsigned char *C=&(S->NewNF[S->NewNB]);
513 PPEnt *ppe=&S->PPE[*spos]; static int AddListLength;
514 if(0==AddListLength) AddListLength=ADD_LIST_LENGTH;
515 #ifdef USE_UNIT_ENCODE
516 if((*uc % 8) > 3) S->NC ++;
517 #endif
518 if(*nuc+2+S->NewNB > S->ANB)
519 { printf("increase CperR_MAX or write/read %s\n",S->iname); exit(0);
520 }
521 while(*spos < (l--)) S->PPE[l+1]=S->PPE[l];
522 ppe->n=*mpos; ppe->pe.n=*lpos; ppe->pe.c=S->NewNB;
523 C[0]=*nv; C[1]=*nuc; C=&C[2]; for(l=0;l<*nuc;l++) C[l]=uc[l];
524 S->NewNB += *nuc+2; S->NP++; if(*C % 4) S->peNM++; else S->peSM++;
525 if(S->PPEN==AddListLength) Insert_PPent_into_Pent(S);
526 }
SearchPPEntList(int * nv,int * nuc,unsigned char * uc,int * Ppos,int * Mpos,int * Lpos,NF_List * S)527 int SearchPPEntList(int *nv,int *nuc,unsigned char *uc, /* new :: 0 */
528 int *Ppos, int *Mpos,int *Lpos, NF_List *S) /* exists :: 1 */
529 { int i, n0, n1;
530 *Ppos=0; if(!S->PPEN) return 0;
531 if((i = PPEntComp(Lpos,Mpos,nv,nuc,uc,&S->PPE[n0=0],S)))
532 { if(i>0) { *Ppos=0; return 0;}}
533 else return 1;
534 if((i = PPEntComp(Lpos,Mpos,nv,nuc,uc,&S->PPE[n1=S->PPEN-1],S)))
535 { if(i<0) { *Ppos=S->PPEN; return 0;}}
536 else {*Ppos=n1; return 1;}
537 while(n1>n0+1)
538 { *Ppos=(n0+n1)/2; i=PPEntComp(Lpos,Mpos,nv,nuc,uc,&S->PPE[*Ppos],S);
539 if(i) if(i>0) n1=*Ppos; else n0=*Ppos;
540 else return 1;
541 }
542 *Ppos=n1; return 0; /* exists: return 1; else Ppos=where it would go */
543 }
SearchFIList(int * v,int * nu,unsigned char * uc,int * pos,FInfoList * FI)544 int SearchFIList(int *v,int *nu,unsigned char *uc, int *pos,FInfoList *FI)
545 { int i, n0, n1; *pos=0; if(!FI->nNF) return 0;
546 if(!FI->NFnum[*v][*nu]) return 0;
547 if((i = RIGHTminusLEFT(uc,&FI->NF[*v][*nu][n0=0],nu)))
548 { if(i>0) return 0;}
549 else return 1;
550 if((i = RIGHTminusLEFT(uc,
551 &FI->NF[*v][*nu][(*nu)*(n1=FI->NFnum[*v][*nu]-1)],nu)))
552 { if(i<0) { *pos=n1+1; return 0;}}
553 else {*pos=n1; return 1;}
554 while(n1>n0+1)
555 { *pos=(n0+n1)/2; i=RIGHTminusLEFT(uc,
556 &FI->NF[*v][*nu][(*nu)*(*pos)],nu);
557 if(i) if(i>0) n1=*pos; else n0=*pos;
558 else return 1;
559 }
560 *pos=n1; return 0; /* exists: return 1; else pos=where it would go */
561 }
SL_Comp(int * nv,int * nuc,unsigned char * uc,int * lipos,NF_List * S)562 int SL_Comp(int *nv,int *nuc,unsigned char *uc,int *lipos,NF_List *S)
563 { unsigned char *C=&S->NewNF[*lipos]; /* lipos - uc */
564 int i = (int) C[0] - *nv; if(i) return i;
565 i = (int) C[1] - *nuc; if(i) return i;
566 return RIGHTminusLEFT(uc,&(C[2]),nuc);
567 }
SearchSL_List(int * nv,int * nuc,unsigned char * uc,int * SLnum,NF_List * S)568 int SearchSL_List(int *nv,int *nuc,unsigned char *uc, int *SLnum,NF_List *S)
569 { int i, n0, n1; /* 1 = exists (mirror or straight) */
570 *SLnum=0; if(!S->SLN) return 0;
571 if((i = SL_Comp(nv,nuc,uc,&S->SLp[n0=0],S)))
572 { if(i>0) { return 0;}}
573 else return 1;
574 if((i = SL_Comp(nv,nuc,uc,&S->SLp[n1=S->SLN-1],S)))
575 { if(i<0) { *SLnum=S->SLN; return 0;}}
576 else { *SLnum=n1; return 1;}
577 while(n1>n0+1)
578 { *SLnum=(n0+n1)/2; i=SL_Comp(nv,nuc,uc,&S->SLp[*SLnum],S);
579 if(i) if(i>0) n1=*SLnum; else n0=*SLnum;
580 else return 1;
581 }
582 *SLnum=n1; return 0; /* exists: return 1; else SLnum=where it would go */
583 }
SL_List_Insert(int * nv,int * nuc,unsigned char * uc,int * SLnum,NF_List * S)584 void SL_List_Insert(int *nv,int *nuc,unsigned char *uc,int *SLnum,NF_List *S)
585 { int l=S->SLN++; unsigned char *C=&(S->NewNF[S->NewNB]);
586 if(*nuc+2 + S->NewNB > S->ANB)
587 { printf("increase CperR_MAX or write/read %s\n",S->iname); exit(0);
588 }
589 while(*SLnum < (l--)) S->SLp[l+1]=S->SLp[l];
590 S->SLp[*SLnum]=S->NewNB;
591 C[0]=*nv; C[1]=*nuc; C=&C[2]; for(l=0;l<*nuc;l++) C[l]=uc[l];
592 S->NewNB+=*nuc+2; S->NP++;S->nSLP++; if(*C % 4)S->slNM++; else S->slSM++;
593 if(S->SLN==SL_Nmax) {puts("Increase SL_Nmax");exit(0);}
594 }
SL_List_Remove(int * NV,int * nUC,int * SLnum,NF_List * S)595 void SL_List_Remove(int *NV,int *nUC, /* unsigned char *UC, */
596 int *SLnum, NF_List *S)
597 { int l=*SLnum; unsigned char *C=&S->NewNF[S->SLp[l]];
598 S->RemNB += (2+C[1]); if(C[2] % 4) S->slNM--; else S->slSM--;
599 assert((*NV==C[0]) && (*nUC==C[1]));
600 /* assert(0==RIGHTminusLEFT(UC,&C[2],nUC)); */
601 while ((++l) < S->SLN) S->SLp[l-1]=S->SLp[l];
602 S->SLN--; S->NP--; S->nSLP--;
603 }
604
605 #define MirTest(A,B) ((((A)+(B)) % 4) != 3) /* 1=same, 0=only mirror */
606 /* ucNF_Sort_Add = 1 = continue = V is honest and new ;
607 * 0 = SL or Ex(honest)
608 * Search...List = 1 : 0 : -1 iff honest : new : sublattice
609 */
ucNF_Sort_Add(int * NV,int * nUC,unsigned char * UC,NF_List * S)610 int ucNF_Sort_Add(int *NV, int *nUC, unsigned char *UC, NF_List *S)
611 { int InPos=0, FIpos=0, PEpos=0, NEWpos=0, NewH=0, NewSL=0;
612 #ifdef INCREMENTAL_WRITE
613 if(SearchFIList(NV,nUC,UC,&InPos, &S->In))
614 if(MirTest(*UC,S->In.NF[*NV][*nUC][(*nUC)*InPos]))
615 { if(!S->SL) S->hc++; goto Ret0;}
616 #endif
617 if(*S->dbname) if(Is_in_DB(NV,nUC,UC,S)) {if(!S->SL) S->hc++; goto Ret0;}
618
619 if(SearchFIList(NV,nUC,UC,&FIpos,&S->Aux))
620 { unsigned char *c=&(S->Aux.NF[*NV][*nUC][(*nUC)*FIpos]);
621 if(MirTest(*UC,*c)) goto Ret0; /* if(S->con==S->AP) S->hc++; */
622 else if(!S->SL)
623 { (*c) += ((*UC) % 4); S->Aux.nNM--; S->NP++; NewH=1;
624 } /* Mir(honest) on Aux => not on New */
625 }
626 else if(SearchPEntList (NV,nUC,UC,&PEpos,&FIpos,S))
627 { unsigned char *c=&(S->NewNF[S->PE[PEpos].c+2]);
628 if(MirTest(*UC,*c)) goto Ret0;
629 else if(!S->SL)
630 { (*c) += ((*UC) % 4); S->peNM--; S->NP++; NewH=1;
631 } /* Mir(honest) on PE => not on PPE */
632 }
633 else if(SearchPPEntList(NV,nUC,UC,&NEWpos,&PEpos,&FIpos,S))
634 { unsigned char *c=&(S->NewNF[S->PPE[NEWpos].pe.c+2]);
635 if(MirTest(*UC,*c)) goto Ret0;
636 else if(!S->SL)
637 { (*c) += ((*UC) % 4); S->peNM--; S->NP++; NewH=1;
638 }
639 }
640 else if(!S->SL)
641 { InsertPNFintoPPEntList(&NEWpos,&PEpos,&FIpos,NV,nUC,UC,S); NewH=1;
642 }
643 assert(NewH + S->SL == 1);
644 if(SearchSL_List(NV,nUC,UC,&NEWpos,S))
645 { unsigned char *c=&(S->NewNF[S->SLp[NEWpos]+2]);
646 if(NewH)
647 { if(MirTest(*UC,*c))
648 { if((*c % 4) != 3) SL_List_Remove(NV,nUC, /*UC,*/ &NEWpos,S);
649 else { (*c) -= (*UC % 4); S->slNM++; S->NP--; S->nSLP--;}
650 NewSL=-1;
651 }
652 }
653 else /* SubLattice case */
654 { if(MirTest(*UC,*c)) goto Ret0;
655 else { (*c) += ((*UC) % 4); S->NP++; S->slNM--;S->nSLP++;NewSL=1;}
656 }
657 }
658 else if(S->SL) {assert(NewH==0);
659 NewSL=1; SL_List_Insert(NV,nUC,UC,&NEWpos,S);}
660
661 if(NewH+NewSL) if(S->NP % WATCHREF == 0)
662 { Print_Statistics(S);
663 if(S->NP - S->savedNP > SAVE_INC - WATCHREF)
664 /* printf("NP=%d PEN=%d peNM=%d peSM=%d\n", S->NP,S->PEN,S->peNM,
665 S->peSM);fflush(stdout);
666 if(2 *(S->PEN+S->PPEN) - S->peNM - S->peSM > SAVE_INC - WATCHREF)*/
667 { ReAlloc_SortList(S); S->con=0; }
668 } if(S->NP%1000 == 0) CheckLastSaveTime(S,(19*GOOD_SAVE_TIME)/20);
669 CheckLastSaveTime(S,GOOD_SAVE_TIME); return NewH;
670 Ret0: CheckLastSaveTime(S,FORCE_SAVE_TIME); return 0;
671 }
672 /* ============================================================== */
673
674
675 /* ============================================================== */
Print_Statistics(NF_List * _L)676 void Print_Statistics(NF_List *_L)
677 { clock_t CLOCK=clock(); time_t DATE=time(NULL);
678 int CPUsec=(CLOCK-_L->CLOCK)/CLOCKS_PER_SEC; /* CLOCKS_PER_SEC::10^6 */
679 int REALsec= (int) difftime(DATE,_L->TIME), NFsec; /* int IPperNF; */
680 char bni[2]; int BNI=(_L->nNF>2000000) ? 1000000 : 1000;
681 strcpy(bni,(BNI==1000) ? "k" : "M");
682 printf("%dkR-%d %dMB %d%sIP %d%sNF-%dk ", (int)_L->NP/1000,
683 (int)_L->nSLP,(int) ((_L->NewNB-_L->RemNB+_L->Aux.NB)/1000000), (int)
684 (_L->nIP/BNI),bni,(int)_L->nNF/BNI,bni,(int)_L->nSLNF/1000);
685 printf("%d%s%d v%dr%d f%dr%d %db%d",
686 _L->Nmin, _L->hc ? "-" : "_", _L->Nmax, _L->VN,_L->V, _L->FN,_L->F,
687 _L->Xdif, _L->Xnuc);
688 #ifdef INCREMENTAL_TIME
689 /* if(_L->NF_Time>0)IPperNF=_L->IP_Time/_L->NF_Time; else IPperNF=-1; */
690 NFsec=_L->NF_Time/CLOCKS_PER_SEC;
691 CPUsec=(_L->IP_Time+_L->NF_Time)/CLOCKS_PER_SEC;
692 if(CPUsec<1000)printf(" %ds %du %dn",REALsec,CPUsec,NFsec);
693 else if((CPUsec)<60000)
694 printf(" %dm %du %dn",REALsec/60,CPUsec/60,NFsec/60);
695 else printf(" %dh %du %dn",REALsec/3600,CPUsec/3600,NFsec/3600);
696 #else
697 CPUsec=(CLOCK-_L->CLOCK)/CLOCKS_PER_SEC;
698 if(CPUsec<1000)printf(" %ds %du",REALsec,CPUsec);
699 else if((CPUsec)<60000) printf(" %dm %du",REALsec/60,CPUsec/60);
700 else printf(" %dh %du",REALsec/3600,CPUsec/3600);
701 #endif
702 _L->Nmin=1000; _L->Nmax=0; puts(""); fflush(stdout);
703 }
Print_Expect(FInfoList * L)704 void Print_Expect(FInfoList *L)
705 { int s=L->nSM; Along m=L->nNF-L->nNM-L->nSM, p=L->nNF+m; double x=p;
706 /* if(d<2) {if(m)printf(" p^2/2m=%d",(p*p)/(2*m));
707 if(s)printf(" pp/ss=%d",(p*p)/(s*s));} */
708 /* if(m) {x*=(x+s)/(2*m);printf(" pP/2m=%g",x);} */
709 if(m) {x=p+s; x*=x/(2*m+s);printf(" pp/2m=%g",x);}
710 if(s) {x=p; x/=s; x*=x; printf(" pp/ss=%g",x);}
711 }
712 /* ============================================================== */
713
714
715 /* ============================================================== */
716 /* SearchLongList:= { new(below *pos)::0, found::1 } */
717 /* SearchPEntLists:={ new::0, found::1 } */
718 /* Add_NF_to_List() assumes that _V and _F are assigned
719 *
720 * if(SL) { if(new) "add to dishonest"; return 0; } ::stop
721 * if(exists on honest) return 0; ::stop
722 * if(exists on dishonest) { "move to honest"; return 1; } ::cont
723 * "add to honest"; return 1; ::cont
724 *
725 * Make_Poly_NF (Make_VPM_NF; Make_PNF_from_VNF;);
726 * PNF_Sort_Add (NF[][],*d,*np,*nv,*nf,*S) := { List::+1, new::0, PEnt::-1 }
727 * SearchLongList:= { new(below *pos)::0, found::1 }
728 * SearchPEntLists:={ new::0, found::1 }
729 */
Add_NF_to_List(PolyPointList * _P,VertexNumList * _V,EqList * _E,NF_List * _L)730 int Add_NF_to_List(PolyPointList *_P, VertexNumList *_V, EqList *_E,
731 NF_List *_L)
732 { unsigned char UC[NB_MAX]; int nUC, NV;
733 int NewNF;
734 #ifdef INCREMENTAL_TIME
735 long long cpuT=clock(), incT=cpuT - _L->CLOCK;
736 if(incT>0) _L->IP_Time+=incT;
737 _L->CLOCK=cpuT;
738 #endif
739 if((_E->ne)>_L->F) _L->F=_E->ne; /* check vertex numbers */
740 if((_V->nv)>_L->V) _L->V=_V->nv;
741 if((_E->ne > VERT_Nmax)||(_V->nv > VERT_Nmax))
742 { printf("Increase VERT_Nmax: f=%d v=%d\n",_E->ne,_V->nv);exit(0);
743 }
744 if(_P->np>_L->Nmax) _L->Nmax=_P->np;
745 if(_P->np<_L->Nmin) _L->Nmin=_P->np;
746 _L->nNF++; if(_L->SL)_L->nSLNF++;
747
748 VF_2_ucNF(_P, _V, _E, &NV, &nUC, UC);
749
750 NewNF=ucNF_Sort_Add(&NV, &nUC, UC, _L); /* 1::new::cont. */
751
752 #ifdef INCREMENTAL_TIME
753 cpuT=clock(), incT=cpuT - _L->CLOCK;
754 if(incT>0) _L->NF_Time+=incT;
755 _L->CLOCK=cpuT;
756 #endif
757 return NewNF;
758 }
759
760 /* #R sl hit #IP #NF TIME( r > u > AddNF) */
Print_Weight_Info(CWS * W,NF_List * _L)761 void Print_Weight_Info(CWS *W,NF_List *_L)
762 { int i,j;
763 #if POLY_Dmax>3
764 Print_Statistics(_L);
765 #endif
766 for(i=0;i<W->nw;i++)
767 { fprintf(outFILE,"%ld ",W->d[i]);
768 for(j=0;j<W->N;j++) fprintf(outFILE,"%ld ", W->W[i][j]);
769 }
770 for(i=0;i<W->nz;i++)
771 { fprintf(outFILE,"/Z%d: ",W->m[i]);
772 for(j=0;j<W->N;j++) fprintf(outFILE,"%d ",W->z[i][j]);
773 }
774 fprintf(outFILE,"R=%lld +%dsl hit=%d IP=%d NF=%d (%d)\n",
775 _L->NP-_L->nSLP, _L->nSLP,_L->hc,_L->nIP,_L->nNF,_L->nSLNF);
776 fflush(outFILE);
777 }
778
779 /* =============== compression package =================== */
780
781 #if (INT_MAX != 2147483647) /* UINT_MAX == 4294967295 */
782 #error use other date types /* ULLONG_MAX == 18446744073709551615 */
783 #endif
784 #define UCM 256 /* Unsigned Char Modulo (= max+1) */
785 #define USM 65536 /* Unsigned Short Modulo (= max+1) */
786 #define Nint_XLong (NB_MAX + 1) / 2
787 #define NX(d,v) ((d)*(v)-((d)*(d-1))/2)
788
789 #ifdef USE_UNIT_ENCODE
790 #define UNIT_NX(d,v) ((d)*((v)-(d))+1)
791 #define UNIT_OFF (8)
792 #else
793 #define UNIT_OFF (4)
794 #endif
795
796 #ifndef LL_BASE
797 #define LL_BASE 32767 /* limit for 64-bit BaseGetInt: cf. NF of */
798 #endif /* 3198174 49 1723 74375 456882 1066058 1599087 */
799
800 typedef struct {int n; unsigned short x[Nint_XLong];} UXLong;
801
XPrint(UXLong * X)802 void XPrint(UXLong *X)
803 { int i; printf("X.n=%d:",X->n);for(i=0;i<X->n;i++)printf(" %d",X->x[i]);
804 }
VPrint(int * d,int * v,Long V[POLY_Dmax][VERT_Nmax])805 void VPrint(int *d,int *v,Long V[POLY_Dmax][VERT_Nmax])
806 { int i,j; puts("");for(i=0;i<*d;i++)
807 { for(j=0;j<*v;j++)printf("%4ld ",V[i][j]);puts("");}
808 }
LLBasePutInt(int * base,int * x,UXLong * X)809 void LLBasePutInt(int *base, int *x, UXLong *X) /* X = X * base + x */
810 { int i=0; long long z, ad=*x;
811 while(i<X->n)
812 { z=ad+(*base)*((long long)X->x[i]); X->x[i++]=z% USM; ad=z / USM;
813 }
814 while(ad) { assert(X->n < Nint_XLong); X->x[X->n++]=ad% USM; ad/=USM; }
815 }
LLBaseGetInt(int * base,int * x,UXLong * X)816 void LLBaseGetInt(int *base, int *x, UXLong *X) /* x = X mod b, X /= b */
817 { int i; long long a, zq, zr;
818 if(!X->n) {*x=0; return;} a=X->x[(i=X->n-1)];
819 while(i) { zq=a/(*base); zr=a%(*base); X->x[i]=zq; a=zr*USM+X->x[--i]; }
820 zq=a/(*base); zr=a%(*base); X->x[i]=zq; *x=zr;
821 for(i=X->n-1; 0<=i; i--) if(X->x[i]) return; else X->n--;
822 }
BasePutInt(int * base,int * x,UXLong * X)823 void BasePutInt(int *base, int *x, UXLong *X) /* X = X * base + x */
824 { int i=0, z, ad=*x; if(*base > LL_BASE) {LLBasePutInt(base,x,X); return;}
825 while(i<X->n)
826 { z=ad+(*base)*X->x[i]; X->x[i++]=z% USM; ad=z / USM;
827 }
828 if(ad) { assert(X->n <= (Nint_XLong-1)); X->x[X->n++]=ad; }
829 }
BaseGetInt(int * base,int * x,UXLong * X)830 void BaseGetInt(int *base, int *x, UXLong *X) /* x = X mod b, X /= b */
831 { int i, a; div_t z; if(*base > LL_BASE) {LLBaseGetInt(base,x,X); return;}
832 if(!X->n) {*x=0; return;} a=X->x[(i=X->n-1)];
833 while(i) { z=div(a,*base); X->x[i]=z.quot; a=z.rem*USM+X->x[--i]; }
834 z=div(a,*base); X->x[i]=z.quot; *x=z.rem;
835 for(i=X->n-1; 0<=i; i--) if(X->x[i]) return; else X->n--;
836 }
BminOff(Long V[POLY_Dmax][VERT_Nmax],int * d,int * v,int * off,int * bmin)837 int BminOff(Long V[POLY_Dmax][VERT_Nmax], int *d, int *v, int *off, int *bmin)
838 { int i, j, moff=0, vmax=0;
839 for(i=0;i<*d;i++) for(j=i;j<*v;j++)
840 { Long x=V[i][j]; if(moff>x) moff=x;
841 if(vmax<x) vmax=x;
842 }
843 *off=-moff; *bmin=*off+vmax+1;
844 #ifdef USE_UNIT_ENCODE
845 for(i=0;i<*d;i++) for(j=i;j<*d;j++) if(V[i][j]!=(i==j)) return 0;
846 return 1;
847 #else
848 return 0;
849 #endif
850 }
AuxBase2nUC(int * base,int * nx,int * nuc)851 void AuxBase2nUC(int *base, int *nx, int *nuc) /* b^nx-2=x0+x1*USM+... */
852 { int i, m=*base-1; UXLong X; X.n=0; /* m=X->n-1 => ... + xm*USM^xm */
853 for(i=0;i<*nx-1;i++) BasePutInt(base,&m,&X);
854 m--; BasePutInt(base,&m,&X);
855 i=UNIT_OFF; m=i-1; BasePutInt(&i,&m,&X); *nuc=2*X.n-(X.x[X.n-1]<256);
856 }
AuxNextGoodBase(int * v,int * nx,Base_List * BL)857 void AuxNextGoodBase(int *v,int *nx,Base_List *BL) /* nx= v*d -d*(d-1)/2 */
858 { int *bo=&(BL->base[*v][BL->v[*v]-1]), nuct,
859 *bn=&(BL->base[*v][BL->v[*v]]), *nucn=&(BL->nuc[*v][BL->v[*v]]);
860 *bn=*bo+1; AuxBase2nUC(bn,nx,nucn);
861 do { (*bn)++; AuxBase2nUC(bn,nx,&nuct); } while (nuct==*nucn);
862 (*bn)--; BL->v[*v]++;
863 }
Init_BaseList(Base_List ** bl,int * d)864 void Init_BaseList(Base_List **bl, int *d) /* once: alloc and init BaseList */
865 { static Base_List *BL=NULL; int i; /* always: set *bl=BL=&BaseList */
866 if(BL!=NULL) {*bl=BL; return;}
867 *bl = BL = (Base_List *) malloc(sizeof(Base_List)); assert(BL!=NULL);
868 for(i=*d; i<=VERT_Nmax; i++)
869 { int nx= NX(*d,i), bn=3, nuco, nucn;
870 BL->v[i]=1; AuxBase2nUC(&bn,&nx,&nuco); nucn=BL->nuc[i][0]=nuco;
871 while(nucn==nuco) { bn++; AuxBase2nUC(&bn,&nx,&nucn); }
872 BL->base[i][0]=bn-1;
873 }
874 }
Bmin2BaseUCn(int * d,int * v,int * bmin,int * base,int * nuc)875 void Bmin2BaseUCn(int *d, int *v, int *bmin, int *base, int *nuc)
876 { int i, *n, nx = NX(*d,*v);
877 Base_List *BL; Init_BaseList(&BL,d); n=&(BL->v[*v]); i=*n-1;
878
879 if(*bmin<=BL->base[*v][i])
880 { while(i > 0 && *bmin <= (BL->base[*v][i-1])) --i;
881 }
882 else
883 { do AuxNextGoodBase(v,&nx,BL); while(*bmin > (BL->base[*v][++i]));
884 }
885 *base=BL->base[*v][i]; *nuc=BL->nuc[*v][i];
886 }
NUCtoBase(int * d,int * v,int * nuc,int * Base)887 void NUCtoBase(int *d, int *v, int *nuc, int *Base)
888 { int i, *n, nx = NX(*d,*v);
889 Base_List *BL; Init_BaseList(&BL,d); n=&(BL->v[*v]); i=*n-1;
890 if(*nuc<=BL->nuc[*v][i])
891 { while(i > 0 && *nuc <= (BL->nuc[*v][i-1])) --i;
892 }
893 else
894 { do AuxNextGoodBase(v,&nx,BL); while(*nuc > (BL->nuc[*v][++i]));
895 }
896 assert(*nuc==BL->nuc[*v][i]); *Base=BL->base[*v][i];
897 }
898
899 #ifdef USE_UNIT_ENCODE
UNIT_Init_BaseList(Base_List ** bl,int * d)900 void UNIT_Init_BaseList(Base_List **bl, int *d) /* once: ainit BaseList */
901 { static Base_List *BL=NULL; int i; /* always: set *bl=BL=&BaseList */
902 if(BL!=NULL) {*bl=BL; return;}
903 *bl = BL = (Base_List *) malloc(sizeof(Base_List)); assert(BL!=NULL);
904 for(i=*d; i<=VERT_Nmax; i++)
905 { int nx= UNIT_NX(*d,i), bn=3, nuco, nucn;
906 BL->v[i]=1; AuxBase2nUC(&bn,&nx,&nuco); nucn=BL->nuc[i][0]=nuco;
907 while(nucn==nuco) { bn++; AuxBase2nUC(&bn,&nx,&nucn); }
908 BL->base[i][0]=bn-1;
909 }
910 }
UNIT_Bmin2BaseUCn(int * d,int * v,int * bmin,int * base,int * nuc)911 void UNIT_Bmin2BaseUCn(int *d, int *v, int *bmin, int *base, int *nuc)
912 { int i, *n, nx = UNIT_NX(*d,*v);
913 Base_List *BL; UNIT_Init_BaseList(&BL,d); n=&(BL->v[*v]); i=*n-1;
914 if(*bmin<=BL->base[*v][i])
915 { while(i > 0 && *bmin <= (BL->base[*v][i-1])) --i;
916 }
917 else
918 { do AuxNextGoodBase(v,&nx,BL); while(*bmin > (BL->base[*v][++i]));
919 }
920 *base=BL->base[*v][i]; *nuc=BL->nuc[*v][i];
921 }
UNIT_NUCtoBase(int * d,int * v,int * nuc,int * Base)922 void UNIT_NUCtoBase(int *d, int *v, int *nuc, int *Base)
923 { int i, *n, nx = UNIT_NX(*d,*v);
924 Base_List *BL; UNIT_Init_BaseList(&BL,d); n=&(BL->v[*v]); i=*n-1;
925 if(*nuc<=BL->nuc[*v][i])
926 { while(i > 0 && *nuc <= (BL->nuc[*v][i-1])) --i;
927 }
928 else
929 { do AuxNextGoodBase(v,&nx,BL); while(*nuc > (BL->nuc[*v][++i]));
930 }
931 assert(*nuc==BL->nuc[*v][i]); *Base=BL->base[*v][i];
932 }
933 #endif
934 /* UNIT :: ms+=4 */
AuxVnf2ucNF(Long nf[POLY_Dmax][VERT_Nmax],int * d,int * v,int * off,int * base,int * nuc,int * ms,unsigned char * UC)935 void AuxVnf2ucNF(Long nf[POLY_Dmax][VERT_Nmax], int *d, int *v,
936 int *off, int *base, int *nuc, int *ms, unsigned char *UC)
937 { int i, j, one=(*ms > 3);
938 UXLong X; X.n=0;
939 for(i=0;i<*d;i++) for(j = (one) ? *d : i ;j<*v;j++)
940 { int x=nf[i][j]+(*off); BasePutInt(base,&x,&X);
941 }
942 if(one) BasePutInt(base,off,&X);
943 i=UNIT_OFF; BasePutInt(&i,ms,&X);
944 for(i=X.n;(2*i)<(*nuc);i++) X.x[i]=0;
945 for(i=0;i<*nuc;i++)if(i%2) UC[i]=X.x[i/2] /UCM; else UC[i]=X.x[i/2] %UCM;
946 }
UCnf2vNF(int * d,int * v,int * nuc,unsigned char * uc,Long NF[POLY_Dmax][VERT_Nmax],int * MS)947 void UCnf2vNF(int *d, int *v, int *nuc, unsigned char *uc, /* IN */
948 Long NF[POLY_Dmax][VERT_Nmax], int *MS) /* OUT */
949 { int i, j, off=0, base, one;
950 UXLong X; X.n=*nuc/2; if(*nuc % 2) X.x[X.n++]=uc[*nuc-1];
951 for(i=0;i<*nuc/2;i++) X.x[i]=uc[2*i]+UCM*uc[2*i+1];
952 i=UNIT_OFF;
953 BaseGetInt(&i,MS,&X); one = (*MS > 3);
954 #ifdef USE_UNIT_ENCODE
955 if(one) {UNIT_NUCtoBase(d,v,nuc,&base); BaseGetInt(&base,&off,&X);} else
956 #endif
957 NUCtoBase(d,v,nuc,&base);
958 i=*d; while(i--)
959 { int m=(one) ? *d : i;
960 j=*v;while(m<(j--)) {int Iau;BaseGetInt(&base,&Iau,&X);NF[i][j]=Iau;}
961 j=i; while(j--) NF[i][j]=0;
962 }
963 #ifdef USE_UNIT_ENCODE
964 if(one) { for(i=0;i<*d;i++) for(j=*d;j<*v;j++) NF[i][j]-=off;
965 for(i=0;i<*d;i++) for(j=i;j<*d;j++) NF[i][j]=(i==j);
966 } else
967 #endif
968 { off=NF[0][0]-1; for(i=0;i<*d;i++) for(j=i;j<*v;j++) NF[i][j]-=off;}
969 for(i=0;i<*d;i++) for(j=0;j<i;j++) NF[i][j]=0;
970 }
971
RIGHTminusLEFT(unsigned char * ucL,unsigned char * ucR,int * nuc)972 int RIGHTminusLEFT(unsigned char *ucL, unsigned char *ucR, int *nuc)
973 { int i=*nuc;
974 while(--i) if(ucL[i]!=ucR[i]) return (ucL[i] < ucR[i]) ? 1 : -1;
975 if(((*ucL)/4)!=((*ucR)/4)) return (((*ucL)/4) < ((*ucR)/4)) ? 1 : -1;
976 return 0;
977 }
VF_2_ucNF(PolyPointList * P,VertexNumList * V,EqList * E,int * NV,int * nUC,unsigned char * UC)978 void VF_2_ucNF(PolyPointList *P, VertexNumList *V, EqList *E, /* IN */
979 int *NV, int *nUC, unsigned char *UC) /* OUT */
980 { Long V_NF[POLY_Dmax][VERT_Nmax], F_NF[POLY_Dmax][VERT_Nmax];
981 Long VM[POLY_Dmax][VERT_Nmax], VPM[VERT_Nmax][VERT_Nmax];
982 unsigned char auxUC[POLY_Dmax*VERT_Nmax];
983 int MS=0, vb,vo,vnuc,vbmin, fb,fo,fnuc,fbmin, vone=0,fone=0,MSone;
984 #ifndef SORT_NUC_FIRST
985 if(V->nv<E->ne) MS=1;
986 if(V->nv>E->ne) MS=2;
987 #endif /* else: priority to sorting by nuc */
988 assert(Init_rVM_VPM(P,V,E,&P->n,&V->nv,&E->ne,VM,VPM));/* make ref VPM */
989 if(MS <2)
990 { Eval_Poly_NF(&P->n,&V->nv,&E->ne,VM,VPM,V_NF,0); /* V_NF */
991 vone=BminOff(V_NF,&P->n,&V->nv,&vo,&vbmin);
992 #ifdef USE_UNIT_ENCODE
993 if(vone) UNIT_Bmin2BaseUCn(&P->n,&V->nv,&vbmin,&vb,&vnuc); else
994 #endif
995 Bmin2BaseUCn(&P->n,&V->nv,&vbmin,&vb,&vnuc);
996 }
997 if(MS!=1)
998 { int i, j, mi=(E->ne>V->nv) ? V->nv : E->ne; /* VM=E */
999 for(i=0;i<P->n;i++)for(j=0;j<E->ne;j++) VM[i][j]=E->e[j].a[i];
1000 for(i=0;i<mi;i++)for(j=0;j<i;j++) /* transpose VPM */
1001 { Long a=VPM[i][j]; VPM[i][j]=VPM[j][i]; VPM[j][i]=a;
1002 }
1003 if(E->ne>V->nv)
1004 for(i=mi;i<E->ne;i++)for(j=0;j<V->nv;j++)VPM[j][i]=VPM[i][j];
1005 else for(i=mi;i<V->nv;i++)for(j=0;j<E->ne;j++)VPM[i][j]=VPM[j][i];
1006 Eval_Poly_NF(&P->n,&E->ne,&V->nv,VM,VPM,F_NF,0); /* compute F_NF */
1007 fone=BminOff(F_NF,&P->n,&E->ne,&fo,&fbmin);
1008 #ifdef USE_UNIT_ENCODE
1009 if(fone) UNIT_Bmin2BaseUCn(&P->n,&E->ne,&fbmin,&fb,&fnuc); else
1010 #endif
1011 Bmin2BaseUCn(&P->n,&E->ne,&fbmin,&fb,&fnuc);
1012 }
1013 if(MS==0) { if(vnuc<fnuc) MS=1; if(fnuc<vnuc) MS=2; } /* sort nuc */
1014 if(MS==0) { if(V->nv<E->ne) MS=1; if(E->ne<V->nv) MS=2; } /* sort vert */
1015 #ifdef USE_UNIT_ENCODE
1016 if(MS==0) { if(vone>fone) MS=1; if(vone<fone) MS=2; } /* sort UNIT */
1017 #endif
1018
1019 *nUC = (MS==2) ? fnuc : vnuc; *NV = (MS==2) ? E->ne : V->nv; MSone=MS;
1020
1021 if(AuxNFLptr!=NULL) /* base:byte statistics */
1022 { /* if(AuxNFLptr->Xmin>-vo) AuxNFLptr->Xmin=-vo;
1023 if(AuxNFLptr->Xmax<vbmin-vo-1) AuxNFLptr->Xmax=vbmin-vo-1; */
1024 if(MS==2) {if(AuxNFLptr->Xdif<fbmin) AuxNFLptr->Xdif=fbmin;}
1025 else {if(AuxNFLptr->Xdif<vbmin) AuxNFLptr->Xdif=vbmin;}
1026 if(AuxNFLptr->Xnuc<*nUC) AuxNFLptr->Xnuc=*nUC;
1027 if(((MS==2)&&(fbmin>BASE_MAX))||((MS<2)&&(vbmin>BASE_MAX)))
1028 { printf("WARNING MS=%d v=%d vb=%d vnuc=%d f=%d fb=%d fnuc=%d\n",
1029 MS,V->nv,vbmin,vnuc,E->ne,fbmin,fnuc);
1030 VPrint(&P->n,&V->nv,V_NF);puts("============= V_NF");
1031 VPrint(&P->n,&E->ne,F_NF);puts("============= F_NF");
1032 puts("BASE_MAX exceeded"); exit(0);
1033 }
1034 }
1035
1036 #ifdef USE_UNIT_ENCODE
1037 if(MS <2) MSone += 4*vone; else MSone += 4*fone;
1038 #endif
1039 if(MS <2) AuxVnf2ucNF(V_NF, &P->n, &V->nv, &vo, &vb, &vnuc, &MSone, UC);
1040 else AuxVnf2ucNF(F_NF, &P->n, &E->ne, &fo, &fb, &fnuc, &MSone, UC);
1041
1042 if(MS==0) /* nuc and nv agree => compare UC[] */
1043 { int RmL;
1044 AuxVnf2ucNF(F_NF, &P->n, &E->ne, &fo, &fb, &fnuc, &MSone, auxUC);
1045 RmL=RIGHTminusLEFT(UC,auxUC,nUC);
1046 if(RmL>0) {MS=1; (*UC)++;} /* o.k. since *UC=4*(...) */
1047 if(RmL<0) {int i; for(i=0;i<*nUC;i++)UC[i]=auxUC[i]; (*UC) += (MS=2);}
1048 }
1049 /* 6+8*(2+17(1+17(0+17(1+17(1+17(1+17(2))))))) {23 189 20 198}
1050 */
1051 #ifdef TEST_UCnf /* One=K3[3/F]: 2 1 0 1 1 1 2 */
1052 { Long tNF[POLY_Dmax][VERT_Nmax]; int tMS, i, j, *d=&P->n, err=0;
1053 static int pn, Err;
1054 pn++; UCnf2vNF(d, NV, nUC, UC, tNF, &tMS);
1055 if(MS<2) for(i=0;i<*d;i++) for(j=i;j<*NV;j++)
1056 err+=(tNF[i][j]!=V_NF[i][j]);
1057 if(err) err*=1000;
1058 if(MS-1) for(i=0;i<*d;i++) for(j=i;j<*NV;j++)
1059 err+=(tNF[i][j]!=F_NF[i][j]);
1060 if(err) Err=1;
1061 if(Err)
1062 { int base;
1063 #ifdef USE_UNIT_ENCODE
1064 if(MSone>3) UNIT_NUCtoBase(&P->n,NV,nUC,&base); else
1065 #endif
1066 NUCtoBase(&P->n,NV,nUC,&base);
1067 printf("\n#%d MS=%d 100*e(V)-e(F) = %d\n",pn,MS,err);
1068 for(i=0;i<*nUC;i++)printf("%d ",UC[i]);
1069 printf("= UC[%d] base=%d vo=%d fo=%d\n",*nUC,base,vo,fo);
1070 VPrint(&P->n,&V->nv,V_NF);puts("============= V_NF");
1071 VPrint(&P->n,&E->ne,F_NF);puts("============= F_NF");
1072 VPrint(&P->n,NV,tNF);puts("============= t_NF"); exit(0);
1073 }
1074 }
1075 #endif
1076 assert((*UC % 4)==MS);
1077 }
Test_ucNF(int * d,int * v,int * nuc,unsigned char * uc,PolyPointList * _P)1078 void Test_ucNF(int *d, int *v, int *nuc, unsigned char *uc, PolyPointList *_P)
1079 { Long tNF[POLY_Dmax][VERT_Nmax]; int i, j, tMS, NV, NUC;
1080 unsigned char UC[POLY_Dmax*VERT_Nmax];
1081 VertexNumList V; EqList E;
1082 UCnf2vNF(d, v, nuc, uc, tNF, &tMS); _P->n=*d; _P->np=*v; tMS %= 4;
1083 for(i=0;i<*v;i++)for(j=0;j<*d;j++)_P->x[i][j]=tNF[j][i];
1084 assert(Ref_Check(_P,&V,&E)); VF_2_ucNF(_P,&V,&E,&NV,&NUC,UC);
1085 assert(*v==NV); assert(*nuc==NUC);
1086 assert(RIGHTminusLEFT(uc,UC,nuc)==0);
1087 assert((tMS==0)==((*uc%4)==0));
1088 }
1089
AuxGet_vn_uc(FILE * F,int * v,int * nu,unsigned char * uc)1090 void AuxGet_vn_uc(FILE *F,int *v, int *nu, unsigned char *uc)
1091 { int i; *v=fgetc(F); *nu=fgetc(F); for(i=0;i<*nu;i++) uc[i]=fgetc(F);
1092 }
AuxGet_uc(FILE * F,int * nu,unsigned char * uc)1093 void AuxGet_uc(FILE *F,int *nu, unsigned char *uc)
1094 { int i; for(i=0;i<*nu;i++) uc[i]=fgetc(F);
1095 }
1096
AuxPut_hNF(FILE * F,int * v,int * nu,unsigned char * Huc,FInfoList * Io,int * slNF,int * slSM,int * slNM,int * slNB,unsigned char * ucSL,int * SLp)1097 void AuxPut_hNF(FILE *F,int *v,int *nu,unsigned char *Huc,FInfoList *Io,
1098 int *slNF,int *slSM,int *slNM,int *slNB,unsigned char *ucSL,int *SLp)
1099 { int i, Hms=(*Huc % 4); static int pos; unsigned char *Suc=NULL;
1100 for(i=0;i<*nu;i++) fputc(Huc[i],F);
1101 if(Hms) {if(Hms<3) Io->nNM++;} else Io->nSM++;
1102 while(pos<*slNF)
1103 { int HmSL; Suc=&ucSL[SLp[pos]]; HmSL=*v-*Suc;
1104 if(HmSL>0) {pos++; continue;} else if(HmSL<0) return;
1105 else HmSL=*nu-Suc[1];
1106 if(HmSL>0) {pos++; continue;} else if(HmSL<0) return;
1107 else HmSL=RIGHTminusLEFT(&Suc[2],Huc,nu);
1108 if(HmSL>0) {pos++; continue;} else if(HmSL<0) return; else break;
1109 }
1110 if(pos<*slNF)
1111 { int Sms=Suc[2]%4; pos++; switch(10*Hms+Sms) {
1112 case 00: case 11: case 22: case 31: case 32: case 33: break;
1113 case 12: case 21: return;
1114 case 13: Suc[2]-=1; (*slNM)++; return; /* remove 1st */
1115 case 23: Suc[2]-=2; (*slNM)++; return; /* remove 2nd */
1116 default: puts("inconsistent MS flags in AuxPut_hNF"); exit(0);}
1117
1118 if(Sms) {if(Sms<3) (*slNM)--;} else (*slSM)--; /* remove SL-entry */
1119 for(i=pos--;i<*slNF;i++) SLp[i-1]=SLp[i];
1120 (*slNF)--; (*slNB)-= *nu+2;
1121 }
1122 }
1123 /* Alloc & read SL; go thru pi / pa; write SL + statistics
1124 */
Add_Polya_2_Polyi(char * polyi,char * polya,char * polyo)1125 void Add_Polya_2_Polyi(char *polyi,char *polya,char *polyo)
1126 { FILE *FI=fopen(polyi,"rb"), *FA=fopen(polya,"rb"), *FO;
1127 FInfoList FIi, FIa, FIo; Along Ipos, Apos, HIpos, HApos;
1128 unsigned char ucI[NUC_Nmax], ucA[NUC_Nmax], *ucSL, *uc; int SLp[SL_Nmax];
1129 int IslNF, IslSM, IslNM, vI, nuI, i, d; unsigned Ili; Along IslNB;
1130 int AslNF, AslSM, AslNM, vA=0, nuA=0, a, j; unsigned Ali; Along AslNB;
1131 int slNF=0,slSM=0,slNM=0, slNB=0, slNP=0, Oli=0, v, nu, AmI=0, ms, tnb=0;
1132 Init_FInfoList(&FIi); Init_FInfoList(&FIa); if(!*polyi||!*polyo) {
1133 puts("With -pa you require -pi and -po or -di and -do");exit(0);}
1134 if(NULL==FI) {printf("Cannot open %s\n",polyi);exit(0);}
1135 if(NULL==FA) {printf("Cannot open %s\n",polya);exit(0);}
1136 if(NULL==(FO=fopen(polyo,"wb"))){printf("Cannot open %s",polyo);exit(0);}
1137 ucSL = (unsigned char *) malloc( SL_Nmax * CperR_MAX * sizeof(char) );
1138 assert(ucSL!=NULL); assert(!fgetc(FI)); assert(!fgetc(FA));
1139 Read_Bin_Info(FI,&i,&Ili,&IslNF,&IslSM,&IslNM,&IslNB,&FIi);
1140 Read_Bin_Info(FA,&d,&Ali,&AslNF,&AslSM,&AslNM,&AslNB,&FIa); assert(d==i);
1141 HIpos=FTELL(FI); FSEEK(FI,0,SEEK_END); Ipos=FTELL(FI);
1142 HApos=FTELL(FA); FSEEK(FA,0,SEEK_END); Apos=FTELL(FA);
1143 printf("Data on %s: %lld+%dsl %lldb (%dd)\n", polyi,
1144 (FIi.nNF-FIi.nSM)+(FIi.nNF-FIi.nNM),
1145 /* Islp= */ 2*IslNF-IslNM-IslSM,Ipos,i);
1146 printf("Data on %s: %lld+%dsl %lldb (%dd)\n", polya,
1147 (FIa.nNF-FIa.nSM)+(FIa.nNF-FIa.nNM),
1148 /* Aslp= */ 2*AslNF-AslNM-AslSM,Apos,d);
1149 /* printf("HIpos=%lld NB=%lld sl=%lld pos=%lld\n",HIpos,FIi.NB,IslNB,Ipos); */
1150 assert(HIpos+FIi.NB+IslNB==Ipos); FSEEK(FI,-IslNB,SEEK_CUR);
1151 /* printf("HApos=%lld NB=%lld sl=%lld pos=%lld\n",HApos,FIa.NB,AslNB,Apos); */
1152 assert(HApos+FIa.NB+AslNB==Apos); FSEEK(FA,-AslNB,SEEK_CUR);
1153 a=0; if(a<AslNF) AuxGet_vn_uc(FA,&vA,&nuA,ucA);
1154 for(i=0;i<IslNF;i++)
1155 { AuxGet_vn_uc(FI,&vI,&nuI,ucI); uc = &ucSL[ 2 + (SLp[slNF++]=slNB) ];
1156 while(a<AslNF)
1157 { if(!(AmI=vA-vI)) if(!(AmI=nuA-nuI))
1158 AmI=RIGHTminusLEFT(ucI,ucA,&nuI);
1159 if(AmI<0) /* put A */
1160 { uc[-2]=vA; uc[-1]=nuA; for(j=0;j<nuA;j++)uc[j]=ucA[j];
1161 slNB+=2+nuA; if((ms=(*uc%4))) {if(ms<3) slNM++;} else slSM++;
1162 if((++a)<AslNF) AuxGet_vn_uc(FA,&vA,&nuA,ucA);
1163 uc = &ucSL[ 2 + (SLp[slNF++]=slNB) ];
1164 }
1165 else break;
1166 }
1167 if((a<AslNF)&&(AmI==0)) /* put I==A */
1168 { uc[-2]=vI; uc[-1]=nuI; for(j=0;j<nuI;j++) uc[j]=ucI[j];
1169 if((*ucI%4)!=(*ucA%4)) *uc = 3+4*(*uc/4);
1170 if((ms=(*uc%4))) {if(ms<3) slNM++;} else slSM++;
1171 if((++a)<AslNF) AuxGet_vn_uc(FA,&vA,&nuA,ucA);
1172 tnb+=2+nuI;
1173 }
1174 else /* put I */
1175 { uc[-2]=vI; uc[-1]=nuI; for(j=0;j<nuI;j++) uc[j]=ucI[j];
1176 if((ms=(*uc%4))) {if(ms<3) slNM++;} else slSM++;
1177 } slNB+=2+nuI;
1178 }
1179 while(a<AslNF) /* put A */
1180 { uc = &ucSL[ 2 + (SLp[slNF++]=slNB) ]; slNB+=2+nuA;
1181 uc[-2]=vA; uc[-1]=nuA; for(j=0;j<nuA;j++)uc[j]=ucA[j];
1182 if((ms=(*uc%4))) {if(ms<3) slNM++;} else slSM++;
1183 if((++a)<AslNF) AuxGet_vn_uc(FA,&vA,&nuA,ucA);
1184 } assert(tnb+slNB==IslNB+AslNB); /* SL done */
1185
1186 printf("SL: %dnf %dsm %dnm %db -> ",slNF,slSM,slNM,slNB);
1187
1188 FSEEK(FI,HIpos,SEEK_SET); FSEEK(FA,HApos,SEEK_SET); Init_FInfoList(&FIo);
1189 FIo.nVmax=max(FIi.nVmax,FIa.nVmax);
1190 FIo.NUCmax=max(FIi.NUCmax,FIa.NUCmax); tnb=0;
1191 for(v=d+1;v<=FIo.nVmax;v++) for(nu=1;nu<=FIo.NUCmax;nu++)
1192 if((FIo.NFnum[v][nu]=FIi.NFnum[v][nu]+FIa.NFnum[v][nu]))
1193 { FIo.nNUC[v]++; Oli++; tnb+=FIo.NFnum[v][nu];
1194 } FIo.nV=0; for(v=d+1;v<=FIo.nVmax;v++) if(FIo.nNUC[v]) FIo.nV++;
1195 fputc(0,FO); fputc(d,FO); fputc(FIo.nV,FO); fputc(FIo.nVmax,FO);
1196 fputc(FIo.NUCmax,FO); fputUI(Oli,FO);
1197 j=InfoSize(0,Oli,&FIo)-5-sizeof(int); for(i=0;i<j;i++) fputc(0,FO);
1198
1199 for(v=d+1;v<=FIo.nVmax;v++) if(FIo.nNUC[v])
1200 for(nu=1;nu<=FIo.NUCmax;nu++) if(FIo.NFnum[v][nu])
1201 { int I_NF=FIi.NFnum[v][nu], A_NF=FIa.NFnum[v][nu], O_NF=0;
1202 int neq=0, peq=0, pi=0, pa=0, po=0;
1203 a=0; if(a<A_NF) {AuxGet_uc(FA,&nu,ucA); pa += 1+(((*ucA)%4)==3);}
1204 for(i=0;i<I_NF;i++)
1205 { AuxGet_uc(FI,&nu,ucI); pi += 1+(((*ucI)%4)==3);
1206 while(a<A_NF)
1207 { AmI=RIGHTminusLEFT(ucI,ucA,&nu);
1208 if(AmI<0) /* put A */
1209 { po += 1+(((*ucA)%4)==3);
1210 AuxPut_hNF(FO,&v,&nu,ucA,&FIo,
1211 &slNF,&slSM,&slNM,&slNB,ucSL,SLp); O_NF++; a++;
1212 if(a<A_NF){AuxGet_uc(FA,&nu,ucA);pa += 1+(((*ucA)%4)==3);}
1213 }
1214 else break;
1215 }
1216 if((a<A_NF)&&(AmI==0)) /* put I==A */
1217 { int mm=10*(*ucI%4) + (*ucA%4); switch(mm) {
1218 case 00: peq++; break; case 11: peq++; break;
1219 case 12: *ucI+=2; break; case 13: *ucI+=2;peq++; break;
1220 case 21: *ucI+=1; break; case 22: peq++; break;
1221 case 23: *ucI+=1;peq++; break; case 31: peq++; break;
1222 case 32: peq++; break; case 33: peq+=2; break;
1223 default: puts("inconsistens mirror flags");exit(0);}
1224 AuxPut_hNF(FO,&v,&nu,ucI,&FIo,&slNF,&slSM,&slNM,&slNB,ucSL,
1225 SLp); po += 1+(((*ucI)%4)==3); neq++; a++;
1226 if(a<A_NF) {AuxGet_uc(FA,&nu,ucA); pa += 1+(((*ucA)%4)==3);}
1227 }
1228 else
1229 { AuxPut_hNF(FO,&v,&nu,ucI,&FIo,&slNF,
1230 &slSM,&slNM,&slNB,ucSL,SLp); po += 1+(((*ucI)%4)==3);
1231 } O_NF++;
1232 }
1233 while(a<A_NF) /* put A */
1234 { O_NF++; po += 1+(((*ucA)%4)==3);
1235 AuxPut_hNF(FO,&v,&nu,ucA,&FIo,&slNF,&slSM,&slNM,&slNB,ucSL,SLp);
1236 ++a; if(a<A_NF) {AuxGet_uc(FA,&nu,ucA); pa += 1+(((*ucA)%4)==3);}
1237 } assert(pi+pa==peq+po); /* checksum(v,nu) */
1238 FIo.NFnum[v][nu]=O_NF; FIo.nNF+=O_NF; FIo.NB+=O_NF*nu;
1239 assert(O_NF+neq==I_NF+A_NF);
1240 /*
1241 {static int list;printf("#%d v=%d nu=%d Inf=%d Anf=%d ",++list,v,nu,
1242 I_NF,A_NF);printf("Onf=%d pi=%d pa=%d po=%d\n",O_NF,pi,pa,po);
1243 */
1244 } tnb=0;
1245 for(i=0;i<slNF;i++) /* write SL */
1246 { uc=&ucSL[SLp[i]+2]; v=uc[-2]; nu=uc[-1]; tnb+=nu+2;
1247 /* printf("#%d: SLp=%d v=%d nu=%d\n",i,SLp[i],v,nu); */
1248 assert(uc[-2]<VERT_Nmax); fputc(uc[-2],FO); slNP+=1+(((*uc)%4)==3);
1249 fputc(nu,FO); for(j=0;j<nu;j++) fputc(uc[j],FO);
1250 } assert(tnb==slNB);
1251 assert(slNP==2*slNF-slNM-slSM);
1252 printf("\nd=%d v%d v<=%d n<=%d vn%d %lld %d %lld %lld %d %d %d %d\n",
1253 d, FIo.nV, FIo.nVmax, FIo.NUCmax,Oli,
1254 FIo.nNF,FIo.nSM,FIo.nNM,FIo.NB, slNF, slSM, slNM, slNB);
1255 /* for(v=d+1;v<=FIo.nVmax;v++) if(FIo.nNUC[v])
1256 { i=0; printf("%d %d\n",v,FIo.nNUC[v]);
1257 for(nu=1;nu<=FIo.NUCmax;nu++) if(FIo.NFnum[v][nu])
1258 { printf("%d %d%s",nu,FIo.NFnum[v][nu],
1259 (++i<FIo.nNUC[v]) ? " " : "\n");
1260 }
1261 } */
1262
1263 FSEEK(FO,5+sizeof(int),SEEK_SET);
1264 fputUI(FIo.nNF,FO); /* write info */
1265 fputUI(FIo.nSM,FO); fputUI(FIo.nNM,FO); fputUI(FIo.NB,FO);
1266 fputUI(slNF,FO); fputUI(slSM,FO); fputUI(slNM,FO); fputUI(slNB,FO);
1267 for(v=d+1;v<=FIo.nVmax;v++) if(FIo.nNUC[v])
1268 { i=0; for(nu=1;nu<=FIo.NUCmax;nu++) if(FIo.NFnum[v][nu]) i++;
1269 fputc(v,FO); fputc(i,FO); /* v #nuc(v) */
1270 for(nu=1;nu<=FIo.NUCmax;nu++) if(FIo.NFnum[v][nu])
1271 { fputc(nu,FO); fputUI(FIo.NFnum[v][nu],FO); /* nuc #NF */
1272 }
1273 }
1274
1275 printf("Writing %s: %lld+%dsl %lldm+%ds %lldb",polyo,2*FIo.nNF-FIo.nNM
1276 -FIo.nSM,slNP,/*tnb=*/ FIo.nNF-FIo.nNM-FIo.nSM,FIo.nSM,FIo.NB+slNB);
1277 /* if(tnb>99)
1278 { long long tnp=(2*FIo.nNF-FIo.nNM-FIo.nSM)/1000; tnp*=tnp;
1279 tnp/=(2*tnb); printf(" [p^2/2m=%ldM]",tnp);
1280 }*/
1281 Print_Expect(&FIo); puts("");
1282 assert(ferror(FI)==0); fclose(FI);
1283 assert(ferror(FA)==0); fclose(FA);
1284 assert(ferror(FO)==0); fclose(FO);
1285 }
1286
1287
1288 /* ================== AffineNF modifications ================== */
1289
UCnf_2_ANF(int * d,int * v,int * nuc,unsigned char * uc,Long NF[POLY_Dmax][VERT_Nmax],int * MS)1290 void UCnf_2_ANF(int *d, int *v, int *nuc, unsigned char *uc, /* IN */
1291 Long NF[POLY_Dmax][VERT_Nmax], int *MS) /* OUT */
1292 { int i,off;UCnf2vNF(d,v,nuc,uc,NF,MS);off=NF[0][*v-1]; /* offset not 1 */
1293 for(i=1;i<*d;i++)assert(NF[i][*v-1]==off); /* but last column ! */
1294
1295 for(i=0;i<*d;i++) {int c; for(c=0;c<i;c++) assert(NF[i][c]==0);
1296 for(c=i;c<*v;c++) NF[i][c]-=off;} assert(*MS%4==1);
1297 }
ANF_2_ucNF(PolyPointList * P,VertexNumList * V,EqList * E,int * NV,int * nUC,unsigned char * UC)1298 void ANF_2_ucNF(PolyPointList *P, VertexNumList *V, EqList *E, /* IN */
1299 int *NV, int *nUC, unsigned char *UC) /* OUT */
1300 { int vb,vo,vnuc,vbmin, vone=0, MSone; Long V_NF[POLY_Dmax][VERT_Nmax];
1301 Make_ANF(P,V,E,V_NF); /* Affine V_NF */
1302 vone=BminOff(V_NF,&P->n,&V->nv,&vo,&vbmin);
1303 #ifdef USE_UNIT_ENCODE
1304 if(vone) UNIT_Bmin2BaseUCn(&P->n,&V->nv,&vbmin,&vb,&vnuc); else
1305 #endif
1306 Bmin2BaseUCn(&P->n,&V->nv,&vbmin,&vb,&vnuc);
1307
1308 *nUC = vnuc; *NV = V->nv; MSone=1;
1309
1310 if(AuxNFLptr!=NULL) /* base:byte statistics */
1311 { if(AuxNFLptr->Xdif<vbmin) AuxNFLptr->Xdif=vbmin;
1312 if(AuxNFLptr->Xnuc<*nUC) AuxNFLptr->Xnuc=*nUC;
1313 if(vbmin>BASE_MAX)
1314 { printf("WARNING MS=%d v=%d vb=%d vnuc=%d\n",
1315 1,V->nv,vbmin,vnuc);
1316 VPrint(&P->n,&V->nv,V_NF);puts("============= V_NF");
1317 puts("BASE_MAX exceeded"); exit(0);
1318 }
1319 }
1320 #ifdef USE_UNIT_ENCODE
1321 MSone += 4*vone;
1322 #endif
1323 AuxVnf2ucNF(V_NF, &P->n, &V->nv, &vo, &vb, &vnuc, &MSone, UC);
1324
1325 #ifdef TEST_UCnf /* One=K3[3/F]: 2 1 0 1 1 1 2 */
1326 { Long tNF[POLY_Dmax][VERT_Nmax]; int tMS, i, j, *d=&P->n, err=0;
1327 static int pn;
1328 pn++; UCnf_2_ANF(d, NV, nUC, UC, tNF, &tMS);
1329 for(i=0;i<*d;i++) for(j=0;j<*NV;j++)
1330 err+=(tNF[i][j]!=V_NF[i][j]);
1331 if(err) err*=1000;
1332 if(err)
1333 { int base;
1334 #ifdef USE_UNIT_ENCODE
1335 if(MSone>3) UNIT_NUCtoBase(&P->n,NV,nUC,&base); else
1336 #endif
1337 NUCtoBase(&P->n,NV,nUC,&base);
1338 printf("\n#%d MS=%d 100*e(V)-e(F) = %d\n",pn,1,err);
1339 for(i=0;i<*nUC;i++)printf("%d ",UC[i]);
1340 printf("= UC[%d] base=%d vo=%d\n",*nUC,base,vo);
1341 VPrint(&P->n,&V->nv,V_NF);puts("============= V_NF");
1342 VPrint(&P->n,NV,tNF);puts("============= t_NF"); exit(0);
1343 }
1344 }
1345 #endif
1346 assert((*UC % 4)==1);
1347 }
1348
Add_ANF_to_List(PolyPointList * _P,VertexNumList * _V,EqList * _E,NF_List * _L)1349 int Add_ANF_to_List(PolyPointList *_P, VertexNumList *_V, EqList *_E,
1350 NF_List *_L)
1351 { unsigned char UC[NB_MAX]; int nUC, NV, NewNF;
1352 #ifdef INCREMENTAL_TIME
1353 long long cpuT=clock(), incT=cpuT - _L->CLOCK;
1354 if(incT>0) _L->IP_Time+=incT;
1355 _L->CLOCK=cpuT;
1356 #endif
1357 if((_E->ne)>_L->F) _L->F=_E->ne; /* check vertex numbers */
1358 if((_V->nv)>_L->V) _L->V=_V->nv;
1359 if((_E->ne > VERT_Nmax)||(_V->nv > VERT_Nmax))
1360 { printf("Increase VERT_Nmax: f=%d v=%d\n",_E->ne,_V->nv);exit(0);
1361 }
1362 if(_P->np>_L->Nmax) _L->Nmax=_P->np;
1363 if(_P->np<_L->Nmin) _L->Nmin=_P->np;
1364 _L->nNF++; if(_L->SL)_L->nSLNF++;
1365
1366 ANF_2_ucNF(_P, _V, _E, &NV, &nUC, UC);
1367
1368 NewNF=ucNF_Sort_Add(&NV, &nUC, UC, _L); /* 1::new::cont. */
1369
1370 #ifdef INCREMENTAL_TIME
1371 cpuT=clock(), incT=cpuT - _L->CLOCK;
1372 if(incT>0) _L->NF_Time+=incT;
1373 _L->CLOCK=cpuT;
1374 #endif
1375 return NewNF;
1376 }
1377
Gen_Ascii_to_Binary(CWS * W,PolyPointList * P,char * dbin,char * polyi,char * polyo)1378 void Gen_Ascii_to_Binary(CWS *W, PolyPointList *P,
1379 char *dbin, char *polyi, char *polyo)
1380 { NF_List *_NFL=(NF_List *) malloc(sizeof(NF_List));
1381 VertexNumList V; EqList F; assert(_NFL!=NULL);
1382 if(!(*polyo)) {
1383 puts("You have to specify an output file via -po in -a-mode!\n");
1384 printf("For more help use option '-h'\n");
1385 exit(0);}
1386 _NFL->of=0; _NFL->rf=0;
1387 _NFL->iname=polyi; _NFL->oname=polyo; _NFL->dbname=dbin;
1388 Init_NF_List(_NFL);
1389 _NFL->SL=0;
1390
1391 while(Read_CWS_PP(W,P)) {
1392 _NFL->hc = 0;
1393 _NFL->V= _NFL->F= _NFL->VN= _NFL->FN = _NFL->Xnuc = _NFL->Xdif = 0;
1394 _NFL->Nmin=P->np; _NFL->Nmax=0;
1395 if(_NFL->d==0) _NFL->d=P->n;
1396 else if(_NFL->d-P->n) {puts("different dim!"); exit(0);}
1397
1398 Find_Equations(P,&V,&F);
1399 if (Add_ANF_to_List(P,&V,&F,_NFL)) if (outFILE!=stdout){
1400 int i, j;
1401 for(i=0;i<W->nw;i++) {
1402 fprintf(outFILE,"%ld ",W->d[i]);
1403 for(j=0;j<W->N;j++) fprintf(outFILE,"%ld ",W->W[i][j]);
1404 if(i+1<W->nw) fprintf(outFILE," "); else fprintf(outFILE,"\n"); }
1405 fflush(0);} }
1406 Write_List_2_File(polyo,_NFL);
1407 free(_NFL);
1408 }
1409