Hi.
Here's the long-awaited description of the numerical algorithm for the Fourier-continuum-sum-based tetration.
It actually took a LOT more to describe than I thought it would...
Periodic approximation
As was mentioned in the posts on the subject, the continuum sum for a holomorphic function is constructed as a limit of continuum sums of periodic approximations of it, that is, a sequence of periodic functions \( f_n(z) \) which converge pointwise to the target function \( f(z) \).
The idea behind the numerical algorithm, then, is to use a periodic function with a very large period as an approximation to a general function in the vicinity of zero.
That is, we work with an approximation to the desired function given by
\( f_{\mathrm{approx}}(z) = \sum_{n=-\infty}^{\infty} a_n e^{\frac{2\pi}{P} i nz} \)
(which in a computer would be truncated to a sum from \( -N \) to \( N \) for some \( N \))
where \( P \) is the period, and large.
Because certain operations in the tetration continuum-sum iteration destroy periodicity (especially the multiplication by \( \log(b)^w \) and the final integral), we must have a way to restore it. We do this via what is called a periodizing function. A periodizing function \( p_P(z) \) generates a periodic approximation \( f_{\mathrm{approx}}(z) \) of a given function \( f(z) \) by composition, i.e. \( f_{\mathrm{approx}}(z) = f(p_P(z)) \), with period \( P \). This means the function should, at bare minimum, satisfy
1. periodic with period \( P \): \( p_P(z + P) = p_P(z) \)
2. \( \forall z \in \mathbb{C}, \lim_{|P| \rightarrow \infty} p_P(z) = z \)
Other desirable properties are
3. entire
4. for any compact subset \( K \) of \( \mathbb{C} \), and a given \( \epsilon > 0 \), there exists an \( M > 0 \) such that \( |p_P(z) - z| < \epsilon \) whenever \( |P| > M \) for all \( z \in K \).
An example of such a function is \( p_P(z) = \frac{P}{2\pi} \sin\left(\frac{2\pi}{P} z\right) \). Property 4 is sort of like compact convergence of a sequence of functions but with a complex index, and is nice because it means that past a certain point, we can just increase \( P \)'s magnitude to make \( p_P \) a better approximation of the identity around 0.
Implementing the tetration formula
We can now discuss how we'd implement the iteration. We work with an approximation of the form
\( F_{\mathrm{approx}}(z) = \sum_{n=-H}^{H} a_n e^{\frac{2\pi}{P} i nz} \).
(the reason for using \( H \) instead of \( N \) will be clear later)
and perform the following operations.
Continuum sum
The continuum sum, \( \Sigma \), is achieved by setting
\( b_n = \frac{a_n}{e^{\frac{2\pi}{P} i n} - 1} \), \( n \ne 0 \)
\( b_0 = -\sum_{n=-H}^{H} \frac{\delta_{n0} a_n}{e^{\frac{2\pi}{P} i n - (1 - \delta_{n0})} \).
(from \( \sum_{n=0}^{z-1} e^{un} = \frac{e^{un} - 1}{e^u - 1} \).)
And so we have \( \sum_{n=0}^{z-1} F_{\mathrm{approx}}(n) = a_0 z + \sum_{n=-H}^{H} b_n e^{\frac{2\pi}{P} i nz} \). We can store \( a_0 \) separately for the next operation.
Exponential
The exponential of the continuum sum is done by treating the series as a pair of Taylor series in \( e^{\frac{2\pi}{P} i z} \) and \( e^{-\frac{2\pi}{P} i z} \), and then applying Faà di Bruno's formula and multiplying the results via the convolution, however this method is not very fast and seems to not be as stable as we might like (though keeping more "slop" terms beyond -N and N index seems to help, but that just makes it even slower).
Though with a good program like Pari/GP, the formal exponential of a power series is already implemented and so there's no need to implement all that ourselves.
However, we can do a lot better.
Fourier transform
A Fourier series, such as the one above, can be thought of as a sort of "Fourier interpolation" between points in the "spatial" domain. It is possible to transform from one to the other and back via the discrete Fourier transform, or DFT. Let's see how that works.
We first transform the function so it has period \( N = 2H + 1 \) (note now that "H" is for "half" as in half of N, rounded down), giving
\( F^{*}_{\mathrm{approx}}(z) = \sum_{n=-H}^{H} b_n e^{\frac{2\pi}{N} i nz} \).
Now consider the "discrete sampling" of this at integer inputs, giving a sequence indexed from \( -H \) to \( H \):
\( x_k = \sum_{n=-H}^{H} b_n e^{\frac{2\pi}{N} i nk} \).
This bears a great similarity to the formula for the inverse discrete Fourier transform (IDFT) of a sequence \( X_n \):
\( x_k = \frac{1}{N} \sum_{n=0}^{N-1} X_n e^{\frac{2\pi}{N} i nk} \).
The only difference is the frequency elements involved, and the presence of the normalization factor \( \frac{1}{N} \). This means that our Fourier series can be thought of as an interpolation of a regularly-spaced grid of points along one period interval. Our sampling at discrete points then creates a list of all those points, to which we can then directly apply transcendental functions like exponential. This method may not be as accurate as the first (note the limited amount of terms), but it can be much faster and also seems numerically more stable and with enough terms to begin with, it can still be quite accurate.
If we can relate the formula to the IDFT, then we can also give the transform to get back to Fourier space, and, even more powerfully, we can harness the high-speed Fast Fourier Transform (FFT) algorithm to rapidly perform this computation, blowing the other method out of the water.
So how do we relate this to the IDFT? Consider the following process: First off, we renumber the elements of \( b_n \) from \( -H...H \) to \( 0...N-1 \), giving:
\( x_k = \sum_{n=0}^{N-1} b^{*}_n e^{\frac{2\pi}{N} i (n - H) k} \).
where \( a^{*}_n \) is the renumbered sequence.
Note the exponent in the exponential -- the displacement by H. Expanding this out, we see that we have
\( x_k = \sum_{n=0}^{N-1} b^{*}_n e^{-\frac{2\pi}{N} i H k} e^{\frac{2\pi}{N} i nk} \).
\( x_k e^{\frac{2\pi}{N} i H k} = \sum_{n=0}^{N-1} b^{*}_n e^{\frac{2\pi}{N} i nk} \).
and thus a regular IDFT (sans the normalization factor) on the right.
However, by the Fourier shift theorem, the multiplication of the spatial-domain data by \( e^{\frac{2\pi}{N} i H k} \) is equivalent a cyclical shift of \( b^{*}_n \) by H. Namely, it corresponds to transforming \( b^{*}_{(n - H) \mathrm{mod} N} \). We can cancel the effects of this shift by cycle-shifting \( b^{*}_n \) the other way, forming \( b^{**}_n = b^{*}_{(n + H) \mathrm{mod} N} \) and taking the inverse DFT of it multiplied by \( N \) (to account for the normalization factor).
There's one final point: the given formula was for \( x_k \) numbered from \( -H \) to \( H \). However, the (I)DFT is usually done with the numbering from \( 0 \) to \( N-1 \). This is not a problem: if we use a regular (I)DFT algorithm, we'll get \( x_0 \) through \( x_H \) followed by \( x_{-H} \) through \( x_{-1} \). A simple right cyclic shift by \( H \) resolves this and gives the data in the correct order, though such is not really needed for evaluating transcendental functions.
Then to transform from the spatial domain back to the frequency domain, we do all that in reverse: left cycle-shift by \( H \), preform forward DFT, divide by \( N \), then cycle-shift the result right by \( H \) and renumber.
Log-power of base and final integral
In the given program, I lump these two steps together. The end result of the last step should be a series of coefficients \( c_n \) giving an approximation to the exponential of the continuum sum of \( F_{\mathrm{approx}} \):
\( \exp_b\left(\sum_{n=0}^{z-1} F_{\mathrm{approx}}(n)\right) \approx b^{a_0 z} \sum_{n=-H}^{H} c_n e^{\frac{2\pi}{P} i nz} \).
We can now multiply by \( \log(b)^z \), giving
\( \begin{align}
\log(b)^z \exp_b\left(\sum_{n=0}^{z-1} F_{\mathrm{approx}}(n)\right) &\approx \log(b)^z b^{a_0 z} \sum_{n=-H}^{H} c_n e^{\frac{2\pi}{P} i nz} \\
&= \sum_{n=-H}^{H} c_n e^{(\frac{2\pi}{P} i n + \log(\log(b)) + \log(b) a_0) z}
\end{align} \)
and integrate, to get
\( \int_{-1}^{z} \log(b)^z \exp_b\left(\sum_{n=0}^{z-1} F_{\mathrm{approx}}(n)\right) dz \approx \sum_{n=-H}^{H} \frac{c_n}{\frac{2\pi}{P} i n + \log(\log(b)) + \log(b) a_0} \left(e^{(\frac{2\pi}{P} i n + \log(\log(b)) + \log(b) a_0) z} - e^{-\frac{2\pi}{P} i n - \log(\log(b)) - \log(b) a_0}\right) \).
(YES, I know I'm abusing the integral notation with the "z" in both the limits and the variable, but it's to try to make it clearer and more consistent.)
We can then composite this with the periodizing function to get a periodic function out the end with the period we desired, and a pointwise evaluation and DFT will give us the new Fourier series for the next iteration cycle. The normalization step (for the derivative-at-0 bit) is real easy, just evaluate at 0 and divide the thing by the result.
Here's the long-awaited description of the numerical algorithm for the Fourier-continuum-sum-based tetration.
It actually took a LOT more to describe than I thought it would...Periodic approximation
As was mentioned in the posts on the subject, the continuum sum for a holomorphic function is constructed as a limit of continuum sums of periodic approximations of it, that is, a sequence of periodic functions \( f_n(z) \) which converge pointwise to the target function \( f(z) \).
The idea behind the numerical algorithm, then, is to use a periodic function with a very large period as an approximation to a general function in the vicinity of zero.
That is, we work with an approximation to the desired function given by
\( f_{\mathrm{approx}}(z) = \sum_{n=-\infty}^{\infty} a_n e^{\frac{2\pi}{P} i nz} \)
(which in a computer would be truncated to a sum from \( -N \) to \( N \) for some \( N \))
where \( P \) is the period, and large.
Because certain operations in the tetration continuum-sum iteration destroy periodicity (especially the multiplication by \( \log(b)^w \) and the final integral), we must have a way to restore it. We do this via what is called a periodizing function. A periodizing function \( p_P(z) \) generates a periodic approximation \( f_{\mathrm{approx}}(z) \) of a given function \( f(z) \) by composition, i.e. \( f_{\mathrm{approx}}(z) = f(p_P(z)) \), with period \( P \). This means the function should, at bare minimum, satisfy
1. periodic with period \( P \): \( p_P(z + P) = p_P(z) \)
2. \( \forall z \in \mathbb{C}, \lim_{|P| \rightarrow \infty} p_P(z) = z \)
Other desirable properties are
3. entire
4. for any compact subset \( K \) of \( \mathbb{C} \), and a given \( \epsilon > 0 \), there exists an \( M > 0 \) such that \( |p_P(z) - z| < \epsilon \) whenever \( |P| > M \) for all \( z \in K \).
An example of such a function is \( p_P(z) = \frac{P}{2\pi} \sin\left(\frac{2\pi}{P} z\right) \). Property 4 is sort of like compact convergence of a sequence of functions but with a complex index, and is nice because it means that past a certain point, we can just increase \( P \)'s magnitude to make \( p_P \) a better approximation of the identity around 0.
Implementing the tetration formula
We can now discuss how we'd implement the iteration. We work with an approximation of the form
\( F_{\mathrm{approx}}(z) = \sum_{n=-H}^{H} a_n e^{\frac{2\pi}{P} i nz} \).
(the reason for using \( H \) instead of \( N \) will be clear later)
and perform the following operations.
Continuum sum
The continuum sum, \( \Sigma \), is achieved by setting
\( b_n = \frac{a_n}{e^{\frac{2\pi}{P} i n} - 1} \), \( n \ne 0 \)
\( b_0 = -\sum_{n=-H}^{H} \frac{\delta_{n0} a_n}{e^{\frac{2\pi}{P} i n - (1 - \delta_{n0})} \).
(from \( \sum_{n=0}^{z-1} e^{un} = \frac{e^{un} - 1}{e^u - 1} \).)
And so we have \( \sum_{n=0}^{z-1} F_{\mathrm{approx}}(n) = a_0 z + \sum_{n=-H}^{H} b_n e^{\frac{2\pi}{P} i nz} \). We can store \( a_0 \) separately for the next operation.
Exponential
The exponential of the continuum sum is done by treating the series as a pair of Taylor series in \( e^{\frac{2\pi}{P} i z} \) and \( e^{-\frac{2\pi}{P} i z} \), and then applying Faà di Bruno's formula and multiplying the results via the convolution, however this method is not very fast and seems to not be as stable as we might like (though keeping more "slop" terms beyond -N and N index seems to help, but that just makes it even slower).
Though with a good program like Pari/GP, the formal exponential of a power series is already implemented and so there's no need to implement all that ourselves.
However, we can do a lot better.
Fourier transform
A Fourier series, such as the one above, can be thought of as a sort of "Fourier interpolation" between points in the "spatial" domain. It is possible to transform from one to the other and back via the discrete Fourier transform, or DFT. Let's see how that works.
We first transform the function so it has period \( N = 2H + 1 \) (note now that "H" is for "half" as in half of N, rounded down), giving
\( F^{*}_{\mathrm{approx}}(z) = \sum_{n=-H}^{H} b_n e^{\frac{2\pi}{N} i nz} \).
Now consider the "discrete sampling" of this at integer inputs, giving a sequence indexed from \( -H \) to \( H \):
\( x_k = \sum_{n=-H}^{H} b_n e^{\frac{2\pi}{N} i nk} \).
This bears a great similarity to the formula for the inverse discrete Fourier transform (IDFT) of a sequence \( X_n \):
\( x_k = \frac{1}{N} \sum_{n=0}^{N-1} X_n e^{\frac{2\pi}{N} i nk} \).
The only difference is the frequency elements involved, and the presence of the normalization factor \( \frac{1}{N} \). This means that our Fourier series can be thought of as an interpolation of a regularly-spaced grid of points along one period interval. Our sampling at discrete points then creates a list of all those points, to which we can then directly apply transcendental functions like exponential. This method may not be as accurate as the first (note the limited amount of terms), but it can be much faster and also seems numerically more stable and with enough terms to begin with, it can still be quite accurate.
If we can relate the formula to the IDFT, then we can also give the transform to get back to Fourier space, and, even more powerfully, we can harness the high-speed Fast Fourier Transform (FFT) algorithm to rapidly perform this computation, blowing the other method out of the water.
So how do we relate this to the IDFT? Consider the following process: First off, we renumber the elements of \( b_n \) from \( -H...H \) to \( 0...N-1 \), giving:
\( x_k = \sum_{n=0}^{N-1} b^{*}_n e^{\frac{2\pi}{N} i (n - H) k} \).
where \( a^{*}_n \) is the renumbered sequence.
Note the exponent in the exponential -- the displacement by H. Expanding this out, we see that we have
\( x_k = \sum_{n=0}^{N-1} b^{*}_n e^{-\frac{2\pi}{N} i H k} e^{\frac{2\pi}{N} i nk} \).
\( x_k e^{\frac{2\pi}{N} i H k} = \sum_{n=0}^{N-1} b^{*}_n e^{\frac{2\pi}{N} i nk} \).
and thus a regular IDFT (sans the normalization factor) on the right.
However, by the Fourier shift theorem, the multiplication of the spatial-domain data by \( e^{\frac{2\pi}{N} i H k} \) is equivalent a cyclical shift of \( b^{*}_n \) by H. Namely, it corresponds to transforming \( b^{*}_{(n - H) \mathrm{mod} N} \). We can cancel the effects of this shift by cycle-shifting \( b^{*}_n \) the other way, forming \( b^{**}_n = b^{*}_{(n + H) \mathrm{mod} N} \) and taking the inverse DFT of it multiplied by \( N \) (to account for the normalization factor).
There's one final point: the given formula was for \( x_k \) numbered from \( -H \) to \( H \). However, the (I)DFT is usually done with the numbering from \( 0 \) to \( N-1 \). This is not a problem: if we use a regular (I)DFT algorithm, we'll get \( x_0 \) through \( x_H \) followed by \( x_{-H} \) through \( x_{-1} \). A simple right cyclic shift by \( H \) resolves this and gives the data in the correct order, though such is not really needed for evaluating transcendental functions.
Then to transform from the spatial domain back to the frequency domain, we do all that in reverse: left cycle-shift by \( H \), preform forward DFT, divide by \( N \), then cycle-shift the result right by \( H \) and renumber.
Log-power of base and final integral
In the given program, I lump these two steps together. The end result of the last step should be a series of coefficients \( c_n \) giving an approximation to the exponential of the continuum sum of \( F_{\mathrm{approx}} \):
\( \exp_b\left(\sum_{n=0}^{z-1} F_{\mathrm{approx}}(n)\right) \approx b^{a_0 z} \sum_{n=-H}^{H} c_n e^{\frac{2\pi}{P} i nz} \).
We can now multiply by \( \log(b)^z \), giving
\( \begin{align}
\log(b)^z \exp_b\left(\sum_{n=0}^{z-1} F_{\mathrm{approx}}(n)\right) &\approx \log(b)^z b^{a_0 z} \sum_{n=-H}^{H} c_n e^{\frac{2\pi}{P} i nz} \\
&= \sum_{n=-H}^{H} c_n e^{(\frac{2\pi}{P} i n + \log(\log(b)) + \log(b) a_0) z}
\end{align} \)
and integrate, to get
\( \int_{-1}^{z} \log(b)^z \exp_b\left(\sum_{n=0}^{z-1} F_{\mathrm{approx}}(n)\right) dz \approx \sum_{n=-H}^{H} \frac{c_n}{\frac{2\pi}{P} i n + \log(\log(b)) + \log(b) a_0} \left(e^{(\frac{2\pi}{P} i n + \log(\log(b)) + \log(b) a_0) z} - e^{-\frac{2\pi}{P} i n - \log(\log(b)) - \log(b) a_0}\right) \).
(YES, I know I'm abusing the integral notation with the "z" in both the limits and the variable, but it's to try to make it clearer and more consistent.)
We can then composite this with the periodizing function to get a periodic function out the end with the period we desired, and a pointwise evaluation and DFT will give us the new Fourier series for the next iteration cycle. The normalization step (for the derivative-at-0 bit) is real easy, just evaluate at 0 and divide the thing by the result.