Thread Rating:
• 0 Vote(s) - 0 Average
• 1
• 2
• 3
• 4
• 5
 Numerical algorithm for Fourier continuum sum tetration theory mike3 Long Time Fellow Posts: 368 Threads: 44 Joined: Sep 2009 09/15/2010, 11:06 AM (This post was last modified: 09/15/2010, 11:21 AM by mike3.) Conclusion Putting this all together, we get the following PARI/GP program. It is probably not the clearest thing out there, but it was assembled more "ad hoc" as these ideas came to me, instead of starting with a "slick" plan like the one above. Code:```/* Radix-3 DIT FFT. */ /* Convert to reversed trinary. */ trinreverse(n, maxtrits) = {         local(nn = n);         local(rv = 0);         for(i=0,maxtrits-1,             rv = rv * 3;             rv = rv + (nn%3);             nn = floor(nn/3);            );         return(rv); } /* Trit reversal. Assumes v has size 3^n. */ DoTritReversal(v) = {         local(N = matsize(v)[2]);         local(Nlg3 = floor(log(N)/log(3) + 0.5));         local(reversed = vector(N));         local(Nrev = 0);         for(n=0,N-1,             Nrev = trinreverse(n, Nlg3);             /*print(Nrev);*/             reversed[Nrev+1] = v[n+1];            );         return(reversed); } /* Do FFT (radix 3, DIT) */ FFTCore(v, voffs, Npow, isign) = {         local(N = 3^Npow);         local(r = vector(3));         local(rr = vector(3));         if(N==1, return([v[voffs]]));         r[1] = FFTCore(v, voffs, Npow-1, isign);         r[2] = FFTCore(v, voffs + N/3, Npow-1, isign);         r[3] = FFTCore(v, voffs + 2*N/3, Npow-1, isign);         for(j=0,2,             for(k=0,N/3-1,                 r[j+1][k+1] = r[j+1][k+1] * exp(-isign*2*Pi*I*j*k/N);                );            );         for(j=0,2,             rr[j+1] = sum(k=0,2,exp(-isign*2*Pi*I*k*j/3)*r[k+1]);            );         return(concat(concat(rr[1], rr[2]), rr[3])); } /* Perform discrete Fourier transform (DFT). */ DFT(v) = {         return(FFTCore(DoTritReversal(v), 1, floor(log(matsize(v)[2])/log(3) + 0.5), 1)); } /* Perform inverse discrete Fourier transform (IDFT). */ IDFT(v) = {         return((1/matsize(v)[2]) * FFTCore(DoTritReversal(v), 1, floor(log(matsize(v)[2])/log(3) + 0.5), -1)); } /* 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))); } /* Increase the number of interpolating points in an approximation */ IncreaseInterp(P, v, new) = {         local(xcs = GetAbscissas(new, P));         return(ConvToFourier(vector(new, i, DoInterpolationPer(P, v, xcs[i])))); } /* 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); pp(x) = sin(x) + (sin(x)^3)/6 + (3*sin(x)^5)/40; Periodizer(P, x) = P/(2*Pi) * pp((2*Pi)/P * x); /*InvPeriodizer(P, x) = P/(2*Pi) * asin((2*Pi)/P * x);*/ /* Multiply continuum-sum exp array by log(base)^x and integrate from -1 to x. */ IntegrationPhase(P, v) = {         /* This is the most complex part of the algorithm.          * We need to first multiply the series by log(base)^x, then          * integrate, and finally periodize that integral back to          * being a periodic function.          *          * Each term is of the form a_n base^(qx) exp(2piinx/P) =          * a_n exp((2piin/P + q log(base)) x). Multiplying          * this by log(base)^x gives a_n exp((2piin/P + q log(base) + log(log(base)))x).          * The integral of these terms from -1 to x is          * a_n/(2piin/P + q log(base) + log(log(base))) *          *               (exp((2piin/P + q log(base) + log(log(base)))x) -          *                exp(-(2piin/P + q log(base) + log(log(base))))).          *          * So we create an auxiliary array of a_n/(2piin/P + q log(base) + log(log(base)))          * and then we evaluate it at the periodized abscissas. There may be a          * faster algorithm for this that is like the FFT.          */         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 this, first create an array of nodes for an initial guess with Code:`vt = ConvToFourier(vector(, i, 1))` This starts us off with a function that is 1 everywhere. Then choose a period for the approximation, say an imaginary one, like PP = 128*I. Now you can iterate with, say, Code:`vt = NormTetIter(PP, , vt)` and to use the approximation, try TetApprox(PP, vt, , x) for a given value x to evaluate at. It seems that with a large number of nodes, the convergence with this initial guess only works for small bases like base 2 -- not sure what's going on there, but that can be used as an initial guess for a wide variety of other bases. The biggest problem with this implementation is the final integration/periodization phase, which is as dog-slow as a plain DFT and has the same crappy $O(N^2)$ growth. I have no idea if it could be accelerated in some manner akin to the Fast Fourier Transform. EDIT: I should mention that I use a different -- perhaps "better", at least for the values of the tetrational near the real axis, periodizing function in the program. « Next Oldest | Next Newest »

 Messages In This Thread Numerical algorithm for Fourier continuum sum tetration theory - by mike3 - 09/15/2010, 11:00 AM RE: Numerical algorithm for Fourier continuum sum tetration theory - by mike3 - 09/15/2010, 11:06 AM RE: Numerical algorithm for Fourier continuum sum tetration theory - by nuninho1980 - 09/15/2010, 01:28 PM RE: Numerical algorithm for Fourier continuum sum tetration theory - by mike3 - 09/15/2010, 07:17 PM RE: Numerical algorithm for Fourier continuum sum tetration theory - by nuninho1980 - 09/16/2010, 01:18 AM RE: Numerical algorithm for Fourier continuum sum tetration theory - by mike3 - 09/16/2010, 08:00 PM RE: Numerical algorithm for Fourier continuum sum tetration theory - by nuninho1980 - 09/17/2010, 12:00 AM RE: Numerical algorithm for Fourier continuum sum tetration theory - by nuninho1980 - 09/17/2010, 12:16 AM RE: Numerical algorithm for Fourier continuum sum tetration theory - by mike3 - 09/17/2010, 10:10 AM RE: Numerical algorithm for Fourier continuum sum tetration theory - by sheldonison - 09/17/2010, 10:32 PM RE: Numerical algorithm for Fourier continuum sum tetration theory - by mike3 - 09/18/2010, 04:11 AM RE: Numerical algorithm for Fourier continuum sum tetration theory - by Ansus - 09/18/2010, 01:38 AM RE: Numerical algorithm for Fourier continuum sum tetration theory - by mike3 - 09/18/2010, 04:12 AM

 Possibly Related Threads... Thread Author Replies Views Last Post fast accurate Kneser sexp algorithm sheldonison 38 107,526 01/14/2016, 05:05 AM Last Post: sheldonison Attempt to make own implementation of "Kneser" algorithm: trouble mike3 9 22,800 06/16/2011, 11:48 AM Last Post: mike3 Road testing Ansus' continuum product formula mike3 33 50,469 09/23/2009, 08:26 PM Last Post: mike3 Software for numerical calculation ala qbasic Ivars 3 8,281 02/22/2008, 04:40 PM Last Post: Ivars numerical examples for approximation for half-iterates for exp(x)-1 Gottfried 0 3,777 08/14/2007, 11:57 AM Last Post: Gottfried

Users browsing this thread: 1 Guest(s)