global(DEBUG:small,GTH1:small,GTH2:small,GTH3:small,GTHMOD:small,SIZ:small,LSIZ,GOOD:small,ONE:small,ZERO:small,INFINIT:small,BAD:int,DEG:small,HODGE,TT:small,WT:small,DISC:int,POLGA,POLGAONE:small,POLDEG:small,VAL,VBE,LFUN,LKEEP,GF,GH,PSQ:int,LIMTAB:small,GPV,TABGAM,TABGAM3,INDP,OFFSETL0:small,LIMMODP2:small,VPOLGA,COUNTP:small,VECA,VECC);
global(PF:small,PFM1:small,VECCPN,VL1,VL0,WPOL,PFL0LIM:small,PRIME:small);
global(DEG2:small,VALNUM,VALDEN,VBENUM,VBEDEN,DIRONE,OFFMPOL:small,PFP1OV2:small);
global(ESA,ESB);
global(w3r4,w3r6);

default(realprecision,38);

/*\r tabgen.gp */
/* TABGAM=readvec("tabgam.gp");
TABGAM3=readvec("tabgam3.gp"); */
/* read("indlim.gp"); */
TABGAM=eval("TABGAM");
TABGAM3=eval("TABGAM3");
INDP=eval("INDP");
LIMTAB=eval("LIMTAB"):small;

\r hyperg-list.txt
w3r4=eval("w3r4");
w3r6=eval("w3r6");

DEBUG=0;

ZERO=0;ONE=1;INFINIT=1000;GOOD=7;

/* Conversion back and forth al,be to ga */

