Consistency in the composition of iterations

[Daniel]

Edit: Changed \( O(n) \) to \( O(z^n) \).

I primarily focus on iterated functions instead of tetration because a general theory of complexly iterated functions serves as a strong foundation for researching the hyperoperators. I study complex tetration, but not real tetration because I don't believe that it can pass the following test. Please let me know if I am wrong.

Let

\( f(0)=0, f^0(z)=z, f^1(z)=f(z) \)

with

\( a,b,z\in \mathbb{C} \) and

\( n\in \mathbb{Z^+} \) where the Taylor series for \( f(z) \) obeys

\( f^a(f^b(z)) - f^{(a+b)}(z) = O(z^n) \).

I consider this the gold standard in extending the iteration of functions from the natural numbers.

I have derived and proven symbolically using Mathematica that the Taylors series for \( f(z) \) satisfies

\( f^a(f^b(z)) - f^{(a+b)}(z) = O(z^{19}) \)

and that iterated sin function

\( \sin^a(\sin^b(z)) - sin^{(a+b)}(z) = O(z^{39}) \)

Has anyone else done this sort of consistency test? If so, what were your results for \( o(z^n) \)?

[JmsNxn]

I'm confused...

What does \(O(n)\) mean, do you mean \(O(z^n)\)? Or do you mean the taylor coefficients grow like \(O(n)\), but then the \(O(constant)\) wouldn't make sense. I'm confused. Could you clarify?

Because:

\[
f^a(f^b(z)) - f^{a+b}(z) = 0\\
\]

And it equals this across all the taylor coefficients. So that:

\[
f^a(f^b(z)) = \sum_{j=0}^\infty c_j z^j\\
\]

\[
f^{a+b}(z) = \sum_{j=0}^\infty c_j z^j\\
\]

They're the same series...

Could you clarify what you mean?

EDIT:

Oh I just realized you meant digit precision! (At least I think you did)

I can usually cap out at about 200 digits using the beta method on the Shell-Thron region, it can vary a bit. Outside the Shell-thron region, the beta method isn't holomorphic (just C^\infty), so the taylor series doesn't converge, but I can still get about 200 digits accuracy for this functional equation. Using the standard Schroder method I can get desired accuracy, but it eats the hell out of my ram.

Using pari-gp basically lets you set precision, and if you program carefully it'll usually agree to the precision minus some odd digits. So If you want 300 digit precision, set about 350, and you should get the Taylor series converging to about 300. But then you should also have series precision really large...

EDIT2:

Oh so it means Daniel was referring to series precision of the Taylor series.

I'll answer again.

This is a bit more difficult using the beta method, it'll start to lag about 100 terms, but I can usually get about \(O(z^{100})\) and that's usually what I shoot for to double check everything. And to do that you have to set the series precision usually larger (some terms start turning to junk). With Schroder, it's again, desired accuracy. So technically you could get it to \(O(z^{1000})\) but be prepared to eat like 10 gbs of ram, lol.

I usually set my series precision to get at like 50-100, depending on what I'm running, then double check that it's continuing to work at something ridiculous like 200-300. But that's the final stages because it's super slow and eatttsss ram like the cookie monster.

[Daniel]

...
EDIT2:

Oh so it means Daniel was referring to series precision of the Taylor series.

I'll answer again.

This is a bit more difficult using the beta method, it'll start to lag about 100 terms, but I can usually get about \(O(z^{100})\) and that's usually what I shoot for to double check everything. And to do that you have to set the series precision usually larger (some terms start turning to junk). With Schroder, it's again, desired accuracy. So technically you could get it to \(O(z^{1000})\) but be prepared to eat like 10 gbs of ram, lol.

I usually set my series precision to get at like 50-100, depending on what I'm running, then double check that it's continuing to work at something ridiculous like 200-300. But that's the final stages because it's super slow and eatttsss ram like the cookie monster.

My precision and accuracy are guaranteed to be prefect because I only use symbolic and rational numbers in my calculations. Are you using floating point? Consider \( f(0)=0 \), then the second derivative of \( f^m(z) \) is always is true regardless of whether Schoeder or Abel is used.

