fixed point formula

[sheldonison]

I'm in the process of updating Mike's merged fixed point code to work for a generic base in the complex plane. But first I needed a way to figure out both primary fixed points for a generic complex tetration base, and then to figure out if the base was inside or outside the Shell Thron region, which tells whether or not one of the fixed points is attracting. The formula to generate the base from the fixed point, \( b=L^{\frac{1}{L}} \) is pretty well known, but what about the function's inverse, if you don't have access to the Lambert W function? What if you want to get the fixed point \( b^L=L \) from the base, using something other than brute force? For this Taylor series formula, for a base b, z is the input to the Taylor series.
\( z=i\sqrt{(\log(\log(b))+1)} \). In this formula, z=0 corresponds to base \( \eta=\exp(1/e) \). Then the two primary fixed points are calculated by putting both the positive and negative square roots in the taylor series, which is accurate to >32 decimal digits for 1.05<=b<=25 or so, and gives the equation for the two fixed points.

Since I don't have access to the Lambert W function, (it is not a part of pari-gp), I think this Taylor series is very helpful. This Taylor series equation was inspired by the observation that you need to loop around \( \eta \), twice to return to your starting point, where the fixed point is slowly modified as you loop around \( \eta \). The log(log(b)) part of the equation was to maximize the radius of convergence, since there is a singularity the fixed point for b=1. The log(log(b)) moves the singularity at b=1 to infinity, so that helps this taylor series converge over a very large range of bases, and probably for any real base>1. I haven't had time to figure out the total range of convergence yet, but it is surprisingly large, and includes negative bases. For any base b, with abs(b)<4000, except in the immediate vicinity of b=1, and b=0, the Taylor series gives double precision accurate results for the fixed point.
- Sheldon
\( z=i\sqrt{(\log(\log(b))+1)} \)
\( L(z)=\sum_{n=0}^{\infty} a_n\times z^n \)

Code:

