Consistency in the composition of iterations

[JmsNxn]

Very cool, given you are getting out to \( z^{100} \) is impressive.

Honestly, all the heavy lifting is done by pari-gp. Granted it's not as user friendly as mathematica, or matlab, but it's a far more powerful language. Especially when you want to iterate exponentials, very large values are much better handled in pari-gp.

See Schroeder Summations for how I derive the Taylors series for an iterated function.

Oh yes, the dreaded matrix coefficient approach. I apologize but I refuse to do that on principal. The majority of my code is recursive in nature, I come from a hard C programming background--only use recursion and basic tools. So much of my code never, and I repeat never uses matrices. I think this is very important just on principle. If you use Sheldon's fatou.gp for example. It's an absolutely beautiful program, and works unbelievably well. But it doesn't account for Taylor series within Taylor series. This is because it uses the matrix approach to grabbing taylor series (essentially, Sheldon has a fuck ton of finesse to it). I don't like that.

When I calculate, for example \(\text{Iexp(s,z,y)} = \exp^{\circ s}_{y^{1/y}}(z)\), I want it to be evaluated recursively, and solely through taylor series. This allows us to add polynomials as input, and it doesn't lag at all. It runs at a reasonable speed.

Code:

Iexp(0.5+s,10,3)

%29 = 15.934117569230230385092587415316582681867959761026 + 19.638343545515097342*s + 24.578667286270856846*s^2 + 27.805807903067047476*s^3 + 29.353331414772487260*s^4 + 29.422879958436546130*s^5 + 28.310941611212828932*s^6 + 26.344062339323627480*s^7 + 23.833324573210982062*s^8 + 21.047238901998545260*s^9 + 18.199431760248257943*s^10 + 15.446854174840686297*s^11 + 12.894590695952666079*s^12 + 10.604128353043934001*s^13 + 8.602820474382437578*s^14 + 6.893074202766851433*s^15 + 5.460427999971121245*s^16 + 4.280150706607455280*s^17 + 3.322303285598587328*s^18 + 2.555388531180141137*s^19 + 1.9488058172990350241*s^20 + 1.4743572297809474250*s^21 + 1.1070422690004214264*s^22 + 0.8253485866395010970*s^23 + 0.6112081756842763786*s^24 + 0.4497497681688004363*s^25 + 0.3289432274340784669*s^26 + 0.23920242703878278126*s^27 + 0.17298996605216733209*s^28 + 0.12444969269431402824*s^29 + 0.08908057983965139875*s^30 + 0.06345705790081507970*s^31 + 0.04499552904418781962*s^32 + 0.03176363976661710682*s^33 + 0.022327301476023239018*s^34 + 0.015629894989807078116*s^35 + 0.010898183551020451669*s^36 + 0.007569911408994115700*s^37 + 0.005238691584274143002*s^38 + 0.003612464697276398120*s^39 + O(s^40)

Iexp(0.5,10+q,3)

%27 = 15.934117569230230427624883880746444576390287061258 + 2.7999731880342584288987705396231613259151125660196*q + 0.15929307687008597870758915334367542451218973925945*q^2 + 0.0035429056491437648784091813679715440313180399013116*q^3 + 3.2717950563800581641208783437659022653046137007582 E-5*q^4 + 1.4571989939034812415993316358928634742977974470800 E-7*q^5 - 1.4466726941573874659675040432308772794562975226253 E-9*q^6 + 1.8050385846879821670404797097163885632417992338526 E-10*q^7 - 1.6935532881258561009117698761420420144296804650632 E-11*q^8 + 1.4910468444524221193344497904266317864710197759836 E-12*q^9 - 1.2596085327034237280797737658677188456604434635678 E-13*q^10 + 1.0318273868462477456938834717873582056391994757241 E-14*q^11 - 8.2262694159076744974636043343064249631702964322443 E-16*q^12 + 6.3774861359175543513936022397691716010904685457160 E-17*q^13 - 4.7833951063140716523179710782841047788609352545588 E-18*q^14 + 3.4335856144645207547218129101385524741887599617793 E-19*q^15 - 2.3081768910860619894741241163812341529387132184486 E-20*q^16 + 1.3841086481882662760726690018792833312651955613400 E-21*q^17 - 6.3822973790649131629338457059963855270295086623230 E-23*q^18 + 4.8717554637806841458878796505377410923291157859884 E-25*q^19 + 4.0436514954431147717534415397104478834384856350286 E-25*q^20 - 7.3894971947386081611749144827388430747026042506157 E-26*q^21 + 9.7105164155666903492872703543172685711577988304886 E-27*q^22 - 1.1149521878709878783414651461706446137737077112424 E-27*q^23 + 1.1834097757486175675340261705553734772448740110044 E-28*q^24 - 1.1878698210223379602964732976691317573771693436787 E-29*q^25 + 1.1386582855437756925572682483039696346678072378017 E-30*q^26 - 1.0451528178485135210062164097957694913707338216443 E-31*q^27 + 9.1592987798608827764512265280826456211330618246676 E-33*q^28 - 7.5888844171910392218120856177312582812125981285382 E-34*q^29 + 5.8135184930724904804464414779359913800500073474853 E-35*q^30 - 3.90149681507602391668780539679898147167[+++]

