• 0 Vote(s) - 0 Average
• 1
• 2
• 3
• 4
• 5
 Single-exp series computation code mike3 Long Time Fellow Posts: 368 Threads: 44 Joined: Sep 2009 04/20/2010, 08:59 PM tommy1729 wanted the details of the method using the single-exp series. I'd post it anyway, but here you go. Periodic interpolation The idea is this. If we denote by $x_n$, a finite sequence of $N$ points ($n = 0 ... N-1$) to be interpolated, then we can find a periodic interpolation by taking their discrete Fourier transform (DFT): $X_k = \sum_{n=0}^{N-1} x_n e^{-2\pi i n k/N}$ We can find the points interpolated by taking the inverse discrete Fourier transform: $x_n = \frac{1}{N} \sum_{k=0}^{N-1} x_k e^{2\pi i n k/N}$ If $N$ is odd, then the interpolation is constructed by $f(z) =\ \sum_{k=-(N-1)/2}^{(N-1)/2} \frac{1}{N} X_k e^{2\pi i k z/N}$. where $X_k$ for $k$ outside $k = 0...N-1$ is obtained by periodic repetition, i.e. $X_k = X_{k\ \mathrm{mod}\ N}$. The values of $f(z)$ will equals $x_k$ when $x = k$ for integers $k$. If we want period $P$, $f(z) =\ \sum_{k=-(N-1)/2}^{(N-1)/2} \frac{1}{N} X_k e^{\frac{2\pi}{P} i k z/N}$. In any case, $f(0)$ will be $x_0$. If we want our nodes $x_k$ to represent values on either side of 0, with the middle node being 0, then we just rotate them by the appropriate factor ($\frac{N-1}{2}$) before interpolating. The discrete Fourier transform can be done much faster using the Fast Fourier Transform (FFT) algorithm. The current program uses only the "brute force" method, since it was just an initial test to see if this route was any good. Periodizing functions The next thing we need is a "periodizing function". The tetration is not periodic, so we must approximate it with a periodic function of long period. A "periodizing function" is a function $P_r(z)$ such that 1. $P_r(z)$ has period $r$. 2. $\lim_{r \rightarrow \infty} P_r(z) = z$. The idea is that given some non-periodic function $f(z)$, $f(P_r(z))$ is a periodic function approximately equal to $f$ for $z$ near 0. Increasing the period $r$ makes the approximation better. For example, $P_r(z) = \frac{r}{2\pi} \sin \left(\frac{2\pi}{r} z\right)$ is one possibility. For the program, we use the slightly better $P_r(z) = \frac{r}{2\pi} 1.25 \sin\left(\frac{2\pi}{r} z\right) - \frac{r}{2\pi} 0.125 \sin\left(\frac{2\pi}{r} 2z\right)$. The algorithm We can put these two together to get the following algorithm: 1. Create a Fourier series array with some odd number of nodes and an initial guess. 2. Compute the continuum sum of the array. Store the linear term separately from the new exp/Fourier series. 3. Compute the exponential to the base $b$ by converting to function space (via IDFT), exponentiating each node, and converting back to Fourier space (via DFT). This may not be the most precise way, and is not very precise for very small numbers of nodes. However it can be done using the FFT if that is available. (Question: Are there any better methods, perhaps with more precision?). Exponentiate the linear term separately (implicit, just interpret its separetly-value $a_0$ from before not as $a_0 z$ but $b^{a_0 z}$ and as a multiplier instead of an addend). 4. Multiply by $\log(b)^z$. This is just "implicit" (just treat the array after this point as if it has a $\log(b)^z$ multiplier). 5. Integrate with lower bound -1. The integral of the Fourier terms $a_n e^{(a_0 \log(b) + \log(\log(b)) + nu)z}$ is $\frac{a_n}{a_0 \log(b) + \log(\log(b)) + nu} \left(e^{(a_0 \log(b) + \log(\log(b)) + nu)z} - e^{-(a_0 \log(b) + \log(\log(b)) + nu)}\right)$. 6. Periodize the result by substituting the periodizing function into the newly-constructed interpolation, evaluating the result at the interpolation nodes' absicssas, and converting it back to Fourier space. (This can also be done using the FFT algorithm, but we should use the DIT/Cooley-Tukey method.) 7. Normalize. Divide the result by its evaluation at 0. 8. Repeat steps 2-7 until convergence is attained. The program The implementation of the algorithm is given below, in PARI/GP. Code:/* Perform discrete Fourier transform (DFT). */ /* We should rewrite this to use FFT */ DFT(v) = {         local(N = matsize(v)[2]);         local(rv = v);         for(k=0,N-1,             rv[k+1] = sum(n=0,N-1,v[n+1]*exp(-2*Pi*I*n*k/N));            );         return(rv); } /* Perform inverse discrete Fourier transform (IDFT). */ /* We should rewrite this to use FFT */ IDFT(v) = {         local(N = matsize(v)[2]);         local(rv = v);         for(k=0,N-1,             rv[k+1] = (1/N)*sum(n=0,N-1,v[n+1]*exp(2*Pi*I*n*k/N));            );         return(rv); } /* Convert a node array from real space to Fourier space. Odd sized array only. */ ConvToFourier(v) = {         local(N = matsize(v)[2]);         local(vsh = vector(N));         local(tmp = vector(N));         local(rv = vector(N));         /*print(v);*/         for(n=0,N-1,vsh[n+1] = v[(n+((N-1)/2))%N + 1]);         /*print(vsh);*/         tmp = (1/N)*DFT(vsh);         /*print(tmp);*/         for(n=-(N-1)/2,(N-1)/2,rv[n+((N-1)/2)+1] = tmp[n%N+1]);         /*print(rv);*/         return(rv); } /* Convert a node array from Fourier space to real space. Odd sized array only. */ ConvToReal(v) = {         local(N = matsize(v)[2]);         local(vsh = vector(N));         local(tmp = vector(N));         local(rv = vector(N));         for(n=0,N-1,vsh[n+1] = v[(n+((N-1)/2))%N+1]);         tmp = IDFT(N*vsh);         for(n=0,N-1,rv[n+1] = tmp[(n-((N-1)/2))%N+1]);         return(rv); } /* Use a node array in Fourier space to do the interpolation. Odd sized array only. */ DoInterpolation(v, x) = {         local(N = matsize(v)[2]);         return(sum(n=-((N-1)/2),(N-1)/2,v[n+((N-1)/2)+1]*exp(2*Pi*I*n*x/N))); } /* Same as above, but with specified period */ DoInterpolationPer(P, v, x) = {         local(N = matsize(v)[2]);         return(sum(n=-((N-1)/2),(N-1)/2,v[n+((N-1)/2)+1]*exp(2*Pi*I*n*x/P))); } /* Get the list of abscissas for the nodes with a given period */ GetAbscissas(N, P) = {         return(vector(N, n, (-((N-1)/2) + (n-1)) * (P/N))); } /* Make the continuum sum of a Fourier-space array with a given period. */ MakeContSum(P, v) = {         local(N = matsize(v)[2]);         local(rv = [0, 0]);         rv[2] = vector(N);         for(n=-((N-1)/2),(N-1)/2,             if(n != 0,                rv[2][n+((N-1)/2)+1] = v[n+((N-1)/2)+1]/(exp(2*Pi*I*n/P) - 1);,                rv[2][n+((N-1)/2)+1] =                  -sum(k=-(N-1)/2,(N-1)/2,if(k==0,0,v[k+((N-1)/2)+1]/(exp(2*Pi*I*k/P) - 1)));               );            );         rv[1] = v[((N-1)/2)+1];         return(rv); } /* Evaluate continuum sum array */ EvalContSumArr(P, v, x) = v[1]*x + DoInterpolationPer(P, v[2], x); /* Perform exponential of a continuum sum array with a given period. */ ExpContSum(P, base, v) = [base, v[1], ConvToFourier(base^ConvToReal(v[2]))]; /* Evaluate continuum sum exp array */ EvalContSumExpArr(P, v, x) = v[1]^(v[2]*x) * DoInterpolationPer(P, v[3], x); /* The periodizing function */ Periodizer(P, x) = P/(2*Pi) * 1.25 * sin((2*Pi)/P * x) - P/(2*Pi) * 0.125 * sin((2*Pi)/P * 2 * x); /* Multiply continuum-sum exp array by log(base)^x and integrate from -1 to x. */ IntegrationPhase(P, v) = {         local(N = matsize(v[3])[2]);         local(vaux = v[3]);         local(q = v[2]);         local(base = v[1]);         local(abscissas = GetAbscissas(N, P));         local(rv = vector(N));         local(expf = vector(N));         for(n=-(N-1)/2,(N-1)/2,expf[n+((N-1)/2)+1] = (2*Pi*I*n/P) + (q*log(base)) + log(log(base)););         for(n=-(N-1)/2,(N-1)/2,             vaux[n+((N-1)/2)+1] = vaux[n+((N-1)/2)+1]/expf[n+((N-1)/2)+1];            );         for(k=1,N,             rv[k] = sum(n=-(N-1)/2,(N-1)/2,                         vaux[n+((N-1)/2)+1]*(exp(expf[n+((N-1)/2)+1]*Periodizer(P, abscissas[k])) -                                              exp(-expf[n+((N-1)/2)+1]))                        );            );         return(rv); } /* Debug: dump log(base)^x times v function */ PreintDbg(P, v, x) = {         local(N = matsize(v[3])[2]);         local(vaux = v[3]);         local(q = v[2]);         local(base = v[1]);         return(sum(n=-(N-1)/2,(N-1)/2,                         vaux[n+((N-1)/2)+1]*exp((2*Pi*I*n/P + q*log(base) + log(log(base)))*x))); } /* Do tetration iteration over Fourier array */ NormTetIter(P, base, v) = {         local(rv);         rv = ConvToFourier(IntegrationPhase(P, ExpContSum(P, base, MakeContSum(P, v))));         rv = rv / DoInterpolationPer(P, rv, 0);         return(rv); } /* Tetration routines */ TetCrit(P, v, X) = DoInterpolationPer(P, v, X); TetApprox(P, v, base, X) = {         if(real(X) < -0.5, return(log(TetApprox(P, v, base, X+1))/log(base)));         if(real(X) > 0.5, return(base^TetApprox(P, v, base, X-1)));         return(TetCrit(P, v, X)); } To use, create an initial guess with Code:tv = ConvToFourier(vector(N, i, 1)) where N should be odd, say N = 129. Decide on some period. I use pure imaginary periods. Assign it to P, say P = 48*I. Then just iterate away with Code:tv = NormTetIter(P, base, tv) where base is the desired base, like 2, or exp(1). When you want to use the approximation to evaluate tetration, use TetApprox(P, tv, base, x) with x being the height. For complex bases, a flat initial guess does not seem to yield convergence. However, starting with a real base's tetration, and then using that as the initial guess for a base somewhat away from it in the imaginary direction, and then using that as the guess for another, and so on, seems to work better. I haven't yet been able to push this into the left halfplane, not sure why. « Next Oldest | Next Newest »

 Possibly Related Threads... Thread Author Replies Views Last Post C++ code for tet, ate and hexp MorgothV8 0 4,794 07/10/2014, 04:24 PM Last Post: MorgothV8 "Kneser"/Riemann mapping method code for *complex* bases mike3 2 10,250 08/15/2011, 03:14 PM Last Post: Gottfried A note on computation of the slog Gottfried 6 15,710 07/12/2010, 10:24 AM Last Post: Gottfried Computations with the double-exp series mike3 0 3,879 04/20/2010, 07:32 PM Last Post: mike3 SAGE code for computing flow matrix for exp(z)-1 jaydfox 4 13,635 08/21/2009, 05:32 PM Last Post: jaydfox A nice series for b^^h , base sqrt(2), by diagonalization Gottfried 19 35,722 06/11/2009, 08:36 PM Last Post: Gottfried An error estimate for fixed point computation of b^x bo198214 0 4,188 05/31/2008, 04:11 PM Last Post: bo198214 SAGE code implementing slog with acceleration jaydfox 4 11,182 10/22/2007, 12:59 AM Last Post: jaydfox Complete SAGE code for tetration eyu100 4 10,043 10/18/2007, 03:55 AM Last Post: jaydfox

Users browsing this thread: 1 Guest(s)