a0=   2.7182818284590452353602874713526624977572470937000
a1=   3.8442310281591168248636716374262768779881984009975
a2=   4.5304697140984087256004791189211041629287451561666
a3=   4.0577994186124010929116533950610700378764316454973
a4=   3.1310579579657891414705533466321408770463105412618
a5=   2.1392433777070640849472839389751781515471363324069
a6=   1.3201485117519207930224182716545764193987429888921
a7=   0.75026349567941987403067580894396345224122700450743
a8=   0.39664914883721664871255223561747527187781503373434
a9=   0.19675048132092445683854734270973946155931655263755
a10=  0.092284535541707553042271049603653939199684877326343
a11=  0.041128857617644107022183948771661806665481904754276
a12=  0.017499926566873175768104436013779089511797807188588
a13=  0.0071383045088424740927613278903873113685388260534036
a14=  0.0027990895495591130046356892963041542988072136991061
a15=  0.0010584397247571801837924827314069442927340627561410
a16=  0.00038690511378324604701554539024015099607150307841871
a17=  0.00013695921477906978645151135090639314281168212780761
a18=  0.000047060925749355490378155571865461380815063050655184
a19=  0.000015718459720078638139982218287502395265847645167202
a20=  0.0000051098833189178642859504435914542633052729886905741
a21=  0.0000016200756562979959213494113078792872671206842607515
a22=  0.00000050114395897125416867359285070262528766885907119212
a23=  0.00000015147122619645160295342350422242820377511300274821
a24=  0.000000044807348486894284474811959719674830325499646189215
a25=  0.000000012962878987582433364248157574392737421319847483067
a26=  0.0000000036777490474599150634033635884535264575305342994411
a27=  0.0000000010235317972261857886939462803228104812154269796464
a28=  0.00000000027892975104614435108979535736539187977244285193364
a29=  7.4913196693434785382673670840314908951403127366391 E-11
a30=  1.9717541191970651018708628895408641932956867632982 E-11
a31=  5.0915368080421762446029435219708682084557437930980 E-12
a32=  1.3066486276084097356358660994233585807052473233683 E-12
a33=  3.2467635216891226343648988598478613134554037958645 E-13
a34=  8.0328830300969462955204663241426866159060111914746 E-14
a35=  1.9884030816356147519332706021476848341643270874558 E-14
a36=  4.5604212355453915576497473836641686140815097764966 E-15
a37=  1.1375475027302582619728186577908772014830927317863 E-15
a38=  2.6388661035640136079598667240173987557818441629236 E-16
a39=  5.2684658313198763292443793305904997392068250988020 E-17
a40=  1.6456895069992339719420263709083815254122876574657 E-17
a41=  2.4709091962079383454024177450388535492756309345447 E-18
a42=  5.3707283652335457424977440370257556179749838766151 E-19
a43=  3.0349615536019343269887649835177688859124766896609 E-19
a44= -3.2787954527686945286108924002948599290667213646190 E-20
a45=  1.8345013123066661922492136389266941966903253270420 E-20
a46=  4.8481796161714217259736691210135868007908205669597 E-21
a47= -3.2824665245176595846850267392135417222393135087407 E-21
a48=  1.4956645179800024112898138808465573719639979782662 E-21
a49= -2.2810205766236867792442627645380203056636207166072 E-22
a50= -8.2025371766299993417377211085132693281802664744834 E-23
a51=  8.2820898235652527449736419417013483792769309534803 E-23
a52= -3.1987649031517273110151312741635882659324832087531 E-23
a53=  5.1028759500570603817123734231141374241956872771981 E-24
a54=  2.0745436935733458976804542834312653952415531741435 E-24
a55= -1.8974009421475881728609023489783737853747016759452 E-24
a56=  7.2399705291174268127829280413850875632584902836527 E-25
a57= -1.0843673887979929422153305212584762102272054654005 E-25
a58= -5.0277375699737374239824380947818496605390245813200 E-26
a59=  4.3954005288785835431790225413101260484756443385679 E-26
a60= -1.6441622742940879676195309996736139807832614792602 E-26
a61=  2.3377978549045549610283334095071378953325737820424 E-27
a62=  1.2153823729634690431272005530504305362371379491894 E-27
a63= -1.0217831791605172358442907077539785681280031170764 E-27
a64=  3.7583342703760100931146730323471004826516979556089 E-28
a65= -5.0842172513431152493662523617793691855778116872791 E-29
a66= -2.9298747678164593371749001842727339027440198115827 E-29
a67=  2.3843736042453943922739674466720939816392076010523 E-29
a68= -8.6404128743033872572631375344081227432768174762623 E-30
a69=  1.1147745631909792790933165843951598884838146816889 E-30
a70=  7.0532998709238197459676161152886849409067455445289 E-31
a71= -5.5840347009750660670277203195942046868244697211916 E-31
a72=  1.9968202962515529653768275377958257984435847587123 E-31
a73= -2.4624153585967391507506022808723761030685594008346 E-32
a74= -1.6971043513573951058377996551655952956127911075482 E-32
a75=  1.3121636252191637725547321104357974814481747681435 E-32
a76= -4.6366491360531875428354175436972903880388681629591 E-33
a77=  5.4759261055169523922054804602053975292206196019485 E-34
a78=  4.0835907618102687341427597190193171384178018179741 E-34
a79= -3.0931489399194580627859122578595983870151745833582 E-34
a80=  1.0813047412434062333080300695597773394885783093044 E-34
a81= -1.2252362103483261754655465170167431865554218135280 E-35
a82= -9.8300622452978664179128765338170742574844836034023 E-36
a83=  7.3130149232712161576150411896903892349851012529721 E-36
a84= -2.5316784946273245634403427410920911852929637185860 E-36
a85=  2.7569030159398000068528639302431024861323649994624 E-37
a86=  2.3678752475291946660830153577567179420330740924808 E-37
a87= -1.7337617902358690948312875490328268568214283564291 E-37
a88=  5.9489649580056122715011899267403080309530471951160 E-38
a89= -6.2353624456438702836611492343755661445986256939609 E-39
a90= -5.7084806755390469129435066349248018249886305456145 E-39
a91=  4.1209734225367270686731246993915008297939802940460 E-39
a92= -1.4025498788164094982161672209129834192788136568375 E-39
a93=  1.4169741817716908448099015810155102528257145161939 E-40
a94=  1.3774829427232759282536454610089933691524288647861 E-40
a95= -9.8186882484920920105275699610740814969572189953598 E-41
a96=  3.3168418862324804650808696747831173926792793354185 E-41
a97= -3.2341718090336638957257547112586099667893533476259 E-42
a98= -3.3272289382182763758362704190111486565549876399064 E-42
a99=  2.3446712211671867041330301250221977724670989098562 E-42
a100= -7.8660813856791455405739446872012737404587511037516 E-43
a101=  7.4117874092723636167895499876920034206925805479953 E-44
a102=  8.0449630071992383377823253016230949067068370627592 E-44
a103= -5.6107673309024071413165294113175917244659515740451 E-44
a104=  1.8703704290344439531751733113502708390740673057030 E-44
a105= -1.7049595279518816756078819029211146676550591619926 E-45
a106= -1.9472226441848835239296603130698560902214718601605 E-45
a107=  1.3452918010782273565438586849286267873298215727367 E-45
a108= -4.4581019582286464605966679709280753360940853885077 E-46
a109=  3.9356988280727884495539422587080135663067916265244 E-47
a110=  4.7179835992286895246829227989155589972139121713110 E-47
a111= -3.2315545147777352407953961572361820900127456511476 E-47
a112=  1.0650043716082824895276880145990182899014721900301 E-47
a113= -9.1147127786496282723800711785279226816331317383552 E-49
a114= -1.1443035082690186137554650391441932946474634683831 E-48
a115=  7.7760383631901205864928578474869898677279380698392 E-49
a116= -2.5495578122471133824914480360868584888184556180062 E-49
a117=  2.1173080587506716879335848811562106778592441836394 E-50
a118=  2.7781889012354101797726076421387624610649014703232 E-50
a119= -1.8741844847909908017849435440876945806324719777844 E-50
a120=  6.1154697933929833224163136341333661094234392463681 E-51
a121= -4.9324234160781771043738905360427488911054147505257 E-52
a122= -6.7516170974169804169671782560329059979919324533539 E-52
a123=  4.5240983802350759934936686549613478193076549996979 E-52
a124= -1.4695734330098409037700447100436584191843150894889 E-52
a125=  1.1521125191958365706864738454301151565496232238249 E-53
a126=  1.6423578046256356048202813386114874500427364236761 E-53
a127= -1.0936528703475943968388872126503033575039972975248 E-53
a128=  3.5375228294839848475909167953331987436497172454718 E-54
a129= -2.6978482304943505310735652292701489001943829219215 E-55
a130= -3.9987970359659040999742140405232720136297398130454 E-55
a131=  2.6474006812280438194473769206724269657102547304092 E-55
a132= -8.5291961019295574135442472443611745856312183189474 E-56
a133=  6.3323381046099492231654979714328684659496610032336 E-57
a134=  9.7449542393413625948030755094469137335603889989868 E-57
a135= -6.4168257129444529858674454277838539098605318969812 E-57
a136=  2.0595681078339135332218823360988545914331945409098 E-57
a137= -1.4896215397166610321145528593568864824078893656240 E-58
a138= -2.3768743298521482002068230707322255247633783728679 E-58
a139=  1.5572247840096992668430577433971873564773219712346 E-58

