Attempt to make own implementation of "Kneser" algorithm: trouble
#1
Hi.

I've been trying to build my own program in Pari/GP to run the Kneser algorithm for tetration, but I can't seem to get it to work so well. The goal is to ultimately be able to generalize it to merge complex bases and bases less than \( \eta \), especially the legendary bases \( b = -1 \) and \( b = 0.04 \).

So far, this is what I've got. Here, I'll use \( RS_{b,L}(z) \) to denote the regular superfunction for base b at the fixed point L, \( R\Psi_{b,L}(z) \) to denote the regular Schroder function, and \( RA_{b,L}(z) \) to denote the regular Abel function. The lower subscripts, then, denote the index of the successive approximations to the given object which is subscripted.

Computation of super-, Schroder, and Abel function

The program uses the following formulas to compute the regular super-, Schroder, and Abel functions of the exponential.

Regular superfunction:

\( RS_{b,L}(z) = \sum_{n=0}^{\infty} r_n \left(L \log(b)\right)^{nz} \)

where the coefficients are recursively generated by

\( r_0 = L \), \( r_1 = 1 \)
\( r_n = \frac{B_n(1! \log(b) r_1, 2! \log(b) r_2, \cdots, (n-1)! \log(b) r_{n-1}, 0)}{n! \left(L^{n-1} \log(b) - \log(b)\right)} \).

(see Bell polynomials.)

Accuracy is best for \( z \) with large negative real part. Then we can use the recurrence

\( RS_{b,L}(z+1) = b^{RS_{b,L}(z)} \)

to extend to other values of \( z \).

Regular Schroder function:

This is generated by taking the inverse of the above when the coefficients are used as a power series, via the "reversion of series" formulae. Pari/GP already implements this so there is no need to do so explicitly in the program. The identity

\( R\Psi_{b,L}(b^z) = L \log(b) R\Psi_{b,L}(z) \)

then is used to extend to the upper halfplane.

Regular Abel function:

\( RA_{b,L}(z) = \log_{L \log(b)}\left(R\Psi_{b,L}(z)\right) \).

The algorithm

From what I could tell, the Kneser algorithm works as follows. Let \( F_n(z) \) be an approximation Taylor series to the tetrational, and \( \theta_n(z) \) be an approximation to the Kneser mapping. sheldonison, did I get this one right?

Here, I will denote by \( {f_m}_n \) the mth Taylor coefficient of \( F_n \), and \( {t_m}_n \) the mth Fourier coefficient of \( \theta_n(z) \).

1. Set up initial coefficients.
We set \( {f_0}_0 = 1 \), \( {f_1}_0 = 1 \), and all the rest to zero, giving \( F_0(z) = 1 + z \). We also set the Fourier coefficients all to zero.

2. Compute Fourier integral.
Now, we compute the Fourier integral to update the theta coefficients.

\( {t_m}_{n+1} = \int_{-1}^{0} \left(RA_{b,L}\left(F_n(x)\right) - x\right) e^{-2\pi i m x} dx \).

3. Compute Cauchy integral.
We now form the Taylor coefficients for the next tetrational approximation via the Cauchy integral.

\( {f_m}_{n+1} = \frac{1}{2\pi i} \left(\int_{\gamma_1} \frac{RS_{b,L}\left(z + \theta_{n+1}(z)\right)}{z^{m+1}} dz + \int_{\gamma_2} \frac{\bar{RS_{b,L}\left(\bar{z} + \theta_{n+1}(\bar{z})\right)}}{z^{m+1}} dz\right) \).

where the curves \( \gamma_1 \) and \( \gamma_2 \) are the upper and lower halves, respectively, of the unit circle, moving counterclockwise, starting at the positive real axis.

4. Normalize Taylor series.
Divide the Taylor series by its leading coefficient, so it acquires the value 1 at \( z = 0 \).

5. Repeat.
Repeat steps 2-4 until maximum precision is reached.

Problems: What's wrong?

The problem is that I can't seem to get this to go. For base \( e \), with 32 terms on the super/etc. series and 100 terms on the theta mapping, the residual (i.e. \( \left|F_n(0.5) - b^{F_n(-0.5)}\right| \), and is used to estimate the precision of the approximation) never drops below 0.00005. After 12 iterations, the residual flattens out at around 0.00005086. Using 64 terms on the super/etc. series and 800 terms on the theta mapping still only yields a residual of 0.00001258. I want a lot more precision than that! Smile I suspect the problem has to do with the fickle nature of the theta mapping, but I'm not sure how to overcome this.

Code

Here is the program. I hope that this is clear enough to be readable.

Code:
/* Compute the sum of a Fourier series:
* f(z) = sum_{n=0...inf} a_n e^(unz)
*/
FourierSum(coef, u, z) = sum(n=0,matsize(coef)[2]-1, coef[n+1]*exp(u*n*z));

/* Compute the sum of a Taylor series centered at u:
* f(z) = sum_{n=0...inf} a_n (z - u)^n
*/
TaylorSum(coef, u, z) = sum(n=0,matsize(coef)[2]-1, coef[n+1]*(z - u)^n);

