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 \( x_n \), a finite sequence of \( N \) points (\( n = 0 ... N-1 \)) to be interpolated, then we can find a periodic interpolation by taking their discrete Fourier transform (DFT):

\( X_k = \sum_{n=0}^{N-1} x_n e^{-2\pi i n k/N} \)

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

\( x_n = \frac{1}{N} \sum_{k=0}^{N-1} x_k e^{2\pi i n k/N} \)

If \( N \) is odd, then the interpolation is constructed by

\( f(z) =\ \sum_{k=-(N-1)/2}^{(N-1)/2} \frac{1}{N} X_k e^{2\pi i k z/N} \).

where \( X_k \) for \( k \) outside \( k = 0...N-1 \) is obtained by periodic repetition, i.e. \( X_k = X_{k\ \mathrm{mod}\ N} \). The values of \( f(z) \) will equals \( x_k \) when \( x = k \) for integers \( k \).

If we want period \( P \),

\( f(z) =\ \sum_{k=-(N-1)/2}^{(N-1)/2} \frac{1}{N} X_k e^{\frac{2\pi}{P} i k z/N} \).

In any case, \( f(0) \) will be \( x_0 \). If we want our nodes \( x_k \) to represent values on either side of 0, with the middle node being 0, then we just rotate them by the appropriate factor (\( \frac{N-1}{2} \)) 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 \( P_r(z) \) such that

1. \( P_r(z) \) has period \( r \).
2. \( \lim_{r \rightarrow \infty} P_r(z) = z \).

The idea is that given some non-periodic function \( f(z) \), \( f(P_r(z)) \) is a periodic function approximately equal to \( f \) for \( z \) near 0. Increasing the period \( r \) makes the approximation better.

For example, \( P_r(z) = \frac{r}{2\pi} \sin \left(\frac{2\pi}{r} z\right) \) is one possibility. For the program, we use the slightly better \( P_r(z) = \frac{r}{2\pi} 1.25 \sin\left(\frac{2\pi}{r} z\right) - \frac{r}{2\pi} 0.125 \sin\left(\frac{2\pi}{r} 2z\right) \).

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 \( b \) 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 \( a_0 \) from before not as \( a_0 z \) but \( b^{a_0 z} \) and as a multiplier instead of an addend).

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

5. Integrate with lower bound -1. The integral of the Fourier terms \( a_n e^{(a_0 \log(b) + \log(\log(b)) + nu)z} \) is \( \frac{a_n}{a_0 \log(b) + \log(\log(b)) + nu} \left(e^{(a_0 \log(b) + \log(\log(b)) + nu)z} - e^{-(a_0 \log(b) + \log(\log(b)) + nu)}\right) \).

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.


Possibly Related Threads…
Thread Author Replies Views Last Post
  Code to calculate tetration using my method Shanghai46 10 2,143 10/19/2023, 09:15 PM
Last Post: marracco
  Writing Kneser's super logarithm using values of Kneser at a single point JmsNxn 1 843 04/21/2023, 04:26 AM
Last Post: JmsNxn
  Terse Schroeder & Abel function code Daniel 1 1,048 10/16/2022, 07:03 AM
Last Post: Daniel
  Mathematica program for tetration based on the series with q-binomial coefficients Vladimir Reshetnikov 0 5,465 01/13/2017, 10:51 PM
Last Post: Vladimir Reshetnikov
  C++ code for tet, ate and hexp MorgothV8 0 5,939 07/10/2014, 04:24 PM
Last Post: MorgothV8
  "Kneser"/Riemann mapping method code for *complex* bases mike3 2 12,669 08/15/2011, 03:14 PM
Last Post: Gottfried
  A note on computation of the slog Gottfried 6 19,433 07/12/2010, 10:24 AM
Last Post: Gottfried
  Computations with the double-exp series mike3 0 4,794 04/20/2010, 07:32 PM
Last Post: mike3
  SAGE code for computing flow matrix for exp(z)-1 jaydfox 4 16,853 08/21/2009, 05:32 PM
Last Post: jaydfox
  A nice series for b^^h , base sqrt(2), by diagonalization Gottfried 19 45,361 06/11/2009, 08:36 PM
Last Post: Gottfried



Users browsing this thread: 1 Guest(s)