Revisting my accelerated slog solution using Abel matrix inversion

[jaydfox]

Note that for Kneser's slog, all of the "c_n" terms in \( \theta_s \) are all zero, which is not the case when approximating \( \theta_s \) from Jay's slog.  So, we need to calculate these undesirable c_n error terms for Jay's slog, and then force all of the c_n terms to zero by using \( \theta_{Rj} \).  Since Jay's slog is an approximation, we'll probably only force the first few terms to zero.  For example, for a 120 term compuation of JaySlog, we might force the first 10 terms c_1...c_10 to zero which should give accuracy>50 decimal digits at the real axis; more in a later post.

We start with the simplest approximation where \( \theta_{Rj} \) has one unknown term, where we approximate \( a_1\mid\;c_1=0 \)
\( \theta_{Rj}(z) = a_1 \exp(2\pi i z) + \overline{a_1} \exp(-2\pi i z) - 2\Re(a_1) \)
Unfortunately, we actually need two equations in two unknown, which are complex conjugates of each other for the upper/lower halves of the complex plane are as follows:
\( c_1 \) in terms of \( a_1;\;\overline{a_1};\;\;\; \)  upper half of complex plane; \( \theta_s \)
\( \overline{c_1} \) in terms of \( \overline{a_1};\;a_1;\;\;\;\; \) lower half of complex plane; \( \overline{\theta_s} \)

Equivalently, one could have two equations in two unknowns, one equation for the real part of a_1, and a second equation for the imaginary part of a_1, and then we use those two equations to force the real and imaginary parts of c_1 to zero.

I've been trying to figure out on my own how to eliminate those coefficients for terms with negative exponents, and my progress has been frustratingly slow. I was re-reading your post more carefully, and I see you already have found a way. I will definitely take the time to understand this, because it seems to be key to solving the Riemann mapping numerically. If you already have a sexp function with unlimited precision and accuracy, the DFT (Cauchy Integral) should return 0 for all the negative exponents. If you use low-accuracy approximation of sexp, even one with very high precision internal consistency (using the reversion of my jslog as just one example; there are others), then you get non-zero terms with negative exponents when you calculate the DFT. In theory, a proper 1-cyclic mapping (what you call a theta function) would allow you to compensate for the negative exponents and effectively cancel them out. I'm still rebuilding my code library, so I'm not at a point where I can test this, but I am very curious to see your next results.

If you're curious what I'm working on, I have a few mini projects going. I'm recoding my LU decomposition library, so that I can solve arbitrary system sizes up to several thousand terms. I'll be able to save the LU decomposition to disk, so I can do forward and back substitution of any system size up to the size of the full system. For example, I could do another 4096x4096 system (which would take days), then I could decide that I want the 3072-term solution. I could solve that in about 10 minutes, give or take.

I'm also working on code to generate the coefficients of the Schroeder and inverse Schroeder functions, in terms of ln(b) and a given fixed point c, for arbitrary base b and fixed point c. I'm not sure it works for attracting or parabolic fixed points, but it should work for repelling fixed points. I'm trying to make the code faster, as the calculation time seems to be increasing at about the 8th power of n, for n terms. For example, it takes about 16 seconds to calculate 20 terms, but over 6 minutes to calculate 30 terms. Once the terms are calculated, the Schroeder and inverse Schroeder functions for an arbitrary base b and fixed point c can be calculated in mere seconds. The solutions will be exact, to machine precision, unlike the DFT method I was using previously, which required performing hundreds (or thousands) of iterated logarithms of hundreds (or thousands) of sample points, with considerable excess precision during the calculations in order to get a desired precision.

I'm also playing around with my FFT library, trying to figure out how to cancel out negative terms. I had made some progress in understanding the problem, but hadn't yet made progress in canceling the terms. I'm going to go through your post one more time, and see if I can figure it out.