Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
Attempt to make own implementation of "Kneser" algorithm: trouble
#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
Reply


Messages In This Thread
RE: Attempt to make own implementation of "Kneser" algorithm: trouble - by mike3 - 06/16/2011, 11:48 AM

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



Users browsing this thread: 1 Guest(s)