Permalink
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
245 lines (200 sloc) 9.73 KB
default(realprecision, 30);
default(parisizemax,3*10^9);
printm(M)=apply(s->print(s),Col(M));print1();
factorial_divide(val,n)={
while(n>1,val=val/n;n=n-1);
return(val);
}
bt(n,t)=return(n^((t/4)*log(n)));
alpha1(s) = 1/(2*s) + 1/(s-1) + (1/2)*log(s/(2*Pi));
alpha1prime(s) = -1/(2*s^2) - 1/(s-1)^2 + 1/(2*s);
H01(s) = (1/2)*s*(s-1)*Pi^(-s/2)*sqrt(2*Pi)*exp((s/2-1/2)*log(s/2)-s/2);
C0(p) = (exp(Pi*I*(p^2/2 + 3/8)) - I*sqrt(2)*cos(Pi*p/2))/(2*cos(Pi*p));
B0_eff(x,y,t) = (1/8)*exp((t/4)*alpha1((1+y-I*x)/2)^2)*H01((1+y-I*x)/2);
abbeff(x,y,t) = {
s = (1 - y + I*x)/2;
N = floor(sqrt((x+Pi*t/4)/(4*Pi)));
alph1 = alpha1(s);
alph2 = alpha1(1-s);
A0 = (1/8)*exp((t/4)*alph1^2)*H01(s);
B0 = (1/8)*exp((t/4)*alph2^2)*H01(1-s);
A_sum = sum(n=1,N,n^((t/4.0)*log(n) - (t/2.0)*alph1 - s));
B_sum = sum(n=1,N,n^((t/4.0)*log(n) - (t/2.0)*alph2 - (1-s)));
A = A0 * A_sum;
B = B0 * B_sum;
return ((A+B)/B0);
}
\\taylorvec(xval,nterms) = vector(nterms,n,xval^(n-1)/factorial(n-1));
taylorvec(xval,nterms) = {
resvec = vector(nterms);
resvec[1]=1; for(ctr=1,nterms-1,resvec[ctr+1] = resvec[ctr]*xval/ctr);
return(resvec);
}
\\powervec(val,nterms) = vectorv(nterms,i,val^(i-1));
powervec(val,nterms) = {
resvec = vectorv(nterms);
resvec[1]=1; for(ctr=1,nterms-1,resvec[ctr+1] = resvec[ctr]*val);
return(resvec);
}
\\powermat(valvec,nterms) = matrix(nterms,#valvec,i,j,valvec[j]^(i-1));
powermat(valvec,nterms) = {
resmat = matrix(nterms,#valvec);
for(j=1,#valvec,resmat[1,j]=1;for(i=2,nterms,resmat[i,j] = resmat[i-1,j]*valvec[j]));
return(resmat);
}
matet(xval,expterms,ttaylorterms) = mattranspose(taylorvec(xval,expterms))*taylorvec(xval^2/4,ttaylorterms);
matsba(vecab,mat) = [vecab[1]*mat,vecab[2]*mat];
mattovect(mat,tval) = vector(#mat[,1],i,mat[i,]*powervec(tval,#mat));
vectoval(vctr,val) = vctr*powervec(val,#vctr);
horner(coeff,val) = {res=0;for(i=1,#coeff,res=res*val+coeff[i]);return(res);};
vechorner(coeff,valvec) = {n=#valvec;res=vector(n,j,0);for(i=1,#coeff,res=vecmul(res,valvec)+vector(n,j,coeff[i]));return(res);};
vecmul(v1,v2) = vector(#v1,i,v1[i]*v2[i]);
vecdiv(v1,v2) = vector(#v1,i,v1[i]/v2[i]);
adds(vctr,sclr) = vctr + vector(#vctr,i,sclr);
ddzboundvec(x,y,tmax=0.22)={
N = floor(sqrt(x/4/Pi));
modgamma = exp(0.02*y)*(x/(4*Pi))^(-y/2);
afac = modgamma*N^(tmax*y/(2*(x-6)));
const1 = log(x/(4*Pi))/4 - max(0,1 - 3*y + 4*y*(1+y)/x^2)/(2*x^2);
k=0; est1=1;
while(abs(est1)>0.00001,\
est1 = N*intnum(u=1/N,1,(((tmax/4)*log(N*u)/(x-6) + log(N*u)/2) + afac*(N*u)^y*((tmax/4)*log(N*u)/(x-6) + (log(abs(1+y+I*x)/(4*Pi)) + Pi + 3/x)*((tmax/4)/(x-6) + 1/2)))*(N*u)^(-(1+y)/2)*(tmax*log(N*u)*(log(N*u)/4 - const1))^k/factorial(k),3);\
\\print([k,est1]);\
k=k+1;\
);
taylorterms=k;
boundvec = sum(n=1,N,(((tmax/4)*log(n)/(x-6) + log(n)/2) + afac*(n^y)*((tmax/4)*log(n)/(x-6) + (log(abs(1+y+I*x)/(4*Pi)) + Pi + 3/x)*((tmax/4)/(x-6) + 1/2)))*n^(-(1+y)/2)*taylorvec(log(n)*(log(n)/4 - const1),taylorterms));
return(boundvec);
}
ddtboundvec(x,y,tmax=0.22)={
N = floor(sqrt(x/4/Pi));
modgamma = exp(0.02*y)*(x/(4*Pi))^(-y/2);
afac = modgamma*N^(tmax*y/(2*(x-6)));
const1 = log(x/(4*Pi))/4 - max(0,1 - 3*y + 4*y*(1+y)/x^2)/(2*x^2);
k=0; est1=1;
while(abs(est1)>0.00001,\
est1 = N*intnum(u=1/N,1,(((1/4)*log(N*u)*log(x/(4*Pi*N*u)) + (Pi/8)*log(N*u) + 2*log(N*u)/(x-6)) + afac*(N*u)^y*((1/4)*log(N*u)*log(x/(4*Pi*N*u)) + (Pi/8)*log(N*u) + 2*log(N*u)/(x-6) + (1/4)*(Pi/2 + 8/(x-6))*(log(x/(4*Pi)) + 8/(x-6))))*(N*u)^(-(1+y)/2)*(tmax*log(N*u)*(log(N*u)/4 - const1))^k/factorial(k),3);\
\\print([k,est1]);\
k=k+1;\
);
taylorterms=k;
boundvec = sum(n=1,N,(((1/4)*log(n)*log(x/(4*Pi*n)) + (Pi/8)*log(n) + 2*log(n)/(x-6)) + afac*(n^y)*((1/4)*log(n)*log(x/(4*Pi*n)) + (Pi/8)*log(n) + 2*log(n)/(x-6) + (1/4)*(Pi/2 + 8/(x-6))*(log(x/(4*Pi)) + 8/(x-6))))*n^(-(1+y)/2)*taylorvec(log(n)*(log(n)/4 - const1),taylorterms));
return(boundvec);
}
ddz_abbeff_Nbound(x,y,t)={
N = floor(sqrt(x/4/Pi));
sig = (1+y)/2 + (t/4)*log(x/(4*Pi)) - (t/(2*x^2))*max(0,1 - 3*y + 4*y*(1+y)/x^2);
modK = t*y/(2*(x-6));
modgamma = exp(0.02*y)*(x/(4*Pi))^(-y/2);
ddxbound_sumb = sum(n=1,N,bt(n,t)*((t/4)*log(n)/(x-6) + log(n)/2)/n^sig);
ddxbound_suma = sum(n=1,N,bt(n,t)*(n^y)*((t/4)*log(n)/(x-6) + (log(abs(1+y+I*x)/(4*Pi)) + Pi + 3/x)*((t/4)/(x-6) + 1/2))/n^sig);
ddxbound = ddxbound_sumb + modgamma*(N^modK)*ddxbound_suma;
return(ddxbound);
}
ddt_abbeff_Nbound(x,y,t)={
N=sqrt(x/4/Pi);
sig = (1+y)/2 + (t/4)*log(x/(4*Pi)) - (t/(2*x^2))*max(0,1 - 3*y + 4*y*(1+y)/x^2);
modK = t*y/(2*(x-6));
modgamma = exp(0.02*y)*(x/(4*Pi))^(-y/2);
ddtb_sumb = sum(n=1,N,bt(n,t)*((1/4)*log(n)*log(x/(4*Pi*n)) + (Pi/8)*log(n) + 2*log(n)/(x-6))/n^sig);
ddtb_suma = sum(n=1,N,bt(n,t)*(n^y)*((1/4)*log(n)*log(x/(4*Pi*n)) + (Pi/8)*log(n) + 2*log(n)/(x-6) + (1/4)*(Pi/2 + 8/(x-6))*(log(x/(4*Pi)) + 8/(x-6)))/n^sig);
ddtbound = ddtb_sumb + modgamma*(N^modK)*ddtb_suma;
return(ddtbound);
}
ddzboundint(x,y,t)={
N = floor(sqrt(x/4/Pi));
afac = exp(0.02*y)*(x/(4*Pi))^(-y/2)*N^(t*y/(2*(x-6)));
const1 = log(x/(4*Pi))/4 - max(0,1 - 3*y + 4*y*(1+y)/x^2)/(2*x^2);
estinit = afac*(log(abs(1+y+I*x)/(4*Pi)) + Pi + 3/x)*((t/4)/(x-6) + 1/2);
intest = estinit + N*intnum(u=1/N,1,(log(N*u)*(((t/4)/(x-6))*(1+afac*(N*u)^y) + 1/2) + (N*u)^y*estinit)*(N*u)^(-(1+y)/2 + t*(log(N*u)/4 - const1)),1);\
return(intest);
}
ddtboundint(x,y,t)={
N = floor(sqrt(x/4/Pi));
afac = exp(0.02*y)*(x/(4*Pi))^(-y/2)*N^(t*y/(2*(x-6)));
const1 = log(x/4/Pi)/4 - max(0,1 - 3*y + 4*y*(1+y)/x^2)/(2*x^2);
const2 = (1/4)*log(x/4/Pi) + (Pi/8) + 2/(x-6);
estinit = afac*(1/4)*(Pi/2 + 8/(x-6))*(log(x/(4*Pi)) + 8/(x-6));
intest = estinit + N*intnum(u=1/N,1,(log(N*u)*(const2 - (1/4)*log(N*u))*(1+afac*(N*u)^y) + estinit*(N*u)^y)*(N*u)^(-(1+y)/2 + t*(log(N*u)/4 - const1)),1);\
return(intest);
}
ttermmagnitude3(X,imax,jmax) = {
N=floor(sqrt(X/4/Pi)); n0=N/2;
print([floor(log(n0)),floor(log(n0)^2/4)]);
preterm = (n0^0.5);
\\preterm = 1;
estinit1 = vecsum(vector(jmax,j,(abs(log(1/n0))^(imax-1)/factorial(imax-1))*((abs(log(1/n0))^2/4)^(j-1)/factorial(j-1))*(0.66^(imax-1))*(0.2^(j-1))));
estinit2 = vecsum(vector(imax,i,(abs(log(1/n0))^(i-1)/factorial(i-1))*((abs(log(1/n0))^2/4)^(jmax-1)/factorial(jmax-1))*(0.66^(i-1))*(0.2^(jmax-1))));
est1 = estinit1 + N*intnum(u=1/N,1,(N*u)^(-1/2)*vecsum(vector(jmax,j,(abs(log(N*u/n0))^(imax-1)/factorial(imax-1))*((abs(log(N*u/n0))^2/4)^(j-1)/factorial(j-1))*(0.66^(imax-1))*(0.2^(j-1)))));
est2 = estinit2 + N*intnum(u=1/N,1,(N*u)^(-1/2)*vecsum(vector(imax,i,(abs(log(N*u/n0))^(i-1)/factorial(i-1))*((abs(log(N*u/n0))^2/4)^(jmax-1)/factorial(jmax-1))*(0.66^(i-1))*(0.2^(jmax-1)))));
print(preterm*(est1 + est2));
}
storedsums(X,taylorterms=50) = {
N = floor(sqrt(X/(4*Pi)));
H=N; n0 = H/2;
X=X+1/2;
expo = (-1+I*X)/2;
finalmat = sum(h=1,H,h^expo*matet(log(h/n0),taylorterms,taylorterms));
X=X-1/2;
return([n0,finalmat]);
}
abbeff_multieval_symmetric_rectangle(X,y,t,num=5,storedsumsdata)={
X=X+1/2;
N = floor(sqrt((X+Pi*t/4)/(4*Pi)));
thtarr = I*vector(num-1,v,-1/2 + (v-1)/(num-1))/2;
zarr = vector(num-1,v,y+(1-y)*(v-1)/(num-1))/2;
s = (1+I*X)/2; ss_expo = conj(s);
\\x lower constant, y upper constant, x upper constant and output to be attached in reverse order, y lower constant and output to be attached in reverse order
sarr = concat([adds(-zarr,s-I/4),adds(thtarr,s-1/2),adds(zarr,s-(1+y)/2+I/4),adds(-thtarr,s-y/2)]);
npoints = #sarr;
[n0,matb]=storedsumsdata;
logtn0 = (t/2)*log(n0); logtn0 = vector(npoints,j,logtn0);
vecbt = Vecrev(mattovect(matb,t));
bexpo = vector(npoints,v,-(t/2)*alpha1(1-sarr[v]) - (1-sarr[v]) + ss_expo);
aexpo = vector(npoints,v,-(t/2)*alpha1(conj(sarr[v])) - conj(sarr[v]) + ss_expo);
afac = vector(npoints,v,exp((t/4)*(alpha1(sarr[v])^2 - alpha1(1-sarr[v])^2))*H01(sarr[v])/H01(1-sarr[v]));
bsums = vecmul(n0^(bexpo+logtn0/2),vechorner(vecbt,bexpo+logtn0));
asums = conj(vecmul(n0^(aexpo+logtn0/2),vechorner(vecbt,aexpo+logtn0)));
ests = bsums+vecmul(afac,asums);
wind_num = (vecsum(arg(vecdiv(ests[1..npoints-1],ests[2..npoints])))+arg(ests[npoints]/ests[1]))/2/Pi;
minmodabb = vecmin(abs(ests));
detailed = matconcat([-I*adds(2.0*sarr,-1.0);ests])~;
X=X-1/2;
return([wind_num,minmodabb,detailed]);
}
\\X=10^3;y=0.2;t=0.2;num=3;taylorterms=60;
\\ttermmagnitude3(X,taylorterms,taylorterms);
\\s=getwalltime();storedsumsdata = storedsums(X,taylorterms);e=getwalltime();print(e-s);
\\[wind_num,minmodabb,detailed] = abbeff_multieval_symmetric_rectangle(X,y,t,num,storedsumsdata);
\\printm(detailed);
\\print(abbeff(X,y,t));\\print(abbeff(X,y,t)-detailed[1][3]);
\\print(abbeff(X+1,1,t));\\print(abbeff(X+1,1,t)-detailed[5][3]);
\\print(abbeff(X,y,t)-detailed[1,2]); print(abbeff(X+1,1,t)-detailed[5,2]); print(abbeff(X,1,t)-detailed[3,2]); print(abbeff(X+1,y,t)-detailed[7,2]);
detailsfile = "~/X_6_10pow10_plus_155019_y0_0.2_t_0.2_detailed_barrier_output_t_agnostic.csv";
summaryfile = "~/X_6_10pow10_plus_155019_y0_0.2_t_0.2_barrier_summary_t_agnostic.csv";
tstart=0.0; tend=0.2; X=6*10^10+155019; y0=0.2;
taylorterms = 50;
st=getwalltime();
storedsumsdata = storedsums(X,taylorterms);
\\ddzvec = ddzboundvec(X+1,y0); ddtvec = ddtboundvec(X+1,y0);
en=getwalltime();
print(en-st);
st=getwalltime();
N = floor(sqrt(X/(4*Pi)));
t=tstart;
counter=0;
while(t<=tend,{
counter=counter+1;
\\Dtabb = ceil(ddt_abbeff_Nbound(N,y0,t)); Dxabb = ceil(ddx_abbeff_Nbound(N,y0,t));
\\Dzabb = ceil(vectoval(ddzvec,t)); Dtabb = ceil(vectoval(ddtvec,t));
Dzabb = ceil(ddzboundint(X+1,y0,t)); Dtabb = ceil(ddtboundint(X+1,y0,t));
sidemesh = Dzabb;
rectmesh = 4*sidemesh - 4;
[wind_num,minmodabb,detailed] = abbeff_multieval_symmetric_rectangle(X,y0,t,sidemesh,storedsumsdata);
\\for(i=1,rectmesh,write(detailsfile,detailed[i,]));
t_summary = [counter,t,Dtabb,Dzabb,wind_num,minmodabb,rectmesh];
print(t_summary); write(summaryfile,t_summary);
t = t + (minmodabb-0.5)/Dtabb;
});
en=getwalltime();
print(en-st);