/* Compute Bell polynomial */
bellcomplete(coefs) = {
            local(n=matsize(coefs)[2]);
            local(M=matrix(n,n));
            for(i=1,n,\
                for(j=1,n,\
                    M[i, j] = 0;
                    if(j-i+1 > 0, M[i, j] = binomial(n-i, j-i)*coefs[j-i+1]);
                    if(j-i+1 == 0, M[i, j] = -1);
                );
            );
            return(matdet(M));
}

/* Compute the coefficients of the regular superfunction of the exponential
* with and at fixed point L and a given base (L must be a fixed point of the base).
*/
RegularCoeffs(L, base, N) = {
            local(a = vector(N));
            local(bellcoeffs);
            local(logbase = log(base));

            a[0+1] = L;
            a[1+1] = 1;

            /* Recurrence formula for Fourier coefficients:
             *   a_n = B_n(1! log(b) a_1, ..., (n-1)! log(b) a_(n-1), 0)/(n! (L^(n-1) log(b)^n - log(b)))
             */
            for(n=2,N-1,
                bellcoeffs = vector(n, k, if(k==n, 0, logbase*k!*a[k+1]));
                a[n+1] = bellcomplete(bellcoeffs)/(n! * (L^(n-1) * logbase^n - logbase));
               );

            return(a);
}

/* Compute the coefficients of the regular Schroder function of the exponential
* with and at fixed point L.
*/
MakeSchroderCoeffs(RegCoeffs) = {
            local(invpoly = 0);
            local(RegCoeffsTrunc = RegCoeffs);
            local(InvCoeffs = RegCoeffs);

            RegCoeffsTrunc[1] = 0; /* truncate constant term */

            /* Do series reversion. This gives the inverse of (InvSchroder(x) - L) */
            invpoly = serreverse(TaylorSum(RegCoeffsTrunc, 0, x) + O(x^matsize(RegCoeffsTrunc)[2]));

            /* Extract coefficients */
            for(i=0,matsize(RegCoeffsTrunc)[2]-1,
                InvCoeffs[i+1] = polcoeff(invpoly, i);
               );

            /* Done! */
            return(InvCoeffs);
}

/* Compute the regular superfunction at a given fixed point L of a base b. */
RegularSuperfunction(L, base, SuperCoeffs, z) = {
            if(real(z) < -5,
               return(FourierSum(SuperCoeffs, log(L*log(base)), z));,
               return(base^RegularSuperfunction(L, base, SuperCoeffs, z-1));
              );
}

/* Compute the regular Schroder function at a given fixed point L of a base b. */
RegularSchroderFunction(L, base, SchroderCoeffs, z) = {
            if(abs(z - L) > 0.5,
               return(L*log(base)*RegularSchroderFunction(L, base, SchroderCoeffs, log(z)/log(base)));,
               return(TaylorSum(SchroderCoeffs, L, z));
              );
}

/* Compute the regular Abel function at a given fixd point L. */
RegularAbelFunction(L, base, SchroderCoeffs, z) = \
            log(RegularSchroderFunction(L, base, SchroderCoeffs, z))/log(L*log(base));

/* --- Quadrature routines --- */

/* Do quadrature using function samples in samp with step delta. */
/* This is an inferior, Simpson's rule quadrature, implemented for testing only. In the future, we will
* drop quadrature and upgrade to FFT for high performance. Note that samp must have an odd number of
* elements.
*/
QuadratureCore(samp, delta) = {
            local(ssum=0);

            ssum = samp[1] + sum(n=1,matsize(samp)[2]-2,(4 + ((-1)^(n-1) - 1))*samp[n+1]) + samp[matsize(samp)[2]];

            return((delta/3)*ssum);
}

/* Compute some regularly-spaced nodes over an interval to use for quadrature. N must be odd. */
QuadratureNodes(a, b, N) = {
            local(nodes);

            nodes = vector(N, i, a + (i-1)*((b-a)/(N-1)));

            return(nodes);
}

/* Compute the delta value for a given interval and node count */
QuadratureDelta(a, b, N) = (b - a)/(N-1);

/* --- End quadrature routines --- */

/* Set up the tetration system. */
Initialize() = {
         /* Initialize the fixed points and quadrature. */
         /* Base and fixed points */
         base = exp(1);
         upperfix = 0.318 + 1.337*I;
         for(i=1,16,upperfix = upperfix - (base^upperfix - upperfix)/(log(base)*base^upperfix - 1));

         /* Quadrature setup */
         NumQuadNodes = 1001;
         ImagOffset = 0.0000001*I; /* offset off the real axis to do Fourier integration */
         ThetaQuadratureNodesUpper = QuadratureNodes(-1.0 + ImagOffset, 0.0 + ImagOffset, NumQuadNodes);
         ThetaQuadratureDeltaUpper = QuadratureDelta(-1.0 + ImagOffset, 0.0 + ImagOffset, NumQuadNodes);
         CauchyCircleQuadratureNodesUpper = QuadratureNodes(0, Pi, NumQuadNodes);
         CauchyCircleQuadratureDeltaUpper = QuadratureDelta(0, Pi, NumQuadNodes);
         CauchyCircleQuadratureNodesLower = QuadratureNodes(Pi, 2*Pi, NumQuadNodes);
         CauchyCircleQuadratureDeltaLower = QuadratureDelta(Pi, 2*Pi, NumQuadNodes);

         /* Term counts */
         NumRegTerms = 32;
         NumTaylorTerms = 32;
         NumThetaTerms = 100;

         /* Generate series */
         RegCoeffsUpper = RegularCoeffs(upperfix, base, NumRegTerms);
         RegSchrCoeffsUpper = MakeSchroderCoeffs(RegCoeffsUpper);

         /* Set up initial Taylor series */
         TaylorCoeffs = vector(NumTaylorTerms);
         TaylorCoeffs[0+1] = 1; /* coefficient of x^0 = 1 */
         TaylorCoeffs[1+1] = 1; /* coefficient of x^1 = 1 */

         /* Set up initial Riemann mapping series */
         RiemannCoeffsUpper = vector(NumThetaTerms);
}

