• 0 Vote(s) - 0 Average
• 1
• 2
• 3
• 4
• 5
 "Kneser"/Riemann mapping method code for *complex* bases mike3 Long Time Fellow Posts: 368 Threads: 44 Joined: Sep 2009 08/15/2011, 03:44 AM (This post was last modified: 08/15/2011, 12:03 PM by mike3.) 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$: 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()". 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*I``` 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. sheldonison Long Time Fellow Posts: 626 Threads: 22 Joined: Oct 2008 08/15/2011, 02:59 PM (08/15/2011, 03:44 AM)mike3 Wrote: 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$: .... 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. Mike, Very nice! For base=3+i, I can calculate two repelling fixed points, each with different periods. Presumably, the upper half of the complex plane converges towards the first fixed point, and the lower half of the complex plane converges towards the second fixed point, while the merged function allows you to maintaining the definition of sexp(0)=1, and sexp(1)=base. L1=0.324536386411256 + 1.00148180609593*I L2=-0.0976763924712228 - 1.39768129114470*I I assume the algorithm works when both bases are repelling? - Sheldon Gottfried Ultimate Fellow Posts: 757 Threads: 116 Joined: Aug 2007 08/15/2011, 03:14 PM Very nice, thank you! I tried it and the Pari/GP code worked immediately. I've nothing more to say at the moment, I think I'll play a bit with it to get a feel. When I find anything concerning your question I'll replay again, may be this will take some days. Gottfried Gottfried Helms, Kassel « Next Oldest | Next Newest »

 Possibly Related Threads... Thread Author Replies Views Last Post Natural complex tetration program + video MorgothV8 1 1,313 04/27/2018, 07:54 PM Last Post: MorgothV8 complex base tetration program sheldonison 23 35,716 10/26/2016, 10:02 AM Last Post: Gottfried fast accurate Kneser sexp algorithm sheldonison 38 67,191 01/14/2016, 05:05 AM Last Post: sheldonison C++ program for generatin complex map in EPS format MorgothV8 0 2,330 09/17/2014, 04:14 PM Last Post: MorgothV8 C++ code for tet, ate and hexp MorgothV8 0 2,711 07/10/2014, 04:24 PM Last Post: MorgothV8 Green Eggs and HAM: Tetration for ALL bases, real and complex, now possible? mike3 27 29,595 07/02/2014, 10:13 PM Last Post: tommy1729 Which method is currently "the best"? MorgothV8 2 4,165 11/15/2013, 03:42 PM Last Post: MorgothV8 Attempt to make own implementation of "Kneser" algorithm: trouble mike3 9 14,226 06/16/2011, 11:48 AM Last Post: mike3 An incremental method to compute (Abel) matrix inverses bo198214 3 8,556 07/20/2010, 12:13 PM Last Post: Gottfried Single-exp series computation code mike3 0 2,737 04/20/2010, 08:59 PM Last Post: mike3

Users browsing this thread: 1 Guest(s)