"Kneser"/Riemann mapping method code for *complex* bases
#1
Hi.

It's been a while but I figured I could post my attempt to tetrate complex bases with the "Kneser"/Riemann mapping algorithm.

This is the code I've got, which will tetrate a complex base, here base \( 3 + i \):

Code:
/* Independent implementation of Kneser mapping algorithm. */

/* 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));
              );
}

rsf(L, base, SchroderCoeffs, z) = {
            if(abs(imag(z)) < 0.02, return(3));
            return(RegularSchroderFunction(L, base, SchroderCoeffs, 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 = 3 + I;
         upperfix = 0.325 + 1.001*I;
         lowerfix = 0.098 - 1.398*I;
         for(i=1,16,upperfix = upperfix - (base^upperfix - upperfix)/(log(base)*base^upperfix - 1));
         for(i=1,16,lowerfix = lowerfix - (base^lowerfix - lowerfix)/(log(base)*base^lowerfix - 1));

         /* Quadrature setup */
         NumQuadNodes = 1001;
         ImagOffset = 0.12*I; /* offset off the real axis to do Fourier integration */
         ThetaQuadratureNodesUpper = QuadratureNodes(-0.5 + ImagOffset, 0.5 + ImagOffset, NumQuadNodes);
         ThetaQuadratureDeltaUpper = QuadratureDelta(-0.5 + ImagOffset, 0.5 + ImagOffset, NumQuadNodes);
         ThetaQuadratureNodesLower = QuadratureNodes(-0.5 - ImagOffset, 0.5 - ImagOffset, NumQuadNodes);
         ThetaQuadratureDeltaLower = QuadratureDelta(-0.5 - ImagOffset, 0.5 - 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);
         RegCoeffsLower = RegularCoeffs(lowerfix, base, NumRegTerms);
         RegSchrCoeffsLower = MakeSchroderCoeffs(RegCoeffsLower);

         /* 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);
         RiemannCoeffsLower = 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);

/* Compute lower regular superfunction */
RegularSuperLower(z) = RegularSuperfunction(lowerfix, base, RegCoeffsLower, z);

/* Compute lower regular Abel function */
RegularAbelLower(z) = RegularAbelFunction(lowerfix, base, RegSchrCoeffsLower, z);

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

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

         /* 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, SamplesUpper[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);

             FourierSamp = vector(NumQuadNodes, i, SamplesLower[i]*exp(2*Pi*I*n*real(ThetaQuadratureNodesLower[i])));
             RiemannCoeffsLower[n+1] = QuadratureCore(FourierSamp, ThetaQuadratureDeltaLower);
             RiemannCoeffsLower[n+1] = RiemannCoeffsLower[n+1] * exp(-2*Pi*I*n*ImagOffset);
            );
}

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

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

/* Calculate the lower "warped" superfunction. */
WarpSuperLower(z) = RegularSuperLower(z + ThetaMappingLower(z));

