The gamma function extends the factorial function. This can be proved inductively using integration-by-parts.
But if you didn't already know that the gamma function had this property, and you wanted to construct a continuation of factorial a priori, how could you construct the gamma function? How could you motivate an exploration which would lead to this answer?
Euler was the first to discover/invent the Gamma function. His intuition was as follows: one complicated way to calculate $n!$ is as the infinite product $$\prod_{k=1}^{\infty} \frac{k}{k+n}.$$ The idea is that as $k$ tends to infinity in the denominator we get "all terms from $1+n$ till $\infty$" and in the numerator we get "all terms from $1$ till $\infty$" so that after cancellations we are left with $n!$. The beauty is that this product also makes sense for non-integer $n$.
There is a little problem with convergence though: the above product diverges to $0$. To see this, consider the logarithm of the product which gives us the series $$\sum_{k=1}^{\infty} \log\left(\frac{k}{k+n}\right)$$ and while the general term does converge to $0$, the series diverges to $-\infty$, which follows since from $\log x \sim x-1$ we see that the series should behave roughly like the (negative of the) harmonic series. To summarize, the series diverges to $-\infty$ so the product diverges to $0$.
Euler's second brilliant idea was to re-write the product instead as $$\prod_{k=1}^{\infty} \frac{k^{1-n}(k+1)^n}{k+n}=\frac{1\cdot 2^n}{1+n}\cdot\frac{2^{1-n}\cdot 3^n}{2+n}\cdot\frac{3^{1-n}\cdot4^n}{3+n}\cdots $$ where the idea is that if we multiply the numerators up to the $k$th term we get $k!\cdot (k+1)^n$, and the $(k+1)^n$ term is again going to cancel with the term after it and so on ad infinitum, so over all in the numerator we are left with "the product of all terms from $1$ to $\infty$" which is the same as we had in the initial product. This product is the first expression Euler gives for the Gamma function, in a letter to Goldbach in 1729. The advantage of this form is simply that the infinite product now converges.
Here we already have an extension of the factorial function to non-integer values. For instance, Euler is able to find $(\frac{1}{2})!=\frac{\sqrt{\pi}}{2}$ by substituting $n=\frac{1}{2}$ into the product and utilizing the formula $$\frac{\pi}{4}=\frac{2\cdot 4}{3\cdot 3}\cdot\frac{4\cdot 6}{5\cdot 5}\cdot\frac{6\cdot 8}{7\cdot 7}\cdots$$ which was found by Wallis some decades earlier.
It is speculated that this value is what made Euler think that it should be possible to express this infinite product in terms of an integral (as $\pi$ is related to areas, and is ubiquitous in the result of many definite integrals).
For this purpose Euler started by looking at $\int_0^1 x^a (1-x)^n dx$ where $n$ is an integer and $a$ is an arbitrary constant. Expanding $(1-x)^n$ by the binomial theorem gives that the integral equals $\frac{n!}{\prod_{k=1}^{n+1} a+k}$. Special cases of this integral were discussed by Wallis, by Newton, and by Stirling, so it stands to reason that Euler was familiar with this integral and chose it as a starting point since evaluating it gives an expression which has $n!$ in it. Some ingenious manipulations on this integral allowed Euler to express $n!$ itself as an integral, the final result being that $n! = \int_0^1 (-\log x)^n dx$. From here it is a simple substitution to get to the familiar modern definition of the Gamma function.
This answer is mostly based on the great How Euler Did It - Gamma the function by Ed Sandifer and the references there, specifically the reference [D] which elaborates on the manipulations required to get to $\int_0^1 (-\log x)^n dx$ (pages 853-855). For a lot more detail check out that article and the references.
The presentation of the Gamma function which I really liked is the one from Ahlfors. Now, the factorial is something defined for non-negative integers only, so if we were to look for nice extensions, we ought to be looking for analytic functions (I'm no history expert, but I believe at the time of "creation" of the Gamma function, people were mostly interested in analytic functions).
Now, a common way of defining analytic functions is by infinite series and infinite products. Factorials are defined as repeated products, so it is only natural that an analytic extension of the factorial involves infinite products in its definition. Now, for $k\geq 1$, $k!$ is the product of all positive integers up to $k$, so we can write it naively as an infinite product of all integers $1$ to infinity divided by all integers from $n+1$ to infinity: \begin{align} k!\,\, &=\prod_{n=1}^{\infty}\frac{n}{n+k}=\frac{1}{\prod_{n=1}^{\infty}\frac{n+k}{n}}=\frac{1}{\prod_{n=1}^{\infty}\left(1+\frac{k}{n}\right)}\tag{$*$} \end{align} All these equal signs should be treated heuristically, since none of the infinite products converge. Now, the reason I wrote the reciprocal is because this now motivates us to consider the expression $\prod_{n=1}^{\infty}\left(1+\frac{z}{n}\right)$ for all $z\in\Bbb{C}$ (again set aside issues of convergence). Therefore, factorials are closely related to (reciprocals of) complex functions which vanish precisely at negative integers.
Thus far we have been very sloppy with convergence, but with the theory of infinite products, there's a standard way to get nice results; this is the theory of elementary Weierstrass factors and the Weierstrass factorization theorem. So, in order to produce an honest analytic function which has simple zeroes precisely at the negative integers, we consider \begin{align} G(z)&:=\prod_{n=1}^{\infty}\left(1+\frac{z}{n}\right)e^{-z/n}; \end{align} this converges absolutely and uniformly on all compact subsets of $\Bbb{C}$, and therefore defines an entire function which has simple zeros precisely at the negative integers (the extra factor of $e^{-z/n}$ is what forces the convergence). Now, $z\mapsto G(z-1)$ is an analytic function which vanishes precisely on the non-positive integers $\Bbb{Z}_{\leq 0}$. Therefore, the quotient $\frac{G(z-1)}{zG(z)}$ is well-defined on $\Bbb{C}\setminus\Bbb{Z}_{\leq 0}$, is non-vanishing there, and is actually bounded near the zeros of the denominator, i.e near $\Bbb{Z}_{\leq 0}$. Therefore, it can be extended analytically (by Riemann's removable singularity theorem) to an entire function which is nowhere vanishing. Such functions can be written as $e^{\gamma(z)}$ for some entire $\gamma$. Therefore, $G$ satisfies the functional equation \begin{align} G(z-1)&=ze^{\gamma(z)}G(z) \qquad(z\in\Bbb{C}). \end{align} Since the infinite product converges absolutely and uniformly on compact subsets of the plane, we can perform term-by-term logarithmic differentiation to get that (for all $z$ in a punctured neighborhood of the origin) \begin{align} \sum_{n=1}^{\infty}\frac{1/n}{1+(z-1)/n}-\frac{1}{n}&= \frac{1}{z} +\gamma'(z)+\sum_{n=1}^{\infty}\frac{1/n}{1+z/n}-\frac{1}{n} \end{align} Some simple index juggling shows that $\gamma'(z)=0$. Since this is true for all $z$ in a punctured neighborhood of the origin, and thus by uniqueness of analytic continuation, $\gamma'=0$ identically, so that $\gamma$ is a constant function. We shall denote this constant value as $\gamma$. Therefore, there is a constant $\gamma\in\Bbb{C}$ such that for all $z\in\Bbb{C}$, we have \begin{align} G(z-1)&=e^{\gamma}zG(z)\tag{$**$}. \end{align} By the way, if you plug in $z=1$ into this equation, you'll see that $G(0)=1=e^{\gamma}G(1)$, and hence $e^{-\gamma}=\prod_{n=1}^{\infty}\left(1+\frac{1}{n}\right)e^{-1/n}$, so that actually $\gamma\in\Bbb{R}$, and this is actually the usual Euler-Mascheroni constant. Now, the functional equation $(**)$ is not so nice. If we define $H(z)=e^{\gamma z}G(z)$, then we can easily verify using $(**)$ that \begin{align} H(z-1)&=zH(z). \end{align} This is a slightly nicer functional equation; it's almost the same as the factorial equation, except it has $z-1$ instead of $z+1$. This is easily fixed: we now define $\Gamma:\Bbb{C}\setminus\Bbb{Z}_{\leq 0}\to\Bbb{C}$ as \begin{align} \Gamma(z)&:=\frac{1}{zH(z)}=\frac{1}{ze^{\gamma z}G(z)} =\frac{e^{-\gamma z}}{z}\prod_{n=1}^{\infty}\left(1+\frac{z}{n}\right)^{-1}e^{z/n}.\tag{$***$} \end{align} Based on the functional equation for $H$, it immediately follows that $\Gamma(z+1)=z\Gamma(z)$, and that $\Gamma(1)=1$. Therefore, $\Gamma(n+1)=n!$ for all integer $n\geq 0$, so that $\Gamma$ is an analytic extension of the factorial to $\Bbb{C}\setminus\Bbb{Z}_{\leq 0}$... i.e it is a meromorphic extension of the factorial to $\Bbb{C}$ having simple poles at the non-positive integers.
Also, $\Gamma$ maps $(0,\infty)$ into $(0,\infty)$, as easily seen from the explicit product formula. Note also that by taking logarithms in $(***)$ and differentiating twice, we get that for all $z\in (0,\infty)$, $(\log \circ \Gamma)''(z)=\sum_{n=0}^{\infty}\frac{1}{(z+n)^2}>0$, so the Gamma function is logarithmically convex on the positive real axis.
To recap: we started off heuristically writing factorials as an infinite product in $(*)$. Through this, we realized that factorials are related to reciprocals of complex functions having zeroes at the negative integers. Using this inspiration, we considered a properly defined function $G$ (it converges nicely because of a simple estimate on the elementary Weierstrass factors $\left(1+\frac{z}{n}\right)e^{-z/n}$). Now, because we introduced the "convergence factor" $e^{-z/n}$ in the definition of $G$, that results in the functional equation $(**)$ being not-so-pretty. But this is easily remedied by introducing $H$, and then finally taking reciprocals, we get $\Gamma$.
This approach also has the added benefit that it easily relates the Gamma function to trigonometry. It is a "standard" fact that using the Weierstrass factorization theorem that $\sin (\pi z)$ has the infinite product expansion \begin{align} \sin(\pi z)&= \pi z \prod_{n\neq 0}\left(1-\frac{z}{n}\right)e^{z/n}=\pi z G(z)G(-z) \end{align} From the definition of $\Gamma$, it thus follows that \begin{align} \sin(\pi z)&=\pi z\frac{1}{ze^{\gamma z}\Gamma(z)}\frac{1}{(-z)e^{-\gamma z}\Gamma(-z)}\\ &=\pi \frac{1}{\Gamma(z)}\frac{1}{\Gamma(1-z)}, \end{align} where we used $(-z)\Gamma(-z)=\Gamma(1-z)$ from the functional equation, and then we cancelled $z$'s and the exponentials. In other words, \begin{align} \Gamma(z)\Gamma(1-z)&=\frac{\pi}{\sin(\pi z)}.\tag{$\ddot{\smile}$} \end{align} In particular, $\Gamma\left(\frac{1}{2}\right)^2=\pi$, and since $\Gamma$ is positive on the positive real axis, we have $\Gamma\left(\frac{1}{2}\right)=\sqrt{\pi}$.
Note that my answer here is approaching from the angle of "how to arrive at a function which extends the factorial", i.e by means of a functional equation, rather than motivating the integral definition of $\Gamma$. I think the integral definition is the most "rabbit out of a hat" if we're looking purely from the perspective of the functional equation.
A-priori, it is not obvious that $\Gamma$ as defined here admits any integral representation. Anyway, what we can do is try to rewrite the definition of $\Gamma$: \begin{align} \Gamma(z)&:=\frac{e^{-\gamma z}}{z}\prod_{n=1}^{\infty}\left(1+\frac{z}{n}\right)^{-1}e^{z/n}=\frac{e^{-\gamma z}}{z}\prod_{n=1}^{\infty}\frac{n}{n+z}e^{z/n}. \end{align} One way of massaging this expression and getting rid of the infinite product is to explicitly plug in the value of the Euler-Mascheroni constant $\gamma=\lim\limits_{n\to\infty}\sum_{k=1}^n\frac{1}{k}-\log n$ (this follows pretty easily from what we mentioned above that $G(0)=1=e^{\gamma}G(1)$). So, we have (by continuity of $\exp$) that for all $z\in\Bbb{C}\setminus\Bbb{Z}_{\leq 0}$, \begin{align} \Gamma(z)&=\lim_{n\to\infty}\frac{e^{-\left(\sum_{k=1}^n\frac{1}{k}-\log n\right)z}}{z}\prod_{k=1}^n\frac{k}{k+z}e^{z/k}\\ &=\lim_{n\to\infty}\frac{e^{z\log n}\prod_{k=1}^ne^{-z/k}}{z}\prod_{k=1}^n\frac{k}{k+z}e^{z/k}\\ &=\lim_{n\to\infty}\frac{n^z}{z}\prod_{k=1}^n\frac{k}{k+z}\\ &=\lim_{n\to\infty}n^z\cdot \frac{n!}{z(z+1)\cdots (z+n)}. \end{align} Up to this step, I think everything is reasonably straightforward (i.e once you know about Weierstrass factorization, it's essentially "follow your nose and grind out the calculations" to arrive at the definition for $\Gamma$ as a function which satisfies the same recurrence as the factorial, and then a simple manipulation yields this other representation).
In my opinion, it is at this stage that one needs a bit of brilliance. The fraction can actually be expressed as an integral if $\text{Re}(z)>0$. How do we see this intuitively? Well, by pattern recognition, we should know that $\int_0^1s^z\,ds=\left[\frac{s^{z+1}}{z+1}\right]_0^1=\frac{1}{z+1}$. So, if we want to introduce many such divisions, then a little trial and error shows that \begin{align} \int_0^1s^{z-1}(1-s)^n\,ds&=\frac{n!}{z(z+1)\cdots (z+n)}. \end{align} Hence, for $\text{Re}(z)>0$, \begin{align} \Gamma(z)&=\lim_{n\to\infty}n^z\frac{n!}{z(z+1)\cdots (z+n)}\\ &=\lim_{n\to\infty}n^z\int_0^1s^{z-1}(1-s)^n\,ds\\ &=\lim_{n\to\infty}\int_0^1(ns)^{z-1}(1-s)^n\,n\,ds\\ &=\lim_{n\to\infty}\int_0^nt^{z-1}\left(1-\frac{t}{n}\right)^n\,dt\tag{put $t=ns$}\\ &=\int_0^{\infty}t^{z-1}e^{-t}\,dt, \end{align} where the final equal sign requires a careful dominated convergence theorem argument to justify why the limit and integral can be swapped.
Remarks.
Note that there is a certain amount of non-uniqueness in the definition of the Gamma function if we only impose that it satisfy the functional equation (pretty much because the Weierstrass factorization of an entire function is not unique).
More precisely, let $\mathcal{F}$ be the set of all meromorphic functions $\Phi$ on $\Bbb{C}$ such that they are non-vanishing and have simple poles precisely at $\Bbb{Z}_{\leq 0}$, and satisfy $\Phi(1)=1$ together with $\Phi(z+1)=z\Phi(z)$. Next, define $\mathcal{G}$ to be the set of all entire functions $f$ such that $f(0)\in 2\pi i \Bbb{Z}$, and there exists a $k\in\Bbb{Z}$ such that for all $z\in\Bbb{C}$, we have $f(z+1)=f(z)+2\pi i kz$. Then, we can show that \begin{align} \mathcal{F}&=\left\{e^{f(\cdot)} \Gamma(\cdot)\,|\, f\in \mathcal{G}\right\} \end{align} Note that the collection $\mathcal{G}$ is not trivial; it contains functions like $f_{k,l}(z):= e^{2\pi i l z}+ 2\pi i kz$ for all $k,l\in\Bbb{Z}$.
However, out of all the functions in the collection $\mathcal{F}$, the Gamma function is the only one mapping $(0,\infty)$ to $(0,\infty)$ and is logarithmically convex (by the Bohr-Mollerup theorem). Logarithmic convexity is desirable because it provides us with with many nice inequalities. This is probably one of the reasons why the Gamma function is so popular (also, the integral defining $\Gamma$ appears all the time, in non-intuitive places, such as in calculating areas/volumes of higher dimensional balls/spheres).
I think a good strategy for finding a generalization of the factorial is to first find a well known formula in which it appears, and then try to extrapolate the formula. One such formula is the $n$-th moment of the exponential distribution:
Suppose you have an exponentially distributed random variable $X \sim \text{Exp}(1)$. By definition, this means that the distribution of $X$ has density $f(t) = e^{-t}1_{[0, \infty)}$ with respect to Lebesgue measure. It is of interest to calculate the moments $E(X^n) = \int_{0}^{\infty}t^ne^{-t}\,dt$ where $n \geq 1$ is an integer. If $p(t)$ is a polynomial of degree at most $n$ and $\lambda \in \mathbb{C} \setminus \{0\}$, then the function $p(t)e^{\lambda t}$ has antiderivative $$\int p(t)e^{\lambda t}\,dt = \frac{1}{\lambda}\left(\sum_{j = 0}^{n}(-1)^{j}\frac{D^j}{\lambda^j}p(t)\right)e^{\lambda t},$$ where $D$ is the differentiation operator. To prove this, note that if $q(t)$ is a polynomial, then $\frac{d}{dt}(q(t)e^{\lambda t}) = (q'(t) + \lambda q(t))e^{\lambda t}$. Hence $\int p(t)e^{\lambda t}\,dt = (D + \lambda I)^{-1}p$. Expanding $(D + \lambda I)^{-1} = \frac{1}{\lambda}(I + \frac{D}{\lambda})^{-1}$ as a geometric series and using $D^{n + 1}p = 0$ gives the result. Specializing gives $$\int t^ne^{-t} \,dt = -\sum_{j = 0}^{n}n(n - 1)\dots(n - j + 1)t^{n - j}e^{-t}.$$ Since $t^{n - j}e^{-t}|_{0}^{\infty} = 0$ for $j < n$, we get $$E(X^n) = -n!e^{-t}|_{0}^{\infty} = n!.$$ To define our extension of the factorial to real $s$, we take $s! := E(X^s) = \int_{0}^{\infty}t^{s}e^{-t}\,dt$. This is absolutely convergent for $s > -1$. Therefore we can even let $s$ be complex as long as $\Re(s) > -1$. By writing $t^s = \exp(s \log(t))$ and differentiating under the integral sign, it follows that this factorial function is holomorphic on $\{s \in \mathbb{C} : \Re(s) > -1\}$.
Derivation of an Extended Expression
For $n\gt m$,
$$
\begin{align}
\prod_{k=1}^n\frac{k}{k+m}
&=\frac{1\cdot2\cdot3\cdots n}{(m+1)(m+2)(m+3)\cdots(m+n)}\tag{1a}\\
&=\frac{\overbrace{1\cdot2\cdot3\cdots m}^{m!}}{\underbrace{(n+1)(n+2)(n+3)\cdots(n+m)}_{\text{between $(n+1)^m$ and $(n+m)^m$}}}\tag{1b}
\end{align}
$$
Explanation:
$\text{(1a):}$ rewrite the product
$\text{(1b):}$ cancel $(m+1)(m+2)(m+3)\cdots n$ in the numerator and denominator
Noting that the denominator in $\text{(1b)}$ is between $(n+1)^m$ and $(n+m)^m$, multiply $(1)$ by $\frac{(n+1)^m}{m!}$ to get $$ \left(\frac{n+1}{n+m}\right)^m\le\frac{(n+1)^m}{m!}\prod_{k=1}^n\frac{k}{k+m}\le1\tag2 $$ Since $n+1=\prod\limits_{k=1}^n\frac{k+1}k$, inequality $(2)$ can be rewritten as $$ \left(1-\frac{m-1}{n+m}\right)^m\le\frac1{m!}\prod_{k=1}^n\frac{k}{k+m}\left(\frac{k+1}k\right)^m\le1\tag3 $$ Applying the Squeeze Theorem and multiplying by $m!$, we get $$ m!=\prod_{k=1}^\infty\frac{k}{k+m}\left(\frac{k+1}k\right)^m\tag4 $$ $(4)$ is a good extension of factorial; it converges as long as $m$ is not a negative integer.
Extension of the Recursion
$$
\begin{align}
\frac{m!}{(m-1)!}
&=\prod_{k=1}^\infty\frac{k+m-1}{k+m}\frac{k+1}{k}\tag{5a}\\
&=\prod_{k=1}^\infty\color{#C00}{\frac{k+m-1}{k}}\color{#090}{\frac{k+1}{k+m}}\tag{5b}\\
&=\prod_{k=1}^\infty\frac{\color{#C00}{f(k)}}{\color{#090}{f(k+1)}}\tag{5c}\\
&=\lim_{k\to\infty}\frac{f(1)}{f(k)}\tag{5d}\\[9pt]
&=m\tag{5e}
\end{align}
$$
Explanation:
$\text{(5a):}$ apply $(4)$
$\text{(5b):}$ rearrange denominators
$\text{(5c):}$ $f(k)=\frac{k+m-1}{k}$
$\text{(5d):}$ apply telescoping product
$\text{(5e):}$ evaluate limit
Thus, $(4)$ obeys the usual recursion: $m!=m(m-1)!$.
Log Convexity
$$
\begin{align}
\frac{\mathrm{d^2}}{\mathrm{d}m^2}\log(m!)
&=\frac{\mathrm{d^2}}{\mathrm{d}m^2}\sum_{k=1}^\infty\left(m\log\left(\frac{k+1}k\right)-\log\left(\frac{k+m}k\right)\right)\tag{6a}\\
&=\sum_{k=1}^\infty\frac1{(k+m)^2}\tag{6b}\\[6pt]
&\gt0\tag{6c}
\end{align}
$$
Explanation:
$\text{(6a):}$ apply $(4)$
$\text{(6b):}$ take two derivatives
$\text{(6c):}$ sum of squares is positive
Thus, $(4)$ is log-convex.
Apply Bohr-Mollerup
Since it is an infinite product of ones, $0!=1$. Therefore, $(5)$, $(6)$, and Bohr-Mollerup say that $(4)$ gives $m!=\Gamma(m+1)$.