Iexp(0.5,10,3+u)

%28 = 15.934117569230230427624883880746444576390287061258 - 0.52385183110766711998424421569174763942871575798842*u - 0.70690057991520575351571261184030325502313074504686*u^2 + 0.44411531342512746547015180738321828208297351126876*u^3 - 0.19308410582405174722014689099569675587295081984791*u^4 + 0.071848131189198944443170083384767685985892827911353*u^5 - 0.024142723142115509036122031466496879526081913862086*u^6 + 0.0070833706999478756045348076171032712818823059262224*u^7 - 0.012982609172377460357410069142856623324988837371939*u^8 - 0.32308949317899733445883466962532456292188922535443*u^9 - 8.6066998598379194230566413127854317765811984574804*u^10 - 207.69476555102571281432527634312719247142254128793*u^11 - 4585.1950885398103578776438170781623005630942711447*u^12 - 93250.241610458153284283320510407932467313593288547*u^13 - 1757412.4383463088105764127273494128605189377949385*u^14 - 30849265.994666013036946264651122856215060783317666*u^15 - 506631652.14371723099678035464573931008156211216010*u^16 - 7814666233.0757550918746067614459270747112401810738*u^17 - 113605431861.66057497151516947354232252352002098782*u^18 - 1561332556196.3480396911547776487218140444871901243*u^19 - 20342239475345.397913656031820154522243811706935980*u^20 - 251878416787650.03635825647325973119456185565738244*u^21 - 2970658101377771.2412558070944015165720411075651196*u^22 - 33440822398454447.465226779294817079273039718215542*u^23 - 359981668460959850.01665627755110335614113806325874*u^24 - 3712038733315423095.3559122839161601200381142104264*u^25 - 36725154490928027925.612836429095687147692516783169*u^26 - 349117523218631764298.00970614405882970726851681481*u^27 - 3193211369997048482807.8730342424182534985159422099*u^28 - 28137206828461070792438.565325532112357962650121852*u^29 - 239135097093725305912053.06040873845601840962042167*u^30 - 1962422582557521506635791.6086945967593804364547216*u^31 - 15565956967737158639698476.637780430970132400631264*u^32 - 119457270745367753801494199.5873400746996137[+++]

So, inserting a polynomial, spits out a polynomial. This is something that really annoyed me about Sheldon's approach. I just did these guys as an example, and I did them quick, these aren't the most accurate. It would've taken too long, I'm just trying to make a point.

Essentially, I avoid everything and anything which will use a matrix, on principle. This is something I even do in mathematics. I just have like an irrational phobia. So for that reason all my code is recursive. And one of the benefits is being able to pull out Taylor series out of a hat. And do Taylor series while we're doing Taylor series, while we're doing taylor series. The entire algorithm runs off of loops like this. And avoids using, what I'd like to refer to as "newtonian root finding". When you use matrices, somewhere in your code, the program is guessing the value based off of \(Ax = b\), which is then itself, usually defined using newtonian root finding. And root finding will fail, when you start nesting more root findings within (the garbage will start to rise to the top, why you aren't seein O(100) yourself, mathematica is failing you).

Code:

Iexp(0.5+s,10+q,3)

