Question

Or rather, what function can be parametrized with some value t in [0,1] such that f(x, t= 0) = x, f(x, t = 1) = e^x, and f(x, 0 < t < 1) is a principled interpolation between those two, kind of like the gamma function is a principled interpolation for the discrete factorial.

Obviously, many functions fit the bill, like f(x,t) = x^(1 + t*(x/ln(x)-1) ), but that seems kind of arbitrary.

What mathematically useful/elegant function exists?

Answer 1

Igor's answer is probably better but I'll observe that we want to interpolate between $e^{\ln{x}}$ at $t=0$ and $e^x$ at $t=1.$ Furthermore $\int x^{-1}dx=\ln(x)+C$ and $\int x^0 dx=x+C.$

Rearranging a bit, since $\displaystyle \lim_{t \rightarrow 0^{+}}\frac{x^t-1}{t}=\ln{x},$ I am tempted by $$e^{\frac{x^t-1}{t}}=\sqrt[t]{e^{x^t-1}}.$$ That does not quite work at $t=1$ since it becomes $e^{x-1}.$ However it is also true that $\ \displaystyle \lim_{t \rightarrow 0^{+}}\frac{(x+t)^t-1}{t}=\ln{x}$ and that works at both ends as do $\frac{x^t-1+t^2}{t}\ $ and $\frac{x^t-1+t^{2-t}}{t}.$

Answer 2

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