fast accurate Kneser sexp algorithm

[sheldonison]

Later, I will probably post results showing how the Taylor and Fourier series terms converge, and how many terms it takes to get appropriate precision in the theta(z) fourier series so as to generate sufficiently accurate samples for the sexp(z) series......

Here's an example, settings which I like as the next step up from the term term double precision results, which take less than a minute. Type "init(exp(1));loop(14);" and this next example runs for 7-8 minutes on my laptop, and reports that resulting sexp(z) Taylor series is accurate to 108 binary bits, or just a bit over 32 decimal digits. That precision would be within a radius of abs(z)<=0.5.

\p 67;
ctrm = 75;
cvec = 150;
rtrm = 1000;
rvec = 2000;
sfuncprec = 1E-32;

Here, gp's precision is set to 67 digits. The goal is 32 decimal digits of precision.

The sfuncprec (superfunction precision) is set to 1E-32 digits, to give superf/isuperf results accurate to 32 digits. That is, the isuperf/superf iteration count was set to 235 iterations. If you take the , the result will be within 10^-32 of the fixed point of L. That is roughly what it takes to get 32 digits of accuracy from the superfunction. We compare exp(L+1E-32) with the fixed point approximation, L+L*1E-32, and we notice that the exponent function has an x^2 term around 1E-64. Thus, the maximum precision we're going to be left with is just a bit over 32 decimal digits.

67-32 leaves 35 digits of precision for gp to work with, or 3 spare bits for rounding errors, which surprisingly seems sufficient! We're asking for a Taylor series with 75 terms. By the Nyquist sampling criteria, we need at least 150 sample points. So cvec (the number of sample terms used to generate the sexp(z) Taylor series) =150. There will be a small sampling error term, due to higher frequencies in the superfunction being sampled. But, for base e sexp(z), the higher frequency terms roll off pretty fast, with each term being roughly half the magnitude of the term before. We want results accurate to 32 decimal digits, so we need 150 sample points, accurate to 32 decimal digits. The sample points will be in the upper half of the complex plane, centered around z=0. The points closest to the real axis are and . This is at

rtrm = 1000 is driven by the requirement that the Riemann approximation, which is be accurate to 32 decimal digits, at . superf is accurate to 32 decimal digits, so the error term is determined by theta(z), which has a radius of convergence at Im(z)=0. The 1000th term in theta(z) has a magnitude of 2.26E-05. The theta(z)=a_n*exp(z*2pi*i). With the known value for Im(z) closest to the real axis, we get magnitude of the a_n term is 6.02*10E-34, which is smaller than the goal of 10E-32. So, using rtrm~=1000 is ok. Sampling for the function requires 2000 samples by the Nyquist criterion, so we need rvec (sample points for the RiemannCircle/theta(z) function) >= 2000. The high frequency roll off is much more gradual, so we need a little bit of extra precision here.

We also need to show that our 75 term sexp(z) Taylor series supports 32 decimal digits of precision at a radius with abs(z)<=0.5. This is for the unit length we use to generate the RiemannCircle/Theta function from. The unit length will be from -0.5+idelta to +0.5+idelta. We want all terms over the unit length accurate to 32 decimal digits. The 75th term is -7.15 E-25. 7.15E-25*(0.5^75)=1.9E-47, which is much much less then 1E-32. The sexp(z) Taylor series has 32 decimal digits of precision out to a radius of about 0.78, which is greater then 1/2.

So, all of the parameters are sufficient to get around 32 bits of precision. Observationally, when iterating the "loop(14)", we get about 108 binary bits equal in the sexp(z) algorithm versus the riemaprx(z) algorithm. This is about 33 decimal digits. Graphing the error gainst a golden copy with higher precision shows similar results.
- Sheldon