I'm trying to find a reference for the following (analytic) proof of Wilson's Theorem. I would appreciate it immensely if some of you knew and told me the person who first came up with it. Here is the proof:
We start by considering the Maclaurin series of the function $f(x)=\ln \left( \frac{1}{1-x} \right)$, that is, $\displaystyle \ln \left( \frac{1}{1-x} \right)= \sum_{n=1}^{\infty} \frac{x^n}{n} \tag*{}$ for $x \in [-1, 1).$ By simply using the definition of $\ln$ on the equation above, we get $\displaystyle e^{x+\frac{x^2}{2}+\frac{x^3}{3}+ \cdots}=\frac{1}{1-x}=1+x+x^{2}+ \cdots + x^{p} + \cdots \tag{1}$ We can then rewrite the leftmost side of $(1)$ as $\displaystyle \begin{align} e^{x+\frac{x^2}{2}+\frac{x^3}{3}+…} & = e^{x}e^{\frac{x^2}{2}}e^{\frac{x^3}{3}} \cdots \\ & = \prod_{n=1}^{\infty} e^{\frac{x^n}{n}} \\ & = \prod_{n=1}^{\infty} \sum_{k=0}^{\infty} \frac{\left( \frac{x^n}{n} \right)^k}{k!} \\ & = \left( 1+\frac{x}{1!}+\frac{x^2}{2!}+ \cdots \right) \cdot \left( 1+\frac{x^{2}/2}{1!}+\frac{(x^{2}/2)^2}{2!} + \cdots \right) \cdots \left(1+\frac{x^{p}/p}{1!}+\frac{(x^{p}/p)^2}{2!}+ \cdots \right) \cdots \\ & = 1+\frac{x}{1!}+x^{2} \left( \frac{1}{2!} + \frac{1}{2}\right)+x^{3} \left( \frac{1}{3!} + \frac{1}{1!} \frac{1/2}{1!} + \frac{1/3}{1!} \right) + \cdots + x^{p} \left( \frac{1}{p!} + \frac{1}{(p-2)!} \frac{1/2}{1!} + \cdots + \frac{1/p}{1!} \right)+ \cdots \end{align} \tag{2}$ Combining the rightmost sides of $(1)$ and $(2)$ together gives us $\displaystyle 1+x+x^{2}+ x^{3}+ \cdots + x^{p} + \cdots= 1+\frac{x}{1!}+x^{2} \left( \frac{1}{2!} + \frac{1}{2}\right)+x^{3} \left( \frac{1}{3!} + \frac{1}{1!} \frac{1/2}{1!} + \frac{1/3}{1!} \right) + \cdots + x^{p} \left( \frac{1}{p!} + \frac{1}{(p-2)!} \frac{1/2}{1!} + \cdots + \frac{1/p}{1!} \right)+ \cdots , \tag{3}$ from which follows that the coefficient of $x^{p}$ on the RHS of $(3)$ is $1$. That is, $\displaystyle \frac{1}{p!}+\frac{1}{(p-2)!} \frac{1/2}{1!}+ \cdots + \frac{1/p}{1!} = 1. \tag{4}$ Notice that such a coefficient is of the form $\frac{1}{p!}+\frac{r}{s}+\frac{1}{p}$, where $\frac{r}{s}$ is the sum of a finite number of rationals, none of which has the factor $p$ in its denominator. Therefore, if $\frac{r}{s}$ is irreducible, then $p$ doesn’t divide $s$. We then proceed by rewriting $(4)$ as $\displaystyle \frac{1}{p!}+\frac{r}{s}+\frac{1}{p}=1. \tag*{}$ Equivalently, $\displaystyle 1-\frac{r}{s}=\frac{1}{p!}+\frac{1}{p}=\frac{1+(p-1)!}{p!}, \tag*{}$ and $\displaystyle (s-r) \cdot (p-1)! = \frac{s(1+(p-1)!)}{p}. \tag*{}$ Hence, $\frac{s(1+(p-1)!)}{p}$ is an integer. Moreover, since $p$ doesn’t divide $s$, $\frac{1+(p-1)!}{p}$ is also an integer. Equivalently, $(p-1)! \equiv -1 \text{ (mod } p)$, as desired.
I don't have a reference but I would not call this an "analytic" proof; although you are using power series every step of the argument works fine in formal power series, so really you are doing combinatorics with generating functions. The combinatorial content of this argument can be extracted as follows.
The identity $\exp \log \frac{1}{1 - x} = \frac{1}{1 - x}$ expresses the uniqueness of cycle decompositions of permutations; it generalizes to an identity called the permutation form of the exponential formula which expresses the same uniqueness but in more detail. Comparing the coefficient of $x^n$ (by the way, I don't think your version of this expansion has enough terms - there should be one for every partition of $n$) expresses the uniqueness of cycle decompositions of permutations in the symmetric group $S_n$, and can be interpreted as its class equation
$$n! = \sum_{[\pi]} \frac{n!}{|C_{S_n}(\pi)|}$$
which expresses $n!$ as the sum of sizes of all conjugacy classes in $S_n$. In more detail, the correct form of this expansion is
$$n! = \sum_{\sum ic_i = n} \frac{n!}{\prod_{i=1}^p i^{c_i} c_i!}.$$
where the term indexed by $\{ c_i \}$ in the sum corresponds to the conjugacy class of permutations with $c_i$ cycles of length $i$. The denominators here are the sizes of the centralizers of each conjugacy class. In terms of the generating function, the factor of $i^{c_i} c_i!$ in the denominator is contributed by the expansion of $\exp \left( \frac{x^i}{i} \right)$. Note that once you've proven this identity in any way, for example by calculating the sizes of the centralizers of every conjugacy class of permutations, you can throw away the generating function (which is equivalent to it).
Now we specialize $n$ to a prime $p$ and take the class equation $\bmod p$. From the above equation it's not hard to see that the denominator can only be divisible by $p$ if either $c_1 = p$ or $c_p = 1$, corresponding to the conjugacy classes of the identity and the $p$-cycles. This can be proven conceptually as follows. $S_p$ has the property that its only elements of order $p$ are the $p$-cycles, which are all conjugate and which generate self-centralizing subgroups: that is, the centralizer of a $p$-cycle is just the cyclic subgroup it generates. It follows that the only conjugacy classes in $S_p$ whose centralizer has order divisible by $p$ (equivalently, whose centralizer contains a $p$-cycle) are the identity and the $p$-cycles. Taking the class equation $\bmod p$, every term except the terms corresponding to the identity and the $p$-cycles vanishes, and we conclude that the number of $p$-cycles plus $1$ is congruent to $0 \bmod p$. But since we calculated the centralizer already, it follows that there are $(p-1)!$ $p$-cycles, so we conclude that
$$(p-1)! + 1 \equiv 0 \bmod p$$
as desired.
So, stripped of the power series trappings, this is a straightforward counting argument in $S_p$. It is related to counting arguments which can be used to prove the Sylow theorems, especially Sylow III (although one of the proofs of Sylow I also uses the class equation). Sylow III applied to $S_p$ actually gives $(p-2)! \equiv 1 \bmod p$ (this is the number of conjugacy classes of cyclic subgroups of order $p$) which is straightforwardly equivalent to Wilson's theorem.
There is also a simpler version of the counting argument which directly uses the $p$-group fixed point theorem, which is just the easy statement that if $P$ is a finite $p$-group acting on a finite set $X$ then
$$|X^P| \equiv |X| \bmod p.$$
We apply the PGFPT to the action of a $p$-cycle on the set of $p$-cycles in $S_p$. As above, there are $(p-1)!$ of these, and the centralizer of any $p$-cycle is the cyclic subgroup it generates. It follows that this action has exactly $p-1$ fixed points, which immediately gives
$$(p-1)! \equiv (p-1) \bmod p.$$
A similarly quick and easy counting proof can be given for Fermat's little theorem and a few other fun facts; see the linked blog post for more.
Just for fun, here is a formal power series proof that $\exp \log \frac{1}{1 - x} = \frac{1}{1 - x}$, which requires no "actual" analysis. As above, this identity is equivalent to a combinatorial statement about permutations and can be proven that way, but it can also be proven as follows. Namely, it's straightforward to show by induction that $\exp(x)$ is the unique formal power series $f(x)$ satisfying $f(0) = 1, f'(x) = f(x)$ (just differentiate to get that every coefficient of the Taylor series is $1$). On the other hand, $\log \frac{1}{1 - x}$ is the unique formal power series $g(x)$ satisfying $g(0) = 0, g'(x) = \frac{1}{1 - x}$ (just compare terms on both sides). (Both of these facts can be given combinatorial interpretations.)
It follows that $f(g(x)) = h(x)$ satisfies $h(0) = 1$ and
$$\frac{d}{dx} \exp \log \frac{1}{1 - x} = \left( \exp \log \frac{1}{1 - x} \right) \frac{1}{1 - x}$$
and there is a unique formal power series satisfying $h(0) = 1$ and $h'(x) = \frac{h(x)}{1 - x}$ (by induction; this follows from a formal version of the Picard-Lindelof theorem, which is much easier to prove), namely $h(x) = \frac{1}{1 - x}$.