# ************************************************************************************************************* # Numerical computation of the density function using the residue sum formula # ************************************************************************************************************* # The function \Phi_a(z) in the paper. # The following should work whenever Re(z)>0 and Im(z)<=0. Phi:=(a,z)->-exp(-Pi/(2*a)*I)/a*(GAMMA(-1/a)-GAMMA(-1/a,-I*z)); # The following works for z on the negative imaginary axis. Phi_ni:=(a,z)->-exp(-Pi/(2*a)*I)/a*(GAMMA(-1/a)-conjugate(GAMMA(-1/a,-I*z))); # Next we implement the computation of the derivatives of the zeros of \Phi_a(z) w.r.t. the variable a. # For this we use some numerical integration. # Constants for Gaussian quadrature. Digits:=40: X:=[0.99470046749582496629807708672516631371, 0.97228751153661628803899420776730417254, 0.93281560119391587194023394885619656619, 0.87770220417750151694755059742372113417, 0.80893812220132187422333588202439550949, 0.72900838882861369317120972149178878677, 0.64080177538962945661523025073024805324, 0.54750625491881872009265966771247903156, 0.45249374508118127990734033228752096843, 0.35919822461037054338476974926975194675, 0.27099161117138630682879027850821121323, 0.19106187779867812577666411797560449050, 0.12229779582249848305244940257627886582, 0.067184398806084128059766051143803433806, 0.027712488463383711961005792232695827454, 0.0052995325041750337019229132748336862868]: W:=[0.027152459411754094851780572456018103512, 0.062253523938647892862843836994377694275, 0.095158511682492784809925107602246226355, 0.12462897125553387205247628219201642014, 0.14959598881657673208150173054747854897, 0.16915651939500253818931207903035996221, 0.18260341504492358886676366796921993938, 0.18945061045506849628539672320828310514, 0.18945061045506849628539672320828310514, 0.18260341504492358886676366796921993938, 0.16915651939500253818931207903035996221, 0.14959598881657673208150173054747854897, 0.12462897125553387205247628219201642014, 0.095158511682492784809925107602246226355, 0.062253523938647892862843836994377694275, 0.027152459411754094851780572456018103512]: # Compute the integral of f(x) over x in [x1,x2], using Gaussian quadrature with N*16 evaluations. genintegral:=proc(f,x1,x2,N) global X,W; local h,j,k,s; s:=0.0; h:=evalf((x2-x1)/N); for j from 0 to N-1 do for k from 1 to 16 do # print("xpt:",x1+h*(j+X[k])); # print("Value:",f(x1+h*(j+X[k]))); s:=evalf(s+W[k]*f(x1+h*(j+X[k]))); # print("Now s is:",s); end do; end do; evalf(h*0.5*s); end proc: # Return a 'primitive function' to exp(-u)*u^(-1-1/a)*log(u) obtained by integrating by parts N times. # In more precise terms we return a 2-vector F of functions such that F[1]+int(F[2]) is a primitive function. exppowlog_primfcn:=proc(a,N) local F1,F2,F3,n; F1:=0; F2:=u^(-1-1/a)*log(u); F3:=exp(-u); for n from 1 to N do F1:=F1+F2*int(F3,u); F2:=-diff(F2,u); F3:=int(F3,u); end do; return([unapply(F1,u),unapply(F2*F3,u)]); end proc: # Compute the integral of u -> exp(-u)*u^(-1-1/a)*log(u) from -I*z to "infinity" where z satisfies # Re(z)>=0 and Im(z)<=0 (thus Re(-I*z)<=0, Im(-I*z)<=0). Note that as contour we can take any ray starting # at -I*z and moving in a direction exp(I*omega) with -Pi/2<=omega<=0. # We return a vector [i,e] where i is the numerical value of the integral and e is a guesstimate of # the absolute error. # # The global variable eplF should be initialized by a call eplF:=exppowlog_primfcn(a,N) before using # this function. # inclgamder:=proc(a,z) global eplF; local fcn,omega,i1,i2,phsomega,N,NN,notdone,pr,i3,F,zalimit; if Digits<=20 then zalimit:=100; else zalimit:=450; end if; # The latter seems ok for Digits=30. if(abs(z)>=zalimit) then # Compute much quicker using integration by parts; we have checked experimentally that the # integral over eplF[2] is negligible; hence just skip here! # (We note that this method could be utilized quite a bit more, to gain further efficiency.) i3:=-eplF[1](-I*z); return([i3,0.0]); end if; fcn:=u->exp(-u)*u^(-1-1/a)*log(u); pr:=10^(3-Digits); #We take this as our 'worst acceptable precision'. if Digits<=20 then N:=8: NN:=5: else N:=13: NN:=10: end if; # The latter seems ok when Digits=30. notdone:=true; while(notdone) do omega:=-Pi/4; phsomega:=evalf(exp(I*omega)); i1:=genintegral(t->fcn(-I*z+phsomega*t)*phsomega,0,8,NN) + genintegral(t->fcn(-I*z+phsomega*t)*phsomega,8,15*N^0.8,N); omega:=-Pi/5; phsomega:=evalf(exp(I*omega)); i2:=genintegral(t->fcn(-I*z+phsomega*t)*phsomega,0,8,NN) + genintegral(t->fcn(-I*z+phsomega*t)*phsomega,8,15*N^0.8,N); if abs(i1-i2)>pr*abs(i1) then N:=N+5: NN:=NN+5: print(i1); print(i2); printf("rel.err: %g. Now raising to %d,%d.\n",abs(i1-i2)/abs(i1),N,NN); else notdone:=false; end if; end do; return([i1,abs(i1-i2)]); end proc: # For testing purposes: Compute the integral of u -> exp(-u)*u^(-1-1/a) from -I*z to "infinity"; # i.e. completely analogous to the function above but without the "log(u)". inclgam:=proc(a,z) local fcn,omega,i1,i2,phsomega,N,NN,notdone,pr,F,zalimit; fcn:=u->exp(-u)*u^(-1-1/a); pr:=10^(3-Digits); #We take this as our 'worst acceptable precision'. if Digits<=20 then N:=8: NN:=5: else N:=13: NN:=10: end if; # The latter seems ok when Digits=30. notdone:=true; while(notdone) do omega:=-Pi/4; phsomega:=evalf(exp(I*omega)); i1:=genintegral(t->fcn(-I*z+phsomega*t)*phsomega,0,8,NN) + genintegral(t->fcn(-I*z+phsomega*t)*phsomega,8,15*N^0.8,N); omega:=-Pi/5; phsomega:=evalf(exp(I*omega)); i2:=genintegral(t->fcn(-I*z+phsomega*t)*phsomega,0,8,NN) + genintegral(t->fcn(-I*z+phsomega*t)*phsomega,8,15*N^0.8,N); if abs(i1-i2)>pr*abs(i1) then N:=N+5: NN:=NN+5: print(i1); print(i2); printf("rel.err: %g. Now raising to %d,%d.\n",abs(i1-i2)/abs(i1),N,NN); else notdone:=false; end if; end do; return([i1,abs(i1-i2)]); end proc: Phi_test2:=proc(a,z) local v; v:=inclgam(a,z); evalf(-exp(-Pi/(2*a)*I)/a*(GAMMA(-1/a)-v[1])); end proc: # The function zerod(a,z): Assuming that z is a zero of \Phi_a(.), then zerod(a,z) gives # the derivative of this zero wrt a. zerod:=proc(a,z) local v,w,pr; v:=inclgamder(a,z); w:=GAMMA(-1/a)*Psi(-1/a)-v[1]; pr:=10^(3-Digits); #We take this as our 'worst acceptable precision'. if v[2]/abs(w)>pr then error("Too bad precision in zerod(%g,%g): %g.\n",a,z,v[2]/abs(w)); end if; -z^(1+1/a)*exp(-I*z-Pi/(2*a)*I)/a^2*w; end proc: # We use the following alternative method when z is near 0. zerod2:=proc(a,z) local v1,v2,pr; # Compute int(exp(I*z*t)*t^(1-1/a)*log(t),t=0..1), after first substituting t=exp(-u): v1:=-genintegral(u->exp(I*z*exp(-u)+(1/a-2)*u)*u,0,60,40); v2:=-genintegral(u->exp(I*z*exp(-u)+(1/a-2)*u)*u,0,60,45); pr:=10^(3-Digits); #We take this as our 'worst acceptable precision'. if abs(v1-v2)/abs(v2)>pr then error("Too bad precision in zerod2(%g,%g): %g.\n",a,z,abs(v1-v2)/abs(v2)); end if; z/(a-1)*(1-z^2*exp(-I*z)*v2); end proc: # For testing purposes: \Phi_a(z)... # We use the following alternative method when z is near 0. Phi_test:=proc(a,z) local v1,v2,pr; # Compute int(exp(I*z*t)*t^(1-1/a),t=0..1). v1:=genintegral(u->exp(I*z*exp(-u)+(1/a-2)*u),0,60,40); v2:=genintegral(u->exp(I*z*exp(-u)+(1/a-2)*u),0,60,45); pr:=10^(3-Digits); #We take this as our 'worst acceptable precision'. if abs(v1-v2)/abs(v2)>pr then error("Too bad precision in Phi_test:",a,z,abs(v1-v2)/abs(v2)); end if; # The following formula follows from (4.7) in the paper by integration by part. exp(I*z)*z^(-1/a)*(1-I*z*a/(a-1))-a/(a-1)*z^(2-1/a)*v2; end proc: # Functions for locating the poles of \Psi_a(z), i.e. the zeros of \Phi_a(z) # ************************************************************************** # Assuming that z0 is a good first approximation of a zero of \Phi_a(z), zoom in on this zero using Newton's method. # Continue until we get abs(\Phi_a(z))=pr) do # Next z := z-Phi(z)/Phi'(z): s:=evalf(p*a*z^(1+1/a)*exp(-I*z)); if(abs(s)>tol) then error("Too big jump: ",a,z0); end if; if niflag then s:=I*Im(s); end if; z:=z+s; old2p:=old1p; old1p:=abs(p); p:=evalf(myPhi(a,z)); if old2p>0.0 and old2p<10*abs(p) then printf("Problem in zoomzero; don't seem to get below %.4g.\n",abs(p)); return z; end if; end do; z; end proc: # For given a, compute the unique z along the negative imaginary axis satisfying \Phi_a(z)=0. findnegimagzero:=proc(a,pr) local eta,y,y1,y2,ym,etam; # Let eta(a,y)=(-I*y)^(1/a)*Phi(a,-I*y), which is real; cf. the proof of Proposition 3 in the paper. eta:=(a,y)->-Re(y^(1/a)*exp(-Pi*I/a)/a*(GAMMA(-1/a)-conjugate(GAMMA(-1/a,-y)))); y:=1.0; while(evalf(eta(a,y))<0) do y:=0.5*y; end do; if y=1.0 then while(evalf(eta(a,y))>0) do y:=2.0*y; end do; y1:=0.5*y; y2:=y; else y1:=y; y2:=2*y; end if; ym:=0.5*(y1+y2); etam:=evalf(eta(a,ym)); while(y2-y1>0.0001 or abs(etam)>0.01) do if etam>0 then y1:=ym; else y2:=ym; end if; ym:=0.5*(y1+y2); etam:=evalf(eta(a,ym)); end do; zoomzero(a,-I*ym,pr,1.0); end proc: # For a>1 and n>=1: Compute an approximate value of the zero zeta_n, using the asymptotic formula. asymptnzeta:=proc(a,n) local g,Y; g:=log(abs(GAMMA(-a^(-1)))); Y:=fsolve(y-0.5*(1+a^(-1))*log((2*Pi*n)^2+y^2)=g,y,0..infinity); evalf(2*n*Pi+(1-a^(-1))*Pi/2-(1+a^(-1))*arctan(Y/(2*Pi*n))-I*Y); end proc: # Find the first N+1 zeros of \Phi_a(z) with Re(z)>=0. Return a list of these in order of increasing Re(z). # pr is the desired precision. findzeros:=proc(a,N,pr) local zl,z,j; zl:=[findnegimagzero(a,pr)]; for j from 1 to N do if j mod 1000 = 0 then printf("%d\n",j); end if; z:=zoomzero(a,asymptnzeta(a,j), pr,0.5); zl:=[op(zl),z]; end do; zl; end proc: # Functions for computing Prob(sigma_{T_j} > a/2) and the corresponding density # ******************************************************************************* # Compute Prob(sigma_{T_j} > a/2) using the formula (A.19) in the paper. # zl is a list containing the poles which we add over, as returned by findzeros(...). computeP:=proc(a,zl) local s,t,j; s:=a*exp(-2*I*zl[1]); print(s); for j from 2 to nops(zl) do t:=2*a*Re(exp(-2*I*zl[j])); s:=s+t; if (j=nops(zl)) then print(s,t); end if; end do; s; end proc: # Compute the corresponding density, using the formula (A.21) in the paper. # (Normalization: Density w.r.t. dc, where c=a/2.) # Return also a rough heuristic "measure" or "bound" on the error: # The sum of the absolute values of all terms of index j>nops(zl)/2. computedensity:=proc(a,zl,writeflag) global zerod,zerod2,eplF; local zd,s,j,t,er; eplF:=exppowlog_primfcn(a,15): if abs(zl[1])<2 then zd:=I*Im(zerod2(a,zl[1])); else zd:=I*Im(zerod(a,zl[1])); end if; s:=2*exp(-2*I*zl[1])*(2*a*I*zd-1); if writeflag then printf("%.15g\n",s); end if; er:=0.0; for j from 2 to nops(zl) do if j mod 1000 = 0 then printf("%d\n",j); end if; t:=2*Re(2*exp(-2*I*zl[j])*(2*a*I*zerod(a+1e-17,zl[j])-1)); s:=evalf(s+t); if 2*j>nops(zl) then er:=evalf(er+abs(t)); end if; if writeflag then printf("%d : %.15g %.5g\n",j,s,t); end if; end do; [s,er]; end proc: # For each a-value in the given range: Compute the zeros of \Phi_a(.); then compute the density at c=a/2! # The zeros are computed to precision pr. Of course the precision obtained for the density is in general much worse. computedensity_graph:=proc(a_start,a_end,a_step,N,pr) local zl,a,d,t,fd,j; t:=[]; for a from a_start to a_end by a_step do zl:=findzeros(a,N,pr); d:=computedensity(a,zl,false); printf("c=%f: %.12f # 'error': %.2g\n",0.5*a,d[1],d[2]); t:=[op(t),[0.5*a,op(d)]]; end do; fd:=fopen("datafiles/density.dat",WRITE); for j in t do fprintf(fd,"%.6f, %.12f # 'error': %.2g\n",j[1],j[2],j[3]); end do; fclose(fd); t; end proc: # Applications. # ************* # We used the following to generate the data for the graph in Figure 1 in our paper. Digits:=20: computedensity_graph(1.01,5.0,0.01,400,1e-14): # Now test some a near 1 and some large a, and compare with the constants K_1 and K_2 in our Corollary 1. # Recall that we prove K_1=39.47841... in the paper. a:=1.01: zl:=findzeros(a,200,6e-14): d:=computedensity(a,zl,false); printf("Rel error, guesstimate: %g\n",d[2]/(2-d[1])): K1_appr:=(2-d[1])/((a-1)/2)^2; # Gives 38.44... a:=1.001: zl:=findzeros(a,200,6e-14): d:=computedensity(a,zl,false); printf("Rel error, guesstimate: %g\n",d[2]/(2-d[1])): K1_appr:=(2-d[1])/((a-1)/2)^2; # Gives 39.373... Digits:=25: a:=1.0001: zl:=findzeros(a,200,5e-18): d:=computedensity(a,zl,false); printf("Rel error, guesstimate: %g\n",d[2]/(2-d[1])): K1_appr:=(2-d[1])/((a-1)/2)^2; # Gives 39.467... Digits:=25: a:=1.00005: zl:=findzeros(a,500,5e-17): d:=computedensity(a,zl,false); printf("Rel error, guesstimate: %g\n",d[2]/(2-d[1])): K1_appr:=(2-d[1])/((a-1)/2)^2; # Gives 39.4731... Digits:=30: a:=1.00001: zl:=findzeros(a,500,5e-22): d:=computedensity(a,zl,false); printf("Rel error, guesstimate: %g\n",d[2]/(2-d[1])): K1_appr:=(2-d[1])/((a-1)/2)^2; # Gives 39.4773... Digits:=60: a:=1.000005: zl:=findzeros(a,2000,1e-40): Digits:=25: d:=computedensity(a,zl,false); printf("Rel error, guesstimate: %g\n",d[2]/(2-d[1])): K1_appr:=(2-d[1])/((a-1)/2)^2; # Gives 39.47789... Digits:=60: a:=1.000001: zl:=findzeros(a,2000,1e-40): Digits:=25: d:=computedensity(a,zl,false); printf("Rel error, guesstimate: %g\n",d[2]/(2-d[1])): K1_appr:=(2-d[1])/((a-1)/2)^2; # Gives 39.47831... # Recall that we prove K_2=0.822467... in the paper. Digits:=20: a:=35.0: zl:=findzeros(a,10000,1e-15): d:=computedensity(a,zl,false); printf("Rel error, guesstimate: %g\n",d[2]/d[1]): K2_appr:=d[1]*(a/2)^3; # Gives 0.82181... a:=50.0: zl:=findzeros(a,10000,1e-15): d:=computedensity(a,zl,false); printf("Rel error, guesstimate: %g\n",d[2]/d[1]): K2_appr:=d[1]*(a/2)^3; # Gives 0.82216... Digits:=30: a:=100.0: zl:=findzeros(a,50000,1e-15): Digits:=20: d:=computedensity(a,zl,false); printf("Rel error, guesstimate: %g\n",d[2]/d[1]): K2_appr:=d[1]*(a/2)^3; # Gives 0.82239...