[jaydfox]

(...) I needed a way to figure out both primary fixed points for a generic complex tetration base, and then to figure out if the base was inside or outside the Shell Thron region, which tells whether or not one of the fixed points is attracting. The formula to generate the base from the fixed point, \( b=L^{\frac{1}{L}} \) is pretty well known, but what about the function's inverse, if you don't have access to the Lambert W function? What if you want to get the fixed point \( b^L=L \) from the base, using something other than brute force? For this Taylor series formula, for a base b, z is the input to the Taylor series.
\( z=i\sqrt{(\log(\log(b))+1)} \). In this formula, z=0 corresponds to base \( \eta=\exp(1/e) \). Then the two primary fixed points are calculated by putting both the positive and negative square roots in the taylor series, which is accurate to >32 decimal digits for 1.05<=b<=25 or so, and gives the equation for the two fixed points.

Since I don't have access to the Lambert W function, (it is not a part of pari-gp), I think this Taylor series is very helpful. This Taylor series equation was inspired by the observation that you need to loop around \( \eta \), twice to return to your starting point, where the fixed point is slowly modified as you loop around \( \eta \). The log(log(b)) part of the equation was to maximize the radius of convergence, since there is a singularity the fixed point for b=1. The log(log(b)) moves the singularity at b=1 to infinity, so that helps this taylor series converge over a very large range of bases, and probably for any real base>1. I haven't had time to figure out the total range of convergence yet, but it is surprisingly large, and includes negative bases. For any base b, with abs(b)<4000, except in the immediate vicinity of b=1, and b=0, the Taylor series gives double precision accurate results for the fixed point.
- Sheldon
\( z=i\sqrt{(\log(\log(b))+1)} \)
\( L(z)=\sum_{n=0}^{\infty} a_n\times z^n \)

Hey Sheldon, I see we found the same Taylor series, though I couldn't tell what algorithm you used to derive it. I have an iterative fixed point finder with quadratic convergence, which I posted perhaps 6-7 years ago. I didn't have a good method of "seeding" the iterative solver, so about 4-6 months ago I decided to come up with one. My initial plan was just to do a polynomial regression to give me decent estimates, but as I was playing around with it, I noticed the same thing you did, that near base eta, the fixed points follow a sort of sqrt-ish behavior. I derived a method to calculate the Taylor series to arbitrary precision, which I initially coded in SAGE/python.

