1 #include "Global.h"
2 #include "Rat.h"
3 #include "LG.h"
4
5 #define SHOW_b01_TWIST (0)
6 #define Tout (0)
7 #define DET_WARN_ONLY (0) /* continue if group has det!=1 */
8
9 #define COEFF_Nmax (d*D+2*N)
10
11 #define ABBREV_POLY_PRINT (4) /* !=0 => #(leading/trailing terms */
12
13 #define NO_COORD_IMPROVEMENT /* switch off weight permutation */
14 #undef W_PERM_CODE
Is_Gen_CY(int index,PolyPointList * P)15 int Is_Gen_CY(int index, PolyPointList *P)
16 { int i; Long *IP=P->x[P->np]; VertexNumList V; /* IP = IP(index * P) */
17 EqList *E = (EqList *) malloc(sizeof(EqList)); assert(E!=NULL);
18 Find_Equations(P,&V,E); for(i=0;i<E->ne;i++) E->e[i].c *= index;
19 for(i=0;i<E->ne;i++) if(1!=Eval_Eq_on_V(&E->e[i],IP,P->n)) break;
20 free(E); return i==E->ne;
21 }
WZerror(char * c)22 void WZerror(char *c)
23 { printf("Format error %s in Read_WZeight\n",c);exit(0);
24 }
auxString2SInt(char * c,int * n)25 int auxString2SInt(char *c,int *n)
26 { int j=0,neg=0; *n=0; if(c[0]=='-')if(('0'<c[1])&&(c[1]<='9')) neg=j=1;
27 while(('0'<=c[j])&&(c[j]<='9')) *n = 10*(*n)+c[j++]-'0';
28 if(neg)*n *=-1;
29 if(j) while(c[j]==' ') j++;
30 return j;
31 }
Read_WZ_PP(Weight * WZ)32 int Read_WZ_PP(Weight *WZ) /* read "d w_i" [ or "w_i d" if last=max ] */
33 { int i,j,k,a,n,d,shift=1, I[W_Nmax+2], *nz=&WZ->M;
34 int FilterFlag=(inFILE==NULL); char C, c[999],b=' ';
35 Long BM[W_Nmax][W_Nmax], *B[W_Nmax], Wa[POLY_Dmax], Za[POLY_Dmax],
36 F[W_Nmax], G[POLY_Dmax][POLY_Dmax], GI[POLY_Dmax][POLY_Dmax],X;
37 if(FilterFlag) inFILE=stdin; else
38 if(inFILE==stdin) printf("type degree and weights [d w1 w2 ...]: ");
39 C=fgetc(inFILE); if( !IsDigit(C) ) return 0; ungetc(C,inFILE); *nz=0;
40 fscanf(inFILE,"%d",I);
41 for(i=1; i<W_Nmax+2; i++)
42 { while(' '==(C=fgetc(inFILE))); ungetc(C,inFILE);
43 if(IsDigit(C)) fscanf(inFILE,"%d",&I[i]); else break;
44 } WZ->N=i-1; if(WZ->N>W_Nmax) { puts("Increase POLY_Dmax"); exit(0); }
45 for(i=0;i<=WZ->N;i++) assert(I[i]>0);
46 if(I[WZ->N]>I[0]) { WZ->d=I[WZ->N]; shift=0; } else WZ->d=I[0];
47 for(i=0;i<WZ->N;i++) {WZ->w[i]=I[i+shift]; assert(WZ->w[i]<WZ->d);}
48 WZ->r=0; for(i=0;i<WZ->N;i++) WZ->r+=WZ->w[i];
49 if(WZ->r%WZ->d)WZ->r=0; else WZ->r/=WZ->d;
50 for(n=0;n<999;n++) /* read /Z*: * * * */
51 { c[n]=fgetc(inFILE);
52 if(feof(inFILE)) {if(FilterFlag) inFILE=NULL; return 0;}
53 if(c[n]=='\n') break;
54 } if(n==999) {puts("Out of space in Read_WZeight");exit(0);} i=0;
55 while(c[i]==b)i++;
56 if((c[i]=='=')&&(c[i+1]=='d'))i+=2;
57 while(c[i]==b)i++;
58 while(i<n)
59 { int j, k, s; if((c[i]!='/')||(c[i+1]!='Z')) break;
60 i+=2; assert(*nz<POLY_Dmax); WZ->m[*nz]=0;
61 while((i<n)&&IsDigit(c[i])) WZ->m[*nz]=10*WZ->m[*nz]+c[i++]-'0';
62 j=0;
63 if(WZ->m[*nz]<=0) WZerror("Order not positive");
64 if(c[i+j]!=':') WZerror(":"); else {i+=j+1; while(c[i]==b)i++;}
65 for(k=0;k<WZ->N;k++)
66 { if((j=auxString2SInt(&c[i],&WZ->z[*nz][k]))) {if(j) i+=j;
67 else WZerror("missing");}
68 } s=0;
69 for(k=0;k<WZ->N;k++) s+=WZ->z[*nz][k];
70 /* if(s % WZ->m[*nz]) WZerror("det!=1") */;
71 (WZ->M)++;
72 } /* Eq_i::{c=Ai a[]=Bi[]} 0<=A+B*x r=sum(w)/d rI=IP(r*P) n=(r,rI) */
73
74 a=WZ->N; for(i=0;i<a;i++)B[i]=BM[i]; assert( 1 == W_to_GLZ(WZ->w,&a,B) );
75 if(1==WZ->r) for(i=0;i<a;i++) WZ->A[i]=F[i]=1; /* Z.(X-F) \cong 0 */
76 else {for(i=0;i<a;i++) { WZ->A[i]=WZ->d*B[0][i]; F[i]=0;}
77 for(j=0;j<WZ->M;j++) {Long s=0; for(n=0;n<a;n++) s+=WZ->z[j][n];
78 if(s%WZ->m[j]){puts("det=1 required if Gorenstein index r>1:");
79 #if DET_WARN_ONLY
80 WZ->M=0; }
81 #else
82 Write_Weight(WZ); exit(0);}
83 #endif
84 }} d=a-1; n=0; shift=1;
85 for(i=0;i<a;i++)for(j=0;j<d;j++)WZ->B[i][j]=B[j+1][i]; /* LATTICE BASIS */
86 for(i=a-1;i>0;i--) { X = 1 - WZ->r * WZ->A[i];
87 for(j=i;j<d;j++) X -= WZ->rI[j]*WZ->B[i][j];
88 assert(0==(X%WZ->B[i][i-1])); WZ->rI[i-1]=X/WZ->B[i][i-1];}
89 for(i=0;i<*nz;i++) /* divide by symmetries: */
90 { Long g=WZ->m[i]; X=0; for(j=0;j<d;j++){Za[j]=0;
91 for(k=0;k<a;k++)Za[j]+=WZ->z[i][k]*WZ->B[k][j];
92 g=NNgcd(Za[j],g);}
93 if(Tout){printf("WZ->m[%d]=%d g=%ld ",i,WZ->m[i],g);if(g==WZ->m[i])puts("");}
94 if(g==WZ->m[i]) continue;
95 if(g>1) for(k=0;k<a;k++) Za[j]/=g;
96 V_to_G_GI(Za,d,G,GI); for(k=0;k<a;k++)
97 { int l; for(l=0;l<d;l++) Wa[l]=WZ->B[k][l]; for(l=0;l<d;l++)
98 { WZ->B[k][l]=0; for(j=0;j<d;j++) WZ->B[k][l] += Wa[j]*G[l][j];
99 } /* shift A to A'=A+B.x with x=dx=Z_ia*(F_a-A_a) */
100 X+=WZ->z[i][k]*(WZ->A[k]-F[k]);
101 } X%=WZ->m[i];/* IP = r*A+B.rI, B'=B.G^T => rI' = GI^T . rI + r dx */
102 if(X<-WZ->m[i]/2)X+=WZ->m[i]; else if(X>WZ->m[i]/2)X-=WZ->m[i];
103 if(Tout) printf("X=%ld\n",X);
104 for(k=0;k<a;k++) WZ->A[k]-=X*WZ->B[k][0];
105 for(j=0;j<d;j++) {Za[j]=WZ->rI[j]; WZ->rI[j]=0;} WZ->rI[0] =X*WZ->r;
106 for(j=0;j<d;j++) for(k=0;k<d;k++) WZ->rI[j]+=Za[k]*GI[k][j];
107 g=WZ->m[i]/g; for(k=0;k<a;k++) WZ->B[k][0]*=g; /* goto sublattice */
108 assert((WZ->rI[0]%g)==0); WZ->rI[0]/=g; shift*=g;
109 } X=1; for(j=0;j<*nz;j++) X*=WZ->m[j]; assert(X%shift==0); shift=X/shift;
110
111 if(Tout){for(j=0;j<a;j++)printf("%2ld ", WZ->A[j]); printf("=A rI=");
112 for(j=0;j<d;j++) printf("%2ld ", WZ->rI[j]);
113 puts("");
114 for(k=0;k<a;k++){for(j=0;j<d;j++)printf("%2ld ",WZ->B[k][j]);puts("=B");}
115 for(j=0;j<a;j++){X=0;for(k=0;k<d;k++)X+=WZ->B[j][k]*WZ->rI[k];
116 printf("%2ld ",X);}
117 puts("=B.rI");Write_Weight(WZ);}
118
119 { PairMat VPM; VertexNumList V; PolyPointList *P=WZ->P;
120 EqList *E=(EqList *)malloc(sizeof(EqList));assert(E!=NULL);
121 for(i=0;i<a;i++){E->e[i].c=WZ->A[i]; VPM[i][0]=WZ->d/WZ->w[i];
122 for(j=0;j<d;j++)E->e[i].a[j]=WZ->B[i][j];} E->ne=a;
123 WZ->P->n=d; WZ->P->np=0; Complete_Poly(VPM,E,1,WZ->P);
124 /* assert(POINT_Nmax>WZ->P->np); / * Print_PPL(WZ->P,""); obsolete */
125 /* for(i=0;i<d;i++) WZ->P->x[WZ->P->np][i]=WZ->rI[i]; / * obsolete */
126 for(i=0;i<d;i++) if(WZ->rI[i]) break;
127 if(i<d) /* rI!=0 => make P0=0 */
128 { for(i=1;i<P->np;i++) for(j=0;j<d;j++) P->x[i][j]-=P->x[0][j];
129 for(i=0;i<a;i++)for(j=0;j<d;j++) WZ->A[i]+=WZ->B[i][j]*P->x[0][j];
130 for(i=0;i<d;i++) {WZ->rI[i]-=WZ->r*P->x[0][i]; P->x[0][i]=0;}
131 }
132 Find_Equations(WZ->P,&V,E); for(i=0;i<E->ne;i++) E->e[i].c *= WZ->r;
133 for(i=0;i<E->ne;i++) if(1!=Eval_Eq_on_V(&E->e[i],WZ->rI,d)) break;
134 WZ->R = (i==E->ne); free(E);
135 } if(shift!=1) /* normalize sublattice data */
136 { Long *PB[W_Nmax], Z[POLY_Dmax][VERT_Nmax], M[POLY_Dmax], D[POLY_Dmax];
137 for(i=0;i<a;i++) PB[i]=WZ->B[i];
138 Sublattice_Basis(d,a,PB,Z,M,&k,G,D);
139 for(i=0;i<k;i++){X=0;for(j=0;j<a;j++)X+=Z[i][j];assert(0==(X%M[i]));}
140 for(i=0;i<k;i++){WZ->m[i]=M[i]; for(j=0;j<a;j++)WZ->z[i][j]=Z[i][j];}
141 WZ->M=k; /* printf("shift=%d rank=%d \n",shift,k); */
142 } X=0; for(j=0;j<a;j++) X+=WZ->w[j]*WZ->A[j]; /* X=w.A */
143 for(i=0;i<WZ->M;i++) {Long A=0,M=WZ->m[i],cx,cm,g; /* A=Z.A */
144 for(j=0;j<a;j++) A+=WZ->z[i][j]*WZ->A[j];
145 if(0==(A%=M)) continue;
146 g=Egcd(X%M,M,&cx,&cm); assert(0==(A%g)); /* Z_i -= cx * A/g * w_ */
147 cx*=A/g; if(cx>0)cx-=M;
148 for(j=0;j<a;j++) {WZ->z[i][j]-=cx*WZ->w[j]; WZ->z[i][j]%=M;}
149 }
150 X=0;for(j=0;j<a;j++)X+=WZ->w[j]*WZ->A[j];assert(X==WZ->d);/* TEST w.A=d */
151 for(j=0;j<a;j++){
152 X=WZ->r*WZ->A[j];
153 for(k=0;k<d;k++) X+=WZ->B[j][k]*WZ->rI[k];
154 assert(X==1);} /* TEST: r*A+B.rI==IP */
155 for(i=0;i<WZ->M;i++) {X=0; for(j=0;j<a;j++)
156 X+=WZ->z[i][j]*WZ->A[j];assert(0==X%WZ->m[i]);}/* TEST: A*z[]%m[]==0 */
157
158 if(FilterFlag) inFILE=NULL;
159 return 1;
160 }
161
162
163 #undef TEST_PP /* print some diagnostic info */
164 #define TEST_PD /* check Poincare duality of PP */
165
166 extern FILE *inFILE, *outFILE;
167
168 typedef struct {Long x[POINT_Nmax][W_Nmax]; int N, np;} AmbiPointList;
169 typedef struct {Long x[POLY_Dmax][W_Nmax]; int N, n;} AmbiLatticeBasis;
170
171 Long Wperm_to_GLZ(Long *W, int *d, Long **G, int *P);
172 void WeightLatticeBasis(Weight *_w, AmbiLatticeBasis *_B);
173 void WeightMakePoints(Weight *_W , AmbiPointList *_P);
174 int ChangeToTrianBasis(AmbiPointList*, AmbiLatticeBasis *, PolyPointList *);
175
176 /* ========== I/O functions: ========== */
177
IsDigit(char c)178 int IsDigit(char c) { return (('0'<=c) && (c<='9'));}
Write_Weight(Weight * W)179 void Write_Weight(Weight *W)
180 { int n; fprintf(outFILE,"%d",(int) W->d);
181 for(n=0; n < W->N; n++) fprintf(outFILE," %d",(int) W->w[n]);
182 #if WZinput
183 { int i,j; if(W->M>0)printf(" "); for(i=0;i<W->M;i++)
184 { fprintf(outFILE,"/Z%d: ",W->m[i]);
185 for(j=0;j<W->N;j++) fprintf(outFILE,"%d ",W->z[i][j]);
186 }
187 }
188 #endif
189 fprintf(outFILE,"\n"); /* puts(" = d w_1 ... w_N"); */
190 }
Read_Weight(Weight * _W)191 int Read_Weight(Weight *_W) /* read "d w_i" [ or "w_i d" if last=max ] */
192 { char c; int i, shift=1, I[W_Nmax+2], FilterFlag=(inFILE==NULL);
193 if(FilterFlag) inFILE=stdin; else
194 if(inFILE==stdin) printf("type degree and weights [d w1 w2 ...]: ");
195 c=fgetc(inFILE); if( !IsDigit(c) ) return 0; ungetc(c,inFILE);
196 fscanf(inFILE,"%d",I);
197 for(i=1; i<W_Nmax+2; i++)
198 { while(' '==(c=fgetc(inFILE))); ungetc(c,inFILE);
199 if(IsDigit(c)) fscanf(inFILE,"%d",&I[i]); else break;
200 } while(fgetc(inFILE )-'\n') if(feof(inFILE))return 0; /* read to EOL */
201 _W->N=i-1; if(_W->N>W_Nmax) { puts("Increase POLY_Dmax"); exit(0); }
202 for(i=0;i<=_W->N;i++) assert(I[i]>0);
203 if(I[_W->N]>I[0]) { _W->d=I[_W->N]; shift=0; } else _W->d=I[0];
204 for(i=0;i<_W->N;i++) {_W->w[i]=I[i+shift]; assert(_W->w[i]<_W->d);}
205 if(FilterFlag) inFILE=NULL;
206 return 1;
207 }
Write_WH(Weight * _W,BaHo * _BH,VaHo * _VH,int rc,int tc,PolyPointList * _P,VertexNumList * _V,EqList * _E)208 void Write_WH(Weight *_W, BaHo *_BH, VaHo *_VH, int rc, int tc,
209 PolyPointList *_P, VertexNumList *_V, EqList *_E){
210 int i, j;
211 #if WZinput
212 fprintf(outFILE,"%d ",(int)_W->d);
213 for(i=0;i<_W->N;i++) fprintf(outFILE,"%d ",(int)_W->w[i]);
214 for(i=0;i<_W->M;i++){fprintf(outFILE,"/Z%d: ",(int)_W->m[i]);
215 for(j=0;j<_W->N;j++)fprintf(outFILE,"%d ",(int)_W->z[i][j]);}
216 #else
217 for(i=0;i<_W->N;i++) fprintf(outFILE,"%d ",(int)_W->w[i]);
218 fprintf(outFILE,"%d=d ",(int)_W->d);
219 #endif
220 fprintf(outFILE,"M:%d %d ",_P->np,_V->nv); /* PolyData */
221 if(rc) fprintf(outFILE,"N:%d %d ",_BH->np,_E->ne);
222 else fprintf(outFILE,"F:%d ",_E->ne); /* END of PolyData */
223 if(tc&&(_P->n!=1+_VH->D)) {
224 int D=0, d=_W->d; /* LG non-geometric */
225 fprintf(outFILE,"LG: "); for(i=0;i<_W->N;i++) D+=d-2*_W->w[i];
226 if(D%d) {int g=Fgcd(D,d); fprintf(outFILE,"c/3=%d/%d\n",D/g,d/g);}
227 else { int r=_W->N-(D/=d); if((r%2)||(r<=2)) r=0;assert(D==_VH->D);
228 for(i=0;i<=D;i++) {fprintf(outFILE,"%sH%d:",i?" ":"",i);
229 for(j=0;j<=D-i;j++)fprintf(outFILE,"%s%d",j?",":"",(int)_VH->h[i][j]);}
230 #if WZinput
231 /* if(r) fputs(_W->R ? " RefGC": " NOrgc",outFILE); */
232 if(r) { if(_W->R) fprintf(outFILE," RefI%d",_W->r);
233 else fputs(" nonRG",outFILE); }
234 #endif
235 fputs("\n",outFILE); }
236 return ;
237 }
238 if((!tc)&&(!rc)){fprintf(outFILE,"non-transversal\n");return;}
239 if(_P->n>2){
240 if(tc&&rc) for(i=1;i<_P->n-1;i++) if(_VH->h[1][i]!=_BH->h1[i])
241 { puts("Vafa!=Batyrev"); for(i=1;i<_P->n-1;i++) printf(
242 "V[1,%d]=%d B[1,%d]=%d\n",i,_VH->h[1][i],i,_BH->h1[i]);exit(0);}
243 if(tc) {
244 fprintf(outFILE,"V:%d",_VH->h[1][1]);
245 for(i=2;i<_P->n-1;i++) fprintf(outFILE,",%d",_VH->h[1][i]);
246 if(_P->n > 5){
247 for (j=2; 2*j <= _P->n-1; j++){
248 fprintf(outFILE,";%d",_VH->h[j][j]);
249 for(i=j+1;i<_P->n-j;i++) fprintf(outFILE,",%d",_VH->h[j][i]);}}
250 }
251 else if(rc) {
252 fprintf(outFILE,"H:%d",_BH->h1[1]);
253 for(i=2;i<_P->n-1;i++) fprintf(outFILE,",%d",_BH->h1[i]); }
254 if(6<_P->n) fprintf(outFILE," [???]\n");
255 else if(tc||rc) /* Euler number */ {
256 int chi=0, *ho=_BH->h1; if(tc) ho=_VH->h[1];
257 if(_P->n==3) chi=4+ho[1];
258 if(_P->n==4) chi=2*(ho[1]-ho[2]);
259 if(_P->n==5) chi=48+6*(ho[1]-ho[2]+ho[3]);
260 if(_P->n==6) chi=24*(ho[1]-ho[2]+ho[3]-ho[4]);
261 fprintf(outFILE," [%d]\n",chi); }
262 else fprintf(outFILE,"\n");}
263 else fprintf(outFILE,"\n");
264 }
265
TEST_LatticeBasis(AmbiLatticeBasis * _B)266 void TEST_LatticeBasis(AmbiLatticeBasis *_B) /* print AmbiLatticeBasis */
267 { int p,a; for(p=0;p<_B->n;p++){
268 for(a=0;a<_B->N;a++) printf(" %3d",(int) _B->x[p][a]);
269 puts(""); }
270 }
271 #ifdef TEST
TEST_WeightMakePoints(AmbiPointList * _P)272 void TEST_WeightMakePoints(AmbiPointList *_P)
273 { int i,j; static int MaxPoNum; if(_P->np>MaxPoNum) MaxPoNum = _P->np;
274 if(_P->np>20)
275 { for(i=0;i<_P->np;i++)
276 { for(j=0;j<_P->N;j++) printf("%d ",_P->x[i][j]); puts("");}
277 }
278 else
279 { for(i=0;i<_P->N;i++)
280 { for(j=0;j<_P->np;j++) printf("%4d",_P->x[j][i]); puts("");}
281 }
282 printf("PointNum=%d [max=%d]\n",_P->np,MaxPoNum);
283 }
284 #endif
285
Ambi_2_Lattice(Long * A,AmbiLatticeBasis * B,Long * P)286 void Ambi_2_Lattice(Long *A,AmbiLatticeBasis *B,Long *P)
287 { int i, p=B->n; while(p--)
288 { int a=p + B->N - B->n; P[p]=A[a] /* -Ai[a]*/ ;
289 for(i=p+1; i<B->n; i++) P[p] -= P[i]*B->x[i][a];
290 if(P[p] % B->x[p][a]) { puts("Error in BasisChange!"); exit(0);}
291 else P[p] /= B->x[p][a];
292 }
293 }
Make_Poly_Points(Weight * _W_in,PolyPointList * _PP)294 void Make_Poly_Points(Weight *_W_in, PolyPointList *_PP)
295 { AmbiLatticeBasis B; Weight *_W=_W_in; int /* index=0,*/ nip;
296 AmbiPointList * _AP = (AmbiPointList *) malloc(sizeof (AmbiPointList));
297 #ifndef NO_COORD_IMPROVEMENT /* ==== Perm Coord Improvement ==== */
298 Weight Waux=(*_W); _W=&Waux;
299 { int i,n=_W->N,pi[AMBI_Dmax]; Long M[AMBI_Dmax][AMBI_Dmax],
300 *G[AMBI_Dmax]; for(i=0;i<n;i++) G[i]=M[i];
301 Wperm_to_GLZ(_W->w,&n,G,pi);for(i=0;i<n;i++)Waux.w[i]=_W_in->w[pi[i]];
302 }
303 #endif /* = End of Perm Coord Improvement = */
304 if(_AP==NULL) {printf("Unable to allocate space for _AP\n"); exit(0);}
305 WeightLatticeBasis(_W, &B); /* TEST_LatticeBasis(&B); */
306 WeightMakePoints(_W, _AP);
307 #ifdef TEST
308 puts("\nWeights:"); Write_Weight(_W);
309 puts("AmbiPoints:"); TEST_WeightMakePoints(_AP);
310 puts("Basis:"); TEST_LatticeBasis(&B);
311 printf("POINT_Num=%d\n",_AP->np);
312 #endif
313 assert(_AP->np <= POINT_Nmax);
314 nip=ChangeToTrianBasis(_AP, &B, _PP);
315 #ifdef TEST
316 Print_PPL(_PP,"PolyPoints:");
317 #endif
318 #ifdef GEN_CY
319 { int i; for(i=0;i<_W->N;i++) index+=_W->w[i]; if(index%_W->d)index=0;}
320 index/=_W->d; if(index>1)
321 { Long A[AMBI_Dmax];int i;for(i=0;i<B.N;i++)A[i]=1-index*_AP->x[nip][i];
322 assert(_PP->np<POINT_Nmax); /* need one more to store IP for GEN_CY */
323 Ambi_2_Lattice(A,&B,_PP->x[_PP->np]);
324 }
325 #endif
326 free(_AP); return;
327 }
328 #if (WZinput)
Read_W_PP(Weight * W,PolyPointList * P)329 int Read_W_PP(Weight *W, PolyPointList *P){ W->P=P; return Read_WZ_PP(W); }
330 #else
Read_W_PP(Weight * _W,PolyPointList * _PP)331 int Read_W_PP(Weight *_W, PolyPointList *_PP){
332 if (!Read_Weight(_W)) return 0; Make_Poly_Points(_W,_PP); return 1;
333 }
334 #endif
335
336 /* ========== WeightLatticeBasis(Weight *_w, AmbiLatticeBasis *_B)
337 *
338 * M[0]:= (-n1, n0, 0, ...) / gcd(n0,n1);
339 * M[i]:= (0,...,0,g/ng,0,...)- (ni/ng) * egcd.Vout(n0,...,n(i-1),0,...);
340 * with g=gcd(n0,...,n[i-1]); ng=gcd(g,ni);
341 * */
342 #ifdef TEST
NormTriangularBasis(AmbiLatticeBasis * _B)343 void NormTriangularBasis(AmbiLatticeBasis *_B)
344 { int p=_B->n-1, a=_B->N-1; while(p--)
345 { int pi=_B->n; while(0==_B->x[p][--a]);
346 while(--pi > p)
347 { int ai=a+1, r = _B->x[pi][a] / _B->x[p][a];
348 /**/ if((_B->x[pi][a] % _B->x[p][a]) < 0) r--;
349 /**/ if(r) while(ai--) _B->x[pi][ai] -= r * _B->x[p][ai];
350 }
351 }
352 }
OldWeightLatticeBasis(Weight * _w,AmbiLatticeBasis * _B)353 void OldWeightLatticeBasis(Weight *_w, AmbiLatticeBasis *_B)
354 { int i, j;
355 Long *V, Vin[W_Nmax], Vout[W_Nmax];
356 for(i=0; i<_w->N; i++) Vin[i]=_w->w[i];
357 V=_B->x[0]; *Vout=Fgcd(Vin[0],Vin[1]);
358 V[0] = -Vin[1] / *Vout; V[1] = *Vin / *Vout; _B->N=_w->N; _B->n=_w->N-1;
359 for(j=2; j<_w->N; j++) V[j]=0;
360 for(i=2; i<_w->N; i++) {
361 int g=REgcd(Vin, &i, Vout), gn=Fgcd(g,_w->w[i]); V=_B->x[i-1];
362 V[i] = g/gn; g = _w->w[i]/gn;
363 for(j=0;j<i;j++) V[j]=-Vout[j]*g;
364 j++; while(j<_w->N) V[j++]=0;
365 }
366 NormTriangularBasis(_B); /* TEST_LatticeBasis(_B); */
367 }
368 #endif
WeightLatticeBasis(Weight * _w,AmbiLatticeBasis * _B)369 void WeightLatticeBasis(Weight *_w, AmbiLatticeBasis *_B)
370 { Long *B[W_Nmax], E[W_Nmax]; int i; _B->N=_w->N; _B->n=_w->N-1;
371 B[0]=E; for(i=0;i<_B->n;i++) B[i+1]=_B->x[i]; W_to_GLZ(_w->w,&(_w->N),B);
372 }
373
374 /* ========== WeightMakePoints(Weight *_W, AmbiPointList *_P)
375 *
376 * MakePoints: X[i]-IP[i] = x[j] * B[j][i] :: i-1<=j<B.n,
377 * with \sum X_i w_i=d and 0<=X[i].
378 * */
WeightMakePoints(Weight * _W,AmbiPointList * _P)379 void WeightMakePoints(Weight *_W, AmbiPointList *_P)
380 { int i=_W->N-1, X_max[W_Nmax], d_Rest[W_Nmax], X[W_Nmax+1]={0};
381 _P->np=0; _P->N=_W->N; d_Rest[i]=_W->d; X_max[i]=d_Rest[i]/_W->w[i];
382 while(i<_W->N)
383 { if(X[i]>X_max[i]) X[++i]++;
384 else if(i)
385 { d_Rest[i-1]=d_Rest[i]-X[i]*_W->w[i]; X[--i]=0;
386 X_max[i]=d_Rest[i]/_W->w[i];
387 }
388 else
389 { int j;
390 Long *Y;
391 *d_Rest=d_Rest[1]-X[1]*_W->w[1];
392 if( !(*d_Rest % *(_W->w)) )
393 { if(POINT_Nmax <= _P->np){puts("increase POINT_Nmax");exit(0);}
394 Y=_P->x[(_P->np)++];
395 *Y= *d_Rest/ *(_W->w); for(j=1;j<_W->N;j++) Y[j]=X[j];
396 }
397 X[1]++; i=1;
398 }
399 }
400 }
401
402 /* ========== ChangeToTrianBasis(AmbiPointList *, Basis *, PolyPointList *)
403 *
404 * Assuming that AmbiLatticeBasis is of triangular form, i.e.:
405 * _B->x[p][a] = 0 for a < p + CoDim ... a < W_Nmax and p < POLY_Dmax
406 * this solves AP=B.PP for AP (to be precise we must take AP.x-IP,
407 * where IP is the unique interior point.)
408 * DP=AP.x-IP; DP[a]=PP[p]*B[p][a] =>
409 * PP[p]*B[p][p+1] = DP[p+1] - sum_{a>p+1} PP[p]*B[p][p+1]
410 * */
ChangeToTrianBasis(AmbiPointList * _AP,AmbiLatticeBasis * _B,PolyPointList * _PP)411 int ChangeToTrianBasis(AmbiPointList *_AP,
412 AmbiLatticeBasis *_B, PolyPointList *_PP)
413 { int n, ipcount=0, nIP=0;
414 if( _AP->N - _B->N) {puts("Dimensions don't match!"); exit(0);}
415 else { _PP->n = _B->n; _PP->np = _AP->np; }
416 for(n=0; n<_AP->np; n++)
417 { int i, p=1; for(i=0; i < _B->N; i++) if( ! _AP->x[n][i] ) p=0;
418 if(p) { ipcount++; nIP=n; }
419 }
420 #ifdef TEST
421 if(ipcount-1) puts("*** Wrong IP-count: take any point as origin! ***");
422 #endif
423 for(n=0; n<_AP->np; n++)
424 { int i, p=_B->n;
425 while(p--)
426 { int a=p + _B->N - _B->n; _PP->x[n][p]=_AP->x[n][a]-_AP->x[nIP][a];
427 for(i=p+1; i<_B->n; i++) _PP->x[n][p] -= _PP->x[n][i]*_B->x[i][a];
428 if(_PP->x[n][p] % _B->x[p][a])
429 { puts("Error in BasisChange!"); exit(0);}
430 else _PP->x[n][p] /= _B->x[p][a];
431 }
432 } return nIP;
433 }
434 /* ============= make weights =========== */
IfRefWWrite(Weight * W,PolyPointList * P)435 int IfRefWWrite(Weight *W, PolyPointList *P)
436 { VertexNumList V; EqList E; Make_Poly_Points(W, P);
437 if(Ref_Check(P,&V,&E)) {Write_Weight(W); fflush(stdout); return 1;}
438 else return 0;
439 }
Rec_RefWeights(Weight * W,PolyPointList * P,int g,int sum,int * npp,int * nrp,int n)440 void Rec_RefWeights(Weight *W, PolyPointList *P, int g, int sum, int *npp,
441 int *nrp, int n)
442 { int wmax=W->d/(W->N-n+1); wmax=min(wmax,W->w[n+1]);
443 wmax=min(wmax,sum-n);
444 if(n) for(W->w[n]=wmax;(n+1)*W->w[n]>=sum;W->w[n]--)
445 Rec_RefWeights(W,P,Fgcd(g,W->w[n]),sum-W->w[n],npp,nrp,n-1);
446 else if(1==Fgcd(g,W->w[0]=sum)) {(*npp)++;if(IfRefWWrite(W,P))(*nrp)++;};
447 }
MakeRefWeights(int N,int from_d,int to_d)448 void MakeRefWeights(int N, int from_d, int to_d)
449 { int npp=0, nrp=0; Weight W;
450 PolyPointList *P = (PolyPointList *) malloc (sizeof(PolyPointList));
451 W.N=N; assert((N<=W_Nmax)&&(N<POLY_Dmax+2)); assert(P!=NULL);
452 for(W.d=from_d;W.d<=to_d;W.d++)
453 for(W.w[N-1]=W.d/2; W.d <= N*W.w[N-1]; W.w[N-1]--)
454 Rec_RefWeights(&W,P,Fgcd(W.d,W.w[W.N-1]),W.d-W.w[W.N-1],&npp,&nrp,N-2);
455 fprintf(outFILE,"#primepartitions=%d #refpolys=%d\n",npp,nrp);exit(0);
456 }
457
458 /* ============= Landau-Ginzburg-Calculations: =========== */
459
PrintPoCoLi(PoCoLi * P)460 void PrintPoCoLi(PoCoLi *P)
461 { int i; for(i=0;i<P->n;i++){
462 #if ABBREV_POLY_PRINT
463 int a=ABBREV_POLY_PRINT;
464 if((a<=i)&&(i<P->n-a)) {if(i==a) printf("+ ... ");} else
465 #endif
466 { int co=abs(P->c[i]); printf("%s ", i ?
467 ((P->c[i] > 0) ? "+" : "-" ) : "");if(!i&&(P->c[0]<0))printf("- ");
468 if(co!=1) printf("%d ",co);
469 printf("x^%d ",(int)P->e[i]); }}
470 puts(i ? "" : " 0");
471 }
UnitPoly(PoCoLi * P)472 void UnitPoly(PoCoLi *P) { P->n=1; P->e[0]=0; P->c[0]=1; }
Add_Mono_2_Poly(int e,Pint c,PoCoLi * P)473 void Add_Mono_2_Poly(int e, Pint c, PoCoLi *P) /* use bisection */
474 { int m=-1, M=P->n;
475 if(P->n) if(e <= P->e[P->n-1]) while(M-m>1)
476 { int pos=(M+m)/2; if(e > P->e[pos]) m=pos; else
477 if (e < P->e[pos]) M=pos; else { P->c[pos]+=c; return; }/* e exists */
478 } if(P->A<=P->n){printf("n=A=%d\n",P->A); fflush(0);
479 assert(P->n < P->A);} /* check #(coeff.) of Poly. */
480 for(m=P->n++;M<m;m--) { /* insert new exponent at M */
481 P->c[m]=P->c[m-1]; P->e[m]=P->e[m-1]; } P->e[m]=e;P->c[m]=c;
482 #ifdef TEST_PP
483 {static int M=1; if(P->n/1000000 > M) {printf("#c=%dM ",++M);fflush(0);}}
484 #endif
485 }
Init1_xN(PoCoLi * P,int N)486 void Init1_xN(PoCoLi *P,int N) /* 1 - x^N */
487 { UnitPoly(P); Add_Mono_2_Poly(N,-1,P);
488 }
Remove_Zeros(PoCoLi * AB)489 void Remove_Zeros(PoCoLi *AB)
490 { int i, s=0; for(i=0; i<AB->n; i++) /* s=shift: eliminatet zeros */
491 { if(0==AB->c[i]) s++; else
492 if(s) { AB->c[i-s]=AB->c[i]; AB->e[i-s]=AB->e[i]; }
493 } AB->n -= s;
494 }
PolyCopy(PoCoLi * X,PoCoLi * Y)495 void PolyCopy(PoCoLi *X,PoCoLi *Y)
496 { assert(X->n<=Y->A); for(Y->n=0;Y->n<X->n;Y->n++){
497 Y->e[Y->n]=X->e[Y->n]; Y->c[Y->n]=X->c[Y->n];}
498 }
BottomUpQuot(PoCoLi * N,PoCoLi * D,PoCoLi * Q,PoCoLi * R)499 int BottomUpQuot(PoCoLi *N,PoCoLi *D,PoCoLi *Q,PoCoLi *R) /* Q*D = N-R */
500 { int i, c, Npos, e, E, dD, dN; Q->n=R->n=0; /* return R==0 */
501 assert(D->n>0); if(N->n==0) return 1; assert(N->n>0); /* assume D != 0 */
502 dD=D->e[D->n-1]-D->e[0]; dN=N->e[N->n-1]-N->e[0];
503 e=N->e[0]-D->e[0]; if((e<0)||(dD>dN)) { PolyCopy(N,R); return 0; }
504 for(Npos=0;Npos<N->n;Npos++)
505 { if(N->e[Npos]<=e+dD)Add_Mono_2_Poly(N->e[Npos],N->c[Npos],R);
506 else break;
507 } E = N->e[N->n-1] - D->e[D->n-1];
508 while(0==(R->c[0]%D->c[0]))
509 { Add_Mono_2_Poly(e,c=R->c[0]/D->c[0],Q); /* next Q.c[] */
510 for(i=0;i<D->n;i++)Add_Mono_2_Poly(e+D->e[i],-c*D->c[i],R); /* R-cQ */
511 Remove_Zeros(R); if((Npos==N->n)&&(0==R->n)) return 1; /* finished? */
512 e = (R->n) ? R->e[0]-D->e[0] : N->e[Npos]-D->e[0]; /* next Q.e[] */
513 while(Npos<N->n) { if(N->e[Npos]<=e+dD)
514 { Add_Mono_2_Poly(N->e[Npos],N->c[Npos],R); Npos++;} else break;
515 Remove_Zeros(R); if(e > E) return 0;
516 }
517 } while(Npos<N->n){Add_Mono_2_Poly(N->e[Npos],N->c[Npos],R); Npos++;}
518 return 0;
519 }
PolyProd(PoCoLi * A,PoCoLi * B,PoCoLi * AB)520 void PolyProd(PoCoLi *A,PoCoLi *B,PoCoLi *AB) /* AB = A*B */
521 { int i, j, s=A->n+B->n-1; AB->n=0; assert(s<AB->A);
522 for(i=0;i<s;i++)
523 { int m = (i<B->n) ? 0 : i+1-B->n, M = (i<A->n) ? i+1 : A->n;
524 for(j=m;j<M;j++)
525 Add_Mono_2_Poly(A->e[j]+B->e[i-j],A->c[j]*B->c[i-j],AB);
526 } Remove_Zeros(AB);
527 }
Poly_Sum(PoCoLi * A,PoCoLi * B,PoCoLi * S)528 void Poly_Sum(PoCoLi *A,PoCoLi *B,PoCoLi *S) /* S = A+B */
529 { int a, b=0; S->n=0; for(a=0;a<A->n;a++) { int q=1;
530 while(b<B->n) if(A->e[a]<=B->e[b]) break; else {Pint s=B->c[b++];
531 if(s) {if(S->n>=S->A){printf("S.n>%d in S=A+B\n",S->n);exit(0);}
532 S->c[S->n]=s; S->e[S->n++]=B->e[b-1];}}
533 if(b<B->n)if(B->e[b]==A->e[a]){Pint s=A->c[a]+B->c[b++]; q=0;
534 if(s) {assert(S->n<S->A); S->e[S->n]=A->e[a]; S->c[S->n++]=s;} }
535 if(q) if(A->c[a]) {assert(S->n<S->A); S->e[S->n]=A->e[a];
536 S->c[S->n++]=A->c[a];} }
537 while(b<B->n) if(B->c[b]) {assert(S->n<S->A);
538 S->e[S->n]=B->e[b]; S->c[S->n++]=B->c[b++];} else b++;
539 }
Poly_Dif(PoCoLi * A,PoCoLi * B,PoCoLi * D)540 void Poly_Dif(PoCoLi *A,PoCoLi *B,PoCoLi *D) /* D = A-B */
541 { int a, b=0; D->n=0; for(a=0;a<A->n;a++) { int q=1;
542 while(b<B->n) if(A->e[a]<=B->e[b]) break; else {Pint d=-B->c[b++];
543 if(d) {assert(D->n<D->A); D->c[D->n]=d; D->e[D->n++]=B->e[b-1];}}
544 if(b<B->n)if(B->e[b]==A->e[a]){Pint d=A->c[a]-B->c[b++]; q=0;
545 if(d) {assert(D->n<D->A); D->e[D->n]=A->e[a]; D->c[D->n++]=d;} }
546 if(q) if(A->c[a]) {assert(D->n<D->A); D->e[D->n]=A->e[a];
547 D->c[D->n++]=A->c[a];} }
548 while(b<B->n) if(B->c[b]) {assert(D->n<D->A);
549 D->e[D->n]=B->e[b]; D->c[D->n++]=-B->c[b++];} else b++;
550 }
AllocPoCoLi(PoCoLi * P)551 void AllocPoCoLi(PoCoLi *P) /* allocate e[A] and c[A] */
552 { assert(0<P->A); assert( NULL != ( P->e = (int *) malloc( P->A *
553 (sizeof(int)+sizeof(Pint) ) ) )); P->c = (Pint *) & P->e[P->A];
554 // printf("AllocPoCoLi: P->A = %d\n", P->A);
555 }
Free_PoCoLi(PoCoLi * P)556 void Free_PoCoLi(PoCoLi *P) { free(P->e); } /* free P.e and P.c */
PoincarePoly(int N,int * w,int d,PoCoLi * P,PoCoLi * Z,PoCoLi * R)557 void PoincarePoly(int N, int *w, int d, PoCoLi *P, PoCoLi *Z,PoCoLi *R)
558 { int i, e[2]; Pint c[2]; PoCoLi B, *In=Z,*Out=P,*aux;
559 if(N==0) {UnitPoly(P); return;}
560 B.e=e; B.c=c; B.A=B.n=2; e[0]=0; c[0]=1; c[1]=-1; Init1_xN(In,d - w[0]);
561 for(i=1;i<N;i++)
562 { e[1]=d-w[i]; PolyProd(In,&B,Out); aux=Out; Out=In; In=aux;
563 } /* printf("\nN =");PrintPoCoLi(Z); */
564 while(i--)
565 { e[1]=w[i]; assert(BottomUpQuot(In,&B,Out,R)); aux=Out; Out=In; In=aux;
566 } /* printf("Q =");PrintPoCoLi(P); */ assert((R->n)==0);
567 #ifdef TEST_PD
568 { int M=P->n-1, I=(M+1)/2, E=P->e[P->n-1];
569 for(i=0;i<I;i++) {assert(P->e[i]==E-P->e[M-i]);
570 assert(P->c[i]==P->c[M-i]); }}
571 #endif
572 }
573
574
575 /* ====== I N D E X & T R A C E / b01 computation of VaHo ==== */
576
Print_VaHo(VaHo * V)577 void Print_VaHo(VaHo *V)
578 { int i,j,D=V->D; for(i=0;i<=D;i++){fprintf(outFILE,"H%d*: ",i);
579 for(j=0;j<=D;j++) fprintf(outFILE,"%d ",V->h[i][j]);}
580 }
DoHodgeTest(VaHo * V)581 int DoHodgeTest(VaHo *V)/*[holo] Poincare duality, Hodge duality, sum rule */
582 { int i,j,D=V->D,X;for(i=D/2;i<=D;i++){if(V->h[i][0]!=V->h[D-i][0])return 0;
583 for(j=0;j<=D;j++)if(V->h[i][j]!=V->h[D-i][D-j]) return 0;} V->E=X=0;
584 for(i=0;i<=D;i++)for(j=i+1;j<=D;j++)if(V->h[i][j]!=V->h[j][i]) return 0;
585 for(i=0;i<=D;i++)for(j=0;j<=D;j++){int x=2*i-D; x=x*x*V->h[i][j];
586 if((i+j)%2) {V->E -= V->h[i][j];X-=x;} else {V->E += V->h[i][j];X+=x;}}
587 assert(3*X==D*V->E); return (V->h[0][0]==1);
588 }
Hodge_Test(VaHo * V)589 int Hodge_Test(VaHo *V) /* [holo] Poincare duality, Hodge duality, sum rule */
590 { if(DoHodgeTest(V)) return 1; else {Print_VaHo(V); return 0;}
591 }
592 /* do {;} while(Multiloop(N,I,j,J)); <--> do forall 0<=I[j]<N[j], 0<=j<J */
Init_Multiloop(int * N,int * I,int * j,int * J)593 int Init_Multiloop(int *N, int *I, int *j, int *J)
594 { long X=1; for(*j=0;*j<*J;(*j)++) {I[*j]=0;assert(N[*j]>0);} assert(0<=*J);
595 if((*J)==0) {N[0]=1; I[0]=0;} /* safe also for J=0 */
596 for(*j=0;*j<*J;(*j)++) X*=N[*j];
597 *j=0; return X; /* return order=prod(N) */
598 }
Multiloop(int * N,int * I,int * j,int * J)599 int Multiloop(int *N,int *I,int *j,int *J)/* need j=I[]=0 => Init_Multiloop */
600 { assert((*j)==0); if(++(I[0])<N[0]) return 1; while(N[*j]<=I[*j])
601 { I[(*j)++]=0; if(*J<=*j) return 0; I[*j]++;} (*j)=0; return 1;
602 }
603
Count_b01(Weight * W)604 int Count_b01(Weight *W)
605 { int b=0,c=0,l,M[POLY_Dmax],I[POLY_Dmax],eq[W_Nmax],j,J=W->M,N=W->N,vac=0;
606 for(j=0;j<J;j++) M[j]=W->m[j];
607 for(l=0;l<W->d;l++) /* Z_d==GSO */
608 { c=Init_Multiloop(M,I,&j,&J); do { /* BEGIN generate group */
609 int i,can=1, sum=0, k; /* can(didate) for b01-contribution */
610 for(i=0;can&&(i<N);i++){Rat th=rR(l*W->w[i],W->d); /* th=0/q -> eq[] */
611 for(k=0;k<J;k++){th=rS(th,rR(I[k]*W->z[k][i],M[k]));} th.N%=th.D;
612 assert(th.N>=0);
613 if(th.N) {th=rP(th,rR(W->d,W->w[i])); if(th.N!=th.D) can=0;
614 else {eq[i]=1;sum+=W->d-2*W->w[i];}}
615 else eq[i]=0;} /* eq=(th_i==q_i); */
616 if(sum&&(sum != W->d)) can=0; /* i.e. charge=1? */
617 for(k=0;can&&(k<J);k++){int gph=0; /* check inv. under G_k: gph(ase) */
618 for(i=0;i<N;i++) {if(eq[i]) gph+=W->z[k][i]; gph%=M[k];}
619 if(gph) can=0; }
620 if(can) {if(sum) b++; else vac++;}
621 #if (SHOW_b01_TWIST)
622 if(can&&sum){ if(b==1)puts("");printf("b01=%d: j^%d, I[]=",b,l);
623 for(k=0;k<J;k++)printf(" %d",I[k]);puts("");}
624 /* printf("M[%d]=",J);for(i=0;i<J;i++)printf("%d ",M[i]); */
625 /* printf(" I=");for(i=0;i<J;i++)printf("%d ",I[i]);printf(" j=%d\n",j);*/
626 #endif
627 assert(0<(c--));} while(Multiloop(M,I,&j,&J)); assert(c==0); /* END */
628 } assert(vac==1); return b;
629 }
Fast_c9_VaHo(Weight * W,VaHo * V)630 void Fast_c9_VaHo(Weight *W,VaHo *V)/* -->> V.D<=3 <<-- via index and trace */
631 { int i,j,k,l,c, b01=0; /* ... from "proced" in lgotwist.c */
632 int fac, wort, nvar, expo, ns=W->M, N=W->N, R, mask[W_Nmax+1], n=N;
633 Long zsum1=0, zsum2=0, mo=W->d, ng, ngb, kgV=mo, omega[POLY_Dmax];
634 int *wo[3],*woG,*woA,*woS, WM=1<<W_Nmax,I[POLY_Dmax]; Rat prod; V->D=0;
635 for(i=0;i<W->N;i++) V->D+= W->d-2*W->w[i];
636 assert(V->D%W->d==0); V->D/=W->d; assert((0<V->D)&&(V->D<=3));
637 V->h[0][0]=V->h[0][V->D]=V->h[V->D][0]=V->h[V->D][V->D]=1;
638 if(V->D==1){assert(Count_b01(W)==1); return;}
639 if(V->D==2){b01=Count_b01(W);
640 V->h[0][1]=V->h[1][0]=V->h[2][1]=V->h[1][2]=b01;
641 assert((b01==0)||(b01==2)); V->h[1][1]=(b01) ? 4 : 20; return;}
642 wo[0] = woG = (int *) malloc(WM*3*sizeof(int));
643 assert(woG!=NULL); wo[1]=woS=&woG[WM]; wo[2]=woA=&woS[WM];
644 for(i=0;i<=N;i++) mask[i]=1<<i;
645 for(i=0;i<=WM;i++) woG[i]=woS[i]=woA[i]=0;
646 for(k=0;k<ns;k++) {kgV *= W->m[k]/Fgcd(kgV,W->m[k]); mo*=W->m[k];}
647 for(k=0;k<ns;k++) omega[k]=kgV/W->m[k];
648 R=kgV/W->d;
649 for(l=0;l<W->d;l++) /* Z_d==GSO */
650 { c=Init_Multiloop(W->m,I,&j,&ns); do { /* BEGIN generate group */
651 nvar=0; for(i=0;i<N;i++) { expo=l*W->w[i]*R; for(k=0;k<ns;k++)expo+=
652 I[k]*omega[k]*W->z[k][i];if(!(expo%kgV))nvar+=mask[i];}woS[nvar]++;
653 assert(0<(c--));}while(Multiloop(W->m,I,&j,&ns));assert(c==0);/* END */
654 } assert(woS[mask[n]-1]==1); /* over=overcount=wo[mask[n]-1][1]; */
655 for(i=0;i<mask[n];i++) { fac = 1;
656 for(j=0;j<n;j++) if ((mask[j] & i) == mask[j]) fac*=-1;
657 if(woS[i]) for(k=0;k<mask[n];k++) if(woS[k]){ /* woS = sum of contrib */
658 wort = i & k;
659 wo[fac+1][wort]+=woS[i]*woS[k]; } }
660 for(i=mask[n]-1;i>=0;i--) if (woG[i]+woA[i]) {
661 prod=rI(1);
662 for (j=0;j<n;j++) if ((mask[j]&i)==mask[j])
663 prod = rP(prod,rR(W->w[j]-W->d,W->w[j]));
664 if (prod.D != 1){ fprintf(outFILE,
665 "\nDenominator != 1 in Fast_c9_VaHo (LG.c)\n");exit(0);}
666 zsum1+= woG[i]*prod.N; zsum2+= woA[i]*prod.N;}
667 ng = -(zsum1/mo+2)/2; ngb = (zsum2/mo-2)/2; /* /over */
668 assert((zsum1==-2*(ng+1)*mo) && (zsum2 == 2*(ngb+1)*mo)); /* *over */
669 free(woG); /* woS[word] = #group elements :: survivors==word */
670 /* woG / woA = contributions to ng / ngb */
671 if(ng==ngb){b01=Count_b01(W); assert((b01==0)||(b01==1)||(b01==3));}
672 for(i=1;i<=2;i++){
673 V->h[i][i]=ngb-2*b01; V->h[i][3-i]=ng-2*b01;
674 V->h[0][i]=V->h[3][i]=b01; V->h[i][0]=V->h[i][3]=b01;
675 V->h[0][0]=V->h[3][3]=V->h[3][0]=V->h[0][3]=1;
676 }
677 }
678
WIndex_HTrace(Weight * W,int * WI,int * T)679 int WIndex_HTrace(Weight *W, int *WI, int *T)/* T=sum(Hij), return over=H00 */
680 { int i,j,k,l,c,vacnum; /* ... from "proced" in lgotwist.c */
681 int fac, wort, nvar, expo, ns=W->M, N=W->N, R, mask[W_Nmax+1], n=N;
682 Long zsum1=0, zsum2=0, mo=W->d, kgV=mo, omega[POLY_Dmax];
683 int *wo[3],*woG,*woA,*woS, WM=1<<W_Nmax,I[POLY_Dmax]; Rat prod;
684 wo[0] = woG = (int *) malloc(WM*3*sizeof(int));
685 assert(woG!=NULL); wo[1]=woS=&woG[WM]; wo[2]=woA=&woS[WM];
686 for(i=0;i<=N;i++) mask[i]=1<<i;
687 for(i=0;i<=WM;i++) woG[i]=woS[i]=woA[i]=0;
688 for(k=0;k<ns;k++) {kgV *= W->m[k]/Fgcd(kgV,W->m[k]); mo*=W->m[k];}
689 for(k=0;k<ns;k++) omega[k]=kgV/W->m[k];
690 R=kgV/W->d;
691 for(l=0;l<W->d;l++) /* Z_d==GSO */
692 { c=Init_Multiloop(W->m,I,&j,&ns); do { /* BEGIN generate group */
693 nvar=0; for(i=0;i<N;i++) { expo=l*W->w[i]*R; for(k=0;k<ns;k++)expo+=
694 I[k]*omega[k]*W->z[k][i];if(!(expo%kgV))nvar+=mask[i];}woS[nvar]++;
695 assert(0<(c--));}while(Multiloop(W->m,I,&j,&ns));assert(c==0);/* END */
696 } vacnum=woS[mask[n]-1]; /* over=overcount=wo[mask[n]-1][1]; */
697 for(i=0;i<mask[n];i++) { fac = 1;
698 for(j=0;j<n;j++) if ((mask[j] & i) == mask[j]) fac*=-1;
699 if(woS[i]) for(k=0;k<mask[n];k++) if(woS[k]){ /* woS = sum of contrib */
700 wort = i & k;
701 wo[fac+1][wort]+=woS[i]*woS[k]; } }
702 for(i=mask[n]-1;i>=0;i--) if (woG[i]+woA[i]) {
703 prod=rI(1);
704 for (j=0;j<n;j++) if ((mask[j]&i)==mask[j])
705 prod = rP(prod,rR(W->w[j]-W->d,W->w[j]));
706 if (prod.D != 1){ fprintf(outFILE,
707 "\nDenominator != 1 in Fast_c9_VaHo (LG.c)\n");exit(0);}
708 zsum1+= woG[i]*prod.N; zsum2+= woA[i]*prod.N;}
709 free(woG); /* woS[word] = #group elements :: survivors==word */
710 /* woG / woA = contributions to ng / ngb */
711 assert(zsum1%mo==0); assert(zsum2%mo==0); assert(vacnum==1);
712 *WI=(zsum2+zsum1)/mo; *T=(zsum2-zsum1)/mo; return vacnum;
713 }
714 /* Z_d is taken care of by the integral charge condition. *
715 * x[j]=lcm(phase denom.), go over G=\prod M[j], denom=\prod(1-t^(w*x)) *
716 * numerator: go over group G and project: this is a critical part that *
717 * can have extremely many terms *
718 * U=charge unit for th[i], X=phase unit for (non-GSO) group projection *
719 * QL=sum(th_i-q_i)+Q_ring; dQ=QR-QL=sum(1-2th_i); *
720 * proj: g|h> = (-1)^(Kg*Kh) \e(g,h)(\det g_h)/(det g)|h>, thus *
721 * K=\e=0 => ph[j] is the sum of the phases on invariant fields. *
722 */
Test_BottomUpQuot(PoCoLi * Num,PoCoLi * Den,PoCoLi * Quo,PoCoLi * Rem)723 int Test_BottomUpQuot(PoCoLi *Num,PoCoLi *Den,PoCoLi *Quo,PoCoLi *Rem)
724 { int i=BottomUpQuot(Num,Den,Quo,Rem); if(Rem->n)assert(i==0);
725 if(i)return 1;
726 printf("Num="); PrintPoCoLi(Num);
727 printf("Den="); PrintPoCoLi(Den); printf("Quo=");PrintPoCoLi(Quo);
728 printf("Rem="); PrintPoCoLi(Rem);
729 /* printf("N=%d D=%d Q=%d R=%d\n",Num.A,Den.A,Quo.A,Rem.A);exit(0);fflush(0);*/
730 { PoCoLi A; A.A=10000;AllocPoCoLi(&A);PolyProd(Den,Quo,&A);
731 Poly_Dif(&A,Num,Den);Poly_Sum(Den,Rem,&A); assert(A.n==0);
732 Free_PoCoLi(&A); puts("BottomUpQuot: Test o.k.");
733 } return 0;
734 }
735
pff(char * c)736 void pff(char *c){puts(c);fflush(0);}
737
738 typedef struct {int X, n, *d, **mt;} /* mt[i][j]=mobius(j,i) */ MobiusData;
739
MakeMobius(MobiusData * M,int X)740 void MakeMobius(MobiusData *M,int X) /* d[i] divisors, mt[i][j] 0<=j<=i<=n */
741 { int i,j=0,n=0,*d;
742 for(i=1;i<X/i;i++) if(X%i==0)n++;
743 M->n=n=2*n+(X==i*i); /* #(Div(X)) */
744 M->d = d = (int *)malloc(((n*(n+3))/2) * sizeof(int) + n * sizeof(int*));
745 assert(d != NULL); M->mt = (int **) &d[(n*(n+3))/2]; M->mt[0]=&(d[n]);
746 for(i=1;i<n;i++)M->mt[i]=&M->mt[i-1][i];
747 for(i=1;i<=X/i;i++) if(X%i==0){d[j]=i;d[n-(++j)]=X/i;}
748 for(i=0;i<n;i++){int k; M->mt[i][i]=1;
749 for(j=i-1;0<=j;j--) {M->mt[i][j]=0; for(k=j+1;k<=i;k++)
750 if(d[i]%d[k]==0) if(d[k]%d[j]==0) M->mt[i][j] -= M->mt[i][k];}}
751 #ifdef PRINT_MOBIUS_FUNCTION
752 printf("Div(%d)=",X);for(i=0;i<n;i++)printf("%d ",d[i]);
753 for(i=0;i<n;i++){printf("m[%d]=",d[i]);
754 for(j=0;j<=i;j++)printf("%d ",M->mt[i][j]);} ;puts("");
755 #endif
756 }
FreeMobius(MobiusData * M)757 void FreeMobius(MobiusData *M){free(M->d);}
758 void Calc_VaHo(Weight *W,VaHo *V);
759 void PoincarePoly(int N, int *w, int d, PoCoLi *P, PoCoLi *Z,PoCoLi *R);
760 void Aux_Phase_Poly(PoCoLi *P,int w,int d, int r, int s, int x);
Index_Trace_Test(VaHo * V,int WI,int T)761 int Index_Trace_Test(VaHo *V,int WI,int T)
762 { int i,j,D=V->D;for(i=0;i<=D;i++)for(j=0;j<=D;j++){T-=V->h[i][j];
763 if((i+j)%2) WI += V->h[i][j]; else WI -= V->h[i][j];}
764 assert((T==0)&&(WI==0)); return 1;
765 }
766 /* Thru 4-folds it is sufficient to know boundary (i.e. H{0i} and H{di} *) *
767 * contributions, Witten index and trace for the twisted sectors; the *
768 * untwisted sector can be reconstructed from "WIndex_HTrace(W,&WI,&T)". *
769 * goint twice over the group is not too costly if the projection is reduced *
770 * to the effectively acting subgroup for each twisted sector */
LGO_VaHo(Weight * W,VaHo * V)771 void LGO_VaHo(Weight *W,VaHo *V)
772 { int i,d=W->d,D=0; for(i=0;i<W->N;i++) D+=d-2*W->w[i];assert(!(D%d));D/=d;
773 /* printf("In LGO_VaHo: W->N: %d W->M: %d D: %d\n", W->N, W->M, D); */
774 /* From palp-2.0 to palp-2.11 the next lines were
775 * if((V->sts)||(D>3)) {if(W->M==0) W->m[0]=1;}
776 * else {if(D<=3) {Fast_c9_VaHo(W,V);return;}
777 * if(W->M==0){Calc_VaHo(W,V);return;}}
778 * but that makes no sense since it's equivalent to
779 * if((V->sts)||(D>3)) {if(W->M==0) W->m[0]=1;}
780 * else {Fast_c9_VaHo(W,V);return;}
781 * therefore: */
782 if ((V->sts == 0) && (W->M == 0)) { /* V->sts ... lg-flag set to 2,
783 W->M ... # of Z_n symmetries */
784 if (D<=3) {Fast_c9_VaHo(W,V);return;}
785 else {Calc_VaHo(W,V);return;} }
786 if(W->M==0) W->m[0]=1;
787 { Long dQ,QL,q; PoCoLi *S,*SO,*SN;/* dQ=sum(1-2th)=QR-QL; QL=sum(th-w) */
788 int j,k,s, J=W->M, M[POLY_Dmax],I[POLY_Dmax],X=W->m[0],x[POLY_Dmax];
789 int WI,T; Long ph[POLY_Dmax], th[W_Nmax], U=X*d/Fgcd(X,d), rd=U/d, G=1;
790 PoCoLi Den,Num,Quo,Rem; MobiusData MX; MakeMobius(&MX,X);
791 Den.A = (1<<W->N); AllocPoCoLi(&Den);
792 Quo.A = X*d*D+2*W->N; AllocPoCoLi(&Quo);
793 Num.A=Rem.A= Quo.A; AllocPoCoLi(&Num); AllocPoCoLi(&Rem);
794 V->D=D; for(i=0;i<=D;i++) for(j=0;j<=D;j++) V->h[i][j]=0;
795 for(j=0;j<J;j++) {G*=(M[j]=W->m[j]); assert(*M%M[j]==0); x[j]=X/M[j];}
796 S = (PoCoLi *) malloc( 2*X*sizeof(PoCoLi) ); assert(S!=NULL);
797 for(k=0;k<W->d;k++) /* Z_d==GSO */
798 { int v,a=Init_Multiloop(W->m,I,&v,&J); do { /*BEGIN gen TWISTS */
799 int n[W_Nmax], N=0,h[POLY_Dmax+1],hn=0;
800 QL=dQ=0; for(j=0;j<J;j++) ph[j]=0; for(i=0;i<=D;i++) h[i]=0;
801 for(i=0;i<W->N;i++) {th[i]=k*W->w[i]*rd;
802 for(j=0;j<J;j++) th[i]+=I[j]*W->z[j][i]*(U/W->m[j]);
803 assert(th[i]>=0); th[i]%=U;
804 if(th[i]) {QL+=th[i]-W->w[i]*rd; dQ+=U-2*th[i];}
805 else{n[N++]=i;for(j=0;j<J;j++)ph[j]+=x[j]*W->z[j][i];}
806 } assert((dQ%U)==0); assert((QL%rd)==0); assert(QL+dQ>=0);
807 dQ/=U; QL/=rd; for(j=0;j<J;j++)ph[j]%=X; /* Q & ph of |h>_gs */
808
809 /* int pntw=0; if(100*k+10*I[0]+I[1]==202)pntw=1; */
810
811 /* U==lcm(X,d)=d*rd -> th[i]/U, X==lcm(M[])=x[j]*M[j] -> ph[j]/X |g> */
812 /* QL/d+dQ G=\prod M[j]=|group| n[i], i<N = untw.fields */
813 /* twist= gso^k*g_j^I[j] r/oz=phase(\Ph_n) s/oz=phase of vacuum */
814 if(N==0) {i=1; for(j=0;j<J;j++) if(ph[j]) i=0;
815 if(i&&(QL%d==0)) hn=h[QL/d]=1; }/* trivial ring sectors */
816 /* s1+l r1 \in X\IZ, g=ar+bX => l=(l1 X - a1 s1)/g1, */
817 /* A g1+B g2=g, if((s1 a1 g2/g - s2 a2 g1/g) % X ==0) then */
818 /* l=l' X/g +(B X-a1 s1)/g1 = l' X/g - (A X +a2 s2)/g2 else break; */
819 /* r'=g, s'=(A X + a2 s2)/g2, a'=1; assert(s' % g' ==0); */
820 if(N==1) {int w=W->w[*n], r=0,g=X; Long A,B; s=0;
821 for(j=0;j<J;j++) { /* cyclic projection */
822 int rj=(W->z[j][*n]*x[j])%X, sj=ph[j]%X; assert(0<=rj);
823 if(rj==0) { if(sj==0) continue; else {g=0;break;} }
824 g=Egcd(rj,X,&A,&B); sj=(sj*A)%X;
825 assert(rj>0);assert(g>0);assert(sj%g==0);
826 if(sj%g) {g=0;break;}
827 rj=g; g=Egcd(r,rj,&A,&B); assert(g>0);
828 if((s*rj-sj*r)% (g*X)) {g=0;break;}
829 r=g; s=(A*X+sj)/rj; if(s%g) {g=0;break;}
830 } assert((j<J)==(g==0));
831 if(g==0) Quo.n=0; else if(g==X) Init1_xN(&Num,d-w);
832 else {/* assert(g>1); */ Aux_Phase_Poly(&Num,w,d,r,s,X/g);}
833 if(g) {Init1_xN(&Den,w*(X/g));
834 BottomUpQuot(&Num,&Den,&Quo,&Rem);
835 for(i=0;i<Quo.n;i++) if(( q=(QL+Quo.e[i]) )%d==0)
836 {h[q/d]+=Quo.c[i]; hn=1;}
837 }}
838 if(N>1) /* N>1: S[] sector -> SN[] new sector PoCoLi */
839 { int u,b,l,L[POLY_Dmax],w[W_Nmax],o[W_Nmax],io[W_Nmax],mm=1;
840 for(i=0;i<N;i++) for(j=i+1;j<N;j++)/* sort invariant weights */
841 if(W->w[n[j]]<W->w[n[i]]) swap(&n[i],&n[j]);
842 for(i=0;i<N;i++) {w[i]=W->w[n[i]]; o[i]=1;
843 for(j=0;j<J;j++){ u=NNgcd(W->z[j][n[i]],W->m[j]); u=W->m[j]/u;
844 o[i]=u*o[i]/Fgcd(o[i],u);} mm=mm*o[i]/Fgcd(mm,o[i]);}
845 /* phase(g.s.) = gs/mm; mm=lcm(o[i]); phase(X(n(i))=z[i]/o[i]; */
846
847 /* BEGIN GROUP PROJECTION */
848 if(mm>1){int *mo, m, ego=1; /* effective group order */
849 PoCoLi *A,*B,*C,Ax; Ax.A=Quo.A; AllocPoCoLi(&Ax); Num.n=0;
850 for(j=0;j<J;j++){int g=W->z[j][n[0]];
851 for(i=1;i<N;i++) g=NNgcd(g,W->z[j][n[i]]);
852 M[j] = W->m[j] / NNgcd(g,W->m[j]);}
853 for(m=0;m<MX.n;m++)if(mm<=MX.d[m])break;assert(mm==MX.d[m]);
854 mo=MX.mt[m]; l=0; for(j=0;j<m;j++)
855 if(mo[j]==1)l++; else if(mo[j]==-1)l++; else assert(mo[j]==0);
856
857 b=ego=Init_Multiloop(M,L,&u,&J); do { /* BEGIN make group */
858 int gs=0,z[W_Nmax];for(i=0;i<N;i++){z[i]=0; /* g.s. phases */
859 for(j=0;j<J;j++) z[i]+=L[j]*((W->z[j][n[i]]*o[i])/W->m[j]);
860 io[i]=mm/o[i]; gs+=io[i]*z[i];} gs%=mm; assert(gs>=0);
861 SO=S; SN=&S[mm]; /* for each projection group element */
862 for(s=0;s<mm;s++){int ds=gs+s; S[s].A=2*mm; /* init SO */
863 AllocPoCoLi(&S[s]); if( ds % io[0] ) S[s].n=0;
864 else if(*o>1)Aux_Phase_Poly(&S[s],w[0],d,z[0],-ds/io[0],
865 o[0]); else Init1_xN(&S[s],d-w[0]);}
866 for(i=1;i<N;i++){PoCoLi *Saux=SO,*ac; int t; for(t=0;t<mm;t++){
867 A=&Ax; B=&Quo; C=&Rem; C->n=0;
868 for(s=0;s<mm;s++)if(s%io[i]==0){
869 if(o[i]>1) Aux_Phase_Poly(A,w[i],d,z[i],-s/io[i],o[i]);
870 else Init1_xN(A,d-w[i]);
871 PolyProd(&SO[(mm+t-s)%mm],A,B); Poly_Sum(B,C,A);
872 ac=A;A=C;C=ac;}
873 SN[t].A=C->A; AllocPoCoLi(&SN[t]); PolyCopy(C,&SN[t]);
874 } for(t=0;t<mm;t++) Free_PoCoLi(&SO[t]);
875 SO=SN; SN=Saux;
876 } /* S[s] finished */
877
878 /* if(pntw)for(s=0;s<mm;s++){
879 printf("gs=%d S[%d]= ",gs,s);PrintPoCoLi(&SO[s]);} */
880
881 Poly_Sum(&SO[0],&Num,&Rem); B=&Num; A=&Rem;
882 for(j=0;j<m;j++) if(mo[j]==1) {Poly_Sum(A,&SO[MX.d[j]],B);
883 C=A;A=B;B=C;} else if(mo[j]==-1)
884 {Poly_Dif(A,&SO[MX.d[j]],B); C=A;A=B;B=C;}
885 assert((A==&Num)==(l%2)); if(l%2==0) PolyCopy(A,&Num);
886 for(s=0;s<mm;s++){Free_PoCoLi(&SO[s]);}
887 assert(0<(b--));} while(Multiloop(M,L,&u,&J)); assert(b==0);
888
889 for(i=0;i<Num.n;i++){assert(0==Num.c[i]%ego); Num.c[i]/=ego;}
890
891 /* if(pntw){printf("Num[%d:%d,%d]=",k,I[0],I[1]);
892 PrintPoCoLi(&Num);fflush(0);} */
893
894 A=&Num; B=&Quo;
895 if(Num.n) for(i=0;i<N;i++){Init1_xN(&Den,w[i]*o[i]);
896 BottomUpQuot(A,&Den,B,&Rem); C=A;A=B;B=C;}
897 for(i=0;i<A->n;i++) if(( q=(QL+A->e[i]) )%d==0)
898 {h[q/d]+=A->c[i]; hn=1;}
899
900 /* if(pntw){printf("Quo[%d:%d,%d]=",k,I[0],I[1]);
901 PrintPoCoLi(A);fflush(0);} */
902
903 Free_PoCoLi(&Ax); /* END GROUP PROJECTION */
904 }
905 else
906 { PoincarePoly(N,w,d,&Quo,&Num,&Den);
907 for(i=0;i<Quo.n;i++) if(( q=(QL+Quo.e[i]) )%d==0)
908 {h[q/d]+=Quo.c[i]; hn=1;} /*END of no projection */
909 }
910 } /* N>1 CASE FINISHED */
911
912 if(hn) for(i=0;i<=D;i++) if(h[i]) V->h[D-i][i+dQ]+=h[i];
913 if(V->sts) if(hn){fprintf(outFILE,"sec[%d",k);
914 for(j=0;j<J;j++) fprintf(outFILE,"%s%d",j?",":":",I[j]);
915 fputs("]",outFILE);
916 fprintf(outFILE," th=%2ld",th[0]);for(i=1;i<W->N;i++)
917 fprintf(outFILE," %2ld",th[i]); /*fprintf(outFILE,"/%d ",U);*/
918 /*fprintf(outFILE," %d*ph=",U);for(i=0;i<J;i++)printf("%d ",ph[i]); */
919 fprintf(outFILE," QL=%2ld/%d dQ=%2ld ",QL,d,dQ);
920 /*fprintf(outFILE,"N=%d ",N);*/
921 for(i=0;i<=D;i++) if(h[i]) fprintf(outFILE," q%d%ld+=%d",
922 i,i+dQ,h[i]); fputs("\n",outFILE);}
923 /* if(!cont)exit(0); */
924 assert(0<(a--));
925 } while(Multiloop(W->m,I,&v,&J));assert(a==0); /* END gen TWISTS */
926 } free(S); FreeMobius(&MX); Free_PoCoLi(&Rem); Free_PoCoLi(&Quo);
927 assert(Hodge_Test(V)); Free_PoCoLi(&Den); Free_PoCoLi(&Num);
928 assert(1==WIndex_HTrace(W,&WI,&T)); assert(Index_Trace_Test(V,WI,T));
929 if(V->sts) printf("WittenIndex=%d, Trace=%d\n",WI,T);
930 }
931 }
932 /* T=\x*t: (1-T^(d-w))/(1-T^w)=(1-T^(d-w))*(1+T^w+T^2w+...+T^(x-1)w)/(1-T^wx)*
933 * phase(\x)=r/x. All terms in numerator with phase s/x, p=0,...,x-1 => *
934 * exp=p mod x in numerator 1+...+T^(x-1)w-T^(d-w)-T^d-...-T^(d+w(x-2)) */
Aux_Phase_Poly(PoCoLi * P,int w,int d,int r,int s,int x)935 void Aux_Phase_Poly(PoCoLi *P,int w,int d, int r, int s, int x)
936 { int i, nocancel = (d%w) || (x*w<d) || ((r*d/w)%x),
937 negmin = nocancel ? d-w : x*w, negmax=d+x*w-2*w,
938 posmax = nocancel ? (x-1)*w : d-2*w; assert(x>1); P->n=0;
939 for(i=0;i<=posmax;i+=w) if((r*(i/w)-s) % x == 0) Add_Mono_2_Poly(i,1,P);
940 for(i=negmin;i<=negmax;i+=w)if((r*(i-d)/w-s)%x==0)Add_Mono_2_Poly(i,-1,P);
941 }
942
Calc_VaHo(Weight * W,VaHo * V)943 void Calc_VaHo(Weight *W,VaHo *V) {
944 int i, j, k, D=0, w[W_Nmax], N=W->N, d=W->d;
945 PoCoLi A,B,C, *Z=&A, *R=&B, *P=&C;
946 for(i=0;i<N;i++) { w[i]=W->w[i]; D+=d-2*w[i]; }
947 if(D % d) {V->D=0; return;} else V->D = (D/=d);
948 for(i=0;i<N;i++)for(j=i+1;j<N;j++)if(w[j]<w[i]) swap(&w[i],&w[j]);
949 #ifndef COEFF_Nmax
950 #define COEFF_Nmax (d*D+2*N) /* 22999000 for W_Nmax == 6 */
951 #endif
952 P->A=Z->A=COEFF_Nmax; R->A=w[N-1]+1; //printf("%d %d %d", P->A, Z->A, R->A);
953 AllocPoCoLi(P); AllocPoCoLi(Z); AllocPoCoLi(R);
954 PoincarePoly(N,w,d,P,Z,R);
955 assert(D*d==P->e[P->n-1]);
956 { int n=P->n, cM=0; long long sum=0,num=1,den=1; /* check sum = P(0) */
957 for(i=0;i<n;i++) {
958 Pint co=P->c[i]; assert(co>0); if(co>cM) cM=co;
959 sum+=co; }
960 for(i=0;i<W->N;i++){ num*=W->d-W->w[i]; den*=W->w[i];
961 {long long g=LFgcd(num,den); num/=g; den/=g; }}
962 #ifdef TEST_PP
963 if(P->n<99) {printf("PP =");PrintPoCoLi(P);} else printf(
964 "#(Exp,Co)=%d Exp<=%d Coeff<=%d sum=%lld\n",n,P->e[n-1],cM,sum);
965 #endif
966 if(den==1) assert(num==sum);
967 else { printf("sum=%lld test=%lld/%lld \n",sum,num,den);}
968 }
969 for(i=0;i<=D;i++) for(j=0;j<=D;j++) V->h[i][j] = 0;
970 for(i=0;i<P->n;i++)if(P->e[i] % d==0) {j=P->e[i]/d; V->h[D-j][j]=P->c[i];}
971 #ifdef TEST_PP
972 for(i=0;i<=D;i++){for(j=0;j<=D;j++)printf("%6d ",(int)V->h[i][j]);
973 puts("= V-pri");}fflush(0);
974 #endif
975 for(k=1;k<d;k++) /* k-twisted sector */
976 { int th[W_Nmax], Tmt, et=0, eT=0, Deff=0, n=0;
977 for(i=0;i<N;i++)
978 { if((th[i]=((long long) k * W->w[i]) % (long long)d))
979 { et+=th[i]-W->w[i];eT+=W->d-th[i]-W->w[i];}
980 else { Deff+=W->d-2*W->w[i]; w[n++]=W->w[i];} /* effective w */
981 }
982 if(0 == ((Tmt=eT-et) % d))
983 { for(i=0;i<n;i++)for(j=i+1;j<n;j++)if(w[j]<w[i]) swap(&w[j],&w[j]);
984 PoincarePoly(n,w,d,P,Z,R); assert(Deff==P->e[P->n-1]); Tmt/=d;
985 for(i=0;i<P->n;i++)
986 { j=et+P->e[i];
987 if(j % d==0) {j/=d; if(j+Tmt>=0)V->h[D-j][j+Tmt]+=P->c[i];}
988 }
989 }
990 } Free_PoCoLi(P); Free_PoCoLi(Z); Free_PoCoLi(R);
991 #ifdef TEST_PP
992 for(i=0;i<=D;i++){for(j=0;j<=D;j++)printf("%6d ",V->h[i][j]);puts("= V");}
993 #endif
994 }
995
996
997 /* ============= Trans_Check(Weight W) ========== */
998 /* OLD version ignores
999 * orbifold projection */
1000
OLDsymcheck(int sum,int link,int * mon,Weight W,int * mask)1001 int OLDsymcheck(int sum, int link, int *mon, Weight W, int *mask){
1002 /* symcheck returns 1 if there is a monomial mon in
1003 * the variables indicated by link whose total weight is sum and which
1004 * transforms under the k'th symmetry with a phase sum[k];
1005 * if X_i occurs in the monomial then mon[i] is set to 1. */
1006 int i, j, newsum, newlink;
1007
1008 for (i=0;i<W.N;i++) if (link&mask[i]){
1009 if (!(sum % W.w[i])){
1010 *mon=*mon|mask[i];
1011 return 1;};
1012 newlink=link-mask[i];
1013 for (j=0;j*W.w[i]<sum;j++){
1014 newsum=sum-j*W.w[i];
1015 if (OLDsymcheck(newsum,newlink,mon,W,mask)) {
1016 if (j) *mon=*mon|mask[i];
1017 return 1; }; }; };
1018 return 0;
1019 }
OLDInit_Trans_Check(int * mask,int ** targets,int ** mighty)1020 void OLDInit_Trans_Check(int *mask, int **targets, int **mighty)
1021 { int i, j, maxPN=1<<W_Nmax; mask[0]=1; /* maxPN=2^W_Nmax */
1022 for(i=1;i<=W_Nmax;i++)mask[i]=2*mask[i-1]; /* mask={1,2,4,8,16,32,...} */
1023
1024 assert(maxPN==mask[W_Nmax]);/* maximum number of pointers at one point */
1025 assert( NULL != (*targets = (int *) malloc(maxPN * sizeof(int))) );
1026 assert( NULL != (*mighty = (int *) malloc(maxPN * sizeof(int))) );
1027
1028 for (i=0;i<W_Nmax;i++) for (j=mask[i];j<mask[i+1];j++)/* make mighty[] */
1029 (*mighty)[j]=(*mighty)[j-mask[i]]+1;
1030 }
1031
OLDTrans_Check(Weight W)1032 int OLDTrans_Check(Weight W){/* returns 1 if non-degenerate potential exists */
1033 int i, j, k, mon; /* j represents the link!!! */
1034 static int mask[1+W_Nmax], *targets, *mighty; assert(W.N<=W_Nmax);
1035 if(targets==NULL) OLDInit_Trans_Check(mask,&targets,&mighty);
1036 targets[0]=0;
1037 for (i=0;i<W.N;i++) for (j=mask[i];j<mask[i+1];j++){
1038 targets[j]=0;
1039 for (k=0;k<=i;k++) targets[j]=targets[j]|targets[j&~mask[k]];
1040 for (k=0;(k<W.N)&&(mighty[targets[j]]<mighty[j]);k++)
1041 if (!(mask[k]&targets[j])) if (OLDsymcheck(W.d-W.w[k],j,&mon,W,mask))
1042 targets[j]=targets[j]|mask[k];
1043 if (mighty[targets[j]]<mighty[j]) return 0;}
1044 return 1;
1045 }
1046
1047 typedef struct {int d, m[POLY_Dmax];} symlist;
1048
symcheck(symlist sum,int link,Weight W,int * mask)1049 int symcheck(symlist sum, int link, Weight W, int *mask){
1050 /* symcheck returns 1 if there is a monomial in
1051 * the variables indicated by link whose total weight is sum.d and which
1052 * transforms under the k'th symmetry with a phase sum.m[k]; */
1053
1054 int i, j, k, check, expo, newlink;
1055 symlist newsum;
1056
1057 for (i=0;i<W.N;i++) if (link&mask[i]){
1058 if (!(sum.d % W.w[i])){
1059 expo=sum.d/W.w[i];
1060 check=1;
1061 for (j=0;(j<W.M)&✓j++)
1062 if ((expo*W.z[j][i]-sum.m[j])% W.m[j]) check=0;
1063 if (check) return 1;}
1064 newlink=link-mask[i];
1065 for (j=0;j*W.w[i]<sum.d;j++){
1066 newsum.d=sum.d-j*W.w[i];
1067 for (k=0;k<W.M;k++) newsum.m[k]=sum.m[k]-j*W.z[k][i];
1068 if (symcheck(newsum,newlink,W,mask)) return 1; } }
1069 return 0;
1070 }
1071
Init_Trans_Check(int * mask,int ** targets,int ** mighty)1072 void Init_Trans_Check(int *mask, int **targets, int **mighty)
1073 { int i, j, maxPN=1<<W_Nmax; mask[0]=1; /* maxPN=2^W_Nmax */
1074 for(i=1;i<=W_Nmax;i++)mask[i]=2*mask[i-1]; /* mask={1,2,4,8,16,32,...} */
1075
1076 assert(maxPN==mask[W_Nmax]);/* maximum number of pointers at one point */
1077 assert( NULL != (*targets = (int *) malloc(maxPN * sizeof(int))) );
1078 assert( NULL != (*mighty = (int *) malloc(maxPN * sizeof(int))) );
1079
1080 for (i=0;i<W_Nmax;i++) for (j=mask[i];j<mask[i+1];j++)/* make mighty[] */
1081 (*mighty)[j]=(*mighty)[j-mask[i]]+1;
1082 }
1083
1084 /* targets[j] bin-encodes the gradients that can be nonzero with variables
1085 bin-encoded in j; transversality <==> mighty[j] <= mighty[targets[j]] */
1086
Trans_Check(Weight W)1087 int Trans_Check(Weight W){ /* returns 1 if non-degenerate potential exists */
1088 int i, j, k, l; /* j represents the link!!! */
1089 symlist dw;
1090 static int mask[1+W_Nmax], *targets, *mighty;
1091 assert(W.N<=W_Nmax);
1092 if(targets==NULL) Init_Trans_Check(mask,&targets,&mighty);
1093 targets[0]=0;
1094 for (i=0;i<W.N;i++) for (j=mask[i];j<mask[i+1];j++){
1095 targets[j]=0; /* union of targets with one variable less -> */
1096 for (k=0;k<=i;k++) targets[j]=targets[j]|targets[j&~mask[k]];
1097 for (k=0;(k<W.N)&&(mighty[targets[j]]<mighty[j]);k++)
1098 if (!(mask[k]&targets[j])) { /* targets=lower bounds only! */
1099 dw.d=W.d-W.w[k]; /* gradient degrees */
1100 for (l=0;l<W.M;l++) dw.m[l]=W.m[l]-W.z[l][k]; /* gradient phases */
1101 if (symcheck(dw,j,W,mask)) targets[j]=targets[j]|mask[k];}
1102 if (mighty[targets[j]]<mighty[j]) return 0;}
1103 return 1;
1104 }
1105