/* Compute upper regular superfunction */
RegularSuperUpper(z) = RegularSuperfunction(upperfix, base, RegCoeffsUpper, z);

/* Compute upper regular Abel function */
RegularAbelUpper(z) = RegularAbelFunction(upperfix, base, RegSchrCoeffsUpper, z);

/* Construct Fourier coefficients for a Riemann theta(z) mapping from a tetration Taylor series. */
RiemannFromTaylorStep() = {
         local(Samples);
         local(FourierSamp);

         /* Sample the function at the node points. */
         Samples = vector(NumQuadNodes, i, \
                          RegularAbelUpper(TaylorSum(TaylorCoeffs, 0, ThetaQuadratureNodesUpper[i])) -\
                              ThetaQuadratureNodesUpper[i]);
         /*
         plot(X=1,NumQuadNodes,real(Samples[floor(X)]));
         plot(X=1,NumQuadNodes,imag(Samples[floor(X)]));
         */

         /* Do the Fourier integrals */
         /* We multiply by exp(-2 pi i n ImagOffset) to compensate for the offset along the imaginary axis. */
         for(n=0,NumThetaTerms-1,
             FourierSamp = vector(NumQuadNodes, i, Samples[i]*exp(-2*Pi*I*n*real(ThetaQuadratureNodesUpper[i])));
             RiemannCoeffsUpper[n+1] = QuadratureCore(FourierSamp, ThetaQuadratureDeltaUpper);
             RiemannCoeffsUpper[n+1] = RiemannCoeffsUpper[n+1] * exp(-2*Pi*I*n*ImagOffset);
            );
}

/* Calculate the theta mapping. */
ThetaMappingUpper(z) = FourierSum(RiemannCoeffsUpper, 2*Pi*I, z);

/* Calculate the upper "warped" superfunction. */
WarpSuperUpper(z) = RegularSuperUpper(z + ThetaMappingUpper(z));

/* Calculate the lower "warped" superfunction. */
WarpSuperLower(z) = conj(WarpSuperUpper(conj(z))); /* fine for real bases > eta only */

/* Do the Cauchy integral to get the Taylor series. */
CauchyTaylorStep() = {
         local(Samples);
         local(UpperSuperSamp);
         local(LowerSuperSamp);
         local(UpperInt);
         local(LowerInt);

         /* Presampling */
         UpperSuperSamp = vector(NumQuadNodes, i, WarpSuperUpper(exp(I*CauchyCircleQuadratureNodesUpper[i])));
         LowerSuperSamp = vector(NumQuadNodes, i, WarpSuperLower(exp(I*CauchyCircleQuadratureNodesLower[i])));

         /* Do Cauchy integral */
         for(n=0,NumTaylorTerms-1,\
             /* Sample the function around the upper half-circle. */\
             Samples = vector(NumQuadNodes, i, UpperSuperSamp[i]/\
                                               (exp(I*CauchyCircleQuadratureNodesUpper[i])^(n+1)) *\
                                               I*exp(I*CauchyCircleQuadratureNodesUpper[i]));\
             UpperInt = QuadratureCore(Samples, CauchyCircleQuadratureDeltaUpper);\

             /* Do the same around the lower half-circle. */\
             Samples = vector(NumQuadNodes, i, LowerSuperSamp[i]/\
                                               (exp(I*CauchyCircleQuadratureNodesLower[i])^(n+1)) *\
                                               I*exp(I*CauchyCircleQuadratureNodesLower[i]));\
             LowerInt = QuadratureCore(Samples, CauchyCircleQuadratureDeltaLower);\

             /* Get the result */\
             TaylorCoeffs[n+1] = 1/(2*Pi*I) * (UpperInt + LowerInt);\
            );

         /* Renormalize Taylor series */
         TaylorCoeffs = TaylorCoeffs / TaylorSum(TaylorCoeffs, 0, 0);
}

/* Do full iterations */
DoKneserIteration(numiters) = {
         for(i=1,numiters,
             print("Iteration ", i, ":");
             print("Building Kneser mapping...");
             RiemannFromTaylorStep();
             print("Performing Cauchy integral...");
             CauchyTaylorStep();
             print("Done. R = ", abs(base^TaylorSum(TaylorCoeffs, 0, -0.5) - TaylorSum(TaylorCoeffs, 0, 0.5)));
            );
}

