• 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. « Next Oldest | Next Newest »

 Messages In This Thread "Kneser"/Riemann mapping method code for *complex* bases - by mike3 - 08/15/2011, 03:44 AM RE: "Kneser"/Riemann mapping method code for *complex* bases - by sheldonison - 08/15/2011, 02:59 PM RE: "Kneser"/Riemann mapping method code for *complex* bases - by Gottfried - 08/15/2011, 03:14 PM

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

Users browsing this thread: 1 Guest(s)