Equations for Kneser sexp algorithm

[sheldonison]

This is a continuation of the thread in the computation forum, http://math.eretrandre.org/tetrationforu...hp?tid=486

This thead contains some of the mathematical equations I used for the fast Kneser algorithm. In this post, B is the base for the sexp function, and L is the fixed point for base B.

This is the complex valued superfunction, developed from the fixed point L for base B, where B>\( \eta \)
\( c=L\times\ln(B) \) for base e, c=L
\( \operatorname{superf}_{B}(z) = \lim_{n \to \infty}
B^{[n](L + c^{z-n})} =
\lim_{n \to \infty} \exp^{[n]} ((L+c^{z-n})/\ln(B)) \)

This is the complex valued inverse superfunction, developed from the fixed point, which is the inverse of the equation above. The inverse superfunction has the property, that isuperf(B^z)=isuperf(B)+1. This particular equation is normalized, so that it converges to the same value as the limit of n approaches infinity. Both of these two functions are implemented in the pari-GP program I wrote.
\( \operatorname{isuperf}_{B}(z) = \lim_{n \to \infty}
\log_{c} (c^n \times ((\log_{B}^{[n]}(z))-L)) \)

If we started with a pefect sexp(z) function, then this is the 1-cyclic theta function linking the sexp with the superf/isuperf.
\( \theta(z)=\operatorname{isuperf}(\operatorname{sexp}(z))-z \)
\( \operatorname{sexp}(z)=\operatorname{superf}(z+\theta(z)) \)

Theta(z) has a singularity at all integer values of n. Theta(z) is represented by an infinite sequence of fourier terms. The fourier series for theta(z) can be developed from any arbitrary unit length on the real axis of sexp(z), where z>-2. Only terms with positive values of n are included, and all terms a_n for negative values of n are zero.
\( \theta(z)=\sum_{n=0}^{\infty}a_n\times \e^{(2\pi n i z)} \)

Theta(z) is intimately connected to the Riemann unit circle mapping, used by Kneser's construction. The Taylor series for the Riemann unit circle function (I'm not sure of the correct notation here) uses the exactly the same a_n coeffecients as the 1-cyclic theta function! This is something that connects the complex fourier analysis of theta(z) to the theory of complex analytic functions, which is really neat! The RiemannCircle has a singularity at z=1, which corresponds to the singularities at the integer values of theta(z).
\( \operatorname{RiemannCircle}(z) = \theta(\ln(z)/2\pi i) \)
\( \theta(z) = \operatorname{RiemannCircle}(e^{2\pi i z}) =\operatorname{isuperf}(\operatorname{sexp}(z))-z \)
\( \operatorname{RiemannCircle(z)=\sum_{n=0}^{\infty}a_n\times z^n \)

If we had a perfect sexp(z) Taylor series, then we have a function for the values of the Riemann unit circle function, which is generated from the theta function, using the equation above, from the inverse superfunction. Now, we can use Cauchy's integral formula to calculate the Taylor series for the Riemann unit circle function. And this also gives us the coefficients of the 1-cyclic theta(z) function.

Of course, there is still the problem of the singularity on the unit circle, which causes problems due to slow convergence. In later posts, I will try to go into some detail, showing values for the Taylor series results for the Riemann circle function, and how the coefficients slowly decay, with poor convergence on the unit circle.

The program I wrote iterates, calculating approximate values the RiemannCircle Taylor series based on an approximation function for sexp(z). And then uses the Taylor series for the RiemannCircle function to calculate another better approximation for the sexp(z) function.

Many, many, many more details to follow! Be patient. This may take a few days....
- Enjoy, Sheldon

[sheldonison]

Continuing on, one more quick post today. We have the following equation for sexp(z).

\( \operatorname{sexp}(z)=\operatorname{superf}(z+\theta(z)) \)

If we substitude in the equation for \( \theta(z) = \operatorname{RiemannCircle}(e^{2\pi i z}) \), then we get the following equation for sexp(z), which is only defined for imag(z)>=0.

\( \operatorname{sexp}(z)=\operatorname{superf}(z+
\operatorname{RiemannCircle}(e^{2\pi i z})) \)

But supposed we want the Taylor series for sexp(z), centered at z=0? This particular equation for sexp(z) is only valid for imag(z)>=0. By the Schwarz reflection theorem, for imag(z)<0
sexp(conj(z)) = conj(sexp(z))

Now, since we have defined sexp(z) on for the entire complex plane. We can generate sexp(z) for a unit circle, centered around z=0. Then generate the Taylor series for sexp(z) using the Cauchy Integral theorem.

So, in summary, earlier I gave the equation for how to generate the Taylor series for the RiemannCircle(z) from the sexp(z) function. And now I have given a different reverse algorithm, whicch uses the Schwarz reflection property to generate a Taylor series for the sexp(z) function from the RiemannCircle(z) function.

In the fast Kneser sexp code post, the pari-GP subroutine for riemaprx(z) is very close to the algorithm I just gave for generating the sexp(z) for imag(z)>0. The loop(n) routine takes those values from the riemaprx, and uses them to generate the updated Taylor series for the sexp(z) function.
- Sheldon

[Gottfried]

Hi Sheldon -

just a short note since I've limited time slices currently. I second Henryk's "wow" :-)
When I'll be back home I'll try out what you have made avaliable in Pari/GP. By the moment I can't say anything, but when I read Kneser's article in 2008 I always wished to have some operational guidance to understand what he was doing.
So thank you much so far...

