06/16/2011, 11:48 AM
Well, I managed to break through that convergence barrier now, by using the Taylor approximation at the real axis.
Here's the new code:
Hope to get this going for complex bases real soon
Here's the new code:
Code:
/* 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));
);
}
/* 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.12*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);
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);
/* 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]);
/* 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(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));
}Hope to get this going for complex bases real soon