/* Calculate full superfunction */
FullTetApprox(z) = {
         if(real(z) > 0.5, return(base^FullTetApprox(z-1)));
         if(real(z) < -0.5, return(log(FullTetApprox(z+1))/log(base)));

         if(imag(z) > 1, return(WarpSuperUpper(z)));
         if(imag(z) < -1, return(WarpSuperLower(z)));

         return(TaylorSum(TaylorCoeffs, 0, z));
}

To use:
Initialize() to initialize the system
then
DoKneserIteration(12) (or however many loops you want)
then
FullTetApprox(x) to compute \( ^x b \).


Attached Files
.gp   MyKneserImpl.gp (Size: 9.53 KB / Downloads: 791)
#2
how do I change base (from exp(1))??
#3
(06/15/2011, 07:43 AM)mike3 Wrote: Hi.

I've been trying to build my own program in Pari/GP to run the Kneser algorithm for tetration, but I can't seem to get it to work so well. The goal is to ultimately be able to generalize it to merge complex bases and bases less than \( \eta \), especially the legendary bases \( b = -1 \) and \( b = 0.04 \).

Hey Mike,

Take a look at my recent post, which has some useful details, and a nice picture.

The first thing I notice is that your "step 2" is trying to compute the coefficients for theta, via the Fourier integral at the real axis, from -1 to zero. Theta(z) has a really nasty singularity at integer coefficient. Analytic Fourier analysis works really great with very minimal number of terms if you're generating the Fourier analysis of an analytic function, but the closer to the singularity you are, the more of a problem you have. As noted in the post above, I got around this by computing the integral at imag(z)=0.12i. I actually computed theta(z) from -0.5+0.12i to +0.5+0.12i, because that was centered very close to the middle of my sexp(z) approximation, where it was most accurate. 0.12i is far enough away from the singularity that the theta(z) function is very well behaved, and requires minimal terms in the Fourier analysis. fyi, my sexp(z) function is sampled in a unit circle around z=0. But the price you pay is that theta(z) now only converges for imag(z)>=0.12i. I went through many many iterations about how to deal with that, originally using a discreet Cauchy sampling, before I got to the current implementation, that splices together the Kneser/Riemann approximation with log(z)/exp(z), from the previous iterations sexp(z) approximation, which works really nicely. So that would be a change to your step (3) for your Cauchy integral.

If you do all of that, you should get about 4-5 bits precision improvement per iteration loop, with far fewer sampling points and terms in the theta(z) mapping. You can double that to about 8-9 bits precision per loop iteration, by adding some stuff to your step (4), normalize the taylor series, doing something I called renormalization to smooth out the function so that sexp(0.5+0.12i)=B^sexp(-0.5+0.12i), which makes the theta(z) Fourier integral behave better, since it is sampling a continuous function.

As far as the algorithm goes, it was originally designed to handle the Kneser mapping with its inherent singularity for \( \text{superf}^{-1}(0,1,B...) \). For some pairs of complex fixed point merger scenarios, where one of the fixed point is attracting, it may not make sense do the merger near the slog singularity, in which case none of the advice in this post is necessary, since the two Fourier series should have excellent convergence on their own, and the sexp intermediary is unnecessary. Eventually, I may support the complex base case where both fixed points are repelling in my Kneser.gp code.

Hope that helps.
- Sheldon
#4
(06/15/2011, 12:59 PM)nuninho1980 Wrote: how do I change base (from exp(1))??

See where it says "base" in the Initialize() part? This needs to be changed. It should work for real bases greater than \( \eta = e^{1/e} \). Not sure how big you can make the base due to the fixed initial fixpoint approximation, though.
#5
(06/15/2011, 03:18 PM)sheldonison Wrote: The first thing I notice is that your "step 2" is trying to compute the coefficients for theta, via the Fourier integral at the real axis, from -1 to zero. Theta(z) has a really nasty singularity at integer coefficient. Analytic Fourier analysis works really great with very minimal number of terms if you're generating the Fourier analysis of an analytic function, but the closer to the singularity you are, the more of a problem you have. As noted in the post above, I got around this by computing the integral at imag(z)=0.12i. I actually computed theta(z) from -0.5+0.12i to +0.5+0.12i, because that was centered very close to the middle of my sexp(z) approximation, where it was most accurate. 0.12i is far enough away from the singularity that the theta(z) function is very well behaved, and requires minimal terms in the Fourier analysis. fyi, my sexp(z) function is sampled in a unit circle around z=0. But the price you pay is that theta(z) now only converges for imag(z)>=0.12i. I went through many many iterations about how to deal with that, originally using a discreet Cauchy sampling, before I got to the current implementation, that splices together the Kneser/Riemann approximation with log(z)/exp(z), from the previous iterations sexp(z) approximation, which works really nicely. So that would be a change to your step (3) for your Cauchy integral.

So you would use a countour like that of Kouznetsov's Cauchy method, right (i.e. half-circle connected to another with straight lines to form a pill- or oval-like shape), and when \( |\Im(z)| < 0.12 \), you switch to integrating over the Taylor appproximation from the last run? I'm not sure what you mean by "log/exp" of it, though -- if you're integrating with zero at the center, the Taylor already converges best there, so I don't see why you need to take log/exp.

