Hi.
I've been trying to build my own program in Pari/GP to run the Kneser algorithm for tetration, but I can't seem to get it to work so well. The goal is to ultimately be able to generalize it to merge complex bases and bases less than \( \eta \), especially the legendary bases \( b = -1 \) and \( b = 0.04 \).
So far, this is what I've got. Here, I'll use \( RS_{b,L}(z) \) to denote the regular superfunction for base b at the fixed point L, \( R\Psi_{b,L}(z) \) to denote the regular Schroder function, and \( RA_{b,L}(z) \) to denote the regular Abel function. The lower subscripts, then, denote the index of the successive approximations to the given object which is subscripted.
Computation of super-, Schroder, and Abel function
The program uses the following formulas to compute the regular super-, Schroder, and Abel functions of the exponential.
Regular superfunction:
\( RS_{b,L}(z) = \sum_{n=0}^{\infty} r_n \left(L \log(b)\right)^{nz} \)
where the coefficients are recursively generated by
\( r_0 = L \), \( r_1 = 1 \)
\( r_n = \frac{B_n(1! \log(b) r_1, 2! \log(b) r_2, \cdots, (n-1)! \log(b) r_{n-1}, 0)}{n! \left(L^{n-1} \log(b) - \log(b)\right)} \).
(see Bell polynomials.)
Accuracy is best for \( z \) with large negative real part. Then we can use the recurrence
\( RS_{b,L}(z+1) = b^{RS_{b,L}(z)} \)
to extend to other values of \( z \).
Regular Schroder function:
This is generated by taking the inverse of the above when the coefficients are used as a power series, via the "reversion of series" formulae. Pari/GP already implements this so there is no need to do so explicitly in the program. The identity
\( R\Psi_{b,L}(b^z) = L \log(b) R\Psi_{b,L}(z) \)
then is used to extend to the upper halfplane.
Regular Abel function:
\( RA_{b,L}(z) = \log_{L \log(b)}\left(R\Psi_{b,L}(z)\right) \).
The algorithm
From what I could tell, the Kneser algorithm works as follows. Let \( F_n(z) \) be an approximation Taylor series to the tetrational, and \( \theta_n(z) \) be an approximation to the Kneser mapping. sheldonison, did I get this one right?
Here, I will denote by \( {f_m}_n \) the mth Taylor coefficient of \( F_n \), and \( {t_m}_n \) the mth Fourier coefficient of \( \theta_n(z) \).
1. Set up initial coefficients.
We set \( {f_0}_0 = 1 \), \( {f_1}_0 = 1 \), and all the rest to zero, giving \( F_0(z) = 1 + z \). We also set the Fourier coefficients all to zero.
2. Compute Fourier integral.
Now, we compute the Fourier integral to update the theta coefficients.
\( {t_m}_{n+1} = \int_{-1}^{0} \left(RA_{b,L}\left(F_n(x)\right) - x\right) e^{-2\pi i m x} dx \).
3. Compute Cauchy integral.
We now form the Taylor coefficients for the next tetrational approximation via the Cauchy integral.
\( {f_m}_{n+1} = \frac{1}{2\pi i} \left(\int_{\gamma_1} \frac{RS_{b,L}\left(z + \theta_{n+1}(z)\right)}{z^{m+1}} dz + \int_{\gamma_2} \frac{\bar{RS_{b,L}\left(\bar{z} + \theta_{n+1}(\bar{z})\right)}}{z^{m+1}} dz\right) \).
where the curves \( \gamma_1 \) and \( \gamma_2 \) are the upper and lower halves, respectively, of the unit circle, moving counterclockwise, starting at the positive real axis.
4. Normalize Taylor series.
Divide the Taylor series by its leading coefficient, so it acquires the value 1 at \( z = 0 \).
5. Repeat.
Repeat steps 2-4 until maximum precision is reached.
Problems: What's wrong?
The problem is that I can't seem to get this to go. For base \( e \), with 32 terms on the super/etc. series and 100 terms on the theta mapping, the residual (i.e. \( \left|F_n(0.5) - b^{F_n(-0.5)}\right| \), and is used to estimate the precision of the approximation) never drops below 0.00005. After 12 iterations, the residual flattens out at around 0.00005086. Using 64 terms on the super/etc. series and 800 terms on the theta mapping still only yields a residual of 0.00001258. I want a lot more precision than that!
I suspect the problem has to do with the fickle nature of the theta mapping, but I'm not sure how to overcome this.
Code
Here is the program. I hope that this is clear enough to be readable.
To use:
Initialize() to initialize the system
then
DoKneserIteration(12) (or however many loops you want)
then
FullTetApprox(x) to compute \( ^x b \).
I've been trying to build my own program in Pari/GP to run the Kneser algorithm for tetration, but I can't seem to get it to work so well. The goal is to ultimately be able to generalize it to merge complex bases and bases less than \( \eta \), especially the legendary bases \( b = -1 \) and \( b = 0.04 \).
So far, this is what I've got. Here, I'll use \( RS_{b,L}(z) \) to denote the regular superfunction for base b at the fixed point L, \( R\Psi_{b,L}(z) \) to denote the regular Schroder function, and \( RA_{b,L}(z) \) to denote the regular Abel function. The lower subscripts, then, denote the index of the successive approximations to the given object which is subscripted.
Computation of super-, Schroder, and Abel function
The program uses the following formulas to compute the regular super-, Schroder, and Abel functions of the exponential.
Regular superfunction:
\( RS_{b,L}(z) = \sum_{n=0}^{\infty} r_n \left(L \log(b)\right)^{nz} \)
where the coefficients are recursively generated by
\( r_0 = L \), \( r_1 = 1 \)
\( r_n = \frac{B_n(1! \log(b) r_1, 2! \log(b) r_2, \cdots, (n-1)! \log(b) r_{n-1}, 0)}{n! \left(L^{n-1} \log(b) - \log(b)\right)} \).
(see Bell polynomials.)
Accuracy is best for \( z \) with large negative real part. Then we can use the recurrence
\( RS_{b,L}(z+1) = b^{RS_{b,L}(z)} \)
to extend to other values of \( z \).
Regular Schroder function:
This is generated by taking the inverse of the above when the coefficients are used as a power series, via the "reversion of series" formulae. Pari/GP already implements this so there is no need to do so explicitly in the program. The identity
\( R\Psi_{b,L}(b^z) = L \log(b) R\Psi_{b,L}(z) \)
then is used to extend to the upper halfplane.
Regular Abel function:
\( RA_{b,L}(z) = \log_{L \log(b)}\left(R\Psi_{b,L}(z)\right) \).
The algorithm
From what I could tell, the Kneser algorithm works as follows. Let \( F_n(z) \) be an approximation Taylor series to the tetrational, and \( \theta_n(z) \) be an approximation to the Kneser mapping. sheldonison, did I get this one right?
Here, I will denote by \( {f_m}_n \) the mth Taylor coefficient of \( F_n \), and \( {t_m}_n \) the mth Fourier coefficient of \( \theta_n(z) \).
1. Set up initial coefficients.
We set \( {f_0}_0 = 1 \), \( {f_1}_0 = 1 \), and all the rest to zero, giving \( F_0(z) = 1 + z \). We also set the Fourier coefficients all to zero.
2. Compute Fourier integral.
Now, we compute the Fourier integral to update the theta coefficients.
\( {t_m}_{n+1} = \int_{-1}^{0} \left(RA_{b,L}\left(F_n(x)\right) - x\right) e^{-2\pi i m x} dx \).
3. Compute Cauchy integral.
We now form the Taylor coefficients for the next tetrational approximation via the Cauchy integral.
\( {f_m}_{n+1} = \frac{1}{2\pi i} \left(\int_{\gamma_1} \frac{RS_{b,L}\left(z + \theta_{n+1}(z)\right)}{z^{m+1}} dz + \int_{\gamma_2} \frac{\bar{RS_{b,L}\left(\bar{z} + \theta_{n+1}(\bar{z})\right)}}{z^{m+1}} dz\right) \).
where the curves \( \gamma_1 \) and \( \gamma_2 \) are the upper and lower halves, respectively, of the unit circle, moving counterclockwise, starting at the positive real axis.
4. Normalize Taylor series.
Divide the Taylor series by its leading coefficient, so it acquires the value 1 at \( z = 0 \).
5. Repeat.
Repeat steps 2-4 until maximum precision is reached.
Problems: What's wrong?
The problem is that I can't seem to get this to go. For base \( e \), with 32 terms on the super/etc. series and 100 terms on the theta mapping, the residual (i.e. \( \left|F_n(0.5) - b^{F_n(-0.5)}\right| \), and is used to estimate the precision of the approximation) never drops below 0.00005. After 12 iterations, the residual flattens out at around 0.00005086. Using 64 terms on the super/etc. series and 800 terms on the theta mapping still only yields a residual of 0.00001258. I want a lot more precision than that!
I suspect the problem has to do with the fickle nature of the theta mapping, but I'm not sure how to overcome this.Code
Here is the program. I hope that this is clear enough to be readable.
Code:
/* Compute the sum of a Fourier series:
* f(z) = sum_{n=0...inf} a_n e^(unz)
*/
FourierSum(coef, u, z) = sum(n=0,matsize(coef)[2]-1, coef[n+1]*exp(u*n*z));
/* Compute the sum of a Taylor series centered at u:
* f(z) = sum_{n=0...inf} a_n (z - u)^n
*/
TaylorSum(coef, u, z) = sum(n=0,matsize(coef)[2]-1, coef[n+1]*(z - u)^n);
/* Compute Bell polynomial */
bellcomplete(coefs) = {
local(n=matsize(coefs)[2]);
local(M=matrix(n,n));
for(i=1,n,\
for(j=1,n,\
M[i, j] = 0;
if(j-i+1 > 0, M[i, j] = binomial(n-i, j-i)*coefs[j-i+1]);
if(j-i+1 == 0, M[i, j] = -1);
);
);
return(matdet(M));
}
/* Compute the coefficients of the regular superfunction of the exponential
* with and at fixed point L and a given base (L must be a fixed point of the base).
*/
RegularCoeffs(L, base, N) = {
local(a = vector(N));
local(bellcoeffs);
local(logbase = log(base));
a[0+1] = L;
a[1+1] = 1;
/* Recurrence formula for Fourier coefficients:
* a_n = B_n(1! log(b) a_1, ..., (n-1)! log(b) a_(n-1), 0)/(n! (L^(n-1) log(b)^n - log(b)))
*/
for(n=2,N-1,
bellcoeffs = vector(n, k, if(k==n, 0, logbase*k!*a[k+1]));
a[n+1] = bellcomplete(bellcoeffs)/(n! * (L^(n-1) * logbase^n - logbase));
);
return(a);
}
/* Compute the coefficients of the regular Schroder function of the exponential
* with and at fixed point L.
*/
MakeSchroderCoeffs(RegCoeffs) = {
local(invpoly = 0);
local(RegCoeffsTrunc = RegCoeffs);
local(InvCoeffs = RegCoeffs);
RegCoeffsTrunc[1] = 0; /* truncate constant term */
/* Do series reversion. This gives the inverse of (InvSchroder(x) - L) */
invpoly = serreverse(TaylorSum(RegCoeffsTrunc, 0, x) + O(x^matsize(RegCoeffsTrunc)[2]));
/* Extract coefficients */
for(i=0,matsize(RegCoeffsTrunc)[2]-1,
InvCoeffs[i+1] = polcoeff(invpoly, i);
);
/* Done! */
return(InvCoeffs);
}
/* Compute the regular superfunction at a given fixed point L of a base b. */
RegularSuperfunction(L, base, SuperCoeffs, z) = {
if(real(z) < -5,
return(FourierSum(SuperCoeffs, log(L*log(base)), z));,
return(base^RegularSuperfunction(L, base, SuperCoeffs, z-1));
);
}
/* Compute the regular Schroder function at a given fixed point L of a base b. */
RegularSchroderFunction(L, base, SchroderCoeffs, z) = {
if(abs(z - L) > 0.5,
return(L*log(base)*RegularSchroderFunction(L, base, SchroderCoeffs, log(z)/log(base)));,
return(TaylorSum(SchroderCoeffs, L, z));
);
}
/* Compute the regular Abel function at a given fixd point L. */
RegularAbelFunction(L, base, SchroderCoeffs, z) = \
log(RegularSchroderFunction(L, base, SchroderCoeffs, z))/log(L*log(base));
/* --- Quadrature routines --- */
/* Do quadrature using function samples in samp with step delta. */
/* This is an inferior, Simpson's rule quadrature, implemented for testing only. In the future, we will
* drop quadrature and upgrade to FFT for high performance. Note that samp must have an odd number of
* elements.
*/
QuadratureCore(samp, delta) = {
local(ssum=0);
ssum = samp[1] + sum(n=1,matsize(samp)[2]-2,(4 + ((-1)^(n-1) - 1))*samp[n+1]) + samp[matsize(samp)[2]];
return((delta/3)*ssum);
}
/* Compute some regularly-spaced nodes over an interval to use for quadrature. N must be odd. */
QuadratureNodes(a, b, N) = {
local(nodes);
nodes = vector(N, i, a + (i-1)*((b-a)/(N-1)));
return(nodes);
}
/* Compute the delta value for a given interval and node count */
QuadratureDelta(a, b, N) = (b - a)/(N-1);
/* --- End quadrature routines --- */
/* Set up the tetration system. */
Initialize() = {
/* Initialize the fixed points and quadrature. */
/* Base and fixed points */
base = exp(1);
upperfix = 0.318 + 1.337*I;
for(i=1,16,upperfix = upperfix - (base^upperfix - upperfix)/(log(base)*base^upperfix - 1));
/* Quadrature setup */
NumQuadNodes = 1001;
ImagOffset = 0.0000001*I; /* offset off the real axis to do Fourier integration */
ThetaQuadratureNodesUpper = QuadratureNodes(-1.0 + ImagOffset, 0.0 + ImagOffset, NumQuadNodes);
ThetaQuadratureDeltaUpper = QuadratureDelta(-1.0 + ImagOffset, 0.0 + ImagOffset, NumQuadNodes);
CauchyCircleQuadratureNodesUpper = QuadratureNodes(0, Pi, NumQuadNodes);
CauchyCircleQuadratureDeltaUpper = QuadratureDelta(0, Pi, NumQuadNodes);
CauchyCircleQuadratureNodesLower = QuadratureNodes(Pi, 2*Pi, NumQuadNodes);
CauchyCircleQuadratureDeltaLower = QuadratureDelta(Pi, 2*Pi, NumQuadNodes);
/* Term counts */
NumRegTerms = 32;
NumTaylorTerms = 32;
NumThetaTerms = 100;
/* Generate series */
RegCoeffsUpper = RegularCoeffs(upperfix, base, NumRegTerms);
RegSchrCoeffsUpper = MakeSchroderCoeffs(RegCoeffsUpper);
/* Set up initial Taylor series */
TaylorCoeffs = vector(NumTaylorTerms);
TaylorCoeffs[0+1] = 1; /* coefficient of x^0 = 1 */
TaylorCoeffs[1+1] = 1; /* coefficient of x^1 = 1 */
/* Set up initial Riemann mapping series */
RiemannCoeffsUpper = vector(NumThetaTerms);
}
/* Compute upper regular superfunction */
RegularSuperUpper(z) = RegularSuperfunction(upperfix, base, RegCoeffsUpper, z);
/* Compute upper regular Abel function */
RegularAbelUpper(z) = RegularAbelFunction(upperfix, base, RegSchrCoeffsUpper, z);
/* Construct Fourier coefficients for a Riemann theta(z) mapping from a tetration Taylor series. */
RiemannFromTaylorStep() = {
local(Samples);
local(FourierSamp);
/* Sample the function at the node points. */
Samples = vector(NumQuadNodes, i, \
RegularAbelUpper(TaylorSum(TaylorCoeffs, 0, ThetaQuadratureNodesUpper[i])) -\
ThetaQuadratureNodesUpper[i]);
/*
plot(X=1,NumQuadNodes,real(Samples[floor(X)]));
plot(X=1,NumQuadNodes,imag(Samples[floor(X)]));
*/
/* Do the Fourier integrals */
/* We multiply by exp(-2 pi i n ImagOffset) to compensate for the offset along the imaginary axis. */
for(n=0,NumThetaTerms-1,
FourierSamp = vector(NumQuadNodes, i, Samples[i]*exp(-2*Pi*I*n*real(ThetaQuadratureNodesUpper[i])));
RiemannCoeffsUpper[n+1] = QuadratureCore(FourierSamp, ThetaQuadratureDeltaUpper);
RiemannCoeffsUpper[n+1] = RiemannCoeffsUpper[n+1] * exp(-2*Pi*I*n*ImagOffset);
);
}
/* Calculate the theta mapping. */
ThetaMappingUpper(z) = FourierSum(RiemannCoeffsUpper, 2*Pi*I, z);
/* Calculate the upper "warped" superfunction. */
WarpSuperUpper(z) = RegularSuperUpper(z + ThetaMappingUpper(z));
/* Calculate the lower "warped" superfunction. */
WarpSuperLower(z) = conj(WarpSuperUpper(conj(z))); /* fine for real bases > eta only */
/* Do the Cauchy integral to get the Taylor series. */
CauchyTaylorStep() = {
local(Samples);
local(UpperSuperSamp);
local(LowerSuperSamp);
local(UpperInt);
local(LowerInt);
/* Presampling */
UpperSuperSamp = vector(NumQuadNodes, i, WarpSuperUpper(exp(I*CauchyCircleQuadratureNodesUpper[i])));
LowerSuperSamp = vector(NumQuadNodes, i, WarpSuperLower(exp(I*CauchyCircleQuadratureNodesLower[i])));
/* Do Cauchy integral */
for(n=0,NumTaylorTerms-1,\
/* Sample the function around the upper half-circle. */\
Samples = vector(NumQuadNodes, i, UpperSuperSamp[i]/\
(exp(I*CauchyCircleQuadratureNodesUpper[i])^(n+1)) *\
I*exp(I*CauchyCircleQuadratureNodesUpper[i]));\
UpperInt = QuadratureCore(Samples, CauchyCircleQuadratureDeltaUpper);\
/* Do the same around the lower half-circle. */\
Samples = vector(NumQuadNodes, i, LowerSuperSamp[i]/\
(exp(I*CauchyCircleQuadratureNodesLower[i])^(n+1)) *\
I*exp(I*CauchyCircleQuadratureNodesLower[i]));\
LowerInt = QuadratureCore(Samples, CauchyCircleQuadratureDeltaLower);\
/* Get the result */\
TaylorCoeffs[n+1] = 1/(2*Pi*I) * (UpperInt + LowerInt);\
);
/* Renormalize Taylor series */
TaylorCoeffs = TaylorCoeffs / TaylorSum(TaylorCoeffs, 0, 0);
}
/* Do full iterations */
DoKneserIteration(numiters) = {
for(i=1,numiters,
print("Iteration ", i, ":");
print("Building Kneser mapping...");
RiemannFromTaylorStep();
print("Performing Cauchy integral...");
CauchyTaylorStep();
print("Done. R = ", abs(base^TaylorSum(TaylorCoeffs, 0, -0.5) - TaylorSum(TaylorCoeffs, 0, 0.5)));
);
}
/* Calculate full superfunction */
FullTetApprox(z) = {
if(real(z) > 0.5, return(base^FullTetApprox(z-1)));
if(real(z) < -0.5, return(log(FullTetApprox(z+1))/log(base)));
if(imag(z) > 1, return(WarpSuperUpper(z)));
if(imag(z) < -1, return(WarpSuperLower(z)));
return(TaylorSum(TaylorCoeffs, 0, z));
}To use:
Initialize() to initialize the system
then
DoKneserIteration(12) (or however many loops you want)
then
FullTetApprox(x) to compute \( ^x b \).
