Hi.
It's been a while but I figured I could post my attempt to tetrate complex bases with the "Kneser"/Riemann mapping algorithm.
This is the code I've got, which will tetrate a complex base, here base \( 3 + i \):
To use, run "Initialize()", then "DoKneserIteration(<some number of iterations>)". The program will generate the theta mappings to fuse the two regular iterations at the two principal fixed points. With the given Taylor/etc. term counts, 16 iterations yields a residual on the order of \( 10^{-15} \), representing what may be some of the most accurate calculations of a tetrational of a complex base to date. We get
\( ^{1/2} (3 + i) \approx 1.73774339741062 + 0.19778352516442i \)
Graphs of \( ^x (3 + i) \) are shown below.
On the complex plane (scale -10 to +10 on both axes):
As we can see, the merger worked successfully. The two different regular iterations have flowed together nice and holomorphically for \( \Re(z) > -2 \).
We can even do the "pentation" (the \( \uparrow \uparrow \uparrow \) operator) of this base now. The sequence of integer pentates is
Note how the pentation looks to converge to a limit, like tetration for STR bases.
It does not work for all complex bases, though. It seems to work only for those relatively near the real axis. Increasing the ImagOffset (how far off the real axis we sample the Fourier integrals) appears to extend the range, however. When it fails, it looks like the iteration used to analytically continue the Schroder function gets sucked into the wrong fixed point. Any way to resolve this? The aforementioned increasing of the ImagOffset causes losses in accuracy and efficiency, and is not a cure-all that allows for tetrating all complex bases outside the STR with the possible exception of the possibly singular base 0.
It's been a while but I figured I could post my attempt to tetrate complex bases with the "Kneser"/Riemann mapping algorithm.
This is the code I've got, which will tetrate a complex base, here base \( 3 + i \):
Code:
/* Independent implementation of Kneser mapping algorithm. */
/* Compute the sum of a Fourier series:
* f(z) = sum_{n=0...inf} a_n e^(uinz)
*/
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));
);
}
rsf(L, base, SchroderCoeffs, z) = {
if(abs(imag(z)) < 0.02, return(3));
return(RegularSchroderFunction(L, base, SchroderCoeffs, 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 = 3 + I;
upperfix = 0.325 + 1.001*I;
lowerfix = 0.098 - 1.398*I;
for(i=1,16,upperfix = upperfix - (base^upperfix - upperfix)/(log(base)*base^upperfix - 1));
for(i=1,16,lowerfix = lowerfix - (base^lowerfix - lowerfix)/(log(base)*base^lowerfix - 1));
/* Quadrature setup */
NumQuadNodes = 1001;
ImagOffset = 0.12*I; /* offset off the real axis to do Fourier integration */
ThetaQuadratureNodesUpper = QuadratureNodes(-0.5 + ImagOffset, 0.5 + ImagOffset, NumQuadNodes);
ThetaQuadratureDeltaUpper = QuadratureDelta(-0.5 + ImagOffset, 0.5 + ImagOffset, NumQuadNodes);
ThetaQuadratureNodesLower = QuadratureNodes(-0.5 - ImagOffset, 0.5 - ImagOffset, NumQuadNodes);
ThetaQuadratureDeltaLower = QuadratureDelta(-0.5 - ImagOffset, 0.5 - ImagOffset, NumQuadNodes);
CauchyCircleQuadratureNodes = QuadratureNodes(0, 2*Pi, NumQuadNodes);
CauchyCircleQuadratureDelta = QuadratureDelta(0, 2*Pi, NumQuadNodes);
/* Term counts */
NumRegTerms = 32;
NumTaylorTerms = 32;
NumThetaTerms = 100;
/* Generate series */
RegCoeffsUpper = RegularCoeffs(upperfix, base, NumRegTerms);
RegSchrCoeffsUpper = MakeSchroderCoeffs(RegCoeffsUpper);
RegCoeffsLower = RegularCoeffs(lowerfix, base, NumRegTerms);
RegSchrCoeffsLower = MakeSchroderCoeffs(RegCoeffsLower);
/* 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);
RiemannCoeffsLower = 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);
/* Compute lower regular superfunction */
RegularSuperLower(z) = RegularSuperfunction(lowerfix, base, RegCoeffsLower, z);
/* Compute lower regular Abel function */
RegularAbelLower(z) = RegularAbelFunction(lowerfix, base, RegSchrCoeffsLower, z);
/* Construct Fourier coefficients for a Riemann theta(z) mapping from a tetration Taylor series. */
RiemannFromTaylorStep() = {
local(SamplesUpper);
local(SamplesLower);
local(FourierSamp);
/* Sample the function at the node points. */
SamplesUpper = vector(NumQuadNodes, i, \
RegularAbelUpper(TaylorSum(TaylorCoeffs, 0, ThetaQuadratureNodesUpper[i])) -\
ThetaQuadratureNodesUpper[i]);
SamplesLower = vector(NumQuadNodes, i, \
RegularAbelLower(TaylorSum(TaylorCoeffs, 0, ThetaQuadratureNodesLower[i])) -\
ThetaQuadratureNodesLower[i]);
/*
plot(X=1,NumQuadNodes,real(SamplesUpper[floor(X)]));
plot(X=1,NumQuadNodes,imag(SamplesUpper[floor(X)]));
plot(X=1,NumQuadNodes,real(SamplesLower[floor(X)]));
plot(X=1,NumQuadNodes,imag(SamplesLower[floor(X)]));
PublicSamplesUpper = SamplesUpper;
PublicSamplesLower = SamplesLower;
*/
/* 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, SamplesUpper[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);
FourierSamp = vector(NumQuadNodes, i, SamplesLower[i]*exp(2*Pi*I*n*real(ThetaQuadratureNodesLower[i])));
RiemannCoeffsLower[n+1] = QuadratureCore(FourierSamp, ThetaQuadratureDeltaLower);
RiemannCoeffsLower[n+1] = RiemannCoeffsLower[n+1] * exp(-2*Pi*I*n*ImagOffset);
);
}
/* Calculate the theta mapping. */
ThetaMappingUpper(z) = FourierSum(RiemannCoeffsUpper, 2*Pi*I, z);
ThetaMappingLower(z) = FourierSum(RiemannCoeffsLower, -2*Pi*I, z);
/* Calculate the upper "warped" superfunction. */
WarpSuperUpper(z) = RegularSuperUpper(z + ThetaMappingUpper(z));
/* Calculate the lower "warped" superfunction. */
WarpSuperLower(z) = RegularSuperLower(z + ThetaMappingLower(z));
/* Do the Cauchy integral to get the Taylor series. */
CauchyTaylorStep() = {
local(Samples);
local(CauchyizedSamples);
/* The integration proceeds in several phases:
* 1. Integrate above Im(z) = ImagOffset using the upper regular warped superfunction,
* 2. Integrate from Im(z) = ImagOffset to Im(z) = -ImagOffset using the old Taylor series,
* 3. Integrate below Im(z) = -ImagOffset using the lower regular warped superfunction,
* 4. Integrate from Im(z) = -ImagOffset to Im(z) = ImagOffset using the old Taylor series again.
* To maximize precision, we use the exp/log of the Taylor series evaluated near zero. The whole
* process is done counterclockwise.
*/
/* Presampling */
Samples = vector(NumQuadNodes);
CirclePositions = vector(NumQuadNodes, i, exp(I*CauchyCircleQuadratureNodes[i]));
for(n=0,NumQuadNodes-1,\
if(abs(imag(CirclePositions[n+1])) < imag(ImagOffset),\
if(real(CirclePositions[n+1]) > 0,\
/* right halfplane */\
Samples[n+1] = base^TaylorSum(TaylorCoeffs, 0, CirclePositions[n+1]-1);,\
/* left halfplane */\
Samples[n+1] = log(TaylorSum(TaylorCoeffs, 0, CirclePositions[n+1]+1))/log(base);\
);,\
if(imag(CirclePositions[n+1]) > 0,\
/* upper halfplane */\
/* use upper regular warped superfunction */\
Samples[n+1] = WarpSuperUpper(CirclePositions[n+1]);,\
/* lower halfplane */\
/* use lower regular warped superfunction */\
Samples[n+1] = WarpSuperLower(CirclePositions[n+1]);\
);\
);\
);
/* Do Cauchy integral */
for(n=0,NumTaylorTerms-1,\
/* Cauchy setup */\
CauchyizedSamples = vector(NumQuadNodes, i, Samples[i]/(CirclePositions[i]^(n+1)) * \
I*CirclePositions[i]);\
/* Integrate */\
TaylorCoeffs[n+1] = 1/(2*Pi*I) * QuadratureCore(CauchyizedSamples, CauchyCircleQuadratureDelta);\
);
/* 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, run "Initialize()", then "DoKneserIteration(<some number of iterations>)". The program will generate the theta mappings to fuse the two regular iterations at the two principal fixed points. With the given Taylor/etc. term counts, 16 iterations yields a residual on the order of \( 10^{-15} \), representing what may be some of the most accurate calculations of a tetrational of a complex base to date. We get
\( ^{1/2} (3 + i) \approx 1.73774339741062 + 0.19778352516442i \)
Graphs of \( ^x (3 + i) \) are shown below.
On the complex plane (scale -10 to +10 on both axes):
As we can see, the merger worked successfully. The two different regular iterations have flowed together nice and holomorphically for \( \Re(z) > -2 \).
We can even do the "pentation" (the \( \uparrow \uparrow \uparrow \) operator) of this base now. The sequence of integer pentates is
Code:
(3 + i) ^^^ 0 ~ 1
(3 + i) ^^^ 1 ~ 3.0000000000000 + 1.0000000000000*I
(3 + i) ^^^ 2 ~ 0.067505751821009 + 0.37038084586491*I
(3 + i) ^^^ 3 ~ 0.95290527513344 + 0.41965918718415*I
(3 + i) ^^^ 4 ~ 1.6689141108302 + 1.5511207891726*I
(3 + i) ^^^ 5 ~ 0.14427450000007 + 1.2422893053044*I
(3 + i) ^^^ 6 ~ 0.59855234346250 + 0.91546250434381*I
(3 + i) ^^^ 7 ~ 0.86698754296965 + 1.1149406202426*I
(3 + i) ^^^ 8 ~ 0.67898538651189 + 1.3074168077768*I
(3 + i) ^^^ 9 ~ 0.60455661809291 + 1.1563634545945*I
(3 + i) ^^^ 10 ~ 0.69330548523182 + 1.1273955819904*I
(3 + i) ^^^ 11 ~ 0.70596729212537 + 1.1855496560953*I
(3 + i) ^^^ 12 ~ 0.66761843466257 + 1.1866711686455*I
(3 + i) ^^^ 13 ~ 0.67126297874107 + 1.1635087706344*I
(3 + i) ^^^ 14 ~ 0.68499772461207 + 1.1678785414537*I
(3 + i) ^^^ 15 ~ 0.68092916998454 + 1.1759457551948*I
(3 + i) ^^^ 16 ~ 0.67640594691649 + 1.1725959600056*I
(3 + i) ^^^ 17 ~ 0.67891225654672 + 1.1701754208531*I
(3 + i) ^^^ 18 ~ 0.68014536910711 + 1.1719549038205*I
(3 + i) ^^^ 19 ~ 0.67892704487410 + 1.1725305088990*I
(3 + i) ^^^ 20 ~ 0.67869891860093 + 1.1717250339640*I
(3 + i) ^^^ 21 ~ 0.67921563866767 + 1.1716670353204*I
(3 + i) ^^^ 22 ~ 0.67919871696253 + 1.1719897863482*I
(3 + i) ^^^ 23 ~ 0.67900254832419 + 1.1719465366940*I
(3 + i) ^^^ 24 ~ 0.67904902631454 + 1.1718306916743*I
(3 + i) ^^^ 25 ~ 0.67911531745648 + 1.1718709649452*I
(3 + i) ^^^ 26 ~ 0.67908388197865 + 1.1719075100622*I
(3 + i) ^^^ 27 ~ 0.67906467678018 + 1.1718845193510*I
(3 + i) ^^^ 28 ~ 0.67908072504425 + 1.1718750848271*I
(3 + i) ^^^ 29 ~ 0.67908487628565 + 1.1718858831902*I
(3 + i) ^^^ 30 ~ 0.67907783385992 + 1.1718873294663*INote how the pentation looks to converge to a limit, like tetration for STR bases.
It does not work for all complex bases, though. It seems to work only for those relatively near the real axis. Increasing the ImagOffset (how far off the real axis we sample the Fourier integrals) appears to extend the range, however. When it fails, it looks like the iteration used to analytically continue the Schroder function gets sucked into the wrong fixed point. Any way to resolve this? The aforementioned increasing of the ImagOffset causes losses in accuracy and efficiency, and is not a cure-all that allows for tetrating all complex bases outside the STR with the possible exception of the possibly singular base 0.