And I guess that why it doesn't converge below \( |\Im(z)| = 0.12 \) is that the function is not perfectly periodic initially, so when the area integrated over is repeated, it forms something not analytic along the given parallel. But I suppose this goes away as the algorithm progresses, and so eventually the theta mapping converges right up to near the real axis.
#6
(06/15/2011, 10:10 PM)mike3 Wrote:
(06/15/2011, 12:59 PM)nuninho1980 Wrote: how do I change base (from exp(1))??

See where it says "base" in the Initialize() part? This needs to be changed. It should work for real bases greater than \( \eta = e^{1/e} \). Not sure how big you can make the base due to the fixed initial fixpoint approximation, though.
ok. I edited base from exp(1) to 10 in your code. I entered 'Initialize()', 'DoKneserIteration(12)', 'FullTetApprox(0.5)' and can evaluate result. Smile but... if base<eta then it doesn't work, ok?
but after I tested your code, Kneser.gp by Sheldon is much better precision than your code.
#7
(06/15/2011, 10:22 PM)mike3 Wrote:
(06/15/2011, 03:18 PM)sheldonison Wrote: The first thing I notice is that your "step 2" is trying to compute the coefficients for theta, via the Fourier integral at the real axis, from -1 to zero. Theta(z) has a really nasty singularity at integer coefficient. Analytic Fourier analysis works really great with very minimal number of terms if you're generating the Fourier analysis of an analytic function, but the closer to the singularity you are, the more of a problem you have. As noted in the post above, I got around this by computing the integral at imag(z)=0.12i. I actually computed theta(z) from -0.5+0.12i to +0.5+0.12i, because that was centered very close to the middle of my sexp(z) approximation, where it was most accurate. 0.12i is far enough away from the singularity that the theta(z) function is very well behaved, and requires minimal terms in the Fourier analysis. fyi, my sexp(z) function is sampled in a unit circle around z=0. But the price you pay is that theta(z) now only converges for imag(z)>=0.12i. I went through many many iterations about how to deal with that, originally using a discreet Cauchy sampling, before I got to the current implementation, that splices together the Kneser/Riemann approximation with log(z)/exp(z), from the previous iterations sexp(z) approximation, which works really nicely. So that would be a change to your step (3) for your Cauchy integral.

So you would use a countour like that of Kouznetsov's Cauchy method, right (i.e. half-circle connected to another with straight lines to form a pill- or oval-like shape), and when \( |\Im(z)| < 0.12 \), you switch to integrating over the Taylor appproximation from the last run? I'm not sure what you mean by "log/exp" of it, though -- if you're integrating with zero at the center, the Taylor already converges best there, so I don't see why you need to take log/exp.

And I guess that why it doesn't converge below \( |\Im(z)| = 0.12 \) is that the function is not perfectly periodic initially, so when the area integrated over is repeated, it forms something not analytic along the given parallel. But I suppose this goes away as the algorithm progresses, and so eventually the theta mapping converges right up to near the real axis.
Hey Mike,
Take a look at the picture in this post, to see how I build the radius=1 sampling circle for the next sexp(z) function. If you had like a trillion terms in the theta(z) Fourier series, then you still wouldn't get nearly as much precision at the real axis as I can get in the sexp(z) series after 30 seconds (13 iterations), with 109 terms in the sexp(z) taylor series, 32 decimal digits accuracy. That sexp(z) was generated from its precursor (as described in the picture), and the \( \text{superf}(z+\theta(z)) \), where the theta(z) series had 93 terms in it. The sexp(z) series is more than accurate enough to generate that one trillion term theta(z) series at the real axis. By definition, the infinite Fourier series will only converge where is sampled or at imag(z)>imag(sample). Precision for the finite series falls off logarithmically, from the imag(sample), to the real axis. The 93 term theta(z) is accurate to 32 decimal digits for imag(z)>=0.12i, but is only accurate to ~3 decimal digits at the real axis.

The number of Fourier series terms required in theta(z) goes up approximately as 1/imag(z), so if we needed 93 terms as 0.12i, then we would need 93x12=~1100 series terms at 0.01i, to get equivalent accuracy. Although, we don't need an infinite number of terms at the real axis, since the terms do get smaller too.
- Sheldon
#8
(06/15/2011, 11:21 PM)sheldonison Wrote:
(06/15/2011, 10:22 PM)mike3 Wrote:
(06/15/2011, 03:18 PM)sheldonison Wrote: The first thing I notice is that your "step 2" is trying to compute the coefficients for theta, via the Fourier integral at the real axis, from -1 to zero. Theta(z) has a really nasty singularity at integer coefficient. Analytic Fourier analysis works really great with very minimal number of terms if you're generating the Fourier analysis of an analytic function, but the closer to the singularity you are, the more of a problem you have. As noted in the post above, I got around this by computing the integral at imag(z)=0.12i. I actually computed theta(z) from -0.5+0.12i to +0.5+0.12i, because that was centered very close to the middle of my sexp(z) approximation, where it was most accurate. 0.12i is far enough away from the singularity that the theta(z) function is very well behaved, and requires minimal terms in the Fourier analysis. fyi, my sexp(z) function is sampled in a unit circle around z=0. But the price you pay is that theta(z) now only converges for imag(z)>=0.12i. I went through many many iterations about how to deal with that, originally using a discreet Cauchy sampling, before I got to the current implementation, that splices together the Kneser/Riemann approximation with log(z)/exp(z), from the previous iterations sexp(z) approximation, which works really nicely. So that would be a change to your step (3) for your Cauchy integral.

