fast accurate Kneser sexp algorithm

[sheldonison]

This is a small update, which improves performance, and allows more flexibility. Its a little faster (25 seconds for 105-110 binary bit precision, as opposed to 1 minute). I added an slog function, as well as a gentaylorseries function that will generate the sexp taylor series about any point in the complex plane. I cleaned up the code, and made it a little simpler, in that the "n" parameter for loop(n) is now optional, and the default base for init; is base e. You can still "loop(1,2...)" at a time if desired.

  kneser.gp (Size: 24.09 KB / Downloads: 72)
For the most recent code version: go to the Nov 21st, 2011 thread.

Code:

\r kneser.gp
init;loop
   base          2.71828182845904523536029
   fixed point   0.318131505204764135312654 +
.... 25 seconds later, after 13 loop iterations
gp > sexp(0.5)
1.6463542337511945809719240315921
gp > slog(I)
-1.3937809364969141717386716853218 + 1.0171347771806901474061430995874*I
gp > gentaylor(3*I,1)
complex sexp Taylor series centered at 3*I
gp > for (s=1,4,print(tseries[s])); /* prints first 4 series terms */
0.37090658903228120608281774428121 + 1.3368216707889122043183467452060*I
0.018300482687991735872774054492565 + 0.069611076949751462108351480615434*I
-0.042221079601603961951042313227903 + 0.024296334049072988793474267457353*I
-0.015851643810851135752181652414535 - 0.014789535958785875826049470725178*I

There is a much larger performance increase for "\p 134", which now takes about 5 1/2 minutes, as opposed to the previous version which took 32 minutes. "\p 134 gives sexp results accurate to 64 decimal digits. The new algorithm cleverly converges the sexp(z) with the riemaprx(z) at imag(z)=0.12*i, thereby avoiding all of the Riemann mapping singularity problems requiring larger and larger numbers of Riemann approximation Fourier series terms. That's where the speed up comes from. I had to fix a couple of pari-gp precision loss problems to make this new algorithm reliable for all bases, init(1.47) - init(100000). The earlier algorithm is still available via "slowmode=1". One reason for the new kneser.gp code is to allow a secondary program, to generate pentation, which I will be posting separately.

Another code example, showing \p 134, and the taylor series function.

Code:

\r kneser.gp
\p 134
init(exp(1));loop
.... wait 5 1/2 minutes, after 26 iterations ....
sexp(0.5)
gp > sexp(0.5)
1.646354233751194580971924031592114518205311648969041530394879048

Code:

type "morestuff" on the command line to see new functions
new functions:
slog(z,est)      inverse sexp(z), using optional est as initial estimate
gentaylor(w,r)   generates tseries sexp taylor series, centered at w, radius r
taprx(z)         evaulates tseries taylor series from gentaylor at z

This implementation shows that excellent convergence can be had at 0.12*I. Here the Cauchy unit circle which is sampled to generate the sexp(z) taylor series is divided into a region where imag(z) at the unit circle is >= 0.12*I. For that region, the Cauchy unit circle is sampled from the Riemann approximation generated from the 1-cyclic transformation of the super function. But, for points where imag(z)<0.12*I, the sample is taken from the exp or logarithm of previous version of sexp(z) itself. Then these samples are stitched together, to generate the next sexp(z) series. Then the next Riemann approximation Fourier series is generated at 0.12*I. The new algorithm required identifying and fixing the pari-gp precision problems, which results in much more predictable stable behavior, which convinces me that the convergence will continue with more and more iterations.
- Sheldon