04/20/2010, 08:59 PM
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.
To use, create an initial guess with
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
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.
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.