So you would use a countour like that of Kouznetsov's Cauchy method, right (i.e. half-circle connected to another with straight lines to form a pill- or oval-like shape), and when \( |\Im(z)| < 0.12 \), you switch to integrating over the Taylor appproximation from the last run? I'm not sure what you mean by "log/exp" of it, though -- if you're integrating with zero at the center, the Taylor already converges best there, so I don't see why you need to take log/exp.

And I guess that why it doesn't converge below \( |\Im(z)| = 0.12 \) is that the function is not perfectly periodic initially, so when the area integrated over is repeated, it forms something not analytic along the given parallel. But I suppose this goes away as the algorithm progresses, and so eventually the theta mapping converges right up to near the real axis.
Hey Mike,
Take a look at the picture in this post, to see how I build the radius=1 sampling circle for the next sexp(z) function. If you had like a trillion terms in the theta(z) Fourier series, then you still wouldn't get nearly as much precision at the real axis as I can get in the sexp(z) series after 45 seconds (13 iterations), with 109 terms in the sexp(z) taylor series, 32 decimal digits accuracy. That sexp(z) was generated from its precursor (as described in the picture), and the \( \text{superf}(z+\theta(z)) \), where the theta(z) series had 93 terms in it. The sexp(z) series is sufficiently accurate to generate a one trillion term theta(z) series at the real axis. By definition, the infinite Fourier series will only converge where is sampled or at imag(z)>imag(sample).
- Sheldon

Ahh. Now I see -- so using the log/exp maximizes precision.

#9
(06/15/2011, 10:49 PM)nuninho1980 Wrote:
(06/15/2011, 10:10 PM)mike3 Wrote:
(06/15/2011, 12:59 PM)nuninho1980 Wrote: how do I change base (from exp(1))??

See where it says "base" in the Initialize() part? This needs to be changed. It should work for real bases greater than \( \eta = e^{1/e} \). Not sure how big you can make the base due to the fixed initial fixpoint approximation, though.
ok. I edited base from exp(1) to 10 in your code. I entered 'Initialize()', 'DoKneserIteration(12)', 'FullTetApprox(0.5)' and can evaluate result. Smile but... if base<eta then it doesn't work, ok?
but after I tested your code, Kneser.gp by Sheldon is much better precision than your code.

Yes, that's right. But if I can improve it with the just-suggested algorithms, I should be able to get it to go to full precision...
#10
Well, I managed to break through that convergence barrier now, by using the Taylor approximation at the real axis.

Here's the new code:

Code:
/* Compute the sum of a Fourier series:
* f(z) = sum_{n=0...inf} a_n e^(uinz)
*/
FourierSum(coef, u, z) = sum(n=0,matsize(coef)[2]-1, coef[n+1]*exp(u*n*z));

/* Compute the sum of a Taylor series centered at u:
* f(z) = sum_{n=0...inf} a_n (z - u)^n
*/
TaylorSum(coef, u, z) = sum(n=0,matsize(coef)[2]-1, coef[n+1]*(z - u)^n);

/* Compute Bell polynomial */
bellcomplete(coefs) = {
            local(n=matsize(coefs)[2]);
            local(M=matrix(n,n));
            for(i=1,n,\
                for(j=1,n,\
                    M[i, j] = 0;
                    if(j-i+1 > 0, M[i, j] = binomial(n-i, j-i)*coefs[j-i+1]);
                    if(j-i+1 == 0, M[i, j] = -1);
                );
            );
            return(matdet(M));
}

/* Compute the coefficients of the regular superfunction of the exponential
* with and at fixed point L and a given base (L must be a fixed point of the base).
*/
RegularCoeffs(L, base, N) = {
            local(a = vector(N));
            local(bellcoeffs);
            local(logbase = log(base));

            a[0+1] = L;
            a[1+1] = 1;

            /* Recurrence formula for Fourier coefficients:
             *   a_n = B_n(1! log(b) a_1, ..., (n-1)! log(b) a_(n-1), 0)/(n! (L^(n-1) log(b)^n - log(b)))
             */
            for(n=2,N-1,
                bellcoeffs = vector(n, k, if(k==n, 0, logbase*k!*a[k+1]));
                a[n+1] = bellcomplete(bellcoeffs)/(n! * (L^(n-1) * logbase^n - logbase));
               );

            return(a);
}

/* Compute the coefficients of the regular Schroder function of the exponential
* with and at fixed point L.
*/
MakeSchroderCoeffs(RegCoeffs) = {
            local(invpoly = 0);
            local(RegCoeffsTrunc = RegCoeffs);
            local(InvCoeffs = RegCoeffs);

            RegCoeffsTrunc[1] = 0; /* truncate constant term */

            /* Do series reversion. This gives the inverse of (InvSchroder(x) - L) */
            invpoly = serreverse(TaylorSum(RegCoeffsTrunc, 0, x) + O(x^matsize(RegCoeffsTrunc)[2]));

            /* Extract coefficients */
            for(i=0,matsize(RegCoeffsTrunc)[2]-1,
                InvCoeffs[i+1] = polcoeff(invpoly, i);
               );

            /* Done! */
            return(InvCoeffs);
}

