I was wondering what would happen if we took Newton's Interpolation formula (series) and applied it the factorial.
It is given as
$$f(x)=\sum_{k=0}^\infty\binom{x-a}k\Delta^k[f](a)$$
where
$$\Delta^n[f](a)=\sum_{k=0}^n\binom nk(-1)^{n-k}f(a+k)$$
and I want $f(x)=x!$
Putting this all in with $a=0$, I get
$$x!=\sum_{k=0}^\infty\sum_{n=0}^k\binom xk\binom nk(-1)^{n-k}k!$$
which, notably, is capable of taking non-integer values for $x$.
So of course, I am concerned with its radius of convergence and especially if it comes out to equal the gamma function.
As a side question, if this fails, then are there any interpolations of the factorial that come out to the factorial?
Notice that the inner sum can be reduced to
$$ \sum_{k=0}^{n} \binom{n}{k}(-1)^{n-k} k! = n! \sum_{j=0}^{n} \frac{(-1)^j}{j!} =: d_n, $$
which is the $n$-th derangement. It is easy to check that $d_n \sim n!/e$ as $n \to \infty.$ On the other hand, if $x \notin \{0, 1, 2, \cdots \}$, then we have
$$\left| \binom{x}{n} \right| = \left| \frac{\Gamma(n-x)}{n!\Gamma(-x)} \right| \sim \frac{1}{|\Gamma(-x)| n^{\Re(x)+1}}$$
as $n \to \infty$. So it follows that the series
$$ \sum_{n=0}^{\infty} \binom{x}{n} \sum_{k=0}^{n} \binom{n}{k}(-1)^{n-k} k! = \sum_{n=0}^{\infty} \binom{x}{n} d_n \tag{1} $$
diverges unless $x$ is a non-negative integer.
We can possibly assign a value to the formal series $\text{(1)}$ under a suitable renormalization. (Of course, the meaningfulness of this assignment is thus purely one's taste of matter.) To this end, consider the following representation of $d_n$:
$$ d_n = \int_{0}^{\infty} (t-1)^n e^{-t} \, dt. $$
Although totally invalid, we may attempt the following ''computation'':
$$ \sum_{n=0}^{\infty} d_n \ ``\text{=''} \int_{0}^{\infty} \sum_{n=0}^{\infty} \binom{x}{n} (t-1)^n e^{-t} \, dt \ ``\text{=''} \int_{0}^{\infty} t^x e^{-t} \, dt = x!. $$
One may try to give a meaning to this nonsense by considering a suitable renormalization. So we consider the following function
$$ d_n(\epsilon) := e^{-\epsilon n} \int_{0}^{1} (t-1)^n e^{-t} \, dt + \int_{1}^{\infty} (t-1)^n e^{-\epsilon n t} e^{-t} \, dt. $$
For each fixed $n$, we have $\lim_{\epsilon \to 0^+} d_n(\epsilon) = d_n (0)$. Now we consider the following renormalized sum:
$$ f(x, \epsilon) := \sum_{n=0}^{\infty} \binom{x}{n} d_n(\epsilon). $$
Then on the region $\mathcal{D} = \{\epsilon \in \Bbb{C} : \Re(\epsilon) > W_0(\frac{1}{e}) \}$, we have
$$ \sup_{t \in [1, \infty)} |(t-1) e^{-\epsilon t}| \leq \frac{e^{-1-\Re(\epsilon)}}{\Re(\epsilon)} < 1 $$
and thus $d_n(\epsilon)$ decays exponentially. Thus the function $\epsilon \mapsto f(x, \epsilon)$ is holomorphic on $\mathcal{D}$. Moreover, using Fubini's theorem and the binomial series we can write
$$ f(x, \epsilon) = \int_{0}^{1} (1 + (t-1) e^{-\epsilon})^x e^{-t} \, dt + \int_{1}^{\infty} (1 + (t-1) e^{-\epsilon t})^x e^{-t} \, dt $$
It is not hard to prove that the right-hand side defines a holomorphic function near an open neighborhood of $(0, \infty)$. So we can use this expression to extend the function $\epsilon \mapsto f(x, \epsilon)$ and we may define
$$ \sum_{n=0}^{\infty} \binom{x}{n} d_n := \lim_{\epsilon \to 0^+} f(x, \epsilon) = \int_{0}^{\infty} t^x e^{-t} \, dt = x!. $$