Is there an "elegant" non-recursive formula for these coefficients? Also, how can one get proofs of these patterns?

Question

Not sure if this is a "good" question for this forum or if it'll get panned, but here goes anyway...

Consider this problem. I've been trying to find a formula to expand the "regular iteration" of "exp". Regular iteration is a special kind of complex function that is a solution of the equation

$$f(z+1) = \exp(f(z))$$

(or more generally for functions other than $\exp$. It is called "regular" because as a solution it is characterized by the fact the the functional iterates $F^t(z) = f(t + f^{-1}(z))$, with $F$ being the function that is $\exp$ in this case, are "regular", or analytic, at a chosen fixed point of $F$, for all non-integer $t$. There are regular iterations for every fixed point.)

This regular iteration in particular is an entire function. To get it, we take a fixed point $L$ of $\exp$ and expand a solution in powers of $L^z$. The result is to obtain a Fourier series

$$f(z) = \sum_{n=0}^{\infty} a_n L^{nz}$$

where

$$a_0 = L$$ $$a_1 = 1$$ $$a_n = \frac{B_n(1! a_1, 2! a_2, ..., (n-1)! a_{n-1}, 0)}{n!(L^{n-1} - 1)}$$

with $B_n$ being the nth "complete" Bell polynomial. This recursive formula yields the following expansions:

$$a_2 = \frac{1}{2L - 2}$$ $$a_3 = \frac{L + 2}{6L^3 - 6L^2 - 6L + 6}$$ $$a_4 = \frac{L^3 + 5L^2 + 6L + 6}{24L^6 - 24L^5 - 24L^4 + 24L^2 + 24L - 24}$$ $$a_5 = \frac{L^6 + 9L^5 + 24L^4 + 40L^3 + 46L^2 + 36L + 24}{120L^{10} - 120L^9 - 120L^8 + 240L^5 - 120L^2 - 120L + 120}$$ ...

It appears that, by pattern recognition and factoring the denominators,

$$a_n = \frac{\sum_{j=0}^{\frac{(n-1)(n-2)}{2}} mag_{n,j} L^j}{\prod_{j=2}^{n} j(L^{j-1} - 1)}$$

where $\mathrm{mag}_{n,j}$ is a sequence of "magic" numbers (integers) that looks like this (with the leftmost column being $j = 0$):

[update]: remark: for readability I transposed the original table. But I did not adapt the use of "columns" and "rows" in the surrounding text, so that must be translated in mind (Gottfried Helms)

n=1 n=2 n=3 n=4 n=5  n=6  n=7    n=8      n=9       n=10   ...
-------------------------------------------------------------- 
  1   1   2   6  24  120  720   5040    40320     362880   ...
          1   6  36  240 1800  15120   141120    1451520 
              5  46  390 3480  33600   352800    4021920 
              1  40  480 5250  58800   695520    8769600 
                 24  514 7028  91014  1204056   16664760 
                  9  416 8056 124250  1855728   28264320 
                  1  301 8252 155994  2640832   44216040 
                     160 7426 177220  3473156   64324680 
                      64 5979 186810  4277156   88189476 
                      14 4208 181076  4942428  114342744 
                       1 2542 163149  5395818  141184014 
                         1295 134665  5561296  166279080 
                          504 102745  5433412  187614312 
                          139  71070  5021790  202901634 
                           20  44605  4391304  210825718 
                            1  24550  3625896  210403826 
                               11712  2820686  201934358 
                                4543  2056845  186191430 
                                1344  1398299  164980407 
                                 265   879339  140216446 
                                  27   504762  114231817 
                                   1   260613   88934355 
                                       117748   66047166 
                                        45178   46576620 
                                        13845   31071602 
                                         3156   19460271 
                                          461   11365652 
                                           35    6112650 
                                            1    2987358 
                                                 1298181 
                                                  488878 
                                                  153094 
                                                   37692 
                                                    6705 
                                                     749 
                                                      44 
                                                       1 

But what is the simplest (or at least "reasonably" simple) non-recursive formula for these numbers, or perhaps the numerators in general? Like a sum formula, or something like that. Is there some kind of "combinatorical"-like formula here (sums/products, perhaps nested, of factorials and powers and stuff like that, binomial coefficients, special numbers, etc.)? I notice that the first column is factorials... (how can one prove that?)

And regardless of the formula for the "mag", can one prove from the recurrence formula that the $a_n$ have the form given, and if so, how? Especially, how can one prove the numerator has degree $\frac{(n-1)(n-2)}{2}$? Perhaps that might provide insight into how to find the formula for the "mag".

