A not-so-arbitrary solution is the following, which is in fact a polynomial approximation whose maximal possible order depends on the resources of your favorite software and computer-memory.

The idea is to use a matrixoperator/"Carleman-matrix" **B** for the generation of a powerseries (in fact a truncated one, so only a polynomial) for the exponential-function and then by the zero'th power **B^0** have the coefficients of the ID-function $f(x)=x $ , by the first power get the coefficients for $f(x)=exp(x)$, by the second power the coefficients of $f(x)=exp(exp(x))$ and so on. One may speak of iteration and iteration-height h here, where h is in the exponent of the matrix: **B^h** and in the iteration-parameter $f(x,h)=exp^{[h]}(x)$ . $h=0$ gives then the id-function, $h=1$ the exponential function.

Fractional powers of **B** for fractional iteration can then be approximated by fractional powers of **B**. Such fractional powers can be computed by matrix-diagonalization:

$ \qquad \qquad B = W*D*W^{-1} $ and $ \qquad \qquad B^h = W*D^h*W^{-1} $ where th diagonal-matrix **D** can raised to any power since we need only the h'th power of the scalar diagonal-entries.

So if you define the matrix **B** by $ B = \operatorname{matrix}_{(r=0 ... n-1),(c=0...n-1)}( c^r/r!) $ for some finite dimension *n*, say 16 or 32, then you get order 16 or -32 polynomials in two variables: $f_n(x,h)\approx \exp^{[h]}(x) $.

An implementation in Pari/GP is simple; however the required float precision is huge already for small *n*. Using *n=16* I get away with float-precision of 200 decimal digits, for *n=32* I need already 1200 or 1600 digits, *n=64* (which is usually my default order for such problems) was unreachable so far:

```
n=16
B = matrix(n,n,r,c,(c-1)^(r-1)/(r-1)!)
W = mateigen(B)
WI = W^-1
D = diag( WI * B * W ) \\ diag extracts the diagonal into a column-vector
bpow(h) = W * matdiagonal (vector(n,r,D[r]^h)) * WI
f(x,h) = Bh = bpow(h); return( sum(r=1,n, x^(r-1)*Bh[r,2]) )
```

Clearly, this code can/should be optimized. With *n=32* you get quite reasonable approximations for some range for, say $-0.5<h<2.5$, so you even get $f(x,h+1) \approx exp(f(x,h)) $ to precision of 10 or so digits.