For \( f^m(z) \) we have \( D^2 f^m(z)= \(\sum_{k=0}^{m-1}f'(0)^{2m-k-2}\)f''(0) \).

[JmsNxn]

...
EDIT2:

Oh so it means Daniel was referring to series precision of the Taylor series.

I'll answer again.

This is a bit more difficult using the beta method, it'll start to lag about 100 terms, but I can usually get about \(O(z^{100})\) and that's usually what I shoot for to double check everything. And to do that you have to set the series precision usually larger (some terms start turning to junk). With Schroder, it's again, desired accuracy. So technically you could get it to \(O(z^{1000})\) but be prepared to eat like 10 gbs of ram, lol.

I usually set my series precision to get at like 50-100, depending on what I'm running, then double check that it's continuing to work at something ridiculous like 200-300. But that's the final stages because it's super slow and eatttsss ram like the cookie monster.

My precision and accuracy are guaranteed to be prefect because I only use symbolic and rational numbers in my calculations. Are you using floating point?

No, I am not using floating point. I'm using pari, which is symbolic, but you'll still get some garbage which comes out in the 100th digit or so, it's just a consequence of any calculator. 

For example, if I take the Schroder function:

\[
\Psi(z)\\
\]

I can set arbitrary precision, and arbitrary series precision (how many terms to the Taylor series), but it will create some garbage data. This is also because I try to optimize the speed, I could technically make it absolutely 100% lossless, but who has time to watch that compile. Then, you have to series invert this, which is not lossless, despite how good Pari is, this is a polynomial inversion, so it can lose a bit of data. Then when you actually plug in the iteration:

\[
f^{\circ s}(z) = \Psi^{-1}(\lambda^s \Psi(z))\\
\]

You do lose a good amount of data, especially for large \(s\) (negative or positive depending on if \(|\lambda|\) is greater/less than 1) and \(z\), but you can recover it by keeping \(z \approx L\), and iterating the functional equation. But still, you will lose data. I could program my code to make sure it's 100% no errors, but, again, it'd be too slow and impractical.

[Daniel]

...
EDIT2:

Oh so it means Daniel was referring to series precision of the Taylor series.

I'll answer again.

This is a bit more difficult using the beta method, it'll start to lag about 100 terms, but I can usually get about \(O(z^{100})\) and that's usually what I shoot for to double check everything. And to do that you have to set the series precision usually larger (some terms start turning to junk). With Schroder, it's again, desired accuracy. So technically you could get it to \(O(z^{1000})\) but be prepared to eat like 10 gbs of ram, lol.

I usually set my series precision to get at like 50-100, depending on what I'm running, then double check that it's continuing to work at something ridiculous like 200-300. But that's the final stages because it's super slow and eatttsss ram like the cookie monster.

My precision and accuracy are guaranteed to be prefect because I only use symbolic and rational numbers in my calculations. Are you using floating point?

No, I am not using floating point. I'm using pari, which is symbolic, but you'll still get some garbage which comes out in the 100th digit or so, it's just a consequence of any calculator. 

For example, if I take the Schroder function:

\[
\Psi(z)\\
\]

I can set arbitrary precision, and arbitrary series precision (how many terms to the Taylor series), but it will create some garbage data. This is also because I try to optimize the speed, I could technically make it absolutely 100% lossless, but who has time to watch that compile. Then, you have to series invert this, which is not lossless, despite how good Pari is, this is a polynomial inversion, so it can lose a bit of data. Then when you actually plug in the iteration:

\[
f^{\circ s}(z) = \Psi^{-1}(\lambda^s \Psi(z))\\
\]

You do lose a good amount of data, especially for large \(s\) (negative or positive depending on if \(|\lambda|\) is greater/less than 1) and \(z\), but you can recover it by keeping \(z \approx L\), and iterating the functional equation. But still, you will lose data. I could program my code to make sure it's 100% no errors, but, again, it'd be too slow and impractical.

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

[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.

[Daniel]

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.

[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.

[Daniel]

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.
...

No, I don't use matrices, please clarify why you think so. The software is a proof of concept that total partitions are the combinatorial structure that comprises all iterated functions. I enumerate the total partitions, assigning an algebraic expression based on a recursive application of Faa Di Bruno's formula. The algebraic expressions are then assembled into the Taylors series of an iterated function.

[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.
...

No, I don't use matrices, please clarify why you think so. The software is a proof of concept that total partitions are the combinatorial structure that comprises all iterated functions. I enumerate the total partitions, assigning an algebraic expression based on a recursive application of Faa Di Bruno's formula. The algebraic expressions are then assembled into the Taylors series of an iterated function.

Oh I apologize. I didn't clarify.

I didn't mean matrices in the mathematics. I'm well aware you are using Faa Di Bruno's method of constructing iterations. What I meant is the manner you are programming it (whether inadvertently or not) is using a matrix algorithm. It's storing matrix solutions as \(Ax = b\). I don't mean your mathematics is doing that. I mean the code is doing that.

You are solving arrays of solutions, which give you the Taylor values. But then, you are asking to the program to solve \(100\times100\) or however values, of linear solutions. By matrix I solely meant in the manner of programming. Not the mathematical idea. I'm well aware Faa Di Bruno has nothing to do with matrices. But the program you are calling is certainly matrix-y in nature.

I apologize, I have an odd way of talking about things... I'm not trying to insult your program. But it is based off array/nested array calculations, which is "matrixy-y" in my book, lol.

EDIT:

You are looking at a vector of taylor coefficients \(c_j\), and multiplying it by \(Ac_j\)--and this is how you've encoded \(F(z)\) next to \(f(F(z))=F(z+1)\). Albeit how you've done it is very good. Honestly, I always though Faa Di Bruno was too hard to be practical. Hats off to you. But the code is clearly matrix based, and based off of solving matrix equations in the taylor coefficients.

This is very different than how I code. That's what I meant. Everything is a forloopy/looping/recursive polynomial into a polynomial.

I hope this makes sense. I just meant, I hate matrix shit, largely because I have an irrational phobia. Not to disannounce your code. But your code is matrix based.