The ultimate goal here is to try and obtain a series expansion for the "tetration" function $^z e$, more specifically, Kneser's tetrational function, described in Kneser's papers on solutions of $f(f(x)) = \exp(x)$ and related equations (paper is in German, I only saw the translations.). Though this may not be the best way to go, since after constructing this regular iteration function, we then need a special mapping derived from a Riemann mapping to "distort" it so it becomes real-valued at the real axis, and I don't know if there's any good way to construct Riemann mappings even as "non-closed" infinite expansions. But I'm still curious to see if at least a formula for this function is possible.

EDIT: Oh, and for all its worth, apparently

$$\sum_{j=0}^{\frac{(n-1)(n-2)}{2}} \mathrm{mag}_{n,j} = \frac{n!(n-1)!}{2^{n-1}}$$

if that helps any (don't see how it would, and this is not proven, I just got it by looking up the sums on the integer sequences dictionary site.). Perhaps maybe some hints as to why it has that value could help in finding the formula, though...

---

Justification for thinking a formula exists

Why do I think this even exists, when there's no guarantee that this kind of really non-trivial recurrence relation should even have a non-recursive solution in the first place? Well, for one, the fact that so much of it could be put in simple form as given, and also I did manage to come up with an explicit formula from a very roundabout way but this formula is excessively complicated and based on very general techniques.

It is difficult to describe that formula here, but the outline of the process to construct it is this, for all its worth:

1. A general recurrence of the form

$$A_1 = r_{1, 1}$$ $$A_n = \sum_{m=1}^{n-1} r_{n,m} A_m$$

has a non-recursive solution formula. This I found myself, but it is hideous and involves binary bit operations. This kind of recurrence is very general, and it also includes the recurrence for the Bernoulli numbers and other kinds of recurrences.

2. The "regular Schroder function" of $F(z) = e^{uz} - 1$, i.e. the function satisfying $\mathrm{RSF}(F(z)) = K \mathrm{RSF}(z)$ (sometimes called the Schroder functional equation, hence the name) which is "regular" in that it can be turned into the regular iteration of $F$ (as we do next), can be given as a Taylor series

$$\mathrm{RSF}(z) = \sum_{n=1}^{\infty} A_n z^n$$

where $A_n$ is given by the recurrence-solving formula with $r_{1,1} = 1$ and $r_{n, m} = \frac{u^{n-1}}{1 - u^{n-1}} \frac{m!}{n!} S(n, m)$ (here, $S(n, m)$ is a Stirling number of the 2nd kind). This is hideous, involving lots of "binary bit manipulation" stuff such as counting 1 bits and positions of 1 bits, which have not-so-nice formulas (the latter involves a set indicator function, at least in the formulation I found myself...). Not sure at all how this could be simplified. The formulas just don't seem to lend themselves to simplification, at least not any that I know of.

3. Invert the regular Schroder function using the Lagrange inversion theorem. This can be expanded in an explicit "non-recursive" form, but it needs so-called "potential polynomials" and other complexity. Plug the huge $A_n$ formula into this. Horrific!

4. Now $U(z) = \mathrm{RSF}^{-1}(u^z)$ is a "regular iteration" of $e^{uz} - 1$, giveable as a Fourier series, or Taylor series in $u^z$.

5. Apply the topological conjugation to conjugate it to iteration of $e^{vz}$ by taking $v = ue^{-u}$ thus $u = -W(-v)$ (Lambert's W-function). Take $H(z) = e^{-u} z - 1$ then find $H^{-1} o U o H$. This gives a regular iteration of $e^{vz}$, thus set $v = 1$ ($u = -W(-1) = \mathrm{fixed\ point\ of\ exponential}$). Though, there may be a constant displacement of some kind offsetting this regular from the one given by the $a_n$-formula. EDIT: Oops!!!! That should be $H^{-1}(U(U^{-1}(H(U(0))) + z))$, but wait, that's just a constant-shift of $H^{-1} o U$, so just take $H^{-1} o U$ as the regular iteration of $e^{vz}$, probably displaced (in $z$) from the one we're trying to solve for by a constant, but should be structurally identical (and you can try and compute $U^{-1}(H(U(0)))$. Perhaps that is the shift required, but I don't know.).

(EDIT: Apparently the step-numbering above isn't working right for some reason.)

So by this, I think an explicit formula exists (though that constant-shift at the end may be a little problem, but not much, since it is immaterial to the structure of the function). I'm just interested in something simpler than this, preferably something to "fill out" the "mag" formula I gave...

EDIT: Now I'm pratically sure explicit non-recursive solution is possible. Using some numerical tests, I figured the constant shift should be (for $v = 1$, i.e. $u = L$) simply -1, that is, take $H^{-1}(U(z - 1))$ and the coefficients of the Fourier expansion will be equal to $a_n$ in explicit non-recursive form (but atrocious, hence my question, to find something more elegant. This at least evidences that an explicit non-recursive solution is possible, addressing any skeptics' concerns that it isn't and so an elegant one wouldn't exist either. And it is a good bet that if an atrocious formula exists derived from very general principles (note Step 1 above), there may be a more elegant one derived from more specific principles.). So, almost a proof. It could probably be turned into one with a little more work, though that would be much too long to post here.

---

Answer 1

Let $\beta_n$ denote the flag $h$-vector (as defined in EC1, Section 3.13) of the partition lattice $\Pi_n$ (EC1, Example 3.10.4). Then $$ \mathrm{mag}_{n,{n-1\choose 2}-j} = \sum_S \beta_n(S), $$ where $S$ ranges over all subsets of $\lbrace 1,2,\dots,n-2\rbrace$ whose elements sum to $j$. An explicit formula for $\beta_n(S)$ is given by $$ \beta_n(S) = \sum_{T\subseteq S} (-1)^{|S-T|} \alpha_n(T), $$ where if the elements of $T$ are $t_1<\cdots < t_k$, then $$ \alpha_n(T) = S(n,n-t_1)S(n-t_1,n-t_2) S(n-t_2,n-t_3)\cdots S(n-t_{k-1},n-t_k). $$ Here $S(m,j)$ denotes a Stirling number of the second kind.

Addendum. A combinatorial description of the mag numbers is somewhat complicated. Consider all ways to start with the $n$ sets $\lbrace 1 \rbrace,\dots, \lbrace n \rbrace$. At each step we take two of our sets and replace them by their union. After $n-1$ steps we will have the single set $\lbrace 1,2,\dots,n \rbrace$. An example for $n=6$ is (writing a set like $\lbrace 2,3,5\rbrace$ as 235) 1-2-3-4-5-6, 1-2-36-4-5, 14-36-2-5, 14-356-2, 14356-2, 123456. At the $i$th step suppose we take the union of two sets $S$ and $T$. Let $a_i$ be the least integer $j$ such that $j$ belongs to one of the sets $S$ or $T$, and some number smaller than $j$ belongs to the other set. For the example above we get $(a_1,\dots,a_5)=(6,4,5,3,2)$. If $\nu$ denotes this merging process, then let $f(\nu) = \sum i$, summed over all $i$ for which $a_i>a_{i+1}$. For the above example, $f(\nu) = 1+3+4=8$. (The number $f(\nu)$ is called the major index of the sequence $(a_1,\dots,a_{n-1})$.) Then $\mathrm{mag}_{n,{n-1\choose 2}-j}$ is the number of merging processes $\nu$ for which $f(\nu)=j$. This might look completely contrived to the uninitiated, but it is very natural within the theory of flag $h$-vectors.

Answer 2

Mike, possibly you know this all, but for my own pleasure (and for the reader who's not yet much familiar with this all) I can give an explicite description for the mag-numbers. However, this is simply an explication of the terms of some series/arrays which are in fact recursive, but the recursion is so flat that we can resolve it without too much hazzle to direct references on factorials/binomials and the log of the fixpoint only. I employ the notation of (operator-)matrices which are known as Bell-/Carleman-matrices.

The text became too long for the answer-field here, so I link to a pdf-file on my homepage.
(If you don't like pdf there is also a html-version, but automatically generated by word and not perfect formatted)

Since I describe this using the known procedure of diagonalization of a triangular matrix some of your questions concerning the structure of coefficients may be answered or possibly a rigorous answer lays on the hand.

(P.s.: if it is more convenient for the MO-readers I could upload or possibly reformat the text for mathjax, but the latter would be much unwanted work...)

[update]: updated the pdf-file for readability

Answer 3

Bell polynomials

The first step to solving the problem is to deal with the recursion combinatorially. Note the Bell polynomial in $$a_n = \frac{B_n(1! a_1, 2! a_2, ..., (n-1)! a_{n-1}, 0)}{n!(L^{n-1} - 1)}.$$

Bell polynomials are generalizations of the combinatorial structure set partitions or Bell numbers. The Faà di Bruno's formula in terms of the Bell polynomials provides a way to express the Taylor's series of a composite function $f(g(x))$.

From Wikipedia,

$${d^n \over dx^n} f(g(x)) = \sum_{k=1}^n f^{(k)}(g(x))\cdot B_{n,k}\left(g'(x),g''(x),\dots,g^{(n-k+1)}(x)\right).$$

The formula has a "combinatorial" form:

$${d^n \over dx^n} f(g(x))=(f\circ g)^{(n)}(x)=\sum_{\pi\in\Pi} f^{(\left|\pi\right|)}(g(x))\cdot\prod_{B\in\pi}g^{(\left|B\right|)}(x)$$

where

π runs through the set Π of all partitions of the set { 1, ..., ''n'' },

$B ∈ π$ means the variable $B$ runs through the list of all of the "blocks" of the partition π, and

$|A|$ denotes the cardinality of the set $A$ (so that $|π|$ is the number of blocks in the partition $π$ and $|B|$ is the size of the block $B$).

The appearance of Bell polynomials in this problem is not surprising as the problem is an example of fractional iteration where $f(g(x))=f(f^{m-1}(x)).$

Combinatorial Structures

Moving the recursion from the analytic to the combinatoric, given that the set partitions are the combinatorial form of the derivatives of composite functions, what is the combinatorial form associated with the derivatives of $D(f(g(x)))=D^nf^m(x)$, iterated functions? Because there are now an arbitrary number of levels of recursion, there are also an arbitrary number of recursive set partitionings. This is equivalent to the combinatorial problem of the total partitions of $n$ or Schroeder's Fourth Problem.

Analytic

Without loss of generalization assume that $f(0)=g(0)=0$. The following proof removes with restriction that $a_1=1$. It will be seen that this approach handles both the cases of Schroeder's equation and Abel's equation.

The First Derivative

The first derivative of a function at its fixed point $Df(0)=f_1$ is often represented by $\lambda$ and referred to as the multiplier or the Lyapunov characteristic number; its logarithm is known as the Lyapunov exponent. Let $g(z)=f^{m-1}(z)$, then $$Df(g(z))=f'(g(z))g'(z)$$ $$=f'(f^{m-1}(z))Df^{m-1}(z)$$ $$=\prod^{m-1}_{k_1=0}f'(f^{m-k_1-1}(z))$$
$$Df^m(0)=f'(0)^m$$ $$=f_1^m = \lambda^m$$

The Second Derivative

$$D^2f(g(z))=f''(g(z))g'(z)^2+f'(g(z))g''(z)$$ $$=f''(f^{m-1}(z))(Df^{m-1}(z))^2+f'(f^{m-1}(z))D^2f^{m-1}(z)$$

Setting $g(z) = f^{m-1}(z)$ results in $$D^2f^m(0)= f_2 \lambda^{2m-2}+\lambda D^2f^{m-1}(0)$$ When $\lambda \neq 0$, a recurrence equation is formed that is solved as a summation. Note that $|\lambda|\neq1$ is consistent with Schroder's equation, while $\lambda=1$ is consistent with Abel's equation. $$ D^2f^m(0)=f_2\lambda^{2m-2}+\lambda D^2f^{m-1}(0)$$ $$ =\lambda^0f_2 \lambda^{2m-2} +\lambda^1f_2 \lambda^{2m-4} +\cdots +\lambda^{m-2}f_2 \lambda^2 +\lambda^{m-1}f_2 \lambda^0 $$ $$ =f_2\sum_{k_1=0}^{m-1}\lambda^{2m-k_1-2} $$

The Third Derivative

Continuing on with the third derivative,

$$ D^3f(g(z))=f'''(g(z))g'(z)^3+3f''(g(z))g'(z)g''(z)+f'(g(z))g'''(z) =f'''(f^{m-1}(z))(Df^{m-1}(z))^3 +3f''(f^{m-1}(z))Df^{m-1}(z)D^2f^{m-1}(z) +f'(f^{m-1}(z))D^3f^{m-1}(z) $$

$$ D^3f^m(0)=f_3\lambda^{3m-3}+3 f_2^2\sum_{k_1=0}^{m-1}\lambda^{3m-k_1-5}+\lambda D^3f^{m-1}(0)$$ $$ =f_3\sum_{k_1=0}^{m-1}\lambda^{3m-2k_1-3} +3f_2^2 \sum_{k_1=0}^{m-1} \sum_{k_2=0}^{m-k_1-2} \lambda^{3m-2k_1-k_2-5} $$ Note that the index $k_1$ from the second derivative is renamed $k_2$ in the final summation of the third derivative. A certain amount of renumbering is unavoidable in order to use a simple index scheme.

The Higher Derivatives

The fourth derivative is comprised of five summations.

$D^n f^m(0) = \sum \frac{n!(D^k f)(0)}{k_1! \cdots k_{n-1}!} $ $\left(\frac{Df^{m-1}(0)}{1!}\right)^{k_1} \cdots $ $\left(\frac{D^n f^{m-1}(0)}{(n-1)!}\right)^{k_{n-1}} $ $ + \lambda D^n f^{m-1}(0)$

Computation

SchroederSummations.nb and Iterate.m are Mathematica notebooks which are an attempt to systematically generate the combinatorial structure of unlabeled total partitions and then compute the analytic form of each combinatorial structure. Tallying the integer coefficients of the derivatives gives the sequence 1,1,4,26,236,2752,39208,660032,$\ldots$ which is the structure total partitions. The number of summations is given by the structure unlabeled total partitions; 1, 1, 2, 5, 12, 33, 90, 261,$\ldots$.

So the eighth derivative with 660032 terms, 261 with symmetry, can be directly computed without the lesser derivatives. Faà di Bruno's formula is indexed by the set partitions, but the derivatives of an iterated function are indexed by total partitions. Each total partition maps into a summation.

Note: I did spend some time trying to find a combinatorial structure associated with $a_n$. I tried setting L to a few reasonable values; -1, -1/2, 1/2 and 1 to see if I could generates a combinatorial sequence in the OEIS, but I didn't find any matching sequences.

Answer 4

Here I'll show how the MAG-numbers of n=4 are computed. I set up a system of linear equations and solve.
I've simply extended this to n=5 but either I must have a bug in my data or there is a more severe problem which I might have overlooked - the rank of the system of linear equations is defective by 2 ranks. I'm still bug-tracking but I'm confident I'll get it working.

Background For some base $t$, with $\small 1 \lt t \lt e$ for the iterable function $\small f(x)=t^x-1$ let's write its logarithm $\small u=\log(t)$ .

• (for the following see the more detailed article APT-decomposition so far; I'll include the essentials here)

The Carleman-matrix F for $\small f(x)$ is lower triangular with powers of $u$ in the diagonal and has coefficients from the Stirling numbers 2nd kind, similarity scaled by the factorials, so that in the second column the entries gived the coefficients for the powerseries of $\small f(x)$ Here is the top left of that matrix $$ F=\begin{bmatrix} \Tiny \begin{array}{} 1 & . & . & . & . & . \\ 0 & u & . & . & . & . \\ 0 & 1/2 u^2 & u^2 & . & . & . \\ 0 & 1/6 u^3 & u^3 & u^3 & . & . \\ 0 & 1/24 u^4 & 7/12 u^4 & 3/2 u^4 & u^4 & . \\ 0 & 1/120 u^5 & 1/4 u^5 & 5/4 u^5 & 2 u^5 & u^5 & \cdots \\ \vdots &\vdots &\vdots &\vdots &\vdots &\vdots & \ddots \end{array} \end{bmatrix} $$

For fractional iterations we need fractional powers of that matrix. We can do this by diagonalization (which can even be done keeping $u$ symbolical)
The diagonalization gives three matrices such that $\small F = M \cdot D \cdot W$ where also $\small W=M^{-1}$. The columns of $M$ (being the matrix of eigenvectors) can arbitrarily be scaled; I usually norm that such that the diagonal has all units. Then $M$ is also the Carlemanmatrix for the Schröder-function for $f(x)$ - that means, it has the coefficients of its formal power series in its second column. That coefficients can be decoded into polynomials in $u$ in the numerator and factorials and other polynomials in $u$ in the denominator and display precisely your formula involving the MAG-numbers.

$$ M[,1]=\small \begin{bmatrix} \Tiny \begin{array} {rll} 0 & & \\ 1 & \cdot u/1! & /u \\ 1 & \cdot u^2/2! & /u(1-u) \\ 2u+1 & \cdot u^3/3! & /u(1-u)(1-u^2) \\ 6u^3+5u^2+6u+1 & \cdot u^4/4! & /u(1-u)(1-u^2)(1-u^3) \\ 24u^6+26u^5+46u^4+45u^3+24u^2+14u+1 & \cdot u^5/5! & /u(1-u)(1-u^2)(1-u^3)(1-u^4)\\ \vdots \end{array}\end{bmatrix}$$

The inverse $W = M^{-1}$ , which is the Carlemanmatrix for the inverse Schröder-function, has the numbers which you are after:

$$ W[,1]=\small \begin{bmatrix} \Tiny \begin{array} {rll} 0 & & \\ 1 & \cdot u/1! & /u \\ -1 & \cdot u^2/2! & /u(1-u) \\ u+2 & \cdot u^3/3! & /u(1-u)(1-u^2) \\ -1u^3-5u^2-6u-6 & \cdot u^4/4! & /u(1-u)(1-u^2)(1-u^3) \\ 1u^6+9u^5+24u^4+40u^3+46u^2+36u+24 & \cdot u^5/5! & /u(1-u)(1-u^2)(1-u^3)(1-u^4)\\ \vdots \end{array}\end{bmatrix}$$

(I hope I did not make sign errors here, the reproduction from the Pari/GP-output is a bit messy)

The problem is, that the diagonalization is a recursive procedure, so the coefficients are defined by recursion which is what you did not want.
In my linked essay on the full symbolic decomposition I proceeded to express from $F^h = M \cdot D^h \cdot W$ (which is the Carlemanmatrix for the h'th iterate $f^{\circ h}(x)$) the symbolic representation of that fractional iterate.
Key-topic: I observed, that the Taylor-coefficients of $f^{\circ h}(x)$ are expressible as bivariate polynomials in $k$ (the index of the coefficient) and in $u$ and (separately!) $u^h$ as arguments and I displayed the coefficients of that polynomials $a_k(u,u^h)$ as matrices $A_k$ such that $$a_k(u,u^h) = V(u)\cdot A_k \cdot V^\tau(u^h) $$
$ \qquad \qquad $ (where $V(u)$ means $\text{rowvector}(1,u,u^2,u^3,...)$ up to the appropriate dimension).

For instance the matrix $A_4$

$ \qquad \qquad $

and the matrix $A_5$

$ \qquad \qquad $

$ \qquad \qquad $ $ \qquad \qquad $ The related MAG-coefficients for $n=4$ and $n=5$ are in the last columns

The denominators are the same as in the coefficents for the Schröder-function, shown in the picture of $M[,1]$ and for the inverse Schröderfunction in the picture of $W[,1]$.

The key for the possibility to compute the MAG-coefficients directly is now, that for each integer $h$ the columns are down- or upshifted and the full evaluation of the polynomial at $u^h$ must be a multiple of the denominator.
From this a system of linear equations can be set, determined by such integer multiplies of numerator and denominator-polynomials.
As starting values we find, that in the first and last rows are Stirlingnumbers of first and second kind (the latter scaled by factorials) and the rest of numbers are unknown and shall be found by solving that system of linear equations. Apriori- knowledge is (well, hypothese) that for $h=0$ the complete expression vanishes and for $h=1$ the numerator must exactly be equal to the denominator. This seems to make the problem solvable.

end background

---

[update 2.Mar] When I computed the following I'd misread your question in that I misread that you wanted the coefficients of the Schröderfunction instead of the inverse. The logic for the computation should not be affected, so I leave the following for the time being [end update]

Here is the solution for $k=n=4$.
As stated above, the matrix $A_4$ provides the coefficients for the polynomial in the numerator of the term at $x^4$ (in the Taylorseries for $f^{\circ h}(x)$) . Actually written as a polynomials in $u$ and $u^h$ the first few such polynomials look (with their matrix-arrangement kept):

The red coefficients are unknown when we want to avoid the recursive procedure. To make it more explicit this is the polynomial with the symbolic names for the coefficients where we know that of the first and last row:

The denominator $D_4(u) $ of the 4'th term is the product

Now my hypothese, from which I construct a system of linear equations is, that for the integer iterates the numerator is a polynomial multiple of the denominator - which, btw., explains then that the limit for $u \to 1$ exists for integer iteration heights, but that is doesn't exist for fractional heights (this is completely equivalent to the q-factorials which exist too if $q$ goes to 1 by the limit consideration).

For $h=0, u^h=1$ we have $A_4(u,1)= 0$ and from the hypothese we must have
$$\begin{array} {rrl} & A_4(u,1) &= z \cdot D_4(u) \\ & 0 &= z \cdot D_4(u) \\ \to & z&=0 \\ \end{array} $$ and also follow some more equations by the observation that this means that the rowsums in the matrix must equal zero.

For $h=1, u^h=u$ we have $A_4(u,u)= ?? $ with the highest exponent equalling the highest exponent in $D_4(u)$ and from the hypothese we must have
$$\begin{array} {rrl} & A_4(u,u) &= y \cdot D_4(u) \\ &&\text{also we "know", that }a_1=-1,d_4=1 \\ \to & y&=1 \\ \end{array} $$ From this follow again some more equations for the unknowns.

For $h>2$ (and higher $h$) the polynomial $A_h(u,u^h)$ has higher order, so -if the hypothesis holds- it must be a polynomial multiple of $D_4(u)$ and actually we have $$ A_4(u,u^2) = (x_1+x_2\cdot u + x_3 \cdot u^2 + x_4 \cdot u^3)\cdot D_4(u)$$ This gives 4 new unknowns (where however the first and last are "trivial") and again a new subset of linear equations.

I constructed the following system of linear equations:

(where I marked the known coefficients with green color and put the 2x4 systematically known coefficients to the rhs of the system to solve it for the unknowns $a_2,...,e_3$ and $x_1...x_4, w_1...w_7$)
Gaussian elimination provided the sought coefficients $d_1,d_2,d_3,d_4$ of the inverse Schröder-function which are the "MAG$\,_4$" numbers (and interestingly $a_1,a_2,a_3,a_4$ that of the non-inverse Schröderfunction).

So this solves your problem for the fourth coefficients $n=4$ in your post.

This all looks very suggestive to simply proceed to higher $n$, but currently, with $n=5$ (meaning $A_5(u.u^h)$ ) I get the system with defective rank ( 2 missing determinations) and I've a catch that if the $y$ to $u$ informations are used for the Gaussian algorithm then a set of $8 $ coefficients of $A_5()$ form a (minimal) collinear subset.
I do not yet know how to resolve this and am studying the reasons for it. If I cannot resolve the problem, this method becomes possibly useless (which were a pity...)

---

Answer 5

This is a comment not an answer. Here are a few apparent patterns:

$$\begin{align}mag_{n,0} &= (n-1)! \\\ mag_{n,1} &= (n-1)! \frac{n-2}{2} \\\ mag_{n,{n-1 \choose 2}} &= 1 \\\ mag_{n,{n-1 \choose 2}-1} &= {n\choose2}-1 \\\ mag_{n,{n-1 \choose 2}-2} &= \frac{(n-2)(n-1)n(3n-5)}{24}-1\end{align}$$

Answer 6

added: there is some resemblance to elements in Richard Stanley's answer (which was also accepted), perhaps this here is just reinventing the wheel...
added: Pari/GP-code, nonrecursive, see bottom of post

I've seemingly found a good pattern heuristically for the direct construction of the mag-numbers. I tried and it works for the indexes up to $a_6$.
The idea for this stems from my earlier linked pdf-file, see here and was extended to $a_5$ and $a_6$ but only by guessing of symmetries and trying of patterns so far.

Notation is for Pari/GP,

Legend:

• 1) what you called $L$ I do call $u$ (it's by my tradition, please forgive that)
• 2) s2(row,col) is in fact $\text{Stirling}_2[r,c] \cdot c!/r!$

Top left of $s2(r,c)$ given as matrix, index zero-based

  1      .     .    .  .  .
  0      1     .    .  .  .
  0    1/2     1    .  .  .
  0    1/6     1    1  .  .
  0   1/24  7/12  3/2  1  .
  0  1/120   1/4  5/4  2  1

---

Mag-Coefficients at $a_3$ :

{3!*(    s2(3,1)                   *(u-1)
        +s2(3,2)*s2(2,1)           
     ) }

 %91 = u + 2  \\ coefficients found
------------------------------------
       1   2  \\ mag-nr expected

---

Mag-Coefficients at $a_4$ :

{4!*(    s2(4,1)                   *(u-1)*(u^2-1)
        +s2(4,3)*s2(3,1)                 *(u^2-1)
        +s2(4,2)        *s2(2,1)   *(u-1)
        +s2(4,3)*s2(3,2)*s2(2,1)
     ) }
     %82 = u^3 + 5*u^2 + 6*u + 6       \\ found coefficients         
    -----------------------------------------------------------     
          1      5       6     6       \\ mag-nr expected

---

Mag-Coefficients at $a_5$:

{5!*( s2(5,1)                          *(u-1)*(u^2-1)*(u^3-1)

     +s2(5,4)*s2(4,1)                        *(u^2-1)*(u^3-1)
     +s2(5,3)        *s2(3,1)          *(u-1)        *(u^3-1)        
     +s2(5,2)                *s2(2,1)  *(u-1)*(u^2-1)

     +s2(5,4)*s2(4,3)*s2(3,1)                        *(u^3-1)
     +s2(5,4)*s2(4,2)        *s2(2,1)        *(u^2-1)
     +s2(5,3)        *s2(3,2)*s2(2,1)  *(u-1)

     +s2(5,4)*s2(4,3)*s2(3,2)*s2(2,1)        
   )   }
 %83 = u^6 + 9*u^5 + 24*u^4 + 40*u^3 + 46*u^2 + 36*u + 24   \\ found coefficients
--------------------------------------------------------------------------------
      1      9       24       40       46       36     24   \\ mag-nr expected

---

Mag-Coefficients at $a_6$:

{6!*(    s2(6,1)                                     *(u-1)*(u^2-1)*(u^3-1)*(u^4-1)

        +s2(6,5)*s2(5,1)                                   *(u^2-1)*(u^3-1)*(u^4-1)
        +s2(6,4)        *s2(4,1)                     *(u-1)        *(u^3-1)*(u^4-1)        
        +s2(6,3)                *s2(3,1)             *(u-1)*(u^2-1)        *(u^4-1)
        +s2(6,2)                         *s2(2,1)    *(u-1)*(u^2-1)*(u^3-1)

        +s2(6,5)*s2(5,4)*s2(4,1)                                   *(u^3-1)*(u^4-1)
        +s2(6,5)*s2(5,3)        *s2(3,1)                   *(u^2-1)        *(u^4-1)
        +s2(6,5)*s2(5,2)                *s2(2,1)           *(u^2-1)*(u^3-1)
        +s2(6,4)        *s2(4,3)*s2(3,1)             *(u-1)                *(u^4-1)
        +s2(6,4)        *s2(4,2)        *s2(2,1)     *(u-1)        *(u^3-1)      
        +s2(6,3)                *s2(3,2)*s2(2,1)     *(u-1)*(u^2-1)

        +s2(6,5)*s2(5,4)*s2(4,3)*s2(3,1)                                   *(u^4-1)
        +s2(6,5)*s2(5,4)*s2(4,2)        *s2(2,1)                    *(u^3-1)
        +s2(6,5)*s2(5,3)        *s2(3,2)*s2(2,1)            *(u^2-1)
        +s2(6,4)        *s2(4,3)*s2(3,2)*s2(2,1)     *(u-1)

        +s2(6,5)*s2(5,4)*s2(4,3)*s2(3,2)*s2(2,1)        
   )   }   
 %85 = u^10 + 14*u^9 + 64*u^8 + 160*u^7 + 301*u^6 + 416*u^5 + 514*u^4 + 480*u^3 + 390*u^2 + 240*u + 120  \\ found coefficients
--------------------------------------------------------------------------------------------------------
      1       14       64       160       301       416       514       480       390       240     120    \\ mag-nr expected 

---

added: A Pari/GP-routine to compute the mag-coefficients for any $a_n$ is the following. Note that I use idx for $n$ and have the factorially scaled Stirlingnumbers precomputed in a matrix fS2F of sufficient dimension, see top of the post. I also optimized/simplified the above representation of the mag-coefficients-compositions by some reformulations which became obvious in the second and third view.
[note: s2()is used by s3() is used by s4r() is used by s5() . Everything is non-recursively in for-loops]. beware of high indexes $n=$idx For n=17 it needed 28 secs, growing exponentially with higher $n$

s2(r,c)=fS2F[1+r,1+c] \\ from precomputed matrix fS2F of stirlingnumbers 2nd kind, s2(4,3) = S2[4,3] *3!/4! 

s3(r,c,idx)=local(x,denom);if(r==idx,denom=1,denom=(u^(idx-r)-1)); return( s2(r,c)/denom );

s4r(vk,idx) = local(pr,l1=#vk+2-1); vk=concat([1],concat(vk,[idx])); return(prod(k=1,l1,s3(vk[k+1],vk[k],idx)));

{s5(idx)=local(su,vk,l1=idx-1,l2=idx-2);
   su=s4r([],idx); \\ get s3([idx,1])=s2(idx,1) as first summand
   for(k=1,l2, 
                  vk=vector(k,c,[2,l1]);
         forvec(vk_x = vk, su+=s4r(vk_x,idx),2); \\ "2" meaning "strict ascending seq of ranges"
      );
     return( su * idx! * prod(k=1,l2,(u^k-1))); }

try it:

s5(4)   
 %148 = u^3 + 5*u^2 + 6*u + 6

s5(5)   
 %149 = u^6 + 9*u^5 + 24*u^4 + 40*u^3 + 46*u^2 + 36*u + 24

s5(6)   
 %150 = u^10 + 14*u^9 + 64*u^8 + 160*u^7 + 301*u^6 + 416*u^5 + 514*u^4 + 480*u^3 + 390*u^2 + 240*u + 120