conv(val)=
{
 my(den:small,z,res);
 den=denominator(val):small;
 z=Mod('y,'y^den-1);
 res=lift(prod(i=1,#val,'x-z^(den*val[i])));
 res=lift(subst(res,'y,Mod('y,polcyclo(den,'y))));
 return (simplify(res));
}

checkalbe()=
{
 my(resal,resbe);
 resal=conv(VAL); if (deriv(resal,y)!=0,return(0));
 resbe=conv(VBE); if (deriv(resbe,y)!=0,return(0));
 if (gcd(resal,resbe)!=1,return(0));
 return ([resal,resbe]);
}

/* Computes POLGA */

getga()=
{
 my(resalbe,resal,resbe,maxden,f);
 resalbe=checkalbe(); if (resalbe==0,error("not Q polynomials"));
 resal=resalbe[1]; resbe=resalbe[2];
 maxden=max(vecmax(vector(DEG,i,denominator(VAL[i]))),vecmax(vector(DEG,i,denominator(VBE[i]))))+1;
 f=-log((polrecip(resal)+O(x^maxden))/(polrecip(resbe)+O(x^maxden)));
 return (sum (n=1,maxden-1,sumdiv(n,d,polcoeff(f,d)*moebius(n/d)*d)/n*x^n));
}

/* should be in GP */

isinvec(v,a)=
{
 for (i=1,#v,if(v[i]==a,return(i)));
 return (0);
}

getalbe(polga)=
{
 my (vtop,vbot,co,w,ct,fl,val,vbe);
 vtop=[]; vbot=[];
 for (nu=1,poldegree(polga),
  co=polcoeff(polga,nu);
  if (co>0,
   w=vector(nu,i,i-1)/nu; for (j=1,co,vtop=concat(vtop,w))
  );
  if (co<0,
   w=vector(nu,i,i-1)/nu; for (j=1,-co,vbot=concat(vbot,w))
  )
 );
 ct=0;
 for (n=1,#vbot,
  fl=isinvec(vtop,vbot[n]);
  if (fl,vtop[fl]=-1;vbot[n]=-1;ct++)
 );
 vtop=vecsort(vtop); vbot=vecsort(vbot);
 val=vector(#vtop-ct,i,vtop[i+ct]); vbe=vector(#vbot-ct,i,vbot[i+ct]);
 return ([val,vbe]);
}

call(polga,x)=
{
 my(s,co);
 s=0;
 for (nu=1,#polga,
  co=polcoeff(polga,nu);
  s+=co*frac(nu*x)
 );
 return (-s);
}

/* computes HODGE */

hodge()=
{	
 my(sw,w,r:small,s:small);
 r=DEG;
 w=vector(r,k,[VAL[k],0]);
 w=concat(w,vector(r,k,[1-VBE[k],1]));
 w=vecsort(w,1);
 sw=0;
 s=0;
 for(k=1,2*r,
  if(w[k][2] == 0, sw+=x^s; s++, s--)
 );
 return (sw);
}

compD()=
{
 my(D:int,ct:small);
 POLGA=getga(); POLGAONE=subst(POLGA,x,1):small; POLDEG=poldegree(POLGA);
 VPOLGA=vector(POLDEG,j,polcoeff(POLGA,j));
 D=prod(i=1,POLDEG,i^abs(polcoeff(POLGA,i))):int;
 ct=0;for(i=1,DEG,if(VAL[i]==0||VBE[i]==0,ct++));
/* if (ct%2,error("odd ct in compD")); */
 D=((-1)^(ct\2)*D):int;
 return (coredisc(D));
}

findeme(v)=
{
  my(d,den,es,e_count,k,fl);
  d=#v;
  den=vector(d,i,denominator(v[i]));
  e_count=1;
  for(i=2,d,fl=0;for(j=1,i-1,if(den[i]==den[j],fl=1));if(fl==0,e_count=e_count+1));
  es=matrix(e_count,2);
  es[1,1]=den[1];
  k=1;	
  for(i=2,d,fl=0;for(j=1,i-1,if(den[i]==den[j],fl=1));if(fl==0,k=k+1;es[k,1]=den[i]));
  for(i=1,e_count,for(j=1,d,if(es[i,1]==den[j],es[i,2]=es[i,2]+1)));
  for(i=1,e_count,es[i,2]=es[i,2]/eulerphi(es[i,1]));
  return(es);
}

/* here use ESA or ESB  and k = +-v_p(t) */
/* Returns [v_p(N),[degs of Euler factors],[weights of Euler factors] */

findrn(k)=
{
 my(n,ve,vme,vr,e,s,ws,ns,su);
 if (k>0,ve=ESA[,1];vme=ESA[,2],ve=ESB[,1];vme=ESB[,2]);
 n=#ve;	
 vr=vector(vecmax(vme));
 for(i=1,n,
   e=ve[i];
   if(k%e==0,vr[vme[i]]+=eulerphi(ve[i]));
 );
 s=0;for(i=1,#vr,if(vr[i]>0,s++));			
 ws=vector(s);ns=vector(s);
 c=0; su=0;
 for(i=1,#vr,if(vr[i]>0,c++;ws[c]=WT+1-i;ns[c]=vr[i];su+=ns[c]));
 return([DEG-su,ns,ws]);
}

/* Fundamental initialization of variety */

initalbe(val,vbe)=
{
 VAL=val;VBE=vbe;
 DEG=#VAL; if (#VBE!=DEG,error("length error")); DEG2=DEG\2;
 ESA=findeme(VAL); ESB=findeme(VBE);
 VALNUM=vector(DEG,j,numerator(VAL[j]));
 VALDEN=vector(DEG,j,denominator(VAL[j]));
 VBENUM=vector(DEG,j,numerator(VBE[j]));
 VBEDEN=vector(DEG,j,denominator(VBE[j]));
 BAD=-lcm(denominator(VAL),denominator(VBE)):int;
 OFFSETL0=0;
 for(j=1,DEG,
  if(VBE[j]==0,OFFSETL0++);
  if(VAL[j]==0,OFFSETL0++)
 );
 HODGE=hodge();
 TT=poldegree(denominator(HODGE));
 WT=poldegree(numerator(HODGE));
/* if (WT%2==0 && DEG%2==0,error("weight even when degree even")); */
/* if (WT%2==0,return(0)); */
 if (WT > 3,error("weight > 3 not yet implemented"));
 LIMMODP2=((2*DEG)^2):small;
 DISC=compD():int;
 OFFMPOL=(OFFSETL0-POLGAONE)\2;
}

/* recognize a p-adic number as an integer */

convpadic(n)=
{
 my(res:int,res2:int,v:small);
 if (type(n)=="t_VEC",return (vector(#n,i,convpadic(n[i]))));
 if (n==0,return(0));
 v=valuation(n,PRIME);n/=PRIME^v;
 res=truncate(n):int; res2=truncate(-n):int;
 if (res<res2,return (PRIME^v*res),return(-PRIME^v*res2));
}

/* Programs to compute gamma_p */

/* compute p^ku_{pk+r} for 0\le r\le p-1 and 0\le k< ell */

compu(ell:small,dfp:small)=
{
 my(v,vall,point:small,oldp:small);
 vall=vector(ell+1);
 v=vector(PRIME);v[1]=1+O(PRIME^dfp);
 for(i=1,PRIME-1,v[i+1]=v[i]/i);
 vall[1]=v; point=1; oldp=PRIME;
 for(k:small=1,ell*PRIME,
  v[point]=(v[oldp]+v[point])/(k+PRIME-1);
  if(oldp==PRIME,oldp=1,oldp++);
  if(point==PRIME,
   point=1; vall[k/PRIME+1]=v,
   point++
  )
 );
 return(vall);
}

/* compute all gamma(m/(p^f-1)) for 0\le m<p^f-1 */

precomp(f:small,dfp:small)=
{
 my(vall,vres,x,v:small,vva:small,vpowv:small,vpo:small,q:int,l2,ga12,m0:small,m1:small,mp,tmp,Q:int,ms2:small,mppro:small,papp:small,zerp,tmp2,pf:int,vtmp,vtmp2);
 if(DEBUG>=2,gettime());
 papp=ceil(dfp/(1-1/PRIME)):small;
 vall=compu(papp,dfp);
 zerp=O(PRIME^dfp);
 if (PRIME==2,
  if (f==1,return([1+zerp]));
  vres=vector(2^f-1);
  for (m=0,2^f-2,
   m0=(m%2):small; m1=((m-m0)\PRIME):small;
   x=(m1+m0*(2^(f-1)))/(2^f-1)+zerp;
   vres[m+1]=sum(k:small=0,papp,binomial(x+k-1,k)*k!*PRIME^k*vall[k+1][m0+1])
  );
  return(vres)
 );
 pf=(PRIME^f):int; Q=pf-1; q=Q\2;
 vva=valuation(Q,2); vpowv=(2^vva):small;
 vtmp=vector((q+1)>>1); vtmp2=vector((q+1)>>1);
 forstep (m:small=1,q,2,
  if (f>1,
   m0=(m%PRIME):small; m1=((m-m0)\PRIME):small;
   x=(m1+m0*(pf/PRIME))/Q+zerp,
   m0=m; x=m/Q+zerp
  );
  tmp2=sum(k:small=0,papp,binomial(x+k-1,k)*k!*PRIME^k*vall[k+1][m0+1]);
  vtmp[(m+1)>>1]=tmp2; vtmp2[(m+1)>>1]=(-1)^((-m-1)%PRIME+1)/tmp2
 );
 vres=vector(Q); vres[1]=1+zerp;
 forstep (m:small=1,q,2,
  vres[m+1]=vtmp[(m+1)>>1]; vres[pf-m]=vtmp2[(m+1)>>1]
 );
 if(DEBUG>=2,print("end odd, begin vva, time = ",gettime()));
 vtmp=vector(((q>>vva)+1)>>1); vtmp2=vector(((q>>vva)+1)>>1);
 forstep (m:small=vpowv,q,2*vpowv,
  if (f>1,
   m0=(m%PRIME):small; m1=((m-m0)\PRIME):small;
   x=(m1+m0*(pf/PRIME))/Q+zerp,
   m0=m; x=m/Q+zerp
  );
  tmp2=sum(k:small=0,papp,binomial(x+k-1,k)*k!*PRIME^k*vall[k+1][m0+1]);
  vtmp[((m>>vva)+1)>>1]=tmp2; vtmp2[((m>>vva)+1)>>1]=(-1)^((-m-1)%PRIME+1)/tmp2
 );
 forstep (m:small=vpowv,q,2*vpowv,
  vres[m+1]=vtmp[((m>>vva)+1)>>1]; vres[pf-m]=vtmp2[((m>>vva)+1)>>1]
 );
 v=1; vpo=2; ga12=1/gamma(1/2+zerp); l2=(PRIME-1)*log(2+zerp);
 if(DEBUG>=2,print("end vva, begin vpo, time = ",gettime()));
 while (vpo<Q,
  if (v!=vva,
   vtmp=vector(((q>>v)+1)>>1); vtmp2=vector(((q>>v)+1)>>1);
   forstep (m:small=vpo,q,2*vpo,
    ms2=((m>>1)+1):small;
    tmp=vres[ms2]*vres[ms2+q]*ga12;
    m0=((-m-1)%PRIME):small; 
    mppro=((m-(m0+1)*Q)\PRIME):small; mp=(mppro+zerp)/Q;
    tmp2=2^m0*exp(mp*l2)*tmp;
    vtmp[((m>>v)+1)>>1]=tmp2; vtmp2[((m>>v)+1)>>1]=(-1)^(m0+1)/tmp2
   );
   forstep (m:small=vpo,q,2*vpo,
    vres[m+1]=vtmp[((m>>v)+1)>>1]; vres[pf-m]=vtmp2[((m>>v)+1)>>1]
   );
  );
  if(DEBUG>=2,if (v<10,print(v," done, time = ",gettime())));
  v++; vpo*=2
 );
 if(DEBUG>=2,print("end precomp, time = ",gettime()));
 return (vres);
}

/* Same for f=1 and mod p^2: compute all gamma_p(m/(p-1)) mod p^2 for 
0\le m<p-1, p>2 */

doprecompmodp2()=initgamp();GPV=vector(PRIME-1,m,gap2simp(m-1));

/* fill global var GPV with gamma_p(m/(p^f-1)) */

doprecomp(f:small,dfp:small)=
{
 my(tmp,pad:small);
 if (PRIME==2,pad=16,if(PRIME==3,pad=11,if(PRIME==5,pad=8,pad=6)));
 if (PRIME>LIMTAB || dfp>pad || (f==3 && PRIME>23),GPV=precomp(f,dfp);return());
 if (f==3,GPV=(1+O(PRIME^dfp))*TABGAM3[INDP[PRIME]]; return());
 if (f==2,
  GPV=(1+O(PRIME^dfp))*TABGAM[INDP[PRIME]]; return(),
  tmp=TABGAM[INDP[PRIME]];
  GPV=(1+O(PRIME^dfp))*vector(PRIME-1,m,tmp[(m-1)*(PRIME+1)+1]); return()
 );
}

/* Assume GPV=vector(p^f-1,m,gamma_p((m-1)/(p^f-1))), read directly from
   TABGAM[INDP[p]] if p<=100, computed by precomp if p>100 */

/* P=c(p,N)*gamma_p(N*s)/N^{[N*s-1-(N-1)\p]} with s=m/(p^f-1), c(p,N)
   given in Theorem 11.6.14 */

/* compute c(p,N) */

gapcpn(f:small,N:small)=
{
 my(N2:small,si);
 if (N%2,return(kronecker(-PRIME,N)));
 if (PRIME==2,error("p=2 and N even in gapcpn"));
 N2=N\2; si=(-1)^(N2+1); 
 return (si*kronecker(si*N2,PRIME)*GPV[PFP1OV2]);
}

dopregapcpn(f:small)=
{
 VECCPN=vector(POLDEG);
 for (N:small=2,POLDEG,if (VPOLGA[N],VECCPN[N]=gapcpn(f,N)));
}
  

/* compute gamma_p(frac(N*m*p^e/(p^f-1))), 0 <= m < p^f-1, 0 <= e < f */

gapgap(f:small,m:int,N:small,e:small)=
{
 return (GPV[(N*m*PRIME^e)%PFM1+1]);
}

/* compute N^{-[N*s-1-(N*s-1)\p]}, N*s=a/(p^f-1) with a=N*m*p^e */

gapnpow(f:small,m:int,N:small,e:small,flag=0)=
{
 my(Nm:int,r0:small,r1:int,r2:int,pad:small);
 if (flag,pad=16,if (PRIME==2,pad=16,if(PRIME==3,pad=11,if(PRIME==5,pad=8,pad=6))));
 if (f==3,
  Nm=((-N*m*PRIME^e)%PFM1):int; r0=Nm%PRIME; 
  Nm=(Nm-r0)\PRIME; r1=Nm%PRIME; r2=(Nm-r1)\PRIME;
  return (((N/teichmuller(N+O(PRIME^pad)))^((r1+r2*PRIME+r0*PRIME^2)/(PRIME^2+PRIME+1)))/N^r0)
 );
 if (f==2,
  Nm=((-N*m*PRIME^e)%PFM1):int; r0=Nm%PRIME; r1=(Nm-r0)\PRIME;
  return (((N/teichmuller(N+O(PRIME^pad)))^((r0*PRIME+r1)/(PRIME+1)))/N^r0),
  Nm=((N*m*PRIME^e)%PFM1):int;
  return (teichmuller(N+O(PRIME^pad))^Nm)
 );
}

/* compute $P=\prod_{0<=j<N}\gamma_p(frac(p^e*(m/(p^f-1)+j/N)))$ */

gapall(f:small,m:int,N:small,e:small,flag=0)=
{
 if (N==1,
  return (GPV[(PRIME^e*m)%PFM1+1]),
  return(VECCPN[N]*gapgap(f,m,N,e)*gapnpow(f,m,N,e,flag))
 );
}

/* no need to reduce */

gapallmodp2(m:small,N:small)=
{
 if (N==1,
  return (GPV[m%PFM1+1]),
  return(VECCPN[N]*GPV[(N*m)%PFM1+1]*lift(WPOL[N]^m));
 );
}

/* compute
prod(j=1,DEG,gamma(frac(p^e*(VAL[j]+m/(p^f-1)))+O(p^6))/gamma(frac(p^e*(VBE[j]+m/(p^f-1)))+O(p^6)));
*/

gapprod(f:small,m:int,e:small,flag=0)=
{
 my(res,c:small);
 res=1;
 for(N:small=1,POLDEG,
  c=VPOLGA[N]:small;
  if (c==0,next());
  res*=gapall(f,m,N,e,flag)^c
 );
 return (res);
}

/* Hypergeometric product */

Cp(f,m,flag=0)=
{
 return (prod(e=0,f-1,gapprod(f,m,e,flag)));
}

/* Power of -p */

/* Stickelberger valuation */

L0(f:small,m:int)=
{
 return (OFFMPOL*f-sum(e=0,f-1,sum(N=1,POLDEG,VPOLGA[N]*((PRIME^e*N*m)\PFM1))));
}

L0simp(m:int)=
{
 if (m<PFL0LIM,
  return(OFFMPOL),
  return (OFFMPOL-sum(N=2,POLDEG,VPOLGA[N]*((N*m)\PFM1)))
 );
}

dopreL0L1(f:small)=
{
 my(tmp:small,vperm,vtmp,vcoe,l0sub,l0coe,co,ww,ct);
 PF=(PRIME^f):small; PFM1=PF-1; PFP1OV2=(PF+1)\2;
 VL1=vector(PFM1);
 for (j=1,DEG,
  if(PFM1%VBEDEN[j]==0,
   tmp=((-(PFM1/VBEDEN[j])*VBENUM[j])%PFM1):small+1;
   VL1[tmp]+=f
  )
 );
 PFL0LIM=ceil(PFM1/POLDEG):small;
 if (f>1,return());
 vtmp=[]; vcoe=[];
 for (N=2,POLDEG,
  if (VPOLGA[N],
   vtmp=concat(vtmp,ceil((PFM1/N)*vector(N,q,q-1)));
   vcoe=concat(vcoe,VPOLGA[N]*vector(N,q,1))
  )
 );
 vperm=vecsort(vtmp,,1);
 vtmp=vecextract(vtmp,vperm); vcoe=vecextract(vcoe,vperm);
 ww=vtmp[1]; l0sub=[ww]; l0coe=[vcoe[1]]; ct=1;
 for (j=2,#vtmp,
  if (vtmp[j]==ww,
   l0coe[ct]+=vcoe[j],
   ct++; ww=vtmp[j]; l0sub=concat(l0sub,ww); l0coe=concat(l0coe,vcoe[j])
  )
 );
 VL0=vector(PFM1); co=OFFMPOL;
 for (j=1,#l0sub-1,
  for (m=l0sub[j],l0sub[j+1]-1,VL0[m+1]=co);
  co-=l0coe[j+1]
 );
 for (m=l0sub[#l0sub],PFM1-1,VL0[m+1]=co);
}
 
C(f:small,m:small,flag=0)=
{
 Cp(f,m,flag)*(-PRIME)^L0(f,m)*PRIME^(-VL1[m+1]);
}

/* General H function */

H(f:small,t,flag=0)=
{
 my(dfp:small,res,L,co,lt:small,tt,vecmo);
 if (DEBUG,gettime());
 if (flag,dfp=16,dfp=(clogp(DEG)+ceil((f*WT)/2)+1):small; if (PRIME==2,dfp++));
 doprecomp(f,dfp); dopreL0L1(f); dopregapcpn(f);
 tt=type(t); if (tt!="t_VEC",t=[t]); lt=#t;
 L=vector(lt,j,teichmuller(t[j]+O(PRIME^dfp))); 
 if (flag,
  vecmo=Pol(vector(PF-1,i,C(f,PF-1-i,flag)));
  vecmo/=(1-PF)*polcoeff(vecmo,0);return(vecmo)
 );
 co=C(f,PF-2,flag); res=vector(lt,j,co);
 for (i=1,PF-2,co=C(f,PF-2-i,flag); for (j=1,lt,res[j]=co+L[j]*res[j]));
 if (tt!="t_VEC",res=res[1]);
 res/=((1-PF)*co);
 if (DEBUG,
  if (f==3,
   GTH3+=gettime(),if (f==2,GTH2+=gettime(),GTH1+=gettime())
  )
 );
 return (res);
}

/* Specific programs when one can work modulo p^2 */

/* compute Gamma_p mod p^2 */

initgamp()=
{
 my(wil:small);
 for(i:small=1,PRIME-1, GF[i+1]=(GF[i]*i)%PSQ);
 wil=((GF[PRIME]+1)\PRIME):small; GH[PRIME]=GH[1]=-wil;
 for(i:small=1,(PRIME-1)\2, GH[PRIME-i]=GH[i+1]=(GH[i]+1/i)%PRIME);
 WPOL=vector(POLDEG,N,Mod(N,PSQ)^(PRIME*N));
}

/* e=0,1,>=2: compute gamma_p(x) modulo p^(2-e) */

gap2(x,e:small)=
{
 my(tmp,r0:small,r1:small,res);
 if (e>=2,return(0));
 tmp=x-1; r0=(tmp%PRIME):small; r1=(((tmp-r0)/PRIME)%PRIME):small;
 res=GF[r0+1]; if ((r0+1+r1*(PRIME+1))%2,res=-res);
 if (e==1,return(res%PRIME));
 return ((res*(1+PRIME*r1*GH[r0+1]))%PSQ);
}

/* compute gamma_p(m/(p-1)) modulo p^2, p odd, 0<=m<p-1 */

gap2simp(m:small)=
{
 my(r0:small,tmp:int);
 r0=PRIME-m-1; 
 tmp=((GF[r0+1]*(1+PRIME*r0*GH[r0+1]))%PSQ):int;
 if (r0%2==0,tmp=PSQ-tmp); return (tmp);
}

/* Called only with f=1, so f suppressed, and new replaced by modp2 */

Cmodp2(m:small)=
{
 my(l0:small,l1:small,e:small);
 l0=VL0[m+1]:small; l1=VL1[m+1]:small-TT; e=l0-l1;
 if (e<2,
  if (e==1,
   return (Cpmodp2(m,PRIME)*(-1)^l0*PRIME),
   return (Cpmodp2(m,PSQ)*(-1)^l0)
  ),
  return(0)
 );
}

Cpmodp2(m:small,M:small)=
{
 return (lift(prod(N=1,POLDEG,Mod(gapallmodp2(m,N),M)^VPOLGA[N])));
}

Hmodp2(f:small,t,flag=0)=
{
 my (tmodp2,lim:small,res,lt:small,tmp,cmo,vecmo,tt2:int);
 if ((f>1)||(PRIME<=LIMMODP2),return(convpadic(PRIME^(f*TT)*H(f,t,flag))));
 if (DEBUG,gettime());
 doprecompmodp2(); dopreL0L1(1); dopregapcpn(1); lim=PRIME-2;
 if (type(t)!="t_VEC",
  tmodp2=truncate(teichmuller(1/t+O(PRIME^3)))%PSQ;
  if (flag,vecmo=vector(lim+1,i,Cmodp2(lim+1-i));print(vecmo);return(Pol(vecmo)));
  res=Cmodp2(0); 
  for (i=1,lim,res=(Cmodp2(i)+tmodp2*res)%PSQ);
  res=(res*tmodp2)%PSQ,
  lt=#t;
  tmodp2=vector(lt,j,truncate(teichmuller(t[j]+O(PRIME^3)))%PSQ);
  cmo=Cmodp2(lim); res=vector(lt);
  vecmo=vector(lim,i,Cmodp2(lim-i));
  for (j=1,lt,
   tmp=cmo; tt2=tmodp2[j]:int;
   for (i=1,lim,tmp=(vecmo[i]+tt2*tmp)%PSQ); res[j]=tmp
  )
 );
 res=centerlift(Mod(res,PSQ)/((1-PRIME)*(-1)^VL0[1]*Cpmodp2(0,PSQ)));
 if (DEBUG,GTHMOD+=gettime());
 return(res);
}

Tr(f:small,t,flag=0)=
{
 if(flag,return(PRIME^(f*TT)*H(f,1,flag)));
 if (f==1,return(Hmodp2(f,t,flag)));
 return (convpadic(PRIME^(f*TT)*H(f,t,flag)));
}

\\ Ceiling of log(x)/log(p)
\\ Actually is  + 1 if x is an exact power of p...

clogp(x:int)=
{
 my(k:small);
 x=abs(x);
 k=0;
 until(x==0,x=x\PRIME;k++);
 return (k);
}

/* guess of the linear factor */

gueli()=
{
 return(1-kronecker(DISC,PRIME)*PRIME^((WT-1)/2)*x);
}

/* Valid for GOOD and ONE. No need for frobpol itself, simply put
   LSIZ such that LSIZ>DEG*log(p) */

frobpoltrunc(t)=
{
 my(B:small,mi:small,pol,qol,S,v,tmp,lt:small,tt);
 B=floor(LSIZ/log(PRIME)):small; mi=max(1,min(DEG2,B));
 S=-sum(f=1,mi,Tr(f,t)*x^f/f); tt=type(t);
 if (DEBUG,COUNTP++; if (COUNTP%200==0,print1(PRIME," ")));
 if (tt!="t_VEC",t=[t];S=[S]);
 lt=#t; v=vector(lt);
 for (i=1,lt,
  pol=truncate(exp(S[i]+O(x^(mi+1))));
  if (valuation(t[i]-1,PRIME)==0 || DEG%2,
   if (mi==DEG2,
    qol=sum(j=max(DEG-B,0),(DEG-1)\2,x^(DEG-j)*PRIME^(WT*(DEG-2*j)/2)*polcoeff(pol,j));
    pol+=qol
   );
   v[i]=pol,
   if (mi==DEG2,
    tmp=gueli(); pol=truncate((pol+O(x^DEG2))/tmp);
    qol=sum(j=max(DEG-2-B,0),(DEG-4)/2,x^(DEG-2-j)*PRIME^(WT*(DEG-2-2*j)/2)*polcoeff(pol,j));
    pol=(pol+qol)*tmp
   );
   v[i]=pol
  )
 );
 if (tt!="t_VEC",return (v[1]),return (v));
}
 
cl(t)=
{
 my(v:small,tt,lt:small,res);
 tt=type(t);
 if (tt!="t_VEC",t=[t]); 
 lt=#t;
 if (BAD%PRIME==0,if (tt!="t_VEC",return(1*BAD),return(vector(lt,i,BAD))));
 res=vector(lt);
 for (j=1,lt,
  v=valuation(t[j],PRIME);
  if (v<0,
   res[j]=INFINIT,
   if (v>0,
    res[j]=ZERO,
    if (valuation(t[j]-1,PRIME)>0,res[j]=ONE,res[j]=GOOD)
   )
  )
 );
 if (tt!="t_VEC",res=res[1]);
 return (res);
}

eulfacspec(t,cla:small)=
{
 if (cla==ZERO,return(eulfaczer(t)));
 if (cla==INFINIT,return(eulfacinf(t)));
 if (cla==BAD,return(eulfacbad(t)));
}
 
eulfac(t,p:small)=
{
 my(cla,tt,lt:small,res,ret,tsub,ct:small);
 PRIME=p; PSQ=(p*p):small;
 tt=type(t); cla=cl(t);
 if (tt!="t_VEC",
  if ((cla==GOOD) || (cla==ONE),
   return(frobpoltrunc(t)),
   return(eulfacspec(t,cla))
  ),
  lt=#t; ct=0; tsub=vector(lt);
  for (j=1,lt,if((cla[j]==GOOD) || (cla[j]==ONE),ct++; tsub[ct]=t[j]));
  tsub=vector(ct,j,tsub[j]);
  if (ct,res=frobpoltrunc(tsub),res=[]);
  ret=vector(lt); ct=0;
  for (j=1,lt,
   if ((cla[j]==GOOD) || (cla[j]==ONE),
    ct++; ret[j]=res[ct],
    ret[j]=eulfacspec(t[j],cla[j])
   )
  );
  return (ret)
 );
}

eulfaczer(t)=
{
 return (1);
}

eulfacinf(t)=
{
 return (1);
}

eulfacbad(t)=
{
 return (1);
}

/* computes euler product direuler(p=2,SIZ,v[p][i]). Used instead of
direuler for compilation with gp2c, otherwise can use commented lines
below. */

mydireuler(v,i:small,t)=
{
 my(res,pol,pd:small,va:small);
 res=DIRONE;
 forprime (p=2,SIZ,
  if (i>0,pol=v[p][i],pol=eulfac(t,p)); 
  pd=poldegree(pol);
  forstep (n=p,SIZ,p,
   va=valuation(n,p); 
   if (va<=pd,res[n]=polcoeff(pol,va)*res[n/p^va]);
  )
 );
 return (res);
}

/* Warning: fills only LKEEP; need to do LFUN=LKEEP or LFUN=LKEEP[in] after,
   e.g using modif, or simply LFUN=filleuler(t); */

filleuler(t,siz:small)=
{
 my(lt,v,tt,ct,lkeep);
 SIZ=siz; LSIZ=log(SIZ); GF=vector(SIZ); GF[1]=1; GH=vector(SIZ); 
 DIRONE=vector(SIZ,i,(i==1));
 GTHMOD=GTH1=GTH2=GTH3=0; COUNTP=0; tt=type(t);
 if (DEBUG,ct=0; forprime(p=2,SIZ,ct++); print(ct," primes to be done"));
 if (tt!="t_VEC",
  lkeep=dirdiv(DIRONE,mydireuler(v,0,t)),
  lt=#t; v=vector(SIZ); forprime(p=2,SIZ,v[p]=eulfac(t,p));
  lkeep=vector(lt,i,dirdiv(DIRONE,mydireuler(v,i,t)));
 );
 if (DEBUG, print([GTH1,GTH2,GTH3,GTHMOD]));
 LKEEP=lkeep;
 return (lkeep);
}
 
modif(lkeep,exc,in=0)=
{
 my(tmp,pr:small,pol,tmpvec,lfun,siz:small);
 if (in<=0, lfun=lkeep, lfun=lkeep[in]);
 siz=#lfun;
 for(i=1,#exc,
  tmp=exc[i]; tmpvec=vector(siz); tmpvec[1]=1;
  pr=tmp[1]:small; pol=tmp[2];
  for (j=1,poldegree(pol),if (pr^j<=siz,tmpvec[pr^j]=polcoeff(pol,j)));
  lfun=dirdiv(lfun,tmpvec)
 );
 return (lfun);
}

kpsi2(z)=exp(-2*Pi*z);

ph2(lfun,N:int,x,eps=10.^(-30))=
{
 my(su,sn,c:int,kp,siz:small);
 su=0.; sn=x/sqrt(N); siz=#lfun;
 for (n=1,siz,
  c=lfun[n]:int;
  if(c,
   kp=kpsi2(n*sn); su+=c*kp;
   if (kp<eps,break())
  )
 );
 return (x^2*su);
}

test2(lfun,N:int,s=1.1)=
{
 my(pre:small,res);
 pre=precision(1.):small;
 default(realprecision,38);
 res=ph2(lfun,N,1/s)/ph2(lfun,N,s);
 default(realprecision,pre);
 return (res);
}

test2wt4(lfun,N:int,s=1.1)=test2(lfun,N,s);
test2wt2(lfun,N:int,s=1.1)=test2(lfun,N,s)*s^2;
test2wt1(lfun,N:int,s=1.1)=test2(lfun,N,s)*s^3;

/* up to now generic. Here assume HODGE = 1+T+T^2+T^3 */

kpsi4(z)=besselk(1,4*Pi*sqrt(z))/(Pi*sqrt(z));

ph4(lfun,N:int,x,eps=10.^(-30))=
{
 my(su,sn,c,kp,siz:small);
 su=0.; sn=x/sqrt(N); siz=#lfun;
 for (n=1,siz,
  c=lfun[n];
  if(c,
   kp=kpsi4(n*sn); su+=c*kp;
   if (kp<eps,break())
  )
 );
 return (x^2*su);
}

errest4(lfun,N:int)=
{
 my(m:int,siz:small);
 m=vecmax(abs(lfun)):int; siz=#lfun;
 return (m*max(10.^(-30),sqrt(siz)*kpsi4(siz/(1.1*sqrt(N)))));
}

test4simp(lfun,N:int,s=1.1)=
{
 my(eps,pre:small,res);
 pre=precision(1.):small;
 default(realprecision,19);
 eps=10.^(-15);
 res=ph4(lfun,N,1/s,eps)/ph4(lfun,N,s,eps);
 default(realprecision,pre);
 return(res);
}

test4(lfun,N:int,s=1.1)=
{
 my(pre:small,res);
 pre=precision(1.):small;
 default(realprecision,38);
 res=ph4(lfun,N,1/s)/ph4(lfun,N,s);
 default(realprecision,pre);
 return (res);
}

/* These are only for gp, not gp2c */

tst4s(N:int,s=1.1)=test4simp(LFUN,N,s);

tst4(N:int,s=1.1)=test4(LFUN,N,s);

erst4(N:int)=errest4(LFUN,N);

mdf(exc,in=0)=modif(LKEEP,exc,in);

/* compute the inverse mellin transform of \Gamma(z)\Gamma(z-1)^2 */

/* For the moment no optimization */

/* f(x) approximately 3.6275...*exp(-3*x^(1/3))/x, max term exp(3*x^(1/3)) 
for n=x^(1/3)e+D*log(10)/3 */

invm(x)=
{
 my(D:small,eps,lo,s,n:small,qn,x3,xn,newprec:small);
 D=precision(1.):small; 
 default(realprecision,9);
 x3=x^(1/3); eps=10.^(-D)*exp(-3*x3)/x;
 newprec=ceil(D+2.61*x3):small;
 default(realprecision,newprec);
 lo=log(x); s=-(lo+3*Euler())/x; n=1; qn=1./2; xn=-x;
 until(abs(qn)<eps,
  s+=qn*((lo-VECA[n])^2+VECC[n]); qn*=xn/(n*(n+1)^2); n++
 );
 s=precision(s,D);
 default(realprecision,D);
 return (s);
}

/* kpsi6(z)=invm((2*Pi)^3*z); */

invm2(x1,x2,lodif)=
{
 my(D:small,eps,lo1,lo2,s1,s2,n:small,qn1,qn2,x3,xn1,xn2,newprec:small,nn:int);
 my(lo12,lo22,lo1s,lo2s);
 D=precision(1.):small; 
 default(realprecision,9);
 x3=x1^(1/3); eps=10.^(-D)*exp(-3*x3)/x1;
 newprec=ceil(D+2.61*x3):small;
 default(realprecision,newprec);
 lo1=log(x1); lo2=lo1+lodif;
 s1=-(lo1+3*Euler())/x1; s2=-(lo2+3*Euler())/x2; 
 n=1; qn1=qn2=1./2; xn1=-x1; xn2=-x2;
 lo12=lo1*lo1; lo22=lo2*lo2; lo1s=lo1<<1; lo2s=lo2<<1;
 until(abs(qn1)<eps,
  nn=n*(n+1)^2; 
  s1+=qn1*(lo12-lo1s*VECA[n]+VECC[n]); qn1*=xn1/nn;
  s2+=qn2*(lo22-lo2s*VECA[n]+VECC[n]); qn2*=xn2/nn;
  n++
 );
 s1=precision(s1,D); s2=precision(s2,D);
 default(realprecision,D);
 return ([s1,s2]);
}

ph6simp(lfun,N:int,x,eps=10.^(-30))=
{
 my(su,sn,c:int,kp,siz:small);
 su=0.; sn=(2*Pi)^3*x/sqrt(N); siz=#lfun;
 for (n=1,siz,
  c=lfun[n]:int;
  if(c,
   kp=invm(n*sn); su+=c*kp; if (kp<eps,break())
  )
 );
 return (x^2*su);
}

ph6(lfun,N:int,x1,x2,eps=10.^(-30))=
{
 my(su,sn1,sn2,c:int,kp,siz:small,lodif);
 siz=#lfun;
 default(realprecision,57);
 su=[0.,0.]; sn1=(2*Pi)^3*x1/sqrt(N); sn2=(2*Pi)^3*x2/sqrt(N);
 lodif=log(x2/x1);
 default(realprecision,38);
 for (n=1,siz,
  c=lfun[n]:int;
  if(c,
   kp=invm2(n*sn1,n*sn2,lodif); 
   su+=c*kp; if (kp[1]<eps,break())
  )
 );
 return ([x1^2*su[1],x2^2*su[2]]);
}

test6simp(lfun,N:int,s=11/10)=
{
 my(eps,pre:small,res,vvh,vvh2,siz:small);
 pre=precision(1.):small;
 default(realprecision,38);
 siz=#lfun;
 vvh=vector(siz+1); vvh2=vector(siz+1);
 vvh[1]=-Euler(); vvh2[1]=Pi^2/6; 
 for (n=1,siz,vvh[n+1]=vvh[n]+1/n;vvh2[n+1]=vvh2[n]+1/n^2);
 VECA=vector(siz,n,2*vvh[n+1]+vvh[n]); 
 VECC=vector(siz,n,2*vvh2[n+1]+vvh2[n]);
 eps=10.^(-18);
 res=ph6simp(lfun,N,1/s,eps)/ph6simp(lfun,N,s,eps);
 default(realprecision,pre);
 return (res);
}

/* Attention: here VECC is the preceding VECC + VECA^2 */

test6(lfun,N:int,s=11/10)=
{
 my(pre:small,res,vvh,vvh2,siz:small);
 pre=precision(1.):small;
 default(realprecision,57);
 siz=#lfun;
 vvh=vector(siz+1); vvh2=vector(siz+1);
 vvh[1]=-Euler(); vvh2[1]=Pi^2/6; 
 for (n=1,siz,vvh[n+1]=vvh[n]+1/n;vvh2[n+1]=vvh2[n]+1/n^2);
 VECA=vector(siz,n,2*vvh[n+1]+vvh[n]); 
 VECC=vector(siz,n,2*vvh2[n+1]+vvh2[n]+VECA[n]^2);
 default(realprecision,38);
 res=ph6(lfun,N,1/s,s);
 default(realprecision,pre);
 return (res[1]/res[2]);
}

/* x<=(2*Pi)^3*siz*s/sqrt(N) */

/* find level of mod form associated to w3r4[j] */

findlevel4(j,siz=3000)=
{
 initalbe(w3r4[j][2],w3r4[j][1]);
 print1(j,": ",w3r4[j],": "); 
 LK=filleuler(1,siz);
 LL=dirdiv(LK,vector(siz,n,n*kronecker(DISC*BAD^2,n)));
}

findlevel6(j,siz=3000)=
{
 initalbe(w3r6[j][2],w3r6[j][1]);
 print1(j,": ",w3r6[j],": "); 
 LK=filleuler(1,siz);
 LL=dirdiv(LK,vector(siz,n,n*kronecker(DISC*BAD^2,n)));
}

findlevel2(lfun,aa,bb)=
{
 LM=subst(subst(lfun,a,aa),b,bb);
 lpr=factor(abs(BAD))[,1]; llpr=#lpr;
 vfor=vector(llpr,i,[0,if(lpr[1]<=3,8,2)]);
 forvec(v=vfor,
  N=prod(i=1,llpr,lpr[i]^v[i]); res=test2wt2(LM,N);
  if (abs(abs(res)-1)<0.00001,print([aa,bb,factor(N)],": ",res))
 );
}

case12(t)=
{
 initalbe([1/3,2/3],[1/6,5/6]);
 LK=filleuler(t,10000);
 eps=10.^(-8);
 k2=valuation(t,2); k3=valuation(t,3); r=t/(2^k2*3^k3);
 d=denominator(r);n=numerator(r);
 fa=factor(abs(numerator(t-1)))[,1];
 fa1=[];for(i=1,#fa,if(fa[i]>3,fa1=concat(fa1,fa[i])));
 fa2=factor(d)[,1];
 fa3=factor(abs(n))[,1];
 faall=concat(concat([2,3],fa1),concat(fa2,fa3)); lfa=#faall;
 vfor=vector(lfa,i,[-1,1]);
 vfor2=vector(lfa,i,if(faall[i]<=3,[0,8],[0,2]));
 forvec(v2=vfor2,
  N=prod(i=1,lfa,faall[i]^v2[i]); res=test2wt2(LK,N);
  if (abs(abs(res)-1)<eps,print([0,factor(N)],": ",res);return())
 );
 forvec(v=vfor,
  if (v,
   w=vector(lfa,i,[faall[i],1-v[i]*x]);
   LL=modif(LK,w);
   forvec(v2=vfor2,
    N=prod(i=1,lfa,faall[i]^v2[i]); res=test2wt2(LL,N);
    if (abs(abs(res)-1)<0.00001,print([v,factor(N)],": ",res))
   )
  )
 );
}


case12cond(t)=
{
 initalbe([1/3,2/3],[1/6,5/6]);
 LK=filleuler(t,2000);
 eps=10.^(-8);
 k2=valuation(t,2); k3=valuation(t,3); r=t/(2^k2*3^k3);
 d=denominator(r);n=numerator(r);
 fa=factor(abs(numerator(t-1)))[,1];
 fa1=[];for(i=1,#fa,if(fa[i]>3,fa1=concat(fa1,fa[i])));
 fa2=factor(d)[,1];
 fa3=factor(abs(n))[,1];
 faall=concat(concat([2,3],fa1),concat(fa2,fa3)); lfa=#faall;
 vfor=vector(lfa,i,[-1,1]);
 vfor2=vector(lfa,i,if(faall[i]<=3,[0,8],[0,2]));
 forvec(v2=vfor2,
  N=prod(i=1,lfa,faall[i]^v2[i]); res=test2wt2(LK,N);
  if (abs(abs(res)-1)<eps,return(N))
 );
 forvec(v=vfor,
  if (v,
   w=vector(lfa,i,[faall[i],1-v[i]*x]);
   LL=modif(LK,w);
   forvec(v2=vfor2,
    N=prod(i=1,lfa,faall[i]^v2[i]); res=test2wt2(LL,N);
    if (abs(abs(res)-1)<0.00001,
     if (abs(abs(res)-1)<eps,return(N))
    )
   )
  )
 );
 return(1);
}