Equations for Kneser sexp algorithm

[sheldonison]

On to the third idea, which is that the sexp(z) is about 100-200x more accurate than the RiemannCircle function from which it is generated.

I don't really know why this is so. The pattern seems to hold when going on to higher precision. Each iteration seems to add 7 or 8 bits of precision to the accuracy of the sexp(z) approximation. By the way, my methodology for testing the precision, is to compare the average precision over a unit length of z=-1 to z=0, comparing the number of significant bits in sexp(z+idelta) and riemaprx(z+idelta), that are the same. idelta for a 50 term approximation is 0.016*i. The comparison has to be done at idelta, since at the real axis, the RiemannCircle approximation for sexp(z) is not very accurate. The number of terms in the RiemannCircle Taylor series are chosen to be large enough to give accurate enough values at idelta. When you graph this comparison, the comparison gets ~150x closer going from one iteration to the next, but the comparison is dominated by low frequency terms.

The biggest anomaly I've seen so far, was when I tried to calculate the Taylor series sexp(z) for base 10, accurate to ~90 bit or precision. The results worked fine for base e and base 2, but the base 10 equation results became chaotic around 64 bits of precision, short of the precision limits I was expecting. So I put in some code to check for decreasing Taylor series terms for the RiemannCircle, and to implement Jay's "uniqueness" checks for the sexp(z) having positive odd derivatives. This prevents precision from getting worst on subsequent iterations, but once chaotic results are seen, with higher frequency terms in the comparison, then additional iterations no longer improve precision.

I was able to cure the base 10 sexp(z) chaos, by recentering the sexp_10(z) around z=-0.5 instead of around z=0. Then the results converged to ~90 bits of precision.

So one possible limit to convergence is chaotic results.

Other minor notes that I made while implementing the program include:

• You need sufficient accuracy in the superfunction, inverse superfunction, and you need sufficient precision in gp.
  
  For the double precision results (52 binary bits), I needed 112 iterations on the isuper/superf functions. Iterating the ln(0.5) 112 times gives a value within 10^-15 of L. At the slightly smaller value of 10^-12, the superf/isuperf limits overall precision to around 42 binary bit. The precision value for gp was set to 38 decimal digits, of which nearly half is lost immediately due to the superfunction. You need as accurate a value of "L" as the precision allows. Otherwise, this can limit convergence if it is not accurate enough.
  
  
• Renormalizing the sexp(z) so that the exp(sexp(-0.5))=sexp(0.5).
  
  This improves convergence, since you want a unit length over which the sexp(z) function goes from z to exp(z). But the sexp(z) taylor series is often a little bit off. The sequence still converges without this improvement, but it probably requires twice as many iterations, or more. I generate the Taylor series for sexp(z) over a unit circle, because the Cauchy integral math is simplest that way. I could try Kouznetsov's suggestion of using a vertical strip to +i*infinity, where the vertical strips are centered at -1 and 1, f(z-1)=ln(f(z)), and f(z+1)=exp(f(z)). I also recentered the sexp(z) at each iteration, but this is mostly cosmetic.
  
  
• Defining sexp(z) over the real axis, by iterating ln(z) and exp(z) around the critical section also helps, since the sexp(z) converges best when abs(z)<=0.5.
  

The algorithm slows down quite a bit when attempting to make the next leap, to a sexp(z) Taylor series with 150 terms. Instead of taking less than a minute to loop seven times, it takes about 15 minutes to loop thirteen times, with results consistent and accurate to 92 binary bits.

• Going from pari-gp 38 bits of precision to 86 bits of precision
  
• idelta, from 0.016 to 0.005, means more Taylor series terms in RiemannCircle. This is worth some discussion by itself. idelta is set based on the number of terms in the sexp(z) Taylor series. We are putting the Riemann circle boundary at i=0.016*i instead of i=0. If idelta is set larger, than the Cauchy integral rapidly loses precision at the sexp(z) point closest to the real axis.
  
• sexp(z) Taylor series from 50 terms to 150 terms (overkill)
  
• RiemannCircle Taylor serie from 400 terms to 1800 terms
  
• superfunction iterations, from 112 to 199 iterations.
  

I don't know what it would take to prove convergence, but if it converges, than it has to converge to the Kneser sexp(z) solution, since one of the two algorithms is based on a Riemann mapping of the superfunction
- Sheldon