The identity theorem in complex analysis states that a holomorphic function $f\colon U \subseteq \mathbb{C} \to \mathbb{C}$ defined on a domain $U$ is identically zero if there exists a subset $S \subseteq U$ with accumulation in $U$ such that $f$ vanishes on $S$. The identity theorem can be used to show that the $\Gamma$ function is the only analytic function such that $$\Gamma(n)=(n-1)!$$ for all $n \in \mathbb{N}$. This can be generalised to several dimensions into what is known as the curve selection lemma (see e.g. Milnor's Singular Points of Complex Hypersurfaces). However, this is possible due to the rigid nature of holomorphic derivatives.
I have recently been wondering about the concept of a "$\mathscr{C}^\infty$ continuation". For example, could we find a non-analytic, $\mathscr{C}^\infty$ function $G\colon I \subseteq \mathbb{R} \to \mathbb{R}$ defined on an open neighbourhood of $\mathbb{N}$ such that $G(n)=(n-1)!$? This implies that one cannot find a bi-analytic function $\phi\colon I \to J \subseteq \mathbb{R}$ such that $\Gamma\circ \phi=G$, so $G$ is not equivalent to $\Gamma$ either (as the so-called $\Pi$ function would be).
I presume someone has come up with this topic at one point (although my MathSciNet-fu is not as strong as I would like it to be, so I am not sure how to go about it), and there may be some condition under which every $\mathscr{C}^\infty$ continuation (using the term loosely) needs to be analytic. Otherwise, is there an example of a function which admits two different non-analytic continuations?
Smooth functions are extremely flexible. Essentially, any shape you can think of (with bounded slope, etc…) can be achieved as a graph of a smooth function. An example of why this is true is given by the tool called partition of unity: this allows to “glue” the graph of two different smooth functions, e.g., given two domains $A$, $B$ with positive distance between each other, which are both compactly embedded into a third domain $D$, and given two smooth functions $f$ and $g$ defined on a neighbourhood of $A, B$ respectively, there always exists a smooth function $h$ that is defined on all $C$ and coincides with $f$ on $A$ and with $g$ on $B$.
Another tool is that of mollification/regularization. For any rough function $f$ (e.g., locally integrable, or even a distribution), you can construct a family of smooth functions $f_\varepsilon$ (parametrized by some parameter $\varepsilon>0$) that locally resemble the original function, with a series of nice properties (e.g., if $f$ is continuous on an open subset, you have local uniform convergence of $f_\varepsilon$ to $f$ on that set,…). The regularization $f_\varepsilon$ is constructed by taking the convolution between $f$ and a fixed smooth, non-negative function $\rho_\varepsilon$ with integral $1$ that is supported on a ball of radius $\varepsilon$ around the origin (this tool is available on $\mathbb R^n$ for all $n$).
With this tool you can easily construct what you want: consider the following function defined on $\mathbb C$: $$ f(z)=\sum_{n\in\mathbb N} (n-1)! \cdot \mathbf 1_{B_{1/2}(n)}(z), $$ where the above are indicator functions of the balls of radius $1/2$ centered around the natural numbers. Then, regularizing $f$, we obtain a smooth function $f_\varepsilon$ which locally resembles $f$ for small $\varepsilon$: in particular, if $\varepsilon$ is small enough, $f_\varepsilon$ councides with $f$ in a neighbourhood of $\mathbb N$ (I can try to give you some references if you want to know more about this method).
Edit: I agree with the comments you received about the non-uniqueness of $\Gamma$, but what I discussed here is unrelated from that and applies to a large variety of settings. For instance take $$ g(z)=\sum_{n\in\mathbb N} \mathbf 1_{B_{1/2}(n)}(z).$$ As before, $g_\epsilon$ is smooth, and equals $1$ on a neighbourhood of $\mathbb N$. Then, $h(z):=g_\varepsilon(z)\Gamma(z)$ is smooth and agrees with $\Gamma$ on a neighbourhood of $\mathbb N$, but it vanishes for all $z$ that lie on a slightly larger neighbourhood of $\mathbb N$.