Numerical algorithm for Fourier continuum sum tetration theory
#2
Conclusion

Putting this all together, we get the following PARI/GP program. It is probably not the clearest thing out there, but it was assembled more "ad hoc" as these ideas came to me, instead of starting with a "slick" plan like the one above.

Code:
/* Radix-3 DIT FFT. */
/* Convert to reversed trinary. */
trinreverse(n, maxtrits) = {
        local(nn = n);
        local(rv = 0);

        for(i=0,maxtrits-1,
            rv = rv * 3;
            rv = rv + (nn%3);
            nn = floor(nn/3);
           );

        return(rv);
}

/* Trit reversal. Assumes v has size 3^n. */
DoTritReversal(v) = {
        local(N = matsize(v)[2]);
        local(Nlg3 = floor(log(N)/log(3) + 0.5));
        local(reversed = vector(N));
        local(Nrev = 0);

        for(n=0,N-1,
            Nrev = trinreverse(n, Nlg3);
            /*print(Nrev);*/
            reversed[Nrev+1] = v[n+1];
           );

        return(reversed);
}

/* Do FFT (radix 3, DIT) */
FFTCore(v, voffs, Npow, isign) = {
        local(N = 3^Npow);
        local(r = vector(3));
        local(rr = vector(3));

        if(N==1, return([v[voffs]]));

        r[1] = FFTCore(v, voffs, Npow-1, isign);
        r[2] = FFTCore(v, voffs + N/3, Npow-1, isign);
        r[3] = FFTCore(v, voffs + 2*N/3, Npow-1, isign);

        for(j=0,2,
            for(k=0,N/3-1,
                r[j+1][k+1] = r[j+1][k+1] * exp(-isign*2*Pi*I*j*k/N);
               );
           );

        for(j=0,2,
            rr[j+1] = sum(k=0,2,exp(-isign*2*Pi*I*k*j/3)*r[k+1]);
           );

        return(concat(concat(rr[1], rr[2]), rr[3]));
}

/* Perform discrete Fourier transform (DFT). */
DFT(v) = {
        return(FFTCore(DoTritReversal(v), 1, floor(log(matsize(v)[2])/log(3) + 0.5), 1));
}

/* Perform inverse discrete Fourier transform (IDFT). */
IDFT(v) = {
        return((1/matsize(v)[2]) * FFTCore(DoTritReversal(v), 1, floor(log(matsize(v)[2])/log(3) + 0.5), -1));
}

/* Convert a node array from real space to Fourier space. Odd sized array only. */
ConvToFourier(v) = {
        local(N = matsize(v)[2]);
        local(vsh = vector(N));
        local(tmp = vector(N));
        local(rv = vector(N));

        /*print(v);*/
        for(n=0,N-1,vsh[n+1] = v[(n+((N-1)/2))%N + 1]);
        /*print(vsh);*/

        tmp = (1/N)*DFT(vsh);

        /*print(tmp);*/
        for(n=-(N-1)/2,(N-1)/2,rv[n+((N-1)/2)+1] = tmp[n%N+1]);
        /*print(rv);*/

        return(rv);
}

/* Convert a node array from Fourier space to real space. Odd sized array only. */
ConvToReal(v) = {
        local(N = matsize(v)[2]);
        local(vsh = vector(N));
        local(tmp = vector(N));
        local(rv = vector(N));

        for(n=0,N-1,vsh[n+1] = v[(n+((N-1)/2))%N+1]);

        tmp = IDFT(N*vsh);

        for(n=0,N-1,rv[n+1] = tmp[(n-((N-1)/2))%N+1]);

        return(rv);
}

/* Use a node array in Fourier space to do the interpolation. Odd sized array only. */
DoInterpolation(v, x) = {
        local(N = matsize(v)[2]);

        return(sum(n=-((N-1)/2),(N-1)/2,v[n+((N-1)/2)+1]*exp(2*Pi*I*n*x/N)));
}

/* Same as above, but with specified period */
DoInterpolationPer(P, v, x) = {
        local(N = matsize(v)[2]);

        return(sum(n=-((N-1)/2),(N-1)/2,v[n+((N-1)/2)+1]*exp(2*Pi*I*n*x/P)));
}

/* Get the list of abscissas for the nodes with a given period */
GetAbscissas(N, P) = {
        return(vector(N, n, (-((N-1)/2) + (n-1)) * (P/N)));
}

/* Increase the number of interpolating points in an approximation */
IncreaseInterp(P, v, new) = {
        local(xcs = GetAbscissas(new, P));

        return(ConvToFourier(vector(new, i, DoInterpolationPer(P, v, xcs[i]))));
}

/* Make the continuum sum of a Fourier-space array with a given period. */
MakeContSum(P, v) = {
        local(N = matsize(v)[2]);
        local(rv = [0, 0]);

        rv[2] = vector(N);

        for(n=-((N-1)/2),(N-1)/2,
            if(n != 0,
               rv[2][n+((N-1)/2)+1] = v[n+((N-1)/2)+1]/(exp(2*Pi*I*n/P) - 1);,
               rv[2][n+((N-1)/2)+1] =
                 -sum(k=-(N-1)/2,(N-1)/2,if(k==0,0,v[k+((N-1)/2)+1]/(exp(2*Pi*I*k/P) - 1)));
              );
           );

        rv[1] = v[((N-1)/2)+1];

        return(rv);
}

/* Evaluate continuum sum array */
EvalContSumArr(P, v, x) = v[1]*x + DoInterpolationPer(P, v[2], x);

/* Perform exponential of a continuum sum array with a given period. */
ExpContSum(P, base, v) = [base, v[1], ConvToFourier(base^ConvToReal(v[2]))];