Gottfried

[sheldonison]

There are main three ideas which improve convergence, plus some other ideas as well. I will discuss the first two main ideas in this post.
(1) Most of the terms of the RiemannCircle Taylor series are not computionally significant to the sexp(z) Taylor series.
(2) Generating the RiemannCircle Taylor series at an idelta (imaginary delta) allows much much more accurate/precise results in the RiemannCircle Taylor series since the RiemannCircle singularity is at a radius of 1, and the idelta corresponds to a radius <1.
(3) the sexp(z) seems to be more precise than the RiemannCircle Taylor series it was generated from, by approximately 8 binary bits.

The sexp(z) for imag(z)>=0, generated from the RiemannCircle function converges slowly at the real axis. In theory, if the Taylor series for the RiemannCircle function has an infinite number of terms, it doesn't matter. The RiemannCircle function converges, and the resulting sexp(z) function, at least for imag(z)>=0, will also be infinitely accurate.

However, the number of terms required in the RiemannCircle function is prohibitively large. This chart extrapolates from the patterns I observe from much smaller number of terms. The largest number of terms I've calculated is 12800 Taylor series terms, where the log_2 of the 12800th RiemannCircle Taylor series term was -19.83. If the pattern holds, a Taylor series with around a billion terms would give a bit more than 35 binary bits of precision, or 10.5 decimal digits of precision. Calculating a one billion term Taylor series is not a practical! Henryk, I haven't succeeded in figuring out how to display images uploaded as attachments, apologies! I think I uploaded an image, attachment=728, but couldn't add it to the post. So I'm still using my personal website for images!


The two ideas I gave above uncover a way out of this conundrum. First off, we use the algorithm given for generating the sexp Taylor series from the Riemann mapping Taylor series. The sexp Taylor series is much more well behaved, and for double precision accuracy (52 binary bits), all we need is an sexp(z) Taylor series with 50 terms in it. This requires 100 sampling points, which can be even placed on the upper half of the unit circle centered at z=0, which is the unit circle I used to calculated the sexp(z) Taylor series. None of these 50 terms need be on the real axis.

Consider this set of 50 points:
e^(pi*I*1/200), e^(pi*I*3/200).... e^(pi*I*199/200)

The imaginary value for the first point is 0.0157*I. I call this value idelta. This corresponds to a radius of 0.905 in the RiemannCircle Taylor series. This means the 400th term in the RiemannCircle Taylor series is being multiplied by 7E-18, which means it is not computationally relevant! Likewise for all terms >=400! So we calculated a billion terms in the RiemannCircle Taylor series to get 13 decimal digits of accuracy at the real axis, and only 300 of those terms are relevant when generating the sexp(z) Taylor series, with 50 terms.

