Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
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 , especially the legendary bases and .

So far, this is what I've got. Here, I'll use to denote the regular superfunction for base b at the fixed point L, to denote the regular Schroder function, and 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:



where the coefficients are recursively generated by

,
.

(see Bell polynomials.)

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



to extend to other values of .

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



then is used to extend to the upper halfplane.

Regular Abel function:

.

The algorithm

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

Here, I will denote by the mth Taylor coefficient of , and the mth Fourier coefficient of .

1. Set up initial coefficients.
We set , , and all the rest to zero, giving . We also set the Fourier coefficients all to zero.

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

.

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

.

where the curves and 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 .

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 , with 32 terms on the super/etc. series and 100 terms on the theta mapping, the residual (i.e. , 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 .


Attached Files
.gp   MyKneserImpl.gp (Size: 9.53 KB / Downloads: 596)
Reply


Messages In This Thread
Attempt to make own implementation of "Kneser" algorithm: trouble - by mike3 - 06/15/2011, 07:43 AM

Possibly Related Threads...
Thread Author Replies Views Last Post
  Kneser-iteration on n-periodic-points (base say \sqrt(2)) Gottfried 11 4,316 05/05/2021, 04:53 AM
Last Post: Gottfried
  fast accurate Kneser sexp algorithm sheldonison 38 115,058 01/14/2016, 05:05 AM
Last Post: sheldonison
  Attempt to find a limit point but each step needs doubling the precision... Gottfried 15 33,765 11/09/2014, 10:25 PM
Last Post: tommy1729
  "Kneser"/Riemann mapping method code for *complex* bases mike3 2 10,091 08/15/2011, 03:14 PM
Last Post: Gottfried
  Numerical algorithm for Fourier continuum sum tetration theory mike3 12 30,619 09/18/2010, 04:12 AM
Last Post: mike3
  Attempting to compute the kslog numerically (i.e., Kneser's construction) jaydfox 11 28,917 10/26/2009, 05:56 PM
Last Post: bo198214



Users browsing this thread: 1 Guest(s)