/* Evaluate continuum sum exp array */
EvalContSumExpArr(P, v, x) = v[1]^(v[2]*x) * DoInterpolationPer(P, v[3], x);

/* The periodizing function */
Periodizer(P, x) = P/(2*Pi) * 1.25 * sin((2*Pi)/P * x) - P/(2*Pi) * 0.125 * sin((2*Pi)/P * 2 * x);
pp(x) = sin(x) + (sin(x)^3)/6 + (3*sin(x)^5)/40;
Periodizer(P, x) = P/(2*Pi) * pp((2*Pi)/P * x);
/*InvPeriodizer(P, x) = P/(2*Pi) * asin((2*Pi)/P * x);*/

/* Multiply continuum-sum exp array by log(base)^x and integrate from -1 to x. */
IntegrationPhase(P, v) = {
        /* This is the most complex part of the algorithm.
         * We need to first multiply the series by log(base)^x, then
         * integrate, and finally periodize that integral back to
         * being a periodic function.
         *
         * Each term is of the form a_n base^(qx) exp(2piinx/P) =
         * a_n exp((2piin/P + q log(base)) x). Multiplying
         * this by log(base)^x gives a_n exp((2piin/P + q log(base) + log(log(base)))x).
         * The integral of these terms from -1 to x is
         * a_n/(2piin/P + q log(base) + log(log(base))) *
         *               (exp((2piin/P + q log(base) + log(log(base)))x) -
         *                exp(-(2piin/P + q log(base) + log(log(base))))).
         *
         * So we create an auxiliary array of a_n/(2piin/P + q log(base) + log(log(base)))
         * and then we evaluate it at the periodized abscissas. There may be a
         * faster algorithm for this that is like the FFT.
         */
        local(N = matsize(v[3])[2]);
        local(vaux = v[3]);
        local(q = v[2]);
        local(base = v[1]);
        local(abscissas = GetAbscissas(N, P));
        local(rv = vector(N));
        local(expf = vector(N));

        for(n=-(N-1)/2,(N-1)/2,expf[n+((N-1)/2)+1] = (2*Pi*I*n/P) + (q*log(base)) + log(log(base)););

        for(n=-(N-1)/2,(N-1)/2,
            vaux[n+((N-1)/2)+1] = vaux[n+((N-1)/2)+1]/expf[n+((N-1)/2)+1];
           );

        for(k=1,N,
            rv[k] = sum(n=-(N-1)/2,(N-1)/2,
                        vaux[n+((N-1)/2)+1]*(exp(expf[n+((N-1)/2)+1]*Periodizer(P, abscissas[k])) -
                                             exp(-expf[n+((N-1)/2)+1]))
                       );
           );

        return(rv);
}

/* Debug: dump log(base)^x times v function */
PreintDbg(P, v, x) = {
        local(N = matsize(v[3])[2]);
        local(vaux = v[3]);
        local(q = v[2]);
        local(base = v[1]);

        return(sum(n=-(N-1)/2,(N-1)/2,
                        vaux[n+((N-1)/2)+1]*exp((2*Pi*I*n/P + q*log(base) + log(log(base)))*x)));
}

/* Do tetration iteration over Fourier array */
NormTetIter(P, base, v) = {
        local(rv);

        rv = ConvToFourier(IntegrationPhase(P, ExpContSum(P, base, MakeContSum(P, v))));

        rv = rv / DoInterpolationPer(P, rv, 0);

        return(rv);
}

/* Tetration routines */
TetCrit(P, v, X) = DoInterpolationPer(P, v, X);
TetApprox(P, v, base, X) = {
        if(real(X) < -0.5, return(log(TetApprox(P, v, base, X+1))/log(base)));
        if(real(X) > 0.5, return(base^TetApprox(P, v, base, X-1)));
        return(TetCrit(P, v, X));
}

To use this, first create an array of nodes for an initial guess with

Code:
vt = ConvToFourier(vector(<power-of-3 number of nodes>, i, 1))

This starts us off with a function that is 1 everywhere. Then choose a period for the approximation, say an imaginary one, like PP = 128*I. Now you can iterate with, say,

Code:
vt = NormTetIter(PP, <base>, vt)

and to use the approximation, try TetApprox(PP, vt, <base>, x) for a given value x to evaluate at.

It seems that with a large number of nodes, the convergence with this initial guess only works for small bases like base 2 -- not sure what's going on there, but that can be used as an initial guess for a wide variety of other bases.

The biggest problem with this implementation is the final integration/periodization phase, which is as dog-slow as a plain DFT and has the same crappy \( O(N^2) \) growth. I have no idea if it could be accelerated in some manner akin to the Fast Fourier Transform.

EDIT: I should mention that I use a different -- perhaps "better", at least for the values of the tetrational near the real axis, periodizing function in the program.


Messages In This Thread
RE: Numerical algorithm for Fourier continuum sum tetration theory - by mike3 - 09/15/2010, 11:06 AM

Possibly Related Threads…
Thread Author Replies Views Last Post
  Road testing Ansus' continuum product formula mike3 40 76,343 07/12/2022, 08:57 AM
Last Post: Gottfried
  fast accurate Kneser sexp algorithm sheldonison 40 145,748 07/03/2022, 06:38 AM
Last Post: JmsNxn
  Attempt to make own implementation of "Kneser" algorithm: trouble mike3 9 30,738 06/16/2011, 11:48 AM
Last Post: mike3
  Software for numerical calculation ala qbasic Ivars 3 11,298 02/22/2008, 04:40 PM
Last Post: Ivars
  numerical examples for approximation for half-iterates for exp(x)-1 Gottfried 0 5,082 08/14/2007, 11:57 AM
Last Post: Gottfried



Users browsing this thread: 1 Guest(s)