The other idea, is that generating the Taylor series for the RiemannCircle at a radius corresponding to idelta improves the convergence of the lower frequency terms too. Basically, the singularity is at radius=1. The unit circle of the RiemannCircle function we're generating a Fourier series for is poorly behaved at a radius of 1. But at a radius of 0.905, corresponding to idelta, it is much better behaved. Here is a plot of the imaginary, and real contours of the theta(z) function (corresponds to the RiemannCircle) at a radius of 1, showing the singularity, at integer values of z.



By moving to idelta=~0.016, the singularity is largely gone. And the Fourier series for the theta(z) (RiemannCircle Taylor series) converges very nicely with only 800 sampling points. The first two images below are blow ups of the imaginary and real contours, at idelta. The third image shows the log_2 of the RiemanCircle terms calculated at i=0, and then calculated at i=0.016i. Notice the loss of precision for terms>300 when the Taylor series at idelta terms are less than 2^-57 (I will discuss significance later). By the way, this scheme works quite nicely, since we only calculate the RiemannCircle for values >= idelta. Again, I will comment more later, in my next post. This is taking even longer than I thought!
- enjoy, Sheldon

[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

[sheldonison]

... 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.

Well, actually, I have some theories, which were going through my head before I got caught up in implementing and debugging the pari-gp code. The high frequency terms in the RiemannCircle function only effect points closest to the real axis due to exponential decay as imag(z) increases.
\( \theta(z)=\sum_{n=0}^{\infty}a_n\times \e^{(2\pi n i z)} \)

In fact, exponential decay, means that all terms a_n with n>=1 should improve in this loop:
Riemann(z) -> sexp(z) -> Riemann(z).

The real error term in a_0 gets cancelled out by the recentering. The imaginary error term in a_0 gets cancelled out by the Schwarz refelection.

The real error term in a_1 is not a factor for terms near the imag(z)=1 part of the circle (e^-2Pi is ~= 0.002), so the error term gets factored down by about 1/10th, if you do the average over 100 points. The imaginary error term probably cancels in the Schwarz reflection. So that's a factor of 1/20th. The error term improvement for higher frequency coefficients gets exponentially better, with each term's error contribution impoving twice as much as the term before it.

So there is some hope that a mathematical argument for congruence that can be made, with some effort taken to describe how the coefficients of the RiemannCircle function converge from one iteration to the next. Also, I could graph the convergence of the various terms in the Taylor series for the RiemannCircle function, and verify that most of the error is in the lowest frequency term, unless the computation limits cause chaotic results (which they sometimes do as noted in my earlier post; I should recreate and analyze the chaotic case for base 10 sexp(z) extended precision with sexp(z) series at 150 terms. I also found a setting that causes chaos for base e.)

update, it occurs to me that this is a general principle. Lets say g(z) is a function with real values at the real axis, and f(z) is an approximation, over a unit length, such that f(0)=g(0) and f(1)=g(1). Then there is some 1-cyclic theta fourier series function defining f(z) over the entire real axis.
\( f(z) = g(z+\theta(z) \)
Now, theta(z) is probably not even an analytic function. But we could wrap the real valued theta(z) around the unit circle. It would have a laurent series. We could throw out all of the terms in the laurent series with a_n*z^-n. Now theta_2(z) is still 1-cyclic, but it is complex valued at the real axis.
\( f_2(z) = g(z+\theta_2(z) \)

You might take this new f_2 function, and use it to generate another function, f_3. Generate the Taylor series over a half circle of f_2(z), using the complex conjugate for the other half of the circle where imag(z)<0. At this point, you should see the connection between the sequence of functions, f(z), f_2(z), f_3(z), and the iterative algorithm for the sexp(z) function I have described, where the sequence of f, f_1, f_2 ... converge to the desired g(z) function.
- Sheldon

[sheldonison]

update, it occurs to me that this is a general principle. Lets say g(z) is a function with real values at the real axis, and f(z) is an approximation, over a unit length, such that f(0)=g(0) and f(1)=g(1). Then there is some 1-cyclic theta fourier series function defining f(z) over the entire real axis.
\( f(z) = g(z+\theta(z) \)
Now, theta(z) is probably not even an analytic function. But we could wrap the real valued theta(z) around the unit circle. It would have a laurent series. We could throw out all of the terms in the laurent series with a_n*z^-n. Now theta_2(z) is still 1-cyclic, but it is complex valued at the real axis.
\( f_2(z) = g(z+\theta_2(z) \)

You might take this new f_2 function, and use it to generate another function, f_3. Generate the Taylor series over a half circle of f_2(z), using the complex conjugate for the other half of the circle where imag(z)<0. At this point, you should see the connection between the sequence of functions, f(z), f_2(z), f_3(z), and the iterative algorithm for the sexp(z) function I have described, where the sequence of f, f_1, f_2 ... converge to the desired g(z) function.

Another quick post, with some ideas I had when I was coming up with this algorithm to generate sexp(z). There may be a general case Riemann mapping algorithm where there is a real valued (at the real axis) auxiliary function, f(z+theta(z)). In our case, the sequence of f(z) functions is desired to converge to sexp(z) which is real valued at the real axis for z>-2.

One important note. It only works if the desired f(z) is analytic, and real valued at the real axis. f(z) is the initial approximation, and we're trying to come up with theta(z), which is exactly related to a Riemann mapping, but theta(z) may have a singularity, making convergence even more difficult.
- Sheldon

[tommy1729]

and those singularities are exactly the main problem.

[sheldonison]

... the sexp(z) is about 100-200x more accurate than the RiemannCircle function from which it is generated.

correction The RiemannCircle(z) function is 100-200x more accurate than the sexp(z) from which it is generated -- however this data does include the renormalize update in the code, which improved convergence. In this case, all of the improvement appears to come from the Riemann mapping. But, as you'll see farther down, there is another alternative without the renormalize step where more of the improvement comes from the sexp(z) step!

Here is a table, showing how binary bits precision improves with each iteration. The initial seed has an accuracy of about 4 bits binary. The accuracy is for the lowest term in the Taylor series. The next term is about the same, for both Taylor series, and accuracy for higher order terms is more accurate (at least for awhile). I could plot it. Anyway, the point is the Riemann series is 7 or 8 bits more accurate than the sexp(z) series from which it is generated. The sexp(z) series has roughly the same accuracy as the Riemann function from which it is generated.
Riemann -8.41108211297959 bits
sexp(z) -9.11708328076014 bits
Riemann -15.1121338561072 bits
sexp(z) -15.1316170140021 bits
Riemann -22.4179232976446 bits
sexp(z) -22.4897567784208 bits
Riemann -30.0155349055434 bits
sexp(z) -30.1087094975908 bits
Riemann -37.7126066320302 bits
sexp(z) -37.8144609513923 bits
Riemann -45.4478768923269 bits
sexp(z) -45.5533186457503 bits
Riemann -53.2001314165620 bits
sexp(z) -53.3071463323265 bits
Riemann -60.9582476200045 bits
sexp(z) -61.0659125192666 bits
Riemann -68.7204813782899 bits
sexp(z) -68.8284373037337 bits
Riemann -76.4826669744581 bits
sexp(z) -76.5907565309966 bits
Riemann -84.2469508431293 bits
sexp(z) -84.3550685136293 bits
Riemann -92.0030319102508 bits
sexp(z) -92.1112561892320 bits

And here is a second table, without the renormalization step. The renormalize step is after the sexp(z) series is generated. It multiplies the sexp(z) series by whatever constant is required to force sexp(0.5)=exp(sexp(-0.5)). After 13 iterations, accuracy is at 69 bits, (as opposed to 92 bits). Here's what's interesting to me. This is a more elegant loop, without the somewhat ugly renormalize step. This leads to a more accurate sexp(z) series, but leads to less improvement overall, in the Riemann mapping, with a little singularity since sexp(0.5)<>exp(sexp(-0.5).

Overall, this approach took 17 or 18 iterations until convergence stopped at 92 bits, as opposed to 13 iterations. So, anyway, that's the data. Both routines are working to improve the results together.

Riemann -8.41108211297959
sexp(z) -9.11708328076014
Riemann -13.5881879738804
sexp(z) -17.1061334888054
Riemann -18.8839530119948
sexp(z) -21.4385286084149
Riemann -24.1419288537565
sexp(z) -26.8704129838725
Riemann -29.4092009637620
sexp(z) -32.1372761148881
Riemann -34.6784789286960
sexp(z) -37.4290603754555
Riemann -39.9512192183627
sexp(z) -42.7181894098719
Riemann -45.2270233243075
sexp(z) -48.0041816783103
Riemann -50.5061031443367
sexp(z) -53.3029348163809
Riemann -55.7895507636155
sexp(z) -58.6076934808030
Riemann -61.0771911987248
sexp(z) -63.9108585125742
Riemann -66.3693959721491
sexp(z) -69.2206441976381
Riemann -71.6653883035104
sexp(z) -74.5285659987165
Riemann -76.9672004092254
sexp(z) -79.8616209389427
Riemann -82.2744841556731
sexp(z) -85.1839830576327
Riemann -87.5851241316643
sexp(z) -90.5153257252998
Riemann -92.8193397830429
(convergence stops improving here)

Some other comments. The sexp(z) isn't even a superfunction. Its an approximation, with its accuracy as a superfunction decaying as the radius increases. Its most important property is that sexp(z) is real valued at the real axis.

Contrast that with the RiemannCircle function which is guaranteed by construction to be a true superfunction, with f(z+1)=exp(f(z)) for all z with imag(z)>=idelta. At idelta, it is the best superfunction matching sexp(z). But we only kept terms from the Taylor series, not the Laurent series, so it is not equal to sexp(z). It has guarenteed convergence for all imag(z)>=idelta.
- Sheldon

[bo198214]

I am still a bit confused about the (even basic) steps in your program and about the key relations that make it work. So let me present from my view:

In the Kneser construction there is the regular Abel function at the primary fixpoint (in the upper halfplane). You call this function \( \text{isuperf} \). It maps the upper halfplane H to some area \( \Omega \) which is bounded by the arcs \( \text{isuperf}((0,1))+k \), \( k\in\mathbb{Z} \).

The Riemann mapping in Kneser's construction then is the bijective map that maps \( \Omega \) back to the upper halfplane. So that Kneser then defines:
\( \text{slog}(z) = \rho(\text{isuperf}(z)) \), or in super expressions:
\( \text{sexp}(z) = \text{superf}(\rho^{-1}(z)) \).

What has this function to do with your \( \text{RiemannCircle} \)???
Apart from
\( \rho^{-1}(z)=\text{isuperf}(\text{sexp}(z))=\theta(z)+z \).
Why is it called RiemannCircle?

Ok, but even if only the naming confuses me here, i am still a bit puzzled what makes your algorithm work. I understand your basic steps as follows:

1. from sexp obtain the Fourier-coefficients (with positive indices) of \( \theta(z)= \text{isuperf}(\text{sexp}(z))-z \) and call the Fourier-Series \( \text{RiemannCircle}(e^{2\pi i z})=\theta(z) \).

2. from RiemannCircle obtain the Powerseries-coefficients of \( \text{Riemaprx}(z)=\text{superf}(z+\text{RiemannCircle}(e^{2\pi i z})) \).

And there my questions start. It seems that Riemaprx is not equal to sexp as should:
\( \text{superf}(z+\text{RiemannCircle}(e^{2\pi i z})) = \text{superf}(z+\theta(z))=\text{sexp}(z) \). I guess this is due to truncation of the negative indices in the Fourier-Series. Can you confirm this?

So if it is not the same as sexp, you force it into a real-valued sexp (which I guess Riemaprx isnt) by this conjugation trick. (Which really causes me headache seen from a theoretic side, I guess the so constructed function is not even continuous on the unit circle) Would it be equally possible to just define sexp just as the real part of Riemaprx, or doesnt the whole algorithm converge then?