+- Tetration Forum (https://tetrationforum.org)
+-- Forum: Tetration and Related Topics (https://tetrationforum.org/forumdisplay.php?fid=1)
+--- Forum: Computation (https://tetrationforum.org/forumdisplay.php?fid=8)
+--- Thread: Numerical algorithm for Fourier continuum sum tetration theory (/showthread.php?tid=516)
Pages:
1
2
Numerical algorithm for Fourier continuum sum tetration theory - mike3 - 09/15/2010
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.
RE: Numerical algorithm for Fourier continuum sum tetration theory - mike3 - 09/15/2010
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.
RE: Numerical algorithm for Fourier continuum sum tetration theory - nuninho1980 - 09/15/2010
Code:
(13:19) gp > vt = ConvToFourier(vector(3, i, 1))
%11 = [0.E-29 + 3.365806530 E-29*I, 1.000000000000000000000000000 + 0.E-29*I, 1.
051814540 E-29 + 2.524354897 E-29*I]
(13:19) gp > vt = NormTetIter(128*I, 2., vt)
%12 = [-0.7725460236329082906309061030 - 9.19165959 E-28*I, 0.769145726068130203
5432582902 - 5.202826184 E-28*I, 1.003400297564778087087647813 + 1.456791331 E-2
7*I]
(13:19) gp > TetApprox(PP, vt, 2., 0.5)
*** exp: negative exponent in gexp.therefore, is this variable "x" integer only?
RE: Numerical algorithm for Fourier continuum sum tetration theory - mike3 - 09/15/2010
(09/15/2010, 01:28 PM)nuninho1980 Wrote:Code:(13:19) gp > vt = ConvToFourier(vector(3, i, 1))
%11 = [0.E-29 + 3.365806530 E-29*I, 1.000000000000000000000000000 + 0.E-29*I, 1.
051814540 E-29 + 2.524354897 E-29*I]
(13:19) gp > vt = NormTetIter(128*I, 2., vt)
%12 = [-0.7725460236329082906309061030 - 9.19165959 E-28*I, 0.769145726068130203
5432582902 - 5.202826184 E-28*I, 1.003400297564778087087647813 + 1.456791331 E-2
7*I]
(13:19) gp > TetApprox(PP, vt, 2., 0.5)
*** exp: negative exponent in gexp.
therefore, is this variable "x" integer only?
It should take any real or complex value. "PP" is not set, that should be the period (note in my post I mentioned "PP = 128*I"), in this case 128*I. This is what is causing your problem.
Though I don't think this'll work with just 3 elements and such a huge period because there won't be enough "resolution" in the Fourier series for it to home in on the function. You should either increase the elements or decrease the PP, say 8*I instead. If you want to use 128*I, try at least 81 (3^4) terms, but even there a smaller PP like 32*I would be best. Usually I use 128*I with at least 243 terms.
(I just tried it. 3 terms at 128*I is not enough.)
RE: Numerical algorithm for Fourier continuum sum tetration theory - nuninho1980 - 09/16/2010
(09/15/2010, 07:17 PM)mike3 Wrote: It should take any real or complex value. "PP" is not set, that should be the period (note in my post I mentioned "PP = 128*I"), in this case 128*I. This is what is causing your problem.
Though I don't think this'll work with just 3 elements and such a huge period because there won't be enough "resolution" in the Fourier series for it to home in on the function. You should either increase the elements or decrease the PP, say 8*I instead. If you want to use 128*I, try at least 81 (3^4) terms, but even there a smaller PP like 32*I would be best. Usually I use 128*I with at least 243 terms.
(I just tried it. 3 terms at 128*I is not enough.)
- vt = ConvToFourier(vector(3^4, i, 1))
- vt = NormTetIter(32*I, 2., vt)
- TetApprox(32*I, vt, 2., 0.5) -> 1.636671277861037786181970844 + ...I
2^^0.5 = 1.63667 - it's incorrect but yes 2^^0.5 = 1.45878181603642 (new knesser code by sheldonison)
RE: Numerical algorithm for Fourier continuum sum tetration theory - mike3 - 09/16/2010
(09/16/2010, 01:18 AM)nuninho1980 Wrote:(09/15/2010, 07:17 PM)mike3 Wrote: It should take any real or complex value. "PP" is not set, that should be the period (note in my post I mentioned "PP = 128*I"), in this case 128*I. This is what is causing your problem.
Though I don't think this'll work with just 3 elements and such a huge period because there won't be enough "resolution" in the Fourier series for it to home in on the function. You should either increase the elements or decrease the PP, say 8*I instead. If you want to use 128*I, try at least 81 (3^4) terms, but even there a smaller PP like 32*I would be best. Usually I use 128*I with at least 243 terms.
(I just tried it. 3 terms at 128*I is not enough.)
- vt = ConvToFourier(vector(3^4, i, 1))
- vt = NormTetIter(32*I, 2., vt)
- TetApprox(32*I, vt, 2., 0.5) -> 1.636671277861037786181970844 + ...I
2^^0.5 = 1.63667 - it's incorrect but yes 2^^0.5 = 1.45878181603642 (new knesser code by sheldonison)
Try running NormTetIter a few more times over it. Like
for(i=1,5,vt = NormTetIter(32*I, 2., vt));
RE: Numerical algorithm for Fourier continuum sum tetration theory - nuninho1980 - 09/17/2010
(09/16/2010, 08:00 PM)mike3 Wrote: Try running NormTetIter a few more times over it. Like1.458810483283048116375249186 + 8.53713983 E-27*I
for(i=1,5,vt = NormTetIter(32*I, 2., vt));
thank!
but it's very probably better code by sheldonison.RE: Numerical algorithm for Fourier continuum sum tetration theory - nuninho1980 - 09/17/2010
your code:
Code:
- vt = ConvToFourier(vector(3^4, i, 1));
- for(i=1,5,vt = NormTetIter(32*I, exp(1), vt));
- TetApprox(32*I, vt, exp(1), Pi)code by sheldonison:
Code:
- init(exp(1.))
base 2.71828182845905
fixed point 0.318131505204764 + 1.33723570143069*I
iterations 126 used for superf/isuperf
- loop(7)
2 ctrm(s) out of 50 sexp(z) generates 800 Riemann samples
60 rtrm(s) out of 400 riemaprx(z) generates 100 sexp samples
sexp(-0.5)= 0.498455596179694
8.16051791246692 Riemann/sexp binary precision bits I=0.0157073173118207*I
-0.00163897418270385 recenter/renorm 0.00120350507213262
14 ctrm(s) out of 50 sexp(z) generates 800 Riemann samples
64 rtrm(s) out of 400 riemaprx(z) generates 100 sexp samples
sexp(-0.5)= 0.498564041905756
15.5432727261501 Riemann/sexp binary precision bits I=0.0157073173118207*I
0.0000252463454086706 recenter/renorm -0.00000302943456839241
50 ctrm(s) out of 50 sexp(z) generates 800 Riemann samples
115 rtrm(s) out of 400 riemaprx(z) generates 100 sexp samples
sexp(-0.5)= 0.498563283958465
22.9408135385566 Riemann/sexp binary precision bits I=0.0157073173118207*I
-0.000000152690627849965 recenter/renorm 0.00000000801647392410488
26 ctrm(s) out of 50 sexp(z) generates 800 Riemann samples
164 rtrm(s) out of 400 riemaprx(z) generates 100 sexp samples
sexp(-0.5)= 0.498563287960307
30.5562661278055 Riemann/sexp binary precision bits I=0.0157073173118207*I
0.000000000772404436364769 recenter/renorm -2.53031553890253 E-11
31 ctrm(s) out of 50 sexp(z) generates 800 Riemann samples
213 rtrm(s) out of 400 riemaprx(z) generates 100 sexp samples
sexp(-0.5)= 0.498563287941024
38.2678540771119 Riemann/sexp binary precision bits I=0.0157073173118207*I
-3.66939771466542 E-12 recenter/renorm 9.18690053231789 E-14
42 ctrm(s) out of 50 sexp(z) generates 800 Riemann samples
264 rtrm(s) out of 400 riemaprx(z) generates 100 sexp samples
sexp(-0.5)= 0.498563287941115
46.0107842465053 Riemann/sexp binary precision bits I=0.0157073173118207*I
1.70948369274334 E-14 recenter/renorm -3.82049499757912 E-16
47 ctrm(s) out of 50 sexp(z) generates 800 Riemann samples
315 rtrm(s) out of 400 riemaprx(z) generates 100 sexp samples
sexp(-0.5)= 0.498563287941114
53.8101076400847 Riemann/sexp binary precision bits I=0.0157073173118207*I
-7.64624004450182 E-17 recenter/renorm 7.97835112074203 E-19
- sexp(Pi)your code may be bad!?
RE: Numerical algorithm for Fourier continuum sum tetration theory - mike3 - 09/17/2010
The problem is the initial guess. I mentioned in the main posts:
Quote: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.
I wasn't specific enough -- the base has to be real small, even e seems too large, at this initial guess. Not sure why this is. Try a loop or two on base 2 first, then use the result of that as the initial guess to run base e. On my system, a single loop at base 2, followed by six loops at base e, using the parameters you give, yields tet(pi) ~ 3.716e10, which is closer to the expected value of 3.71498e10. You'll probably need more than 81 nodes and 32*I period to get more accuracy than that.
Yes, the Kneser process seems to yield more precision more quickly, but is restricted only to real bases \( b > e^{1/e} \) (where the conjugate symmetry holds).
ADDENDUM: Actually, it does work in your case, but with more nodes than what you used, it won't. You just need to iterate it a lot more times -- try 30 or so.
RE: Numerical algorithm for Fourier continuum sum tetration theory - sheldonison - 09/17/2010
Mike, thanks for your detailed description. So, am I correct that your approximating the sexp(z) with a long periodic function, to approximate switching from the space domain to the frequency domain. Then, to get a more accurate version of the sexp(z), you take sexp(z+1)=e^sexp(z), using Faà di Bruno's formula? I mean, high level overview, is that more or less correct?
I do remember convolutions, and converting the time domain to the frequency domain, but its been 27 or 28 years.... So much more to learn, and yet so much already forgotten...
So, as to the general applicability to complex domains, and getting other solutions then we're used to seeing -- here's my intuitive feeling. Any analytic solution, especially one with limiting behavior matching the super function, is a 1-cyclic transformation, via theta(z), of the superfunction.
f(z)=regularsuper(z+theta(z)).
To preserve the behavior at either +I*infinity or -I*infinity, then theta(z) must go to zero at +I or -I infinity. Then, I believe it is most likely that there will be a singularity in theta(z), and we usually get around the singularity with a schwarz transformation about the real axis. Maybe something similar applies in the complex domain. The other possibility, is that theta(z), when wrapped around the unit circle, is an annular Laurent series, with convergence between the two radius's of singularities. This hasn't been explored (to my knowledge), but I think that's what would be going on when converting between the regularsuperf_sqrt(2)(z) and the "regular iteration of log" developed from the other fixed point for sqrt(2). That's what I'd like to explore when I have time.
Nininho, I found Andrew's tetration website via the wayback archive, Andrew's Website. I remember those exp^^pi calculations! I will add one more number to the list...
e^^pi=37149801960.556985498745240109389505725500669657852560614478645123...
This one is accurate to 65 decimal digits, with the Taylor series calculated in about 30 minutes, using the latest version of my kneser.gp code. I'll probably update the code online again; the current version is getting almost 2x more precision for the same amount of compute time (100 binary bits precision in under a minute), but it slows way down once the number of iterated logs used to calculate the inverse superf gets out of hand, so I was thinking of updating that part of the algorithm as well. Mostly, I'd like to figure out how to prove it converges...
- Sheldon