Representing a factorial through its pseudo Z transform

JPaquim
Ok, so I was playing around with some Z transforms. I'm sorry about the long derivation, but I'm a bit unsure of the mathematical rigor, and want to make sure every step is clear. I started with the recurrence relation defining the factorial:
$$n!: u_{n+1}=(n+1)u_n=u_n+nu_n$$ $$u_0 = 1$$
I imagined that, even though the series diverges for every z, that I could formally work with it, and so took its Z transform
$$Z\{u_n\}(z)=\sum_{n=0}^\infty u_n z^{-n}=U(z)$$
I was a bit rusty and decided to derive a couple of Z transform properties by hand, resulting in:
$$Z\{u_{n+1}\}(z)=\sum_{n=0}^\infty u_{n+1} z^{-n}= \sum_{n=1}^\infty u_{n} z^{-n+1}= z\sum_{n=1}^\infty u_{n} z^{-n}= z(\sum_{n=0}^\infty u_{n} z^{-n}-u_0)= z(U(z)-u_0)$$ $$\frac{d}{dz}U(z)=\frac{d}{dz}\sum_{n=0}^\infty u_n z^{-n}= \sum_{n=0}^\infty u_n \frac{d}{dz}(z^{-n})= \sum_{n=0}^\infty u_n (-n) z^{-n-1}= -\frac{1}{z}\sum_{n=0}^\infty n\,u_n z^{-n}= -\frac{1}{z}Z\{n\,u_n\}(z)$$ $$\Leftrightarrow Z\{n\,u_n\}(z) = -z\frac{d}{dz}U(z)$$
Applying these two properties to the recurrence relation results in the following differential equation, which I rearranged and solved using an integrating factor:
$$z(U(z)-u_0)=U(z)-z\frac{d}{dz}U(z)$$ $$\Leftrightarrow z\frac{d}{dz}U(z) + (z-1)U(z)=u_0z$$ $$\Leftrightarrow \frac{d}{dz}U(z) + (1-\frac{1}{z})U(z)=u_0$$ $$\Leftrightarrow \frac{d}{dz}(U(z)\frac{e^z}{z})=u_0\frac{e^z}{z}$$ $$\Leftrightarrow U(z)=ze^{-z}(u_0\int \frac{e^z}{z} dz + K)$$
Having the Z transform, I then applied the inversion formula:
$$u_n = Z^{-1}\{U(z)\}(n) = \frac{1}{2\pi i}\oint_\gamma U(z)z^{n-1}dz$$ $$=\frac{1}{2\pi i}\oint_\gamma z^n e^{-z}(u_0\int \frac{e^z}{z} dz + K)dz$$
Only terms with non-zero residue will contribute to the contour integral. Expanding the inner integral as a power series:
$$\int \frac{e^z}{z} dz = \int \sum_{m=0}^\infty \frac{z^{m-1}}{m!}= \sum_{m=0}^\infty \frac{1}{m!}\int z^{m-1}dz = \log z + \sum_{m=1}^\infty \frac{z^m}{m!m} = \log z + f(z),$$ where f(z) is holomorphic, and so doesn't contribute to the contour integral. $$u_n = \frac{u_0}{2\pi i}\oint_\gamma z^n e^{-z}\log z\,dz$$
Substituting $z = e^{i\theta}, d\theta=\frac{dz}{iz}$, as well as $u_0 = 1$ results in:
$$n! = \frac{i}{2\pi}\int_{-\pi}^\pi e^{i(n+1)\theta} e^{-e^{i\theta}}\theta\,d\theta$$
Plugging this into Mathematica, to my astonishment, produces correct results for integer n. Curious, I tried plugging in non-integer values, to see if it matched the Gamma function, or some other extension of the factorial. However, Mathematica doesn't seem to be able to carry out the calculation...
Any thoughts on this sort of computation/representation?
I also tried playing around with a similar recurrence relation:
$$u_{n+1}=nu_n,$$
which corresponds to a somewhat simpler, constant coefficients, differential equation, with an elementary and holomorphic solution, which when plugged into the inversion formula, gives out zero for every n other than zero...