/* Compute the regular superfunction at a given fixed point L of a base b. */
RegularSuperfunction(L, base, SuperCoeffs, z) = {
            if(real(z) < -5,
               return(FourierSum(SuperCoeffs, log(L*log(base)), z));,
               return(base^RegularSuperfunction(L, base, SuperCoeffs, z-1));
              );
}

/* Compute the regular Schroder function at a given fixed point L of a base b. */
RegularSchroderFunction(L, base, SchroderCoeffs, z) = {
            if(abs(z - L) > 0.5,
               return(L*log(base)*RegularSchroderFunction(L, base, SchroderCoeffs, log(z)/log(base)));,
               return(TaylorSum(SchroderCoeffs, L, z));
              );
}

/* Compute the regular Abel function at a given fixd point L. */
RegularAbelFunction(L, base, SchroderCoeffs, z) = \
            log(RegularSchroderFunction(L, base, SchroderCoeffs, z))/log(L*log(base));

/* --- Quadrature routines --- */

/* Do quadrature using function samples in samp with step delta. */
/* This is an inferior, Simpson's rule quadrature, implemented for testing only. In the future, we will
* drop quadrature and upgrade to FFT for high performance. Note that samp must have an odd number of
* elements.
*/
QuadratureCore(samp, delta) = {
            local(ssum=0);

            ssum = samp[1] + sum(n=1,matsize(samp)[2]-2,(4 + ((-1)^(n-1) - 1))*samp[n+1]) + samp[matsize(samp)[2]];

            return((delta/3)*ssum);
}

/* Compute some regularly-spaced nodes over an interval to use for quadrature. N must be odd. */
QuadratureNodes(a, b, N) = {
            local(nodes);

            nodes = vector(N, i, a + (i-1)*((b-a)/(N-1)));

            return(nodes);
}

/* Compute the delta value for a given interval and node count */
QuadratureDelta(a, b, N) = (b - a)/(N-1);

/* --- End quadrature routines --- */

/* Set up the tetration system. */
Initialize() = {
         /* Initialize the fixed points and quadrature. */
         /* Base and fixed points */
         base = exp(1);
         upperfix = 0.318 + 1.337*I;
         for(i=1,16,upperfix = upperfix - (base^upperfix - upperfix)/(log(base)*base^upperfix - 1));

         /* Quadrature setup */
         NumQuadNodes = 1001;
         ImagOffset = 0.12*I; /* offset off the real axis to do Fourier integration */
         ThetaQuadratureNodesUpper = QuadratureNodes(-1.0 + ImagOffset, 0.0 + ImagOffset, NumQuadNodes);
         ThetaQuadratureDeltaUpper = QuadratureDelta(-1.0 + ImagOffset, 0.0 + ImagOffset, NumQuadNodes);
         CauchyCircleQuadratureNodes = QuadratureNodes(0, 2*Pi, NumQuadNodes);
         CauchyCircleQuadratureDelta = QuadratureDelta(0, 2*Pi, NumQuadNodes);

         /* Term counts */
         NumRegTerms = 32;
         NumTaylorTerms = 32;
         NumThetaTerms = 100;

         /* Generate series */
         RegCoeffsUpper = RegularCoeffs(upperfix, base, NumRegTerms);
         RegSchrCoeffsUpper = MakeSchroderCoeffs(RegCoeffsUpper);

         /* Set up initial Taylor series */
         TaylorCoeffs = vector(NumTaylorTerms);
         TaylorCoeffs[0+1] = 1; /* coefficient of x^0 = 1 */
         TaylorCoeffs[1+1] = 1; /* coefficient of x^1 = 1 */

         /* Set up initial Riemann mapping series */
         RiemannCoeffsUpper = vector(NumThetaTerms);
}

/* Compute upper regular superfunction */
RegularSuperUpper(z) = RegularSuperfunction(upperfix, base, RegCoeffsUpper, z);

/* Compute upper regular Abel function */
RegularAbelUpper(z) = RegularAbelFunction(upperfix, base, RegSchrCoeffsUpper, z);

/* Construct Fourier coefficients for a Riemann theta(z) mapping from a tetration Taylor series. */
RiemannFromTaylorStep() = {
         local(Samples);
         local(FourierSamp);

         /* Sample the function at the node points. */
         Samples = vector(NumQuadNodes, i, \
                          RegularAbelUpper(TaylorSum(TaylorCoeffs, 0, ThetaQuadratureNodesUpper[i])) -\
                              ThetaQuadratureNodesUpper[i]);

         /* Do the Fourier integrals */
         /* We multiply by exp(-2 pi i n ImagOffset) to compensate for the offset along the imaginary axis. */
         for(n=0,NumThetaTerms-1,
             FourierSamp = vector(NumQuadNodes, i, Samples[i]*exp(-2*Pi*I*n*real(ThetaQuadratureNodesUpper[i])));
             RiemannCoeffsUpper[n+1] = QuadratureCore(FourierSamp, ThetaQuadratureDeltaUpper);
             RiemannCoeffsUpper[n+1] = RiemannCoeffsUpper[n+1] * exp(-2*Pi*I*n*ImagOffset);
            );
}