%30 = (15.934117569230230427624883880746444576390287061258 + 2.7999731880342584288987705396231613259151125660196*q + 0.15929307687008597870758915334367542451218973925945*q^2 + 0.0035429056491437648784091813679715440313180399013116*q^3 + 3.2717950563800581641208783437659022653046137007582 E-5*q^4 + 1.4571989939034812415993316358928634742977974470800 E-7*q^5 - 1.4466726941573874659675040432308772794562975226253 E-9*q^6 + 1.8050385846879821670404797097163885632417992338526 E-10*q^7 - 1.6935532881258561009117698761420420144296804650632 E-11*q^8 + 1.4910468444524221193344497904266317864710197759836 E-12*q^9 - 1.2596085327034237280797737658677188456604434635678 E-13*q^10 + 1.0318273868462477456938834717873582056391994757241 E-14*q^11 - 8.2262694159076744974636043343064249631702964322443 E-16*q^12 + 6.3774861359175543513936022397691716010904685457160 E-17*q^13 - 4.7833951063140716523179710782841047788609352545587 E-18*q^14 + 3.4335856144645207547218129101385524741887599617793 E-19*q^15 - 2.3081768910860619894741241163812341529387132184486 E-20*q^16 + 1.3841086481882662760726690018792833312651955613398 E-21*q^17 - 6.3822973790649131629338457059963855270295086623259 E-23*q^18 + 4.8717554637806841458878796505377410923291157859711 E-25*q^19 + 4.0436514954431147717534415397104478834384856350279 E-25*q^20 - 7.3894971947386081611749144827388430747026042506152 E-26*q^21 + 9.7105164155666903492872703543172685711577988304814 E-27*q^22 - 1.1149521878709878783414651461706446137737077112450 E-27*q^23 + 1.1834097757486175675340261705553734772448740110015 E-28*q^24 - 1.1878698210223379602964732976691317573771693436790 E-29*q^25 + 1.1386582855437756925572682483039696346678072378045 E-30*q^26 - 1.0451528178485135210062164097957694913707338216367 E-31*q^27 + 9.1592987798608827764512265280826456211330618248063 E-33*q^28 - 7.5888844171910392218120856177312582812125981284386 E-34*q^29 + 5.8135184930724904804464414779359913800500073474744 E-35*q^30 - 3.9014968150760239166878053967989814716[+++]

With this I can grab taylor series in two variables (relatively fast, this part can lag, you don't see the s variable here, but it is being stored). So I just made the taylor series:

\[
\sum_{j=0}^{39} \sum_{k=0}^{39}c_{jk}s^jq^k + O(s^{40}q^{40})= \exp^{\circ 0.5+s}_{3^{1/3}}(10+q)\
\]

But the important part, is that I am nesting the amount of taylor series i want. If I write:

\[
\exp^{\circ \exp(s)}(1)\\
\]

I'd be able to grab taylor series as well.

I find the recursive approach definitely is slower. Matrix is much fucking faster. But recursive gives you much more fluidity in how you call things. I guess, we're just sacrificing time for a more fluid result.

EDIT:

I thought I'd do something funny just to make a point:

Code:

Iexp(sin(s),10,3)
%31 = 9.9999999999999999777022175474883791262524916173650 + 7.013761285086568797*s + 5.979552711854540612*s^2 + 3.553290722224585034*s^3 + 1.5426940905950570248*s^4 + 0.23989164387692366606*s^5 - 0.3213073171905238572*s^6 - 0.4086292120259967000*s^7 - 0.2866996241873069205*s^8 - 0.13488124477623914370*s^9 - 0.030399075560803115114*s^10 + 0.016875734098622677627*s^11 + 0.025841691564653603040*s^12 + 0.018282825751704728937*s^13 + 0.008303706138553343140*s^14 + 0.0016393063394353790689*s^15 - 0.0011743377902937807577*s^16 - 0.0015673197287081274631*s^17 - 0.0010224091808505190733*s^18 - 0.0004138225351810333758*s^19 - 4.691203532047800570 E-5*s^20 + 8.609126394027138524 E-5*s^21 + 8.899171848423584755 E-5*s^22 + 5.073748599652385864 E-5*s^23 + 1.6786599564808138607 E-5*s^24 - 8.418983150066324215 E-7*s^25 - 5.721731733421904335 E-6*s^26 - 4.606370046530561666 E-6*s^27 - 2.205538946145741645 E-6*s^28 - 5.043493261571677544 E-7*s^29 + 2.2191737354157169363 E-7*s^30 + 3.326492645475257261 E-7*s^31 + 2.1306233760445312169 E-7*s^32 + 8.100863061056472969 E-8*s^33 + 5.800713243957596227 E-9*s^34 - 1.8451725514362570687 E-8*s^35 - 1.6872107064662877546 E-8*s^36 - 8.614866386137634214 E-9*s^37 - 2.272557235219667921 E-9*s^38 + 5.804907505273511860 E-10*s^39 + 1.1154809183220092231 E-9*s^40 + O(s^41)

This is \(\exp^{\circ \sin(s)}_{3^{1/3}}(10)\). And it does this very fucking fast. It's only about \(15\) digits accurate (I did less iterations) but still.