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.
To use this, first create an array of nodes for an initial guess with
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,
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.
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.