Permalink
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
154 lines (132 sloc) 6.75 KB
subsets(A, k) = {
my (lst = List());
forvec(v = vector(k, i, [1, #A]), listput(lst, vecextract(A, v)), 2);
return(Vec(lst));
};
position = (elt, array) -> select((x) -> x == elt, array, 1)[1];
bt(n,t) = n^((t/4)*log(n));
at(n,t,y) = n^(y+(t/4)*log(n));
alpha1(s) = 1/(2*s) + 1/(s-1) + (1/2)*log(s/(2*Pi));
divdelta(n,d) = if(n%d==0,1,0);
threshdelta(n,dN) = if(n<=dN,1,0);
combdelta(n,divisor,N) = divdelta(n,divisor)*threshdelta(n,divisor*N);
mollifier(mtype,t,mvariety="normal")={
if(mtype==1, return(List([1,1,[1],[1]])));
if(mtype>=2,\
np = position(mtype,primes(200));\
p = primes(np); mollp = vector(np,i,bt(p[i],t));\
divsubs = subsets(p,1); mollsubs=subsets(mollp,1);\
if(mvariety!="prime",for(i=2,np,divsubs=concat(divsubs,subsets(p,i)); mollsubs=concat(mollsubs,subsets(mollp,i))));\
d = vector(length(divsubs),j,prod(k=1,length(divsubs[j]),divsubs[j][k]));\
moll = vector(length(mollsubs),j,(-1)^length(mollsubs[j])*prod(k=1,length(mollsubs[j]),mollsubs[j][k]));\
d=concat([1],d); moll=concat([1],moll);\
ndivs = length(d); D = d[ndivs];\
return(List([D,ndivs,d,moll]));\
);
}
bound_constants(N,y,t) = {
xN = 4*Pi*N^2 - Pi*t/4;
xNp1 = 4*Pi*(N+1)^2 - Pi*t/4;
sig = (1+y)/2 + (t/4)*log(xN/(4*Pi)) - (t/(2*xN^2))*max(0,1 - 3*y + 4*y*(1+y)/xN^2);
modK = t*y/(2*(xN-6));
modgamma = exp(0.02*y)*(xN/(4*Pi))^(-y/2);
return([sig,modK,modgamma]);
}
abbeff_largex_ep_bound(N,y,t,btype="T",mtype=1,mvariety="normal",tail="no")={
[sig,modK,modgamma] = bound_constants(N,y,t);
[D,ndivs,d,moll] = mollifier(mtype,t,mvariety);
modmoll = sum(nd=1,ndivs,abs(moll[nd])/d[nd]^sig);
if(btype=="T",\
if(tail=="no",\
b = sum(n=2,D*N,abs(sum(nd=1,ndivs,if(combdelta(n,d[nd],N)>0,moll[nd]*bt(n/d[nd],t))))/n^sig);\
a = modgamma*sum(n=2,D*N,abs(sum(nd=1,ndivs,if(combdelta(n,d[nd],N)>0,moll[nd]*at(n/d[nd],t,y))))/n^sig);\
lbound = (1/modmoll)*(1 - modgamma - b - a);\
);\
if(tail=="yes",\
b = sum(n=2,N,abs(sum(nd=1,ndivs,if(combdelta(n,d[nd],N)>0,moll[nd]*bt(n/d[nd],t))))/n^sig);\
b = b + sum(nd=2,ndivs,(abs(moll[nd])*bt(d[nd],t)/N^((t/2)*log(d[nd])))*(N^(1-(1+y)/2-(t/4)*log(N))+(d[nd]*N)^(1-(1+y)/2-(t/4)*log(N/d[nd])))*log(d[nd])/d[nd]/2);\
a = modgamma*sum(n=2,N,abs(sum(nd=1,ndivs,if(combdelta(n,d[nd],N)>0,moll[nd]*at(n/d[nd],t,y))))/n^sig);\
a = a + modgamma*sum(nd=2,ndivs,(abs(moll[nd])*bt(d[nd],t)/(d[nd]^y*N^((t/2)*log(d[nd]))))*(N^(1-(1-y)/2-(t/4)*log(N))+(d[nd]*N)^(1-(1-y)/2-(t/4)*log(N/d[nd])))*log(d[nd])/d[nd]/2);\
lbound = (1/modmoll)*(1 - modgamma - b - a);\
);\
);
if(btype=="L",\
common = (1-modgamma)/(1+modgamma);\
if(tail=="no",\
b = vector(D*N-1,n,sum(nd=1,ndivs,if(combdelta(n+1,d[nd],N)>0,moll[nd]*bt((n+1)/d[nd],t))));\
a = modgamma*vector(D*N-1,n,sum(nd=1,ndivs,if(combdelta(n+1,d[nd],N)>0,moll[nd]*at((n+1)/d[nd],t,y))));\
lbound = (1/modmoll)*(1 - modgamma - sum(n=1,D*N-1,max(abs(b[n]-a[n]),common*abs(b[n]+a[n]))/(n+1)^sig));\
);\
if(tail=="yes",\
b = vector(N-1,n,sum(nd=1,ndivs,if(combdelta(n+1,d[nd],N)>0,moll[nd]*bt((n+1)/d[nd],t))));\
a = modgamma*vector(N-1,n,sum(nd=1,ndivs,if(combdelta(n+1,d[nd],N)>0,moll[nd]*at((n+1)/d[nd],t,y))));\
tailb = sum(nd=2,ndivs,(abs(moll[nd])*bt(d[nd],t)/N^((t/2)*log(d[nd])))*(N^(1-(1+y)/2-(t/4)*log(N))+(d[nd]*N)^(1-(1+y)/2-(t/4)*log(N/d[nd])))*log(d[nd])/d[nd]/2);\
taila = modgamma*sum(nd=2,ndivs,(abs(moll[nd])*bt(d[nd],t)/(d[nd]^y*N^((t/2)*log(d[nd]))))*(N^(1-(1-y)/2-(t/4)*log(N))+(d[nd]*N)^(1-(1-y)/2-(t/4)*log(N/d[nd])))*log(d[nd])/d[nd]/2);\
lbound = (1/modmoll)*(1 - modgamma - sum(n=1,N-1,max(abs(b[n]-a[n]),common*abs(b[n]+a[n]))/(n+1)^sig) - tailb - taila);\
);\
);
lbound = lbound - modgamma*sum(n=1,N,bt(n,t)*(n^modK-1)/n^(sig-y));
return(lbound);
}
ba_lemma(n,d,moll,N,t,y) = sum(nd=1,#d,if(combdelta(n,d[nd],N)>0,moll[nd]*[at(n/d[nd],t,0),at(n/d[nd],t,y)],[0,0]));
lsummand(mollba,modgamma,modgamend,n,sig) = if(mollba==[0,0],0,vecmax([(1-modgamend)/(1+modgamend)*abs(mollba[1]+modgamma*mollba[2]),abs(mollba[1]-modgamma*mollba[2]),abs(mollba[1]-modgamend*mollba[2])])/n^sig);
abbeff_largex_ep_sawtooth_incremental_lbounds(boundfile,initbound,threshold,Nstart,Nend,y,t,mtype)={
modgamend = exp(0.02*y)*((4*Pi*Nend^2 - Pi*t/4)/(4*Pi))^(-y/2);
[D,ndivs,d,moll] = mollifier(mtype,t);
prevbound=initbound;
for(nextN=Nstart,Nend,\
[sig,modK,modgamma] = bound_constants(nextN,y,t);
modmoll = sum(nd=1,ndivs,abs(moll[nd])/d[nd]^sig);\
offset = vector(ndivs,nd,d[nd]*(nextN-1));\
lbound = prevbound - (1/modmoll)*sum(v=1,ndivs,sum(n=1,d[v],lsummand(ba_lemma(n+offset[v],d,moll,nextN,t,y),modgamma,modgamend,n+offset[v],sig))) - modgamma*bt(nextN,t)*(nextN^modK-1)/nextN^(sig-y);\
\\print(nextN);
if(lbound<threshold,return(nextN));\
row = [nextN,lbound];\
if(nextN%10000==0,print(row);write(boundfile,row););\
prevbound=lbound;\
);
return(Nend);
}
abbeff_largex_ep_sawtooth_base_lbound(boundfile,Nmin,Nend,y,t,mtype)={
modgamend = exp(0.02*y)*((4*Pi*Nend^2 - Pi*t/4)/(4*Pi))^(-y/2);
[D,ndivs,d,moll] = mollifier(mtype,t);
[sig,modK,modgamma] = bound_constants(Nmin,y,t);
modmoll = sum(nd=1,ndivs,abs(moll[nd])/d[nd]^sig);
lbound = (1/modmoll)*(1 - modgamma - sum(n=1,D*Nmin-1,lsummand(ba_lemma(n+1,d,moll,Nmin,t,y),modgamma,modgamend,n+1,sig))) - modgamma*sum(n=1,Nmin,bt(n,t)*(n^modK-1)/n^(sig-y));
row = [Nmin,lbound,"Nmin bound"];
print(row);write(boundfile,row);
return(lbound);
}
abbeff_nomoll_analytic_bound(N,y,t,N0=2000)={
\\based on older bound but can provide a bound estimate for arbitrarily large N
N0 = if(N<N0,N,N0);
xN = 4*Pi*N^2 - Pi*t/4;
xNp1 = 4*Pi*(N+1)^2 - Pi*t/4;
delta = Pi*y/(2*(xN - 6 - (14 + 2*y)/Pi)) + 2*y*(7+y)*log(sqrt((1+y)^2+xNp1^2)/(4*Pi))/xN^2;
afac = exp(delta + t*y*log(N)/(2*(xN-6)))/N^y;
sigb = (1+y)/2;
siga = (1-y)/2;
sb = sigb + (t/2)*log(N);
sa = siga + (t/2)*log(N);
bbound = sum(n=1,N0,bt(n,t)/n^sb) + max(N0^(1-sigb-(t/4)*log(N^2/N0)),N^(1-sigb-(t/4)*log(N)))*log(N/N0);
abound = sum(n=1,N0,bt(n,t)/n^sa) + max(N0^(1-siga-(t/4)*log(N^2/N0)),N^(1-siga-(t/4)*log(N)))*log(N/N0);
lbound = 2 - (bbound + afac*abound);
return(lbound);
}
erfi(z) = if(z>0,abs(1-erfc(I*z)),-abs(1-erfc(I*z)));
tailint(llim,ulim,a,b) = (sqrt(Pi)/sqrt(b)/2)*exp(-(a-1)^2/(4*b))*(erfi((2*b*log(ulim)-a+1)/sqrt(b)/2)-erfi((2*b*log(llim)-a+1)/sqrt(b)/2));
abbeff_moll2_analytic_bound(N,y,t,N0=1000)={
S1_main = sum(n=1,N0,(1+n^y/N^y)*n^(-(1+y)/2-(t/4)*log(N^2/n)));
S1_tail = tailint(N0,N,(1+y)/2+(t/2)*log(N),t/4) + (1/N^y)*tailint(N0,N,(1-y)/2+(t/2)*log(N),t/4);
S1 = S1_main + S1_tail;
modb2 = 2^(-(1+y)/2-(t/4)*log(N^2/2));
bound = 2 - (1-modb2)*S1 - 2*modb2*(tailint(N/2,N,(1+y)/2+(t/2)*log(N),t/4) + (1/N^y)*tailint(N/2,N,(1-y)/2+(t/2)*log(N),t/4));
return(bound);
}
boundfile="~/t_0.2_y_0.2_moll5_sawtooth_v2.txt";
Nend=1.5*10^6;threshold=0.03;y=0.2;t=0.2;mtype=5;
Nmin = 80000;
while(Nmin<=Nend,\
initbound = abbeff_largex_ep_sawtooth_base_lbound(boundfile,Nmin,Nend,y,t,mtype);\
Nmin = abbeff_largex_ep_sawtooth_incremental_lbounds(boundfile,initbound,threshold,Nmin+1,Nend,y,t,mtype);\
);