/* 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));
}

To use, run "Initialize()", then "DoKneserIteration(<some number of iterations>)". The program will generate the theta mappings to fuse the two regular iterations at the two principal fixed points. With the given Taylor/etc. term counts, 16 iterations yields a residual on the order of \( 10^{-15} \), representing what may be some of the most accurate calculations of a tetrational of a complex base to date. We get

\( ^{1/2} (3 + i) \approx 1.73774339741062 + 0.19778352516442i \)

Graphs of \( ^x (3 + i) \) are shown below.

   

On the complex plane (scale -10 to +10 on both axes):

   

As we can see, the merger worked successfully. The two different regular iterations have flowed together nice and holomorphically for \( \Re(z) > -2 \).

We can even do the "pentation" (the \( \uparrow \uparrow \uparrow \) operator) of this base now. The sequence of integer pentates is

Code:
(3 + i) ^^^ 0 ~ 1
(3 + i) ^^^ 1 ~ 3.0000000000000 + 1.0000000000000*I
(3 + i) ^^^ 2 ~ 0.067505751821009 + 0.37038084586491*I
(3 + i) ^^^ 3 ~ 0.95290527513344 + 0.41965918718415*I
(3 + i) ^^^ 4 ~ 1.6689141108302 + 1.5511207891726*I
(3 + i) ^^^ 5 ~ 0.14427450000007 + 1.2422893053044*I
(3 + i) ^^^ 6 ~ 0.59855234346250 + 0.91546250434381*I
(3 + i) ^^^ 7 ~ 0.86698754296965 + 1.1149406202426*I
(3 + i) ^^^ 8 ~ 0.67898538651189 + 1.3074168077768*I
(3 + i) ^^^ 9 ~ 0.60455661809291 + 1.1563634545945*I
(3 + i) ^^^ 10 ~ 0.69330548523182 + 1.1273955819904*I
(3 + i) ^^^ 11 ~ 0.70596729212537 + 1.1855496560953*I
(3 + i) ^^^ 12 ~ 0.66761843466257 + 1.1866711686455*I
(3 + i) ^^^ 13 ~ 0.67126297874107 + 1.1635087706344*I
(3 + i) ^^^ 14 ~ 0.68499772461207 + 1.1678785414537*I
(3 + i) ^^^ 15 ~ 0.68092916998454 + 1.1759457551948*I
(3 + i) ^^^ 16 ~ 0.67640594691649 + 1.1725959600056*I
(3 + i) ^^^ 17 ~ 0.67891225654672 + 1.1701754208531*I
(3 + i) ^^^ 18 ~ 0.68014536910711 + 1.1719549038205*I
(3 + i) ^^^ 19 ~ 0.67892704487410 + 1.1725305088990*I
(3 + i) ^^^ 20 ~ 0.67869891860093 + 1.1717250339640*I
(3 + i) ^^^ 21 ~ 0.67921563866767 + 1.1716670353204*I
(3 + i) ^^^ 22 ~ 0.67919871696253 + 1.1719897863482*I
(3 + i) ^^^ 23 ~ 0.67900254832419 + 1.1719465366940*I
(3 + i) ^^^ 24 ~ 0.67904902631454 + 1.1718306916743*I
(3 + i) ^^^ 25 ~ 0.67911531745648 + 1.1718709649452*I
(3 + i) ^^^ 26 ~ 0.67908388197865 + 1.1719075100622*I
(3 + i) ^^^ 27 ~ 0.67906467678018 + 1.1718845193510*I
(3 + i) ^^^ 28 ~ 0.67908072504425 + 1.1718750848271*I
(3 + i) ^^^ 29 ~ 0.67908487628565 + 1.1718858831902*I
(3 + i) ^^^ 30 ~ 0.67907783385992 + 1.1718873294663*I

Note how the pentation looks to converge to a limit, like tetration for STR bases.

It does not work for all complex bases, though. It seems to work only for those relatively near the real axis. Increasing the ImagOffset (how far off the real axis we sample the Fourier integrals) appears to extend the range, however. When it fails, it looks like the iteration used to analytically continue the Schroder function gets sucked into the wrong fixed point. Any way to resolve this? The aforementioned increasing of the ImagOffset causes losses in accuracy and efficiency, and is not a cure-all that allows for tetrating all complex bases outside the STR with the possible exception of the possibly singular base 0.
#2
(08/15/2011, 03:44 AM)mike3 Wrote: Hi.

It's been a while but I figured I could post my attempt to tetrate complex bases with the "Kneser"/Riemann mapping algorithm.

This is the code I've got, which will tetrate a complex base, here base \( 3 + i \):
....
It does not work for all complex bases, though. It seems to work only for those relatively near the real axis. Increasing the ImagOffset (how far off the real axis we sample the Fourier integrals) appears to extend the range, however. When it fails, it looks like the iteration used to analytically continue the Schroder function gets sucked into the wrong fixed point. Any way to resolve this? The aforementioned increasing of the ImagOffset causes losses in accuracy and efficiency, and is not a cure-all that allows for tetrating all complex bases outside the STR with the possible exception of the possibly singular base 0.

Mike,

Very nice!

For base=3+i, I can calculate two repelling fixed points, each with different periods. Presumably, the upper half of the complex plane converges towards the first fixed point, and the lower half of the complex plane converges towards the second fixed point, while the merged function allows you to maintaining the definition of sexp(0)=1, and sexp(1)=base.
L1=0.324536386411256 + 1.00148180609593*I
L2=-0.0976763924712228 - 1.39768129114470*I

I assume the algorithm works when both bases are repelling?
- Sheldon
#3
Very nice, thank you! I tried it and the Pari/GP code worked immediately.

I've nothing more to say at the moment, I think I'll play a bit with it to get a feel. When I find anything concerning your question I'll replay again, may be this will take some days.

Gottfried
Gottfried Helms, Kassel


Possibly Related Threads…
Thread Author Replies Views Last Post
  Code to calculate tetration using my method Shanghai46 10 2,151 10/19/2023, 09:15 PM
Last Post: marracco
  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
  Terse Schroeder & Abel function code Daniel 1 1,053 10/16/2022, 07:03 AM
Last Post: Daniel
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,751 07/03/2022, 06:38 AM
Last Post: JmsNxn
  The beta method program JmsNxn 0 1,398 02/25/2022, 03:05 AM
Last Post: JmsNxn
  Kneser-iteration on n-periodic-points (base say \sqrt(2)) Gottfried 11 11,868 05/05/2021, 04:53 AM
Last Post: Gottfried
  Natural complex tetration program + video MorgothV8 1 6,796 04/27/2018, 07:54 PM
Last Post: MorgothV8
  complex base tetration program sheldonison 23 87,330 10/26/2016, 10:02 AM
Last Post: Gottfried
  C++ program for generatin complex map in EPS format MorgothV8 0 5,737 09/17/2014, 04:14 PM
Last Post: MorgothV8



Users browsing this thread: 1 Guest(s)