Then, after doing all that work, I found this thread.

Anyway, after that, I was working on making improvements to your kneser.gp / tetcomplex.gp code (which I'll get around to posting some day), so I ported the code to PARI/gp. In the process, I checked your coefficients, and the accuracy degrades pretty quickly. By the last few coefficients, they're not accurate to more than a few decimal places.

I've been meaning to post it here for a few months, but I got busy with other projects (most recently, the binary partition sequence). MorgothV8's recent thread on fixed points prompted me to finally post my fixed point code. If you're interested, the code is posted in that thread:
http://math.eretrandre.org/tetrationforu...60#pid7460

By the way, you mentioned that there is a large range of convergence. When I first analyzed the Taylor series, it seemed to have a radius of convergence of about 2.5-2.6, give or take. Keep in mind we're working in units of sqrt(log(log(z))+1), which corresponds to a radius of about 6.0-6.5 in units of log(log(z))+1. That number is really close to 2*pi, so I was initially tempted to assume that it had something to do with a singularity in a nearby branch.

However, when I looked at where the values of the function "explode", it was near large negative values of log(log(z))+1, which corresponds to bases just above 1, e.g., 1.0001.

Well, when I thought about it, this actually makes sense. For bases very close to 1, there will be a fixed point very close to the base, and another fixed point (on the other side of e) that is very large. By very large, we could be talking a thousand, a million, a billion, etc. This is what is constricting the radius. The radius of convergence should be infinite (I think), but the coefficients in the Taylor series decay so slowly that it behaves like a singularity when you approach base 1.

I've only calculated the first 500 or so coefficients of the Taylor series for analysis, but I've found a method to roughly estimate the rate at which the coefficients decay, and it's very slow. In fact, the "effective radius" seems to climb up to about 2.6 before actually coming back down to 2.5 or less, and then eventually, slowly, climbing back up.

Consider that for a base just slightly larger than 1, the upper fixed point will be very large.

For example, for a base with an upper fixed point of 2^12 (4096), the value of abs(sqrt(ln(ln(z))+1)) is about 2.28. For 2^20 (1048576), the value of abs(sqrt(ln(ln(z))+1)) is about 3.20. Think about how many terms it would take, at a radius of 3.20, to be able to encode a value of about a million. A hundred thousand at least, I would think. So even by a hundred thousand terms, the function will behave like it has a radius of convergence of 3.5 or less. At least, that's how I perceive it, intuitively. My complex analysis isn't sharp enough to do a formal analysis.

Anyway, the good thing is, on the top side, a radius of convergence of 2.5 implies that we can find fixed points for bases up to about exp(exp(2.5^2-1)), which is about 10^80 (rounding conservatively). That's provided we use enough terms in the Taylor series. If we conservatively stay within a radius of 2, that still gives us a range up to exp(exp(2^2-1)) = e^e^3 = 5.3 x 10^8, down to about 1.00676.

I haven't tried calculating anything that large yet, so I don't know how many terms are needed. I've played around with the small numbers, and they don't seem to blow up until about 1.001 or lower.

As for my algorithm, I thought I had a pretty TeX-ified description of it somewhere, but I can't find it. I at least had the presence of mind to put comments in my SAGE/python code (4-6 months ago), so I'll post those until I have time to TeX-ify it. This is just the comments: I started modifying the (SAGE/python) code, separate from the comments, so I still need to merge the code and comments back together.

Code:

# The base is developed as a function of the fixed point (fp).
    # base = B(fp)     = fp^(1/fp)
    # ln(ln(base)) + 1 = ln(ln(fp^(1/fp))) + 1
    #                  = ln(ln(fp)) - ln(fp) + 1
    # Next we calculate the derivatives at fp=e to develop a power
    # series in fp. Given fp=e, then ln(ln(fp)) - ln(fp) + 1 = 0.
    # The first derivative is: 1/(fp*ln(fp)) - 1/fp = 1/e - 1/e = 0.
    # Therefore, the first non-zero term is the second order term.

    # Since the first non-zero term is the second order term, we will
    # take the square root of the power series.
    # sqrt(ln(ln(base)) + 1) = sqrt(ln(ln(fp)) - ln(fp) + 1)

    # Reversion of the RH series gives us the fixed point as a function of
    # the base. However, note that the left-hand side was sqrt(ln(ln(base)) + 1)
    # Finally, recall that we developed the original power series at fp=e,
    # so add e back to the final answer.

One final thought. I did an analysis, and I found that precision decays at about 3 bits per term in the Taylor series. Therefore, to calculate 200 terms to 150 bits of precision, you need to start with about 150+3*200, or 750 bits, plus a little more just to be safe. This turns out to be easy in PARI/gp, because precision is measured in decimal digits. Each decimal digit of precision is about 3.3 bits. So I just added an extra digit of precision for each term.

[jaydfox]

By the way, you mentioned that there is a large range of convergence. When I first analyzed the Taylor series, it seemed to have a radius of convergence of about 2.5-2.6, give or take. Keep in mind we're working in units of sqrt(log(log(z))+1), which corresponds to a radius of about 6.0-6.5 in units of log(log(z))+1. That number is really close to 2*pi, so I was initially tempted to assume that it had something to do with a singularity in a nearby branch.

However, when I looked at where the values of the function "explode", it was near large negative values of log(log(z))+1, which corresponds to bases just above 1, e.g., 1.0001.

Well, when I thought about it, this actually makes sense. For bases very close to 1, there will be a fixed point very close to the base, and another fixed point (on the other side of e) that is very large. By very large, we could be talking a thousand, a million, a billion, etc. This is what is constricting the radius. The radius of convergence should be infinite (I think), but the coefficients in the Taylor series decay so slowly that it behaves like a singularity when you approach base 1.

I've only calculated the first 500 or so coefficients of the Taylor series for analysis, but I've found a method to roughly estimate the rate at which the coefficients decay, and it's very slow. In fact, the "effective radius" seems to climb up to about 2.6 before actually coming back down to 2.5 or less, and then eventually, slowly, climbing back up.

Okay, so I've calculated the Taylor series out to 1600 terms. It took almost two full days. The bulk of the time (over 99%) was spent on the series reversion. I'm using SAGE/python to do the series reversion, but I checked the source code, and it's using the PARI library to do reversion. So the SAGE code and the PARI code should take just about as much time, for a given level of precision. And as I mentioned before, the minimum precision during the reversion step is dictated by the number of terms in the series.

Anyway, I should have trusted my gut about the 2*pi thing. The root test is still limited to about 2.51, or really close to sqrt(2*pi). I did an analysis of where the function "blows up", and found something interesting. In terms of absolute value, the function blows up in the direction of theta=pi/2, e.g., sqrt(2*pi)*I = sqrt(-2*pi). However, there is no obvious singularity in that direction.

Instead of looking at where the absolute value is the largest, I looked for anything strange, and found it at odd multiples of pi/4. So at theta=pi*(2*k+1)/4, e.g., sqrt(2*pi)*(sqrt(1/2)+sqrt(1/2)*I), there is a singularity. It's not an "infinite" singularity. Instead, it turns out to be the location of base eta, in other branches. Once I realized this, it was like, duh! Why didn't I think of that! You see, the value is well-defined at base eta, unlike say, base 0 or 1. But it's an algebraic singularity, because there are pairs of fixed points that converge when we get to that base.

Remember, we're working with sqrt(ln(ln(base))+1). Well, consider base eta. The first logarithm is 1/e. The second iterated logarithm is -1. Add 2*pi*I and then add 1, and now you have sqrt(2*pi*I). We could just as easily have added -2*pi*I, and the sqrt has two branches, so this gives us all four possibilities for the sqrt: sqrt(2*pi)*sqrt(1/2)*(+/-1 +/- I) = sqrt(pi) * (+/-1 +/i I)

\(
\sqrt{\pi} (\pm 1 \pm i)
\)

This means that, without analytically continuing the function, we will be limited to a radius of convergence of about 2.5. For reference, sqrt(pi) is about 2.506628... Edit:Oops, I forgot the factor of 2: sqrt(2*pi) = 2.506628... End Edit

Now, I haven't put much thought yet into how to analytically extend the function. However, just to get something to work with, we could calculate a bunch of fixed points for a list of bases, and then use that information to come up with an approximate Taylor series. This would boil down to either polynomial interpolation or least squares regression. You can do both at the same time if you use a discrete fourier transform, or DFT (use a list of equally spaces points on the circumference of a circle). I believe this is what Sheldon refers to as a Cauchy integral. It has the benefit of being a polynomial interpolation and a least squares fit, which is why it's so useful.

Another benefit of a DFT (compared to, say, Lagrange or Newton interpolation), is that it uses considerably less memory, and it's quite a bit faster if you use an FFT (Fast Fourier Transform). I'll see what I can put together over the next few days.

[jaydfox]

Now, I haven't put much thought yet into how to analytically extend the function. However, just to get something to work with, we could calculate a bunch of fixed points for a list of bases, and then use that information to come up with an approximate Taylor series. This would boil down to either polynomial interpolation or least squares regression. You can do both at the same time if you use a discrete fourier transform, or DFT (use a list of equally spaces points on the circumference of a circle). I believe this is what Sheldon refers to as a Cauchy integral. It has the benefit of being a polynomial interpolation and a least squares fit, which is why it's so useful.

Another benefit of a DFT (compared to, say, Lagrange or Newton interpolation), is that it uses considerably less memory, and it's quite a bit faster if you use an FFT (Fast Fourier Transform). I'll see what I can put together over the next few days.

Okay, here's a somewhat stable version of my DFT/FFT library, a slightly updated version of the fixed point finder, and a piece of code that calculates the fixed point Taylor series using a DFT (i.e., a "Cauchy integral").

I was able to calculate 512 terms to 96 digits of precision, in just a couple seconds. It took my other code (closed-form calculation of terms, using reversion of a power series) over 2 minutes to get the same result. I then calculated 1600 terms, again in only a few seconds. It took nearly two days using my other code.

Needless to say, in this specific scenario, the DFT method is extremely fast. And just as accurate.

So, the files. First, the fixed point iterative solver:
  jsloglib_v0.7.gp (Size: 6.3 KB / Downloads: 991)
It also has code for finding the inverse Schroeder function at the fixed point of an arbitrary Taylor series. This is used to find the inverse Schroeder function for the exponential function, but I assume it can be used for arbitrary functions. Edit: I forgot to mention, in this version of the library, I decided to go with Sheldon's convention of applying a factor of sqrt(-1) to the "z" term: z = I*sqrt(ln(ln(base))+1). The benefit of this is that the Taylor series takes on real coefficients only, which makes it a little easier to study. I haven't double checked all my other code yet to make sure I didn't break something in the process, so if you get an unexpected result, check with the previous version of the code. End Edit

The second file is the DFT library:
  jDFTlib_v0.6.gp (Size: 16.45 KB / Downloads: 959)
It has a few comments here and there, but it needs more. I was originally writing it to speed up the "Cauchy integrals" in Sheldon's kneser/tetcomplex code, by converting them from simple DFT's to Fast Fourier Transforms. But FFT's turned out to be so much fun, that I got carried away and built a library, useful for other applications as well.

Finally, here's the code that uses these two libraries to calculate the fixed point formula:
  jdft_fp_v0.2.gp (Size: 6.76 KB / Downloads: 988)
There's a lot of experimental code in here, so I apologize in advance for any bugs, or for any code that's hard to reverse engineer.

Note the performance of the FFT. I can use thousands of samples and still get a result in a few seconds. Heck, I can do tens of thousands of samples and still get a result in a reasonable timeframe. I did 8! (40320) samples in 11 seconds at 67 digits of precision, or 30 seconds at 144 digits.

Update: I suppose I should post a screenshot. The following three lines of code were typed at the prompt. It's hard to tell because of the long help text that displays when you load the file):

Code:

\r jdft_fp_0.2.gp
#
DFT_FPFunc(-30, 2.0, 67)

   

[sheldonison]

...
Another benefit of a DFT (compared to, say, Lagrange or Newton interpolation), is that it uses considerably less memory, and it's quite a bit faster if you use an FFT (Fast Fourier Transform). I'll see what I can put together over the next few days.
....

Okay, here's a somewhat stable version of my DFT/FFT library, a slightly updated version of the fixed point finder, and a piece of code that calculates the fixed point Taylor series using a DFT (i.e., a "Cauchy integral").

I was able to calculate 512 terms to 96 digits of precision, in just a couple seconds. It took my other code (closed-form calculation of terms, using reversion of a power series) over 2 minutes to get the same result. I then calculated 1600 terms, again in only a few seconds. It took nearly two days using my other code.

Needless to say, in this specific scenario, the DFT method is extremely fast. And just as accurate.
.....
Note the performance of the FFT. I can use thousands of samples and still get a result in a few seconds. Heck, I can do tens of thousands of samples and still get a result in a reasonable timeframe. I did 8! (40320) samples in 11 seconds at 67 digits of precision, or 30 seconds at 144 digits.

I missed the posts in the computation section! This looks like great fun. I've never written an FFT, but I would like to. Its just that a generic DFT works pretty good within reasonable radius of convergence scenarios.

Back to the original topic, generating the fixed points from the \( \ln(\ln(b))+1 \). The \( k=\ln(\ln(b))+1 \) form is also very interesting because it means the tetration itreration function switches from iterating \( z \mapsto b^z \) to iterating \( z \mapsto \exp(z)-1+k \). Here, \( \text{sexp}_k(z) \mapsto \text{sexp}_b(z)\cdot \ln(b)+\ln(\ln(b)) \)

Another comment. Sometime after this post, I did succeed in getting a complex base tetration algorithm to work; its posted online here http://math.eretrandre.org/tetrationforu...hp?tid=729. But it doesn't converge for bases near the Shell Thron boundary. Its also dog slow for bases close to \( \eta=\exp(1/e) \). So I'm thinking about a tetration algorithm that iterates generating the slog, centered exactly between the two fixed points, using the "k" form ... and using the very same fixed point Taylor series as from above; well actually the \( L_k \mapsto \ln(\ln(L_b))+1 \). The algorithm I have in mind should have convergences of \( n^{-r} \), where I think r should be the real(pseudo period+1), so a 500 Taylor series terms for base e would be accurate to 15 decimal digits or so. Smaller values of k (base closer to eta) have larger real periods.... An FFT would help because it allow for thousands of terms instead of hundreds of terms.... So I may eventually need your FFT.

[sheldonison]

\( z=i\sqrt{(\ln(\ln(b))+1)} \)
\( L(z)=\sum_{n=0}^{\infty} a_n\times z^n \)
....

This evening, I stumbled across this equivalent formal power series for the same fixed-point function.... Instead of the fixed point function above, the function is the ln(ln(fixed-point)) of the fixed-point function quoted above.
\( z=\sqrt{-2\left(\ln(\ln(b))+1\right)} \)
\( L(b)=\exp(\exp(\text{xfixed}(z))) \)

Here are the first 60 terms of the "xfixed" formal power series; update xfixed is defined as follows: see my answer to this mathstack question
\( g(x)=\sqrt{2(\exp(x)-x-1)};\;\;\;\;\text{xfixed}=g^{-1}(x) \)

Code:

a1=   1
a2=  -1/6
a3=   1/36
a4=  -1/270
a5=   1/4320
a6=   1/17010
a7=  -139/5443200
a8=   1/204120
a9=  -571/2351462400
a10= -281/1515591000
a11=  163879/2172751257600
a12= -5221/354648294000
a13=  5246819/10168475885568000
a14=  5459/7447614174000
a15= -534703531/1830325659402240000
a16=  91207079/1595278956070800000
a17= -4483131259/2987091476144455680000
a18= -2650986803/818378104464320400000
a19=  432261921612371/337123143997663268044800000
a20= -6171801683/24551343133929612000000
a21=  6232523202521089/1189370452023756009662054400000
a22=  4283933145517/279517041579788632620000000
a23= -25834629665134204969/4267461181861236562667451187200000
a24=  11963983648109/10062613496872390774320000000
a25= -1579029138854919086429/76814301273502258128014121369600000000
a26= -208697624924077/2747093484646162681389360000000
a27=  746590869962651602203151/24887833612614731633476575323750400000000
a28= -29320119130515566117/4968750398484053826709377112800000000
a29=  1511513601028097903631961/17321932194379853216899696425330278400000000
a30=  2700231121460756431181/6931406805885255088259581072356000000000
a31= -8849272268392873147705987190261/57510547078560550665428682101739057315840000000000
a32=  10084288256532215186381/332707526682492244236459891473088000000000
a33= -6208770108287283939483943525987/15842905490859115174615483904200811615354880000000000
a34= -6782242429223267933535073/3308776352857385368931593620699860160000000000
a35=  2355444393109967510921431436000087153/2907806873792281999148925915777016963882234675200000000000
a36= -51748587106835353426330148693/323935822497443742389140878653757709384320000000000
a37=  2346608607351903737647919577082115121863/1278155589444135475545901875538945576644075073830912000000000000
a38=  7007277101869903281324331583/636701444219113562626932071836696187410560000000000
a39= -2603072187220373277150999431416562396331667/598176815859855402555482077752226529869427134552866816000000000000
a40=  585302872633292617248814587726421/681334215458873423367080010072448590148040256000000000000
a41= -73239727426811935976967471475430268695630993/8240483815285368025604321103114672675481228205600293257216000000000000
a42= -110855495796575034381969281033555329/1845734389678088103901419747286263230711041053504000000000000
a43=  34856851734234401648335623107688675640839679447003/1466970928797101215918081242776474029689168245160964205649592320000000000000
a44= -18447986573777204063499607563765439/3929628055443671447015925913577205587965442242944000000000000
a45=  909773124599542506852275229422593983242880452145053/20596271840311301071489860648581695376835922162059937447320276172800000000000000
a46=  38650132745379700438031566826935471987259957/116354250486669482021689508823610789591151200820537602492800000000000000
a47= -1527335577854677023023224272800947125313629267269390501/11616297317935573804320281405800076192535460099401804720288635761459200000000000000
a48=  217784448556937372678947372805330071920344629/8377506035040202705561644635299976850562886459078707379481600000000000000
a49= -183856455668177802003316143799518064719008299958634826921/819645938753534087632839055993253376145302064613791341063566139328561152000000000000000
a50= -1167289109751840227800236733417523750884898531/628312952628015202917123347647498263792216484430903053461120000000000000000
a51=  2583312098861137963745902036370496943872138148651712093816393/3511363201620140031419082515875097463406474044805482105116317340883555975168000000000000000
a52= -107748081854646619391722638838074116233224341741059/740272037656801231924925556964805979417351539791645668557356972800000000000000000
a53=  5180134290822682443757710427952467581918233549140896702364013/4466453992460818119965072960193123973453034984992573237707955657603883200413696000000000000000
a54=  27346403208634415483181063970969158506217340077059/2607045002182647816779085657136925405774151074918404311006344121600000000000000000
a55= -527550309097873396592733540579928993424142983691519876840948418433873/126949575182314587341955272240273143385058268074945409878049452860509411817038440038400000000000000000
a56=  377036553764192941179202019520271416437725306277603527/458349795923735677963564607552557129439965048783554498726647372706739200000000000000000
a57= -2114866241537081164613223324215572812504648703648482437460602956015127/347334037698812710967589624849387320301519421453050641426343303026353750731417171945062400000000000000000
a58= -2107144283266473668026539971128155003797327940672559775477/35290642537148028524804656958509136181230108931089778629458214461555384704000000000000000000
a59=  180394412915538782140015777241228025103785450235726235175126981743099027459/7623287459413541380316657086194352905977748262051555478025382814822412121053144089850229555200000000000000000
a60= -907975882295290895046750344009772888231118193554130319911809/193745627528942676601177566702215157634953298031682884675725597393939062024960000000000000000000

The series above calculates the fixed point x for f(x) function, as z varies in the neighborhood of zero.  Here, if x is the desired fixed point L, then f(x)=x
\( f(x)=\exp(x)-1+k \;\mapsto \; f(x) = \exp(x)-1-\frac{ z^2}{2}\;\;\; k \mapsto -\frac{ z^2}{2};\; \) this mapping gives rational coefficients for L(z)
k=0 corresponds to the parabolic case, or \( f(x)=\exp(x)-1\; \) which has a parabolic fixed point of zero.
\( L(z) = \sum_{n=1}^{\infty}a_n z^n\;\; \) The Taylor series from above
\( L(z) = \exp(L(z))-1-\frac{ z^2}{2}\;\; \)plugging L(z) into f(x) as x
For \( \text{sexp}_e; \; k=1; \; f(x)=\exp(x); \; z=\pm\sqrt{-2} \)

With an 80 term formal power series, one gets the fixed point for base(e) accurate to 35 decimal digits.  Smaller bases closer to exp(1/e) are more accurate.  For base(2), the fixed point is accurate to 43 decimal digits.

So for example, for tetration base(e), the fixed point is \( L\approx 0.318132 + 1.337236i \)
To get this value for L and its conjugate, plug \( z=\pm\sqrt{-2(\ln(\ln(e))+1)}=\pm\sqrt{-2} \) into the series above...

\( L = \sqrt{-2}
- \frac{\left(\sqrt{-2}\right)^2}{6}
+ \frac{\left(\sqrt{-2}\right)^3}{36}
- \frac{\left(\sqrt{-2}\right)^4}{270}
+ \frac{\left(\sqrt{-2}\right)^5}{4320}
+ \frac{\left(\sqrt{-2}\right)^6}{17010}
- \frac{139\left(\sqrt{-2}\right)^7}{5443200}+... \approx 0.318132 + 1.337236i \)
For base(e) conveniently \( L=\exp(L) \)  For other bases, \( L=\exp(L)-1+k \)

For tetration base2, we first convert it to the conjugate form
\( b=\ln(\ln(2))+1 \approx 0.63348708   \)
\( f(x)\approx \exp(x)-1+0.63348708 \)
If we have a fixed point L of f(x), than the fixed point for tetration base(2) is \( \exp(\exp(L)) \)
Then the value z for the series above is \( z = \sqrt{-2\left( \ln(\ln(2))+1 \right)} \)

[mike3]

Want to point out that there are far superior methods:

http://math.eretrandre.org/tetrationforu...70#pid7470
http://math.eretrandre.org/tetrationforu...12#pid7512
and below.

This code calculates arbitrary branches of the Lambert W function for arbitrary inputs and is extremely fast (I can get thousands of digits of precision in a fraction of a second on my machine) by using Newton's method.