(* -*-Mathematica-*- code *) (* Time-stamp: *) (* * Code for computation of "sieve functions" described in the book A * HIGHER-DIMENSIONAL SIEVE METHOD WITH WITH PROCEDURES FOR COMPUTING SIEVE * FUNCTIONS, by Harold G. Diamond, H. Halberstam, and William F. Galway. * * Please note that we cannot guarantee the accuracy of any results * computed with this package. Questions about or comments on this package * should be addressed to sievetheorybook@math.uiuc.edu. *) (* :Title: SieveFunctions.m -- utilities for computing functions used in sieve theory. *) (* :Author: William F. Galway *) (* :Copyright: 2009 by William F. Galway *) (* :Package Version: 0.9 *) (* :Mathematica Version: 5.2 *) (********************************************************************************* * KNOWN BUGS AND ISSUES: * -- When kappa is large and u is much larger than kappa, e.g. kappa=10, u=100, * then evaluating UpperSieveFunc[kappa,u], LowerSieveFunc[kappa,u], * generates messages that * "$RecursionLimit::reclim: Recursion depth of 256 exceeded", * at which point Mathematica may crash. *********************************************************************************) BeginPackage["SieveFunctions`"] SieveFunctions::usage = "SieveFunctions is not a function, but Options[SieveFunctions] gives the default options for the functions UpperSieveFunc[...], LowerSieveFunc[...], SieveP[...], SieveQ[...]." AlphaBeta::usage = "AlphaBeta is an option for many of the functions in the SieveFunction package which should have one of the two forms AlphaBeta->{alpha,beta}, where alpha>=beta are real numbers, or the form AlphaBeta->Automatic." ComputationMethod::usage = "ComputationMethod is an option for many of the functions in the SieveFunctions package that specifies the preferred method for computing those functions. Choices include ComputationMethod->IntegralRep; ComputationMethod->DDE; and the default ComputationMethod->Automatic, which specifies, in effect, that a mixture of ComputationMethod->IntegralRep and ComputationMethod->DDE should be used in such a way as to attempt optimally to balance accuracy against speed. See the documentation for IntegralRep, RightIntegralRep, LeftIntegralRep, and DDE for further details. Also, note that the ComputationMethod option will be ignored when there is no corresponding method for computing a given function." IntegralRep::usage = "IntegralRep is a choice for the option ComputationMethod used by several functions in the SieveFunctions package, ComputationMethod->IntegralRep specifies that the function should be computed, when it makes sense to do so, using an integral representation for that function (usually a contour integral). See also the documentation for LeftIntegralRep and RightIntegralRep." LeftIntegralRep::usage = "LeftIntegralRep is a choice for the option ComputationMethod used by several functions in the SieveFunctions package, ComputationMethod->LeftIntegralRep applies to the computation of sigmaFunc and jFunc, and specifies that the function should be computed using a contour integral representation for that function where the contour passes to the left of the origin." RightIntegralRep::usage = "RightIntegralRep is a choice for the option ComputationMethod used by several functions in the SieveFunctions package, ComputationMethod->RightIntegralRep applies to the computation of sigmaFunc and jFunc, and specifies that the function should be computed using a contour integral representation for that function where the contour passes to the right of the origin." DDE::usage = "DDE is a choice for the option ComputationMethod, ComputationMethod->DDE specifies that the function should be computed by numerical integration of the delay differential equation satisfied by that function." Ein::usage = "Ein[z] = Integrate[(1-Exp[-s])/s,{s,0,z}] and is an entire function (an exponential integral). It is implemented as EulerGamma + Log[z] + ExpIntegralE[1,z]." jLaplace::usage = "jLaplace[kappa,z] is the Laplace transform of the Grupp-Richert \"j\" function: jFunc[kappa,u]." jIntegrand::usage = "jIntegrand[kappa,u,z] is the integrand used to compute jFunc[kappa,u] by the inverse Laplace transform integral formula." jSadIter::usage = "jSadIter[kappa,u,z0] refines an approximation z0 to the negative real saddle-point of jIntegrand[kappa,u,z]." jSad::usage = "jSad[kappa,u] returns the negative real saddle-point of jIntegrand[kappa,u,z]. jSad[kappa,u,m] returns the saddle-point with negative real part, and imaginary part near 2*Pi*m." jIntLeftPath::usage = "jIntLeftPath[kappa,u,digits] returns points {z0,z1, ...} defining a path in the complex plane to be used by NIntegrate[jIntegrand[kappa,u,z], {z, z0, z1, ...}] in approximating jFuncInternal[LeftIntegralRep,...] to digits worth of accuracy." jIntRightPath::usage = "jIntRightPath[kappa,u,digits] returns points {z0,z1, ...} defining a path in the complex plane to be used by NIntegrate[jIntegrand[kappa,u,z], {z, z0, z1, ...}] in approximating jFuncInternal[RightIntegralRep,...] to digits worth of accuracy." jFuncSeg0::usage = "jFuncSeg0[kappa,u] returns jFunc[kappa,u], valid for 0= 1; jFunc[...] also takes optional arguments, see Options[jFunc] for further details." sigmaFunc::usage = "sigmaFunc[kappa,u] = jFunc[kappa, u/2]. See the documentation for jFunc. sigmaFunc[...] takes the same optional arguments as jFunc." AdjointSieveP::usage = "AdjointSieveP[kappa,u] returns the value of the \"p\" sieve function of u with parameter kappa>0. When kappa<=1 it is important that u>0." qIntegrand::usage = "qIntegrand[kappa,u,z] is the integrand of a contour-integral used in the computation of AdjointSieveQ[kappa,u]." LoglikeXi::usage = "LoglikeXi[t] is defined to be the solution \"x\" to t=(e^x-1)/x." DLogqIntegrand::usage = "DLogqIntegrand[kappa,u,z] is the logarithmic derivitive with respect to z of qIntegrand[kappa,u,z]." AdjointSieveQ::usage = "AdjointSieveQ[kappa,u] returns the value of the \"q\" sieve function with parameters kappa>0, u>0. When 2*kappa is an integer the result is a polynomial in u, so u <= 0 is acceptable in that case." qPoly::usage = "qPoly[deg,u] returns AdjointSieveQ[(deg+1)/2, u], where deg is a non-negative integer." qPolyCoef::usage = "qPolyCoef[deg,n] is defined so that qPoly[deg,u] = Sum[(-1)^n*Binomial[deg,n]*qPolyCoef[deg,n]*u^(deg-n), {n,0,deg}]." qDDEcheck::usage = "qDDEcheck[kappa] returns D[u*AdjointSieveQ[kappa,u], u] - kappa*(AdjointSieveQ[kappa,u]+AdjointSieveQ[kappa,u+1]), which should yield (or simplify to) zero. (Only works when 2*kappa is an integer.)" qIntegralCheck::usage = "qIntegralCheck[kappa,u] returns (u+1)*AdjointSieveQ[kappa,u+1] - u*AdjointSieveQ[kappa,u] - kappa*Integrate[AdjointSieveQ[kappa,t],{t,u,u+2}], which should return zero (or near zero, when 2*kappa is not an integer)." CriticalRatio::usage = "CriticalRatio is the value of r := u/kappa at which the \"interesting\" saddle point of qIntegrand[kappa,u,z] leaves the real axis. (The saddle point is real for r>CriticalRatio and complex-valued for r 0." qSad::usage = "qSad[u/kappa] gives a good approximation to the \"interesting\" saddle point of qIntegrand[kappa,u,z]." qIntPath::usage = "qIntPath[kappa,u] returns points {z0,z1, ...} defining a path in the complex plane to be used by NIntegrate[qIntegrand[kappa,u,z], {z, z0, z1, ...}] in computing AdjointSieveQ[...]." PiSieveFunc::usage = "Can be called as PiSieveFunc[kappa,u,v], or as PiSieveFunc[kappa,u] := PiSieveFunc[kappa,u,u]." PiTildeSieveFunc::usage = "PiTildeSieveFunc[kappa,u] := PiSieveFunc[kappa,u,u-1]." XiSieveFunc::usage = "Can be called as XiSieveFunc[kappa,u,v], or as XiSieveFunc[kappa,u] := XiSieveFunc[kappa,u,u]." XiTildeSieveFunc::usage = "XiTildeSieveFunc[kappa,u] := XiSieveFunc[kappa,u,u-1]." alphaSearchSieveFunc::usage = "One of the roots of alphaSearchSieveFunc[kappa,u]==0 is u=alphaParam[kappa]; this is the function l(u) in the notation of Diamond and Halberstam." betaSearchSieveFunc::usage = "betaSearchSieveFunc[kappa,alpha,u]:=IwaniecInnerQ[kappa,u,alpha,u]/u^(2*kappa-1), where the denominator is chosen to keep the magnitude of betaSearchSieveFunc[...,u] bounded as u->Infinity. betaSearchSieveFunc[kappa,alphaParam[kapp],betaParam[kappa]] = 0, in theory." ChopUpInterval::usage = "ChopUpInterval[u0,u1,avec] returns a list {u0,...,u1} where the entries are in increasing order and intervening entries are congruent modulo 1 to one of the avec[[k]] and the entries congruent modulo 1 to avec[[k]] are >= avec[[k]]." UpperSieveFunc::usage = "UpperSieveFunc[kappa,u] is the sieve upper bound function F(u), of dimension kappa, satisfying its defining delay differential equation up to the limits of accuracy of the numerical method used." UpperSieveFuncInternal::usage = "Internal function." UpperSieveFuncSegment::usage = "Internal function used to implement UpperSieveFunc." LowerSieveFunc::usage = "LowerSieveFunc[kappa,u] is the sieve lower bound function f(u), of dimension kappa, satisfying its defining delay differential equation up to the limits of accuracy of the numerical method used." LowerSieveFuncInternal::usage = "Internal function." LowerSieveFuncSegment::usage = "Internal function used to implement LowerSieveFunc." SieveP::usage = "SieveP[kappa,u] == UpperSieveFunc[kappa,u] + LowerSieveFunc[kappa,u]. Optional arguments are the same as for UpperSieveFunc[...] and LowerSieveFunc[...]." SieveQ::usage = "SieveQ[kappa,u] == UpperSieveFunc[kappa,u] - LowerSieveFunc[kappa,u]. Optional arguments are the same as for UpperSieveFunc[...] and LowerSieveFunc[...]." IwaniecInnerP::usage = "IwaniecInnerP[kappa,u] is an expression involving an integral that should be a constant as a function of u, for u >= alpha >= beta >= 2. Note that alpha and beta can be passed as optional arguments using the form IwaniecInnerP[kappa,u,AlphaBeta->{alpha,beta}]." IwaniecInnerQ::usage = "IwaniecInnerQ[kappa,u] is an expression involving an integral that should be a constant as a function of u, for u >= alpha >= beta >= 2. Note that alpha and beta can be passed as optional arguments using the form IwaniecInnerQ[kappa,u,AlphaBeta->{alpha,beta}]." IwaniecInnerPintegrand::usage = "IwaniecInnerPintegrand[kappa,u,alpha,beta] is the integrand evaluated by IwaniecInnerP[...]." IwaniecInnerQintegrand::usage = "IwaniecInnerQintegrand[kappa,u,alpha,beta] is the integrand evaluated by IwaniecInnerQ[...]." alphaInterpolatingFunction::usage = "alphaInterpolatingFunction is an InterpolatingFunction for which alphaInterpolatingFunction[kappa] approximates alphaParam[kappa]." betaInterpolatingFunction::usage = "betaInterpolatingFunction is an InterpolatingFunction for which betaInterpolatingFunction[kappa] approximates betaParam[kappa]." OldalphaInterpolatingFunction::usage = "OldalphaInterpolatingFunction is a previously calculated version of alphaInterpolatingFunction, retained for comparison and \"quality control\" tests." OldbetaInterpolatingFunction::usage = "OldbetaInterpolatingFunction is a previously calculated version of betaaInterpolatingFunction, retained for comparison and \"quality control\" tests." AlphaBetaExtrapolate::usage = "AlphaBetaExtrapolate extends the range of alphaInterpolatingFunction and betaInterpolatingFunction by 0.5." AlphaBetaParams::usage = "AlphaBetaParams[kappa] returns the \"proper\" values of {alpha,beta} to pass to UpperSieveFunc[kappa,u,alpha,beta], LowerSieveFunc[kappa,u,alpha,beta], etc." alphaParam::usage = "alphaParam[kappa] returns the \"proper\" value of alpha to pass to UpperSieveFunc[kappa,u], LowerSieveFunc[kappa,u], etc." betaParam::usage = "betaParam[kappa] returns the \"proper\" value of beta to pass to UpperSieveFunc[kappa,u], LowerSieveFunc[kappa,u], etc." Begin["`Private`"] LoglikeXi[t_?NumberQ] := Module[{xi}, (* Our starting guess for the solution is somewhat arbitrarily * split into 3 different cases, based on the size of "t". *) xi/.FindRoot[t==(Exp[xi]-1)/xi, {xi,Piecewise[{{Log[t]+Log[1+Log[t]], t>2}, {-1/t, t<0.5}}, 2*t-1]}][[1]] ] /; t>0; (* * Note that Mathematica is not clever enough to deduce that all * singularities and branch-cuts cancel in this definition, giving an * entire function of z. *) Ein[z_] := EulerGamma + Log[z] + ExpIntegralE[1,z]; jLaplace[kappa_,z_] := Exp[-kappa*Ein[z]]/z; jIntegrand[kappa_?NumberQ,u_?NumberQ,z_?NumberQ] := Exp[u*z]*jLaplace[kappa,z]; jSadIter[kappa_,u_,z_] := Log[kappa/(1+kappa-u*z)]; jSad[kappa_,u_] := jSad[kappa,u,0]; (* * This would be faster if we'd switch to something like * FindRoot[u*z + kappa*Exp[-z] - kappa-1 == 0, ...] after a few * iterations of jSadIter. *) jSad[kappa_,u_,m_] := FixedPoint[(2*Pi*I*m + jSadIter[kappa,u,#])&, -Log[1.0+1/kappa]]; (* * jIntLeftPath[...] returns a list of consecutive saddle-points of jIntegrand[...], * which lie to the left of the origin, truncated at the first saddle point at which the * integrand is bounded by "eps" as calculated below. *) jIntLeftPath[kappa_,u_,digits_] := Module[{m,path = {jSad[kappa,u,0], jSad[kappa,u,1]}, zsad, eps}, eps = 10.0^(-digits)/2; For[m=2, Abs[jIntegrand[kappa,u, (zsad=jSad[kappa,u,m])]] > eps, m++, path = Append[path, zsad]; ]; (* Include the last saddle point for good measure. *) path = Append[path, zsad]; path ]; octahedralBasisPoint = Exp[2*Pi*I/16]; OctahedralPathPoints = {1, octahedralBasisPoint/Re[octahedralBasisPoint], (I/octahedralBasisPoint)/Re[octahedralBasisPoint]}; (* * jIntRightPath[...] returns a path suitable for jIntegrand[...], this * path passes through the real saddle point to the right of the * origin and then moves to the left of the imaginary axis to pass * through one or more additional saddle points, truncated at the * first saddle point at which the integrand is bounded by "eps" as * calculated below. *) jIntRightPath[kappa_,u_,digits_] := Module[{m, zsad, rslt, eps, val}, If[1+kappa > 10*u, (* * If the saddle point is too far to the right (>10, * roughly), then we don't bother to pass through the saddle * point. *) zsad = 10, (* * Otherwise set zsad to a rough approximation to * its ultimate value. *) zsad = (1+kappa)/u; (* Then we give FindRoot[...] an even better starting * approx. to start off with. *) zsad = z/.FindRoot[u*z-(1+kappa)+kappa*Exp[-z]==0, {z,zsad-Exp[-zsad]*kappa/u}][[1]]; (* zsad is the location of the positive real saddle point of * jIntegrand[...]. *) ]; path = zsad*OctahedralPathPoints; eps = 10.0^(-1-digits)*Abs[jIntegrand[kappa,u,zsad]]; For[m=Ceiling[zsad/(2*Pi)], Abs[jIntegrand[kappa,u, (zsad=jSad[kappa,u,m])]] > eps, m++, path = Append[path, zsad]; ]; path = Append[path, zsad]; path ]; jFuncSeg0[kappa_,u_] := Exp[-kappa*EulerGamma]*u^kappa/Gamma[kappa+1]; (* By convention, jFunc[kappa,u] == 0 for u<=0. *) jFuncInternal[method_,kappa_?NumberQ,u_?NumberQ] := 0 /; u<=0; (* This handles most cases cases with 0IntegralRep, u>1, jFunc is implemented as an * inverse Laplace transform choosing a path based on the relative * sizes of u versus kappa, to "encourage" getting the most accurate * possible result. *) jFuncInternal[IntegralRep,kappa_?NumberQ,u_?NumberQ] := If[u<=kappa, jFuncInternal[RightIntegralRep,kappa,u], jFuncInternal[LeftIntegralRep,kappa,u]]; jFuncInternal[LeftIntegralRep,kappa_?NumberQ,u_?NumberQ] := Module[{z,digits = 8, path}, (* Note "z" plays a symbolic role here. *) path = {z, Sequence @@ jIntLeftPath[kappa, u, digits]}; 1 + Im[NIntegrate[jIntegrand[kappa,u,z], Evaluate[path], AccuracyGoal->digits, MaxRecursion->20] ]/Pi ]; jFuncInternal[RightIntegralRep,kappa_?NumberQ,u_?NumberQ] := Module[{z,digits = 8, path}, (* Note "z" plays a symbolic role here. *) path = {z, Sequence @@ jIntRightPath[kappa, u, digits]}; Im[NIntegrate[jIntegrand[kappa,u,z], Evaluate[path], MaxRecursion->20]]/Pi ]; jFuncInternal[Automatic,kappa_?NumberQ,u_?IntegerQ] := jFuncInternal[IntegralRep,kappa,u]; (* * This handles ComputationMethod->DDE and ComputationMethod->Automatic * when u>1, u not an integer. *) jFuncInternal[method_?(MemberQ[{DDE, Automatic},#]&), kappa_?NumberQ, u_?NumberQ] := Module[{u1=Floor[u]}, u1 = If[u>u1,u1,u1-1]; jFuncViaDDESegment[method,kappa,u1,u1+1][u] ] /; u-1 > 0; (* * Return (and store a cached copy of) an InterpolatingFunction object giving * jFunc[kappa,u] for u1 <= u<= u2. *) jFuncViaDDESegment[method_,kappa_?NumberQ,u1_?IntegerQ,u2_?IntegerQ] := jFuncViaDDESegment[method,kappa,u1,u2] = Module[{u,j}, j/.NDSolve[{j[u1] == jFuncInternal[method,kappa,u1], u*D[j[u],u] == kappa*(j[u]-jFuncInternal[method,kappa,u-1])}, j,{u,u1,u2}, AccuracyGoal->Infinity][[1]] ]; Options[jFunc] = {ComputationMethod->Automatic}; jFunc[kappa_?NumberQ,u_?NumberQ, opts___?OptionQ] := jFuncInternal[ComputationMethod/.{opts}/.Options[jFunc],kappa,u]; Options[sigmaFunc] = {ComputationMethod->Automatic}; sigmaFunc[kappa_,u_,opts___?OptionQ] := jFunc[kappa, u/2, Sequence @@ Join[{opts},Options[sigmaFunc]]]; (* * AdjointSieveP[kappa,u] is implemented as a Laplace transform. *) AdjointSieveP[kappa_?NumberQ,u_?NumberQ] := Module[{t}, NIntegrate[Exp[-u*t - kappa*Ein[t]],{t,0,Infinity}]] /; kappa>0; qIntegrand[kappa_,u_,z_] := Exp[u*z + kappa*(Ein[-z] - 2*Log[z])]; DLogqIntegrand[kappa_, u_, z_] := u - kappa*(1+Exp[z])/z; (* * Note: in the Diamond-Halberstam sieve book, the authors observe * that qPolyCoef[deg,n] is a polynomial in kappa=(deg+1)/2. We don't * use that fact here, since it would give an O(deg^2) running time to * compute qPoly[deg,u], while this approach should give an O(deg) * running time (with, of course, a different and perhaps "worse" * O-constant). *) qPolyCoef[deg_Integer,0] = 1; (* Note that (deg+1)/2 denotes "kappa". *) qPolyCoef[deg_Integer,n_Integer] := (qPolyCoef[deg,n] = Sum[(-1)^m*Binomial[n-1,m]*qPolyCoef[deg,n-1-m]/(m+1), {m,0,n-1}]*(deg+1)/2) /; deg>0; (* * We enclose the following two definitions within a Module[...] * to be REALLY sure that no name conflict occurs with "uvar". *) Module[{uvar}, qPolyInternal[deg_Integer] := qPolyInternal[deg] = Sum[(-1)^n*Binomial[deg,n]*qPolyCoef[deg,n]*uvar^(deg-n), {n,0,deg}]; qPoly[deg_Integer,u_] := (qPolyInternal[deg]/.{uvar->u}) /; deg>=0; ]; CriticalRatio = Module[{r}, r /. FindRoot[r*(Log[r]-1)==1, {r,3.59}]]; (* * Map z to a better approximation of the "interesting" saddle point * of qIntegrand[kappa,u,z], where r = u/kappa. *) qSadNewtonIteration[r_,z_] := ((z-1)*Exp[z]-1)/(Exp[z]-r); qSadApproxBigRatio[r_] := 2/r + 2/r^2; qSadApproxSmallRatio[r_] := (Pi*r)^2/2 + (1-r+r^2)*Pi*I; (* * The following has an error of order (r - CriticalRatio), barely * good enough? *) qSadApproxIntermediateRatio[r_] := Conjugate[ Log[CriticalRatio] + (r-CriticalRatio - Sqrt[(r-CriticalRatio)^2 + 2*CriticalRatio*Log[CriticalRatio]*(r-CriticalRatio)])/CriticalRatio ]; (* * The points at which we switch from one approximation to another * were determined by trial-and-error. Piecewise[...] is new to Mathematica 5.1. *) qSadApprox[r_] := Piecewise[{{qSadApproxSmallRatio[r], r < 0.517}, {qSadApproxBigRatio[r], r > 5.3}}, qSadApproxIntermediateRatio[r] (* Default case. *) ]; (* * Up to 4 Newton iterations suffice for an error bound of 10^(-5) or * so, but Newton's method encounters division by zero for * r=CriticalRatio. *) qSad[r_] := If[Abs[r-CriticalRatio] < 10.0^(-12), qSadApproxIntermediateRatio[r], FixedPoint[qSadNewtonIteration[r,#]&,qSadApprox[r],4]]; (* * Note, the integral should run to -Infinity, but Mathematica gets * confused by that in some cases, so we use a kludge to truncate the * path to a finite range. ALSO NOTE that qIntPath[...] is a bit * dubious for u/kappa near CriticalRatio. *) qIntPath[kappa_, u_] := Module[{zs, zdelt, ztrunc}, zs=qSad[u/kappa]; ztrunc = Min[-1,(Log[u] - kappa*(Exp[-1]+EulerGamma) - 20.0)/u]; If[Im[zs]===0, Return[{zs, zs*(1+I),ztrunc + zs*I}], zdelt = Sqrt[-1/Residue[DLogqIntegrand[kappa,u,z]/(z-zs)^2,{z,zs}]]; zdelt = Im[zs]*zdelt/Im[zdelt]; (* Note that Im[zs-zdelt]==0, to within the limits of numerical accuracy. *) Return[{Re[zs-zdelt], zs+zdelt, Min[Re[zs+zdelt], ztrunc] + I*Max[1,Im[zs+zdelt]]} ]; ]; ]; AdjointSieveQ[kappa_,u_] := qPoly[2*kappa-1,u] /; IntegerQ[2*kappa] && kappa>0; (* WHY do we do the {Complex[x_,y_] -> x + I*y} below? *) AdjointSieveQ[kappa_?NumberQ,u_?NumberQ] := Module[{z, path}, path = {z, Sequence @@ qIntPath[kappa,u]}/.{Complex[x_,y_] -> x + I*y}; Gamma[2*kappa]*Im[NIntegrate[qIntegrand[kappa,u,z],Evaluate[path]]]/Pi ] /; !IntegerQ[2*kappa]; qDDEcheck[kappa_] := Simplify[D[u*AdjointSieveQ[kappa,u], u] - kappa*(AdjointSieveQ[kappa,u]+AdjointSieveQ[kappa,u+1])]; qIntegralCheck[kappa_, u_] := ((u+1)*AdjointSieveQ[kappa,u+1] - u*AdjointSieveQ[kappa,u] - kappa*Integrate[AdjointSieveQ[kappa,t],{t,u,u+2}] // Simplify) /; IntegerQ[2*kappa]; qIntegralCheck[kappa_?NumberQ, u_?NumberQ] := Module[{t}, (u+1)*AdjointSieveQ[kappa,u+1] - u*AdjointSieveQ[kappa,u] - kappa*NIntegrate[AdjointSieveQ[kappa,t],{t,u,u+2}] /; !IntegerQ[2*kappa]]; PiSieveFunc[kappa_?NumberQ, u_?NumberQ, v_?NumberQ] := Module[{t}, u*AdjointSieveP[kappa,u]/sigmaFunc[kappa,u] + kappa*NIntegrate[AdjointSieveP[kappa,t+1]/sigmaFunc[kappa,t],{t,v-1,u}]]; PiSieveFunc[kappa_, u_] := PiSieveFunc[kappa,u,u]; PiTildeSieveFunc[kappa_, u_] := PiSieveFunc[kappa,u,u-1]; XiSieveFunc[kappa_?NumberQ, u_?NumberQ, v_?NumberQ] := Module[{t}, u*AdjointSieveQ[kappa,u]/sigmaFunc[kappa,u] - kappa*NIntegrate[AdjointSieveQ[kappa,t+1]/sigmaFunc[kappa,t],{t,v-1,u}]]; XiSieveFunc[kappa_, u_] := XiSieveFunc[kappa,u,u]; XiTildeSieveFunc[kappa_, u_] := XiSieveFunc[kappa,u,u-1]; alphaSearchSieveFunc[kappa_,u_] := AdjointSieveQ[kappa,u-1]*(PiTildeSieveFunc[kappa,u]-2) + AdjointSieveP[kappa,u-1]*XiTildeSieveFunc[kappa,u]; betaSearchSieveFunc[kappa_,alpha_,u_] := IwaniecInnerQ[kappa,u,alpha,u]/u^(2*kappa-1); (* * We somewhat sloppily generate ranges for each member of avec and then use * Union to sort the results and eliminate duplicate entries. * We eliminate entries that are less than u0, but we cannot generate any which * are greater than u1. *) ChopUpInterval[u0_?NumberQ,u1_?NumberQ,avec_List] := Select[Union[{u0,u1}, Flatten[Range[Floor[Max[u0,#]]+Mod[#,1],u1]& /@ avec]], (#-u0 >= 0)&]; (* * The options for both UpperSieveFunc[...] and LowerSieveFunc[...] * are stored as Options[SieveFunctions]. *) Options[SieveFunctions] = {ComputationMethod->Automatic, AlphaBeta->Automatic}; UpperSieveFunc[kappa_?NumberQ,u_?NumberQ,opts___?OptionQ] := Module[{method=ComputationMethod/.{opts}/.Options[SieveFunctions], alphabeta = AlphaBeta/.{opts}/.Options[SieveFunctions]}, UpperSieveFuncInternal[method,alphabeta,kappa,u] ]; UpperSieveFuncInternal[method_,Automatic,kappa_?NumberQ,u_?NumberQ] := UpperSieveFuncInternal[method,AlphaBetaParams[kappa,ComputationMethod->method], kappa,u]; UpperSieveFuncInternal[method_,{alpha_?NumberQ,beta_?NumberQ},kappa_?NumberQ,u_?NumberQ] := 1/sigmaFunc[kappa,u,ComputationMethod->method] /; alpha-beta>=0 && beta-2>=0 && u>0 && u-alpha<=0; UpperSieveFuncInternal[method_,alphabeta:{alpha_?NumberQ,beta_?NumberQ},kappa_?NumberQ,u_?NumberQ] := Module[{u0=Ceiling[u], u1, u2}, u1 = Max[Select[{u0-1,u0-1+Mod[alpha,1], u0-1+Mod[beta,1]}, (#-u < 0)&]]; u2 = Min[Select[{u0-1+Mod[alpha,1], u0-1+Mod[beta,1],u0}, (#-u >= 0)&]]; UpperSieveFuncSegment[method,kappa,alphabeta,u1,u2][u] ] /; alpha-beta >= 0 && beta-2 >= 0 && u>0 && u-alpha > 0; LowerSieveFunc[kappa_?NumberQ,u_?NumberQ,opts___?OptionQ] := Module[{method=ComputationMethod/.{opts}/.Options[SieveFunctions], alphabeta = AlphaBeta/.{opts}/.Options[SieveFunctions]}, LowerSieveFuncInternal[method,alphabeta,kappa,u] ]; LowerSieveFuncInternal[method_,Automatic,kappa_?NumberQ,u_?NumberQ] := LowerSieveFuncInternal[method,AlphaBetaParams[kappa,ComputationMethod->method], kappa,u]; LowerSieveFuncInternal[method_,{alpha_?NumberQ,beta_?NumberQ},kappa_?NumberQ,u_?NumberQ] := 0 /; alpha-beta>=0 && beta-2>=0 && u>0 && u-beta<=0; LowerSieveFuncInternal[method_,alphabeta:{alpha_?NumberQ,beta_?NumberQ},kappa_?NumberQ,u_?NumberQ] := Module[{u0=Ceiling[u], u1, u2}, u1 = Max[Select[{u0-1,u0-1+Mod[alpha,1], u0-1+Mod[beta,1]}, (#-u < 0)&]]; u2 = Min[Select[{u0-1+Mod[alpha,1], u0-1+Mod[beta,1],u0}, (#-u >= 0)&]]; LowerSieveFuncSegment[method,kappa,alphabeta,u1,u2][u] ] /; alpha-beta >= 0 && beta-2 >=0 && u>0 && u-beta > 0; (* * Note that we don't bother with AccuracyGoal->Infinity when we do NDSolve[...] * below since it should be unnecessary -- because UpperSieveFunc[...] never gets * close to zero (at least for alpha=alphaParam[kappa], beta=betaParam[kappa]). *) UpperSieveFuncSegment[method_,kappa_?NumberQ,alphabeta:{alpha_?NumberQ,beta_?NumberQ}, u1_?NumberQ,u2_?NumberQ] := UpperSieveFuncSegment[method,kappa,alphabeta,u1,u2] = Module[{u,F}, F/.NDSolve[{F[u1]==UpperSieveFuncInternal[method,alphabeta,kappa,u1], D[u^kappa*F[u],u] == kappa*u^(kappa-1)*LowerSieveFuncInternal[method,alphabeta,kappa,u-1]}, F,{u,u1,u2}][[1]] ]; (* * We need to do further investigation on the proper choice for AccuracyGoal in * the code for NDSolve[...] below. AccuracyGoal->Infinity would seem to be the * best choice, since the function can be near zero when u is near beta, but it * "confuses" NDSolve when the function is identically zero (for u<=beta). *) LowerSieveFuncSegment[method_,kappa_?NumberQ,alphabeta:{alpha_?NumberQ,beta_?NumberQ}, u1_?NumberQ,u2_?NumberQ] := LowerSieveFuncSegment[method,kappa,alphabeta,u1,u2] = Module[{u,f}, f/.NDSolve[{f[u1]==LowerSieveFuncInternal[method,alphabeta,kappa,u1], D[u^kappa*f[u],u] == kappa*u^(kappa-1)*UpperSieveFuncInternal[method,alphabeta,kappa,u-1]}, f,{u,u1,u2}][[1]] ]; SieveP[kappa_?NumberQ,u_?NumberQ,opts___?OptionQ] := UpperSieveFunc[kappa,u,opts] + LowerSieveFunc[kappa,u,opts]; SieveQ[kappa_?NumberQ,u_?NumberQ,opts___?OptionQ] := UpperSieveFunc[kappa,u,opts] - LowerSieveFunc[kappa,u,opts]; IwaniecInnerP[kappa_?NumberQ,u_?NumberQ,opts___?OptionQ] := Module[{method=ComputationMethod/.{opts}/.Options[SieveFunctions], alphabeta = AlphaBeta/.{opts}/.Options[SieveFunctions]}, IwaniecInnerPInternal[method,alphabeta,kappa,u] ]; IwaniecInnerPInternal[method_,Automatic,kappa_?NumberQ,u_?NumberQ] := IwaniecInnerPInternal[method,AlphaBetaParams[kappa],kappa,u]; IwaniecInnerPInternal[method_,alphabeta:{alpha_?NumberQ,beta_?NumberQ},kappa_?NumberQ,u_?NumberQ] := Module[{t,options = Sequence[ComputationMethod->method,AlphaBeta->alphabeta], range = ChopUpInterval[u-1,u,{1,alpha,beta}]}, u*AdjointSieveP[kappa,u]*SieveP[kappa,u,options] + kappa*NIntegrate[IwaniecInnerPintegrand[kappa,t,options], Evaluate[{t, Sequence @@ range}]]] /; alpha-beta >= 0 && beta-2 >= 0 && u-1 > 0; IwaniecInnerPintegrand[kappa_?NumberQ,u_?NumberQ,opts___?OptionQ] := SieveP[kappa,u,opts]*AdjointSieveP[kappa,u+1]; IwaniecInnerQ[kappa_?NumberQ,u_?NumberQ,opts___?OptionQ] := Module[{method=ComputationMethod/.{opts}/.Options[SieveFunctions], alphabeta = AlphaBeta/.{opts}/.Options[SieveFunctions]}, IwaniecInnerQInternal[method,alphabeta,kappa,u] ]; IwaniecInnerQInternal[method_,Automatic,kappa_?NumberQ,u_?NumberQ] := IwaniecInnerQInternal[method,AlphaBetaParams[kappa],kappa,u]; IwaniecInnerQInternal[method_,alphabeta:{alpha_?NumberQ,beta_?NumberQ},kappa_?NumberQ,u_?NumberQ] := Module[{t,options = Sequence[ComputationMethod->method,AlphaBeta->alphabeta], range = ChopUpInterval[u-1,u,{1,alpha,beta}]}, u*AdjointSieveQ[kappa,u]*SieveQ[kappa,u,options] - kappa*NIntegrate[IwaniecInnerQintegrand[kappa,t,options], Evaluate[{t, Sequence @@ range}]]] /; alpha-beta >= 0 && beta-2 >= 0 && u-1 > 0; IwaniecInnerQintegrand[kappa_?NumberQ,u_?NumberQ,opts___?OptionQ] := SieveQ[kappa,u,opts]*AdjointSieveQ[kappa,u+1]; OldalphaInterpolatingFunction = Interpolation[{{1.0, 2.0}, {1.5, 3.91148053729465}, {2.0, 5.35772744559447}, {2.5, 6.83999829752238}, {3.0, 8.37193124087476}, {3.5, 9.938884278702993}, {4.0, 11.53179887333135}, {4.5, 13.144726313935339}, {5.0, 14.773560564422173}, {5.5, 16.41534952473994}, {6.0, 18.06789882325006}, {6.5, 19.729536743282214}, {7.0, 21.39895238871602}, {7.5, 23.0751034155587}, {8.0, 24.757152857145783}, {8.5, 26.44440066898821}, {9.0, 28.136279888607895}, {9.5, 29.832305732604354}, {10.0, 31.532073832271813}} ]; OldbetaInterpolatingFunction = Interpolation[{{1.0, 2.0}, {1.5, 3.11582102793118}, {2.0, 4.26645028414864}, {2.5, 5.44406767347209}, {3.0, 6.64085945080066}, {3.5, 7.85146280929457}, {4.0, 9.07224832988245}, {4.5, 10.300628297351295}, {5.0, 11.534709043155392}, {5.5, 12.773073718692846}, {6.0, 14.014643713794339}, {6.5, 15.258587779754821}, {7.0, 16.50425792691003}, {7.5, 17.751145703520976}, {8.0, 18.99885243735298}, {8.5, 20.247056225760552}, {9.0, 21.4955099366742}, {9.5, 22.744012660084884}, {10.0, 23.992407668870914}} ]; (* * Values of alpha[kappa],beta[kappa] were initially provided by Diamond (and * Wheeler) up to kappa=2, as given in Table 1 of the paper "A Boundary Value * Problem for a Pair of Differential Delay Equations Related to Sieve Theory, * II", by Diamond, Halberstam and Richert. Values for larger kappa were * computed by repeated calls to AlphaBetaExtrapolate[...]. We then printed * tables of {kappa,alphaParam[kappa]}//InputForm and of * {kappa,betaParam[kappa]}//InputForm over the range of kappa listed below and * cut-and-pasted the results into this file. *) alphaInterpolatingFunction = Interpolation[{{1.0, 2.0}, {1.5, 3.9114803664122135}, {2.0, 5.35772721287637}, {2.5, 6.83999806677035}, {3.0, 8.371931201514588}, {3.5, 9.9388843575464}, {4.0, 11.531798689326964}, {4.5, 13.14472600027376}, {5.0, 14.773560338450814}, {5.5, 16.415349546518847}, {6.0, 18.0678982882347}, {6.5, 19.72953507336269}, {7.0, 21.398950084400397}, {7.5, 23.075102907700888}, {8.0, 24.757153227573674}, {8.5, 26.444400589662646}, {9.0, 28.136279334072093}, {9.5, 29.832306940018704}, {10.0, 31.532073341373458} }]; betaInterpolatingFunction = Interpolation[{{1.0, 2.0}, {1.5, 3.1158210084394864}, {2.0, 4.266450191078886}, {2.5, 5.444067484665609}, {3.0, 6.640859342208907}, {3.5, 7.85146280680935}, {4.0, 9.072248439870867}, {4.5, 10.30062836655699}, {5.0, 11.534709164251854}, {5.5, 12.77307384200808}, {6.0, 14.014643812637978}, {6.5, 15.258588242295511}, {7.0, 16.504258451747017}, {7.5, 17.751146059428116}, {8.0, 18.998852968997024}, {8.5, 20.2470567471739}, {9.0, 21.49551049175982}, {9.5, 22.744013537264475}, {10.0, 23.99240838782658} }]; (* By convention we always set the range of alphaInterpolatingFunction * and betaInterpolatingFunction to lie 0.5 beyond the last datapoint. *) alphaInterpolatingFunction[[1,1,2]] += 0.5; betaInterpolatingFunction[[1,1,2]] += 0.5; (* * In AlphaBetaParams IwaniecInnerQ[...] gets huge for larger values of * "u", so we scale it down to keep it more nearly constant in * magnitude. *) Options[AlphaBetaParams] = {ComputationMethod->Automatic}; AlphaBetaParams[1.0] = N[AlphaBetaParams[1] = {2,2}]; AlphaBetaParams[kappa_?NumberQ, opts___?OptionQ] := AlphaBetaParamsInternal[ComputationMethod/.{opts}/.Options[AlphaBetaParams],kappa]; AlphaBetaParamsInternal[method_,kappa_] := AlphaBetaParamsInternal[method,kappa] = Module[{a,b, aX = alphaInterpolatingFunction[kappa], bX = betaInterpolatingFunction[kappa], abPair, scaleFactor}, scaleFactor = aX^(2*kappa-1); abPair = FindRoot[{IwaniecInnerP[kappa,a, ComputationMethod->method,AlphaBeta->{a,b}]==2, IwaniecInnerQ[kappa,a, ComputationMethod->method,AlphaBeta->{a,b}]/scaleFactor==0}, {a,aX},{b,bX}]; {a,b}/.abPair]; alphaParam[kappa_?NumberQ] := AlphaBetaParams[kappa][[1]]; betaParam[kappa_?NumberQ] := AlphaBetaParams[kappa][[2]]; (* * Note, we force kappa to be an integer multiple of 1/2, as opposed * to a multiple of 0.5, since AdjointSieveQ[kappa,u] is a polynomial in the * former case, and for kappa in the range we're interested in * AdjointSieveQ[...] can be much more quickly evaluated as a polynomial. *) AlphaBetaExtrapolate[] := Module[{kappa=Round[2*alphaInterpolatingFunction[[1,1,2]]]/2, a, b, alphaTable, betaTable}, {a,b} = AlphaBetaParams[kappa]; alphaTable = Append[Table[{k,alphaInterpolatingFunction[k]},{k,1,kappa-0.5,0.5}], {kappa,a}]; betaTable = Append[Table[{k,betaInterpolatingFunction[k]},{k,1,kappa-0.5,0.5}], {kappa,b}]; alphaInterpolatingFunction = Interpolation[alphaTable]; betaInterpolatingFunction = Interpolation[betaTable]; alphaInterpolatingFunction[[1,1,2]] += 0.5; betaInterpolatingFunction[[1,1,2]] += 0.5; ]; End[] EndPackage[]