Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
Single-exp series computation code
#1
tommy1729 wanted the details of the method using the single-exp series. I'd post it anyway, but here you go.

Periodic interpolation
The idea is this. If we denote by , a finite sequence of points () to be interpolated, then we can find a periodic interpolation by taking their discrete Fourier transform (DFT):



We can find the points interpolated by taking the inverse discrete Fourier transform:



If is odd, then the interpolation is constructed by

.

where for outside is obtained by periodic repetition, i.e. . The values of will equals when for integers .

If we want period ,

.

In any case, will be . If we want our nodes to represent values on either side of 0, with the middle node being 0, then we just rotate them by the appropriate factor () before interpolating.

The discrete Fourier transform can be done much faster using the Fast Fourier Transform (FFT) algorithm. The current program uses only the "brute force" method, since it was just an initial test to see if this route was any good.

Periodizing functions

The next thing we need is a "periodizing function". The tetration is not periodic, so we must approximate it with a periodic function of long period. A "periodizing function" is a function such that

1. has period .
2. .

The idea is that given some non-periodic function , is a periodic function approximately equal to for near 0. Increasing the period makes the approximation better.

For example, is one possibility. For the program, we use the slightly better .

The algorithm

We can put these two together to get the following algorithm:

1. Create a Fourier series array with some odd number of nodes and an initial guess.

2. Compute the continuum sum of the array. Store the linear term separately from the new exp/Fourier series.

3. Compute the exponential to the base by converting to function space (via IDFT), exponentiating each node, and converting back to Fourier space (via DFT). This may not be the most precise way, and is not very precise for very small numbers of nodes. However it can be done using the FFT if that is available. (Question: Are there any better methods, perhaps with more precision?). Exponentiate the linear term separately (implicit, just interpret its separetly-value from before not as but and as a multiplier instead of an addend).

4. Multiply by . This is just "implicit" (just treat the array after this point as if it has a multiplier).

5. Integrate with lower bound -1. The integral of the Fourier terms is .

6. Periodize the result by substituting the periodizing function into the newly-constructed interpolation, evaluating the result at the interpolation nodes' absicssas, and converting it back to Fourier space. (This can also be done using the FFT algorithm, but we should use the DIT/Cooley-Tukey method.)

7. Normalize. Divide the result by its evaluation at 0.

8. Repeat steps 2-7 until convergence is attained.

The program

The implementation of the algorithm is given below, in PARI/GP.

Code:
/* Perform discrete Fourier transform (DFT). */
/* We should rewrite this to use FFT */
DFT(v) = {
        local(N = matsize(v)[2]);
        local(rv = v);

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

        return(rv);
}

/* Perform inverse discrete Fourier transform (IDFT). */
/* We should rewrite this to use FFT */
IDFT(v) = {
        local(N = matsize(v)[2]);
        local(rv = v);

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

        return(rv);
}

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

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

/* Multiply continuum-sum exp array by log(base)^x and integrate from -1 to x. */
IntegrationPhase(P, v) = {
        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, create an initial guess with

Code:
tv = ConvToFourier(vector(N, i, 1))

where N should be odd, say N = 129. Decide on some period. I use pure imaginary periods. Assign it to P, say P = 48*I.

Then just iterate away with

Code:
tv = NormTetIter(P, base, tv)

where base is the desired base, like 2, or exp(1).

When you want to use the approximation to evaluate tetration, use TetApprox(P, tv, base, x) with x being the height.

For complex bases, a flat initial guess does not seem to yield convergence. However, starting with a real base's tetration, and then using that as the initial guess for a base somewhat away from it in the imaginary direction, and then using that as the guess for another, and so on, seems to work better. I haven't yet been able to push this into the left halfplane, not sure why.
Reply


Possibly Related Threads...
Thread Author Replies Views Last Post
  C++ code for tet, ate and hexp MorgothV8 0 2,700 07/10/2014, 04:24 PM
Last Post: MorgothV8
  "Kneser"/Riemann mapping method code for *complex* bases mike3 2 6,134 08/15/2011, 03:14 PM
Last Post: Gottfried
  A note on computation of the slog Gottfried 6 8,800 07/12/2010, 10:24 AM
Last Post: Gottfried
  Computations with the double-exp series mike3 0 2,252 04/20/2010, 07:32 PM
Last Post: mike3
  SAGE code for computing flow matrix for exp(z)-1 jaydfox 4 7,542 08/21/2009, 05:32 PM
Last Post: jaydfox
  An error estimate for fixed point computation of b^x bo198214 0 2,348 05/31/2008, 04:11 PM
Last Post: bo198214
  SAGE code implementing slog with acceleration jaydfox 4 6,535 10/22/2007, 12:59 AM
Last Post: jaydfox
  Complete SAGE code for tetration eyu100 4 5,582 10/18/2007, 03:55 AM
Last Post: jaydfox



Users browsing this thread: 1 Guest(s)