/* Calculate the theta mapping. */
ThetaMappingUpper(z) = FourierSum(RiemannCoeffsUpper, 2*Pi*I, z);

/* Calculate the upper "warped" superfunction. */
WarpSuperUpper(z) = RegularSuperUpper(z + ThetaMappingUpper(z));

/* Calculate the lower "warped" superfunction. */
WarpSuperLower(z) = conj(WarpSuperUpper(conj(z))); /* fine for real bases > eta only */

/* Do the Cauchy integral to get the Taylor series. */
CauchyTaylorStep() = {
         local(Samples);
         local(CauchyizedSamples);

         /* The integration proceeds in several phases:
          *   1. Integrate above Im(z) = ImagOffset using the upper regular warped superfunction,
          *   2. Integrate from Im(z) = ImagOffset to Im(z) = -ImagOffset using the old Taylor series,
          *   3. Integrate below Im(z) = -ImagOffset using the lower regular warped superfunction,
          *   4. Integrate from Im(z) = -ImagOffset to Im(z) = ImagOffset using the old Taylor series again.
          * To maximize precision, we use the exp/log of the Taylor series evaluated near zero. The whole
          * process is done counterclockwise.
          */

         /* Presampling */
         Samples = vector(NumQuadNodes);
         CirclePositions = vector(NumQuadNodes, i, exp(I*CauchyCircleQuadratureNodes[i]));
         for(n=0,NumQuadNodes-1,\
             if(abs(imag(CirclePositions[n+1])) < imag(ImagOffset),\
                if(real(CirclePositions[n+1]) > 0,\
                   /* right halfplane */\
                   Samples[n+1] = base^TaylorSum(TaylorCoeffs, 0, CirclePositions[n+1]-1);,\
                   /* left halfplane */\
                   Samples[n+1] = log(TaylorSum(TaylorCoeffs, 0, CirclePositions[n+1]+1))/log(base);\
                  );,\
                if(imag(CirclePositions[n+1]) > 0,\
                   /* upper halfplane */\
                   /* use upper regular warped superfunction */\
                   Samples[n+1] = WarpSuperUpper(CirclePositions[n+1]);,\
                   /* lower halfplane */\
                   /* use lower regular warped superfunction */\
                   Samples[n+1] = WarpSuperLower(CirclePositions[n+1]);\
                  );\
               );\
            );

         /* Do Cauchy integral */
         for(n=0,NumTaylorTerms-1,\
             /* Cauchy setup */\
             CauchyizedSamples = vector(NumQuadNodes, i, Samples[i]/(CirclePositions[i]^(n+1)) * \
                                        I*CirclePositions[i]);\

             /* Integrate */\
             TaylorCoeffs[n+1] = 1/(2*Pi*I) * QuadratureCore(CauchyizedSamples, CauchyCircleQuadratureDelta);\
            );

         /* Renormalize Taylor series */
         TaylorCoeffs = TaylorCoeffs / TaylorSum(TaylorCoeffs, 0, 0);
}

/* Do full iterations */
DoKneserIteration(numiters) = {
         for(i=1,numiters,
             print("Iteration ", i, ":");
             print("Building Kneser mapping...");
             RiemannFromTaylorStep();
             print("Performing Cauchy integral...");
             CauchyTaylorStep();
             print("Done. R = ", abs(base^TaylorSum(TaylorCoeffs, 0, -0.5) - TaylorSum(TaylorCoeffs, 0, 0.5)));
            );
}

/* Calculate full superfunction */
FullTetApprox(z) = {
         if(real(z) > 0.5, return(base^FullTetApprox(z-1)));
         if(real(z) < -0.5, return(log(FullTetApprox(z+1))/log(base)));

         if(imag(z) > 1, return(WarpSuperUpper(z)));
         if(imag(z) < -1, return(WarpSuperLower(z)));

         return(TaylorSum(TaylorCoeffs, 0, z));
}

Hope to get this going for complex bases real soon Smile


Possibly Related Threads…
Thread Author Replies Views Last Post
  Writing Kneser's super logarithm using values of Kneser at a single point JmsNxn 1 846 04/21/2023, 04:26 AM
Last Post: JmsNxn
Question Computing Kneser's Super Logarithm and its Analytic Continuation Catullus 2 1,692 07/10/2022, 04:04 AM
Last Post: Catullus
  fast accurate Kneser sexp algorithm sheldonison 40 145,748 07/03/2022, 06:38 AM
Last Post: JmsNxn
  Kneser-iteration on n-periodic-points (base say \sqrt(2)) Gottfried 11 11,866 05/05/2021, 04:53 AM
Last Post: Gottfried
  Attempt to find a limit point but each step needs doubling the precision... Gottfried 15 42,501 11/09/2014, 10:25 PM
Last Post: tommy1729
  "Kneser"/Riemann mapping method code for *complex* bases mike3 2 12,672 08/15/2011, 03:14 PM
Last Post: Gottfried
  Numerical algorithm for Fourier continuum sum tetration theory mike3 12 38,330 09/18/2010, 04:12 AM
Last Post: mike3
  Attempting to compute the kslog numerically (i.e., Kneser's construction) jaydfox 11 36,416 10/26/2009, 05:56 PM
Last Post: bo198214



Users browsing this thread: 1 Guest(s)