Numerical algorithm for Fourier continuum sum tetration theory  Printable Version + Tetration Forum (https://math.eretrandre.org/tetrationforum) + Forum: Tetration and Related Topics (https://math.eretrandre.org/tetrationforum/forumdisplay.php?fid=1) + Forum: Computation (https://math.eretrandre.org/tetrationforum/forumdisplay.php?fid=8) + Thread: Numerical algorithm for Fourier continuum sum tetration theory (/showthread.php?tid=516) Pages:
1
2

Numerical algorithm for Fourier continuum sum tetration theory  mike3  09/15/2010 Hi. Here's the longawaited description of the numerical algorithm for the Fouriercontinuumsumbased tetration. It actually took a LOT more to describe than I thought it would... Periodic approximation As was mentioned in the posts on the subject, the continuum sum for a holomorphic function is constructed as a limit of continuum sums of periodic approximations of it, that is, a sequence of periodic functions which converge pointwise to the target function . The idea behind the numerical algorithm, then, is to use a periodic function with a very large period as an approximation to a general function in the vicinity of zero. That is, we work with an approximation to the desired function given by (which in a computer would be truncated to a sum from to for some ) where is the period, and large. Because certain operations in the tetration continuumsum iteration destroy periodicity (especially the multiplication by and the final integral), we must have a way to restore it. We do this via what is called a periodizing function. A periodizing function generates a periodic approximation of a given function by composition, i.e. , with period . This means the function should, at bare minimum, satisfy 1. periodic with period : 2. Other desirable properties are 3. entire 4. for any compact subset of , and a given , there exists an such that whenever for all . An example of such a function is . Property 4 is sort of like compact convergence of a sequence of functions but with a complex index, and is nice because it means that past a certain point, we can just increase 's magnitude to make a better approximation of the identity around 0. Implementing the tetration formula We can now discuss how we'd implement the iteration. We work with an approximation of the form . (the reason for using instead of will be clear later) and perform the following operations. Continuum sum The continuum sum, , is achieved by setting , . (from .) And so we have . We can store separately for the next operation. Exponential The exponential of the continuum sum is done by treating the series as a pair of Taylor series in and , and then applying Faà di Bruno's formula and multiplying the results via the convolution, however this method is not very fast and seems to not be as stable as we might like (though keeping more "slop" terms beyond N and N index seems to help, but that just makes it even slower). Though with a good program like Pari/GP, the formal exponential of a power series is already implemented and so there's no need to implement all that ourselves. However, we can do a lot better. Fourier transform A Fourier series, such as the one above, can be thought of as a sort of "Fourier interpolation" between points in the "spatial" domain. It is possible to transform from one to the other and back via the discrete Fourier transform, or DFT. Let's see how that works. We first transform the function so it has period (note now that "H" is for "half" as in half of N, rounded down), giving . Now consider the "discrete sampling" of this at integer inputs, giving a sequence indexed from to : . This bears a great similarity to the formula for the inverse discrete Fourier transform (IDFT) of a sequence : . The only difference is the frequency elements involved, and the presence of the normalization factor . This means that our Fourier series can be thought of as an interpolation of a regularlyspaced grid of points along one period interval. Our sampling at discrete points then creates a list of all those points, to which we can then directly apply transcendental functions like exponential. This method may not be as accurate as the first (note the limited amount of terms), but it can be much faster and also seems numerically more stable and with enough terms to begin with, it can still be quite accurate. If we can relate the formula to the IDFT, then we can also give the transform to get back to Fourier space, and, even more powerfully, we can harness the highspeed Fast Fourier Transform (FFT) algorithm to rapidly perform this computation, blowing the other method out of the water. So how do we relate this to the IDFT? Consider the following process: First off, we renumber the elements of from to , giving: . where is the renumbered sequence. Note the exponent in the exponential  the displacement by H. Expanding this out, we see that we have . . and thus a regular IDFT (sans the normalization factor) on the right. However, by the Fourier shift theorem, the multiplication of the spatialdomain data by is equivalent a cyclical shift of by H. Namely, it corresponds to transforming . We can cancel the effects of this shift by cycleshifting the other way, forming and taking the inverse DFT of it multiplied by (to account for the normalization factor). There's one final point: the given formula was for numbered from to . However, the (I)DFT is usually done with the numbering from to . This is not a problem: if we use a regular (I)DFT algorithm, we'll get through followed by through . A simple right cyclic shift by resolves this and gives the data in the correct order, though such is not really needed for evaluating transcendental functions. Then to transform from the spatial domain back to the frequency domain, we do all that in reverse: left cycleshift by , preform forward DFT, divide by , then cycleshift the result right by and renumber. Logpower of base and final integral In the given program, I lump these two steps together. The end result of the last step should be a series of coefficients giving an approximation to the exponential of the continuum sum of : . We can now multiply by , giving and integrate, to get . (YES, I know I'm abusing the integral notation with the "z" in both the limits and the variable, but it's to try to make it clearer and more consistent.) We can then composite this with the periodizing function to get a periodic function out the end with the period we desired, and a pointwise evaluation and DFT will give us the new Fourier series for the next iteration cycle. The normalization step (for the derivativeat0 bit) is real easy, just evaluate at 0 and divide the thing by the result. RE: Numerical algorithm for Fourier continuum sum tetration theory  mike3  09/15/2010 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: /* Radix3 DIT FFT. */ To use this, first create an array of nodes for an initial guess with Code: vt = ConvToFourier(vector(<powerof3 number of nodes>, 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, <base>, vt) and to use the approximation, try TetApprox(PP, vt, <base>, 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 dogslow as a plain DFT and has the same crappy 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. RE: Numerical algorithm for Fourier continuum sum tetration theory  nuninho1980  09/15/2010 Code: (13:19) gp > vt = ConvToFourier(vector(3, i, 1)) therefore, is this variable "x" integer only? RE: Numerical algorithm for Fourier continuum sum tetration theory  mike3  09/15/2010 (09/15/2010, 01:28 PM)nuninho1980 Wrote: It should take any real or complex value. "PP" is not set, that should be the period (note in my post I mentioned "PP = 128*I"), in this case 128*I. This is what is causing your problem. Though I don't think this'll work with just 3 elements and such a huge period because there won't be enough "resolution" in the Fourier series for it to home in on the function. You should either increase the elements or decrease the PP, say 8*I instead. If you want to use 128*I, try at least 81 (3^4) terms, but even there a smaller PP like 32*I would be best. Usually I use 128*I with at least 243 terms. (I just tried it. 3 terms at 128*I is not enough.) RE: Numerical algorithm for Fourier continuum sum tetration theory  nuninho1980  09/16/2010 (09/15/2010, 07:17 PM)mike3 Wrote: It should take any real or complex value. "PP" is not set, that should be the period (note in my post I mentioned "PP = 128*I"), in this case 128*I. This is what is causing your problem.  vt = ConvToFourier(vector(3^4, i, 1))  vt = NormTetIter(32*I, 2., vt)  TetApprox(32*I, vt, 2., 0.5) > 1.636671277861037786181970844 + ...I 2^^0.5 = 1.63667  it's incorrect but yes 2^^0.5 = 1.45878181603642 (new knesser code by sheldonison) RE: Numerical algorithm for Fourier continuum sum tetration theory  mike3  09/16/2010 (09/16/2010, 01:18 AM)nuninho1980 Wrote:(09/15/2010, 07:17 PM)mike3 Wrote: It should take any real or complex value. "PP" is not set, that should be the period (note in my post I mentioned "PP = 128*I"), in this case 128*I. This is what is causing your problem. Try running NormTetIter a few more times over it. Like for(i=1,5,vt = NormTetIter(32*I, 2., vt)); RE: Numerical algorithm for Fourier continuum sum tetration theory  nuninho1980  09/17/2010 (09/16/2010, 08:00 PM)mike3 Wrote: Try running NormTetIter a few more times over it. Like1.458810483283048116375249186 + 8.53713983 E27*I thank! but it's very probably better code by sheldonison. RE: Numerical algorithm for Fourier continuum sum tetration theory  nuninho1980  09/17/2010 your code: Code:  vt = ConvToFourier(vector(3^4, i, 1)); code by sheldonison: Code:  init(exp(1.)) your code may be bad!? RE: Numerical algorithm for Fourier continuum sum tetration theory  mike3  09/17/2010 The problem is the initial guess. I mentioned in the main posts: Quote: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. I wasn't specific enough  the base has to be real small, even e seems too large, at this initial guess. Not sure why this is. Try a loop or two on base 2 first, then use the result of that as the initial guess to run base e. On my system, a single loop at base 2, followed by six loops at base e, using the parameters you give, yields tet(pi) ~ 3.716e10, which is closer to the expected value of 3.71498e10. You'll probably need more than 81 nodes and 32*I period to get more accuracy than that. Yes, the Kneser process seems to yield more precision more quickly, but is restricted only to real bases (where the conjugate symmetry holds). ADDENDUM: Actually, it does work in your case, but with more nodes than what you used, it won't. You just need to iterate it a lot more times  try 30 or so. RE: Numerical algorithm for Fourier continuum sum tetration theory  sheldonison  09/17/2010 Mike, thanks for your detailed description. So, am I correct that your approximating the sexp(z) with a long periodic function, to approximate switching from the space domain to the frequency domain. Then, to get a more accurate version of the sexp(z), you take sexp(z+1)=e^sexp(z), using Faà di Bruno's formula? I mean, high level overview, is that more or less correct? I do remember convolutions, and converting the time domain to the frequency domain, but its been 27 or 28 years.... So much more to learn, and yet so much already forgotten... So, as to the general applicability to complex domains, and getting other solutions then we're used to seeing  here's my intuitive feeling. Any analytic solution, especially one with limiting behavior matching the super function, is a 1cyclic transformation, via theta(z), of the superfunction. f(z)=regularsuper(z+theta(z)). To preserve the behavior at either +I*infinity or I*infinity, then theta(z) must go to zero at +I or I infinity. Then, I believe it is most likely that there will be a singularity in theta(z), and we usually get around the singularity with a schwarz transformation about the real axis. Maybe something similar applies in the complex domain. The other possibility, is that theta(z), when wrapped around the unit circle, is an annular Laurent series, with convergence between the two radius's of singularities. This hasn't been explored (to my knowledge), but I think that's what would be going on when converting between the regularsuperf_sqrt(2)(z) and the "regular iteration of log" developed from the other fixed point for sqrt(2). That's what I'd like to explore when I have time. Nininho, I found Andrew's tetration website via the wayback archive, Andrew's Website. I remember those exp^^pi calculations! I will add one more number to the list... e^^pi=37149801960.556985498745240109389505725500669657852560614478645123... This one is accurate to 65 decimal digits, with the Taylor series calculated in about 30 minutes, using the latest version of my kneser.gp code. I'll probably update the code online again; the current version is getting almost 2x more precision for the same amount of compute time (100 binary bits precision in under a minute), but it slows way down once the number of iterated logs used to calculate the inverse superf gets out of hand, so I was thinking of updating that part of the algorithm as well. Mostly, I'd like to figure out how to prove it converges...  Sheldon 