While trying to understand a paper of Cayley, he left something unexplained, I managed to show that it is equivalent to the following formula, which I got stuck at:
$$k \cdot (f^k)^{(k-1)} = \sum_{j=0}^{k-1} {{k} \choose {j}} (f^{j})^{(j)}(f^{k-j})^{(k-j-1)}. $$
Can somebody help?
A detailed historical discussion of identities like this one can be found in Warren P. Johnson's paper The Pfaff/Cauchy derivative identities and Hurwitz type extensions, The Ramanujan Journal 13 (2007) pp. 167–201. In particular, his formula (1.3) is $$\frac{d^n\ }{dx^n} \phi^n(x) u(x) v(x) = \sum_{k=0}^n \binom nk \left(\frac{d^k\ }{dx^k} \phi^k(x) u(x)\right) \left(\frac{d^{n-k-1}}{dx^{n-k-1}}\phi^{n-k}(x) v'(x)\right)\tag{1}\label{1}$$ where for $n=k$, $\frac{d^{n-k-1}}{dx^{n-k-1}}\phi^{n-k}(x) v'(x)$ is to be interpreted as $v(x)$. This identity was given by Cauchy in 1826, but is equivalent to an identity given by Pfaff in 1795.
The OP's identity is equivalent to the case $\phi(x)=f$, $u(x)=1$, $v(x) = x$.
There is a related formula in Cayley's paper On the partitions of a polygon, Coll. Math. Papers 13 (1897), 93–113; Proc. London Math. Soc. (1) 22 (1891), 237–262, but I didn't see this formula there.
Formulas like \eqref{1} are discussed in my survey paper Lagrange Inversion, Journal of Combinatorial Theory, Series A 144 (2016), pp. 212–249, section 2.6.
There is also a discussion of a somewhat related formula on Terry Tao's blog at Another problem about power series.
Added November 12, 2022: A new paper with multivariable generalizations of these identities is Wenchang Chu, Multiple derivative inversions and Lagrange-Good expansion formulae (DOI is being registered), Mathematics 10 (2022), 4234.
In addition to Ira Gessel's answer, let me explain how such Pfaff–Cauchy–Hurwitz type identities are connected to polynomial interpolation techniques and to Raney lemma.
Denote $\mathbf k=\{1,\dotsc,k\}$.
We prove a more general (not polynomial in $f$, but multilinear in $f_1,\dotsc,f_k$) identity $$ k(f_1\dotsm f_k)^{(k-1)}=\sum_{\emptyset\ne A\subset\mathbf k}\left(\prod_{i\in A} f_i\right)^{(\lvert A\rvert-1)}\left(\prod_{i\notin A} f_i\right)^{(n-\lvert A\rvert)}. $$ It suffices to prove this identity when all $f_i$'s are power functions $f_i(t)=t^{x_i}$ (by some standard abstract nonsense). In this case we get a polynomial identity in $x_i's$ (where we denote $(x)_n=x(x-1)\dots(x-n+1)$, and also denote $s(A)=\sum_{i\in A}x_i$): $$ k(s(\mathbf k))_{k-1}=\sum_{A\sqcup B=\mathbf k,A\ne \emptyset}(s(A))_{\lvert A\rvert-1}\,(s(B))_{\lvert B\rvert}. $$ Both parts are polynomials of degree at most $k-1$, and it suffices to check that they coincide on the following combinatorial simplex $\Delta$: $x_i$ take non-negative integral values such that $\sum x_i\leqslant k-1$.
Note that for any $(x_1,\dotsc,x_k)\in \Delta$ and any partition $\mathbf k=A\sqcup B$, $A\ne \emptyset$, we have $s(A)+s(B)=s(\mathbf k)\leqslant k-1$, thus either $s(A)<\lvert A\rvert-1$ or $s(B)<\lvert B\rvert$ or $s(A)=\lvert A\rvert-1$, $s(B)=\lvert B\rvert$, $s(\mathbf k)=k-1$. Thus if $s(\mathbf k)<k-1$ both parts vanish, and if $s(\mathbf k)=k-1$, we have to prove that $$ k!=\sum_{A:s(A)=\lvert A\rvert-1} (\lvert A\rvert-1)! (k-\lvert A\rvert)! $$ I know this identity from Andrei Zelevinsky (2008). His manuscript A peculiar algebraic identity contains a clever proof both of combinatorial and falling factorials identities. Andrei asked not to circulate this text, but I think that this request is out of date….
Consider a permutation $\pi=(c_1,\dotsc,c_k)$ of $1,2,\dotsc,k$. Choose minimal $m$ for which $s(\{c_1,\dotsc,c_m\})<m$. Then by minimality $s(\{c_1,\dotsc,c_m\})=m-1$. For how many permutations did we have $\{c_1,\dotsc,c_m\}=A$, where $A$ is a fixed subset of $\mathbf k$ with $s(A)=\lvert A\rvert-1$? The answer is $(\lvert A\rvert-1)!\lvert B\rvert!$, since for any fixed cyclic order $(c_1,\dotsc,c_m)$ of $A$ exactly one out of $m$ permutations $(c_{i+1},c_{i+2},\dotsc,c_i)$ satisfies $s(c_{i+1},\dotsc,c_{i+t})\geqslant t$ for all $t=1,2,\dotsc,m-1$. This is a variant of Raney's lemma. Thus the identity.
We may avoid combinatorial arguments with Raney's lemma and remain entirely in algebra. The following trick is suggested by Vlad Volkov: we instead check the points for which $(x_1,\dotsc,x_{k-1},x_k+1)\in \Delta$. In this case $\sum x_i\in \{-1,0,1,\dotsc,k-2\}$ and LHS equals 0 unless $x_1=\dotsb=x_{k-1}=0$, $x_k=-1$. In this last case LHS equals $(-1)^{k-1}k!$. Check the same for RHS. Always either $s(A)<\lvert A\rvert-1$ or $s(B)<\lvert B\rvert=k-\lvert A\rvert$, thus if the corresponding summand in RHS is non-zero, one of the sets $A$, $B$ has a sum of $x$'s equal to $-1$, i.e., consists of $k$ and several indices $i$ such that $x_i=0$. After some calculations everything reduces to Vandermonde–Chu identity.
Now a bit about multidimensional interpolation. Consider the space $\Pi(n,k)$ of polynomials $f(x_1,\dotsc,x_k)$ with degree at most $n$. Its dimension equals the number of monomials $x_1^{c_1}\dotsm x_k^{c_k}$, where $c_i$ are non-negative integers such that $\sum c_i\leqslant n$, we write this as $(c_1,\dotsc,c_k)\in \Delta_n^k$. Thus if we manage to find $\lvert\Delta_n^k\rvert$ linearly independent linear forms on $\Pi(n,k)$, then the polynomial $f\in \Pi(n,k)$ is uniquely determined by the values of these forms. I claim that the values at points of $\Delta_n^k$ serve as such linear forms. For proving that they are linearly independent we may simply construct a polynomial from $\Pi(n,k)$ which vanishes at all points of $\Delta_n^k$ but the given point $(t_1,\dotsc,t_k)$. Here it comes: $$\prod_{i=1}^k \prod_{j=0}^{t_i-1}(x_i-j)\cdot \prod_{\ell=t_1+\dotsb+t_k+1}^n (x_1+\dotsb+x_k-\ell).$$
Extensions of the Leibnitz rule to fractional calculus were investigated by Osler in the 1970s. See, e.g., "Fractional derivatives and special functions" by Lavoie, Osler, and Tremblay.
One can even make expansions of the type
$$D^{-1} f(x)g(x) = (D_L+D_R)^{-1} f(x)g(x) = D_L^{-1}(1+D_R/D_L)^{-1}f(x)g(x)$$
$$ = \sum_{n \geq 0} \binom{-1}{n} [D^{-n-1}f(x)][D^ng(x)]= \sum_{n \geq 0} (-1)^n [D^{-n-1}f(x)][D^ng(x)]$$
for some classes of functions with, e.g., $D^{-1}f(x) = \int_0^x f(t)dt$ and the subscripts on $D$ designating action on either the Left or Right factor.
For example, let $f(x) = x^{q-r}$ and $g(x) = x^{p+r}$ . Then
$$D^{-1}f(x)g(x) = D^{-1}x^{p+q} = \frac{x^{p+q+1}}{p+q+1}$$
and
$$\sum_{n \geq 0} (-1)^n [D^{-n-1}f(x)][D^ng(x)] = \sum_{n \geq 0} (-1)^n [D^{-n-1}x^{q-r}][D^nx^{p+r}]$$
$$ =x^{p+q+1} \sum_{n \geq 0} (-1)^n \frac{(q-r)!}{(q-r+n+1)!} \frac{(p+r)!}{(p+r-n)!} , $$
so the generalized Leibnitz rule holds in this case with the generalized integration $D^{-1} \frac{x^s}{s!} = \frac{x^{s+1}}{(s+1)!}$ for $x >0$ when
$$\frac{1}{p+q+1} =\sum_{n \geq 0} (-1)^n \frac{(q-r)!}{(q-r+n+1)!} \frac{(p+r)!}{(p+r-n)!} $$
which is valid for $p+q +1 > 0$
(see p. 268 in Section 105 of Die Gammafunktion by Niels Nielsen for convergence of this sum).
