Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
"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 :

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 , representing what may be some of the most accurate calculations of a tetrational of a complex base to date. We get



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

We can even do the "pentation" (the 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.
Reply
#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 :
....
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
Reply
#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
Reply


Possibly Related Threads...
Thread Author Replies Views Last Post
  complex base tetration program sheldonison 23 20,733 10/26/2016, 10:02 AM
Last Post: Gottfried
  fast accurate Kneser sexp algorithm sheldonison 38 45,124 01/14/2016, 05:05 AM
Last Post: sheldonison
  C++ program for generatin complex map in EPS format MorgothV8 0 1,210 09/17/2014, 04:14 PM
Last Post: MorgothV8
  C++ code for tet, ate and hexp MorgothV8 0 1,657 07/10/2014, 04:24 PM
Last Post: MorgothV8
  Green Eggs and HAM: Tetration for ALL bases, real and complex, now possible? mike3 27 14,437 07/02/2014, 10:13 PM
Last Post: tommy1729
  Which method is currently "the best"? MorgothV8 2 2,217 11/15/2013, 03:42 PM
Last Post: MorgothV8
  Attempt to make own implementation of "Kneser" algorithm: trouble mike3 9 9,285 06/16/2011, 11:48 AM
Last Post: mike3
  An incremental method to compute (Abel) matrix inverses bo198214 3 6,414 07/20/2010, 12:13 PM
Last Post: Gottfried
  Single-exp series computation code mike3 0 1,999 04/20/2010, 08:59 PM
Last Post: mike3
  Attempting to compute the kslog numerically (i.e., Kneser's construction) jaydfox 11 11,155 10/26/2009, 05:56 PM
Last Post: bo198214



Users browsing this thread: 1 Guest(s)