Involutions of $S_n$

An involution is any permutation $\sigma$ such that $\sigma^2 = \text{id}$. Hence, the number of involutions of $S_n$ will be the number of permutations $\sigma \in S_n$ such that $\sigma$ has order $2$.

The key idea is that $\sigma \in S_n$ has order $2 \iff \sigma$ is a $2$-cycle or a product of disjoint $2$-cycles. If $\sigma$ is a product of $m$ $2$-cycles, then there are $\displaystyle \binom{n}{2m}$ ways of choosing $2m$ elements, and from there $\displaystyle \frac{(2m)!}{m!2^m}$ ways of pairing the selected elements together into disjoint $2$-cycles. Bearing in mind that $m$ cannot exceed $\displaystyle \frac{n}{2}$, the total number of involutions is thus given by:

$$\sum_{k=0}^{\lfloor n/2 \rfloor}\binom{n}{2k}\frac{(2k)!}{k!2^k}$$

Edit: Depending on your definition of "involution", you may instead want the sum to start at $k=1$. The $k=0$ summand is accounting for the identity permutation.

Another approach because I'm partial to generating functions: Let $a_n$ denote the number of involutions in $S_n,$ and note that $a_0 = a_1=1.$ Suppose we want to construct involutions in $S_n, n\geq 2$ (of which there are $a_n$ of). Such an involution either fixes $n$ or moves it somewhere. If it fixes $n$, then it is just an involution on the remaining $n-1$ elements, and there are $a_{n-1}$ of these. If it moves $n,$ then it is free to move it to any of the other $n-1$ elements, and then the rest of the involution is just an involution on the remaining $n-2$ elements.

So we have the recurrence relation $a_n = a_{n-1} + (n-1)a_{n-2}.$ This is good for a practical computation for small $n.$

Now let $\displaystyle A(z) = \sum_{n=0}^{\infty} a_n \frac{z^n}{n!}$ be the Exponential Generating Function of the sequence. Replace $a_n$ by $a_{n-1}+(n-1)a_{n-2}$ into this sum and differentiate both sides to get the equation $A' = (1+z)A.$ Solving this gives $A(z) = \exp( z + z^2/2).$

$$ A(z) = \sum_{m=0}^{\infty} \frac{ (z+z^2/2)^m }{m!}.$$

The coefficient of $z^n$ in this sum is $a_n/n! \ .$ We see that the terms with $m< \lfloor n/2 \rfloor$ and $m>n$ don't contribute to the $z^n$ coefficient.

For $m=n-k$ with $k=0,1,\ldots, \lfloor n/2 \rfloor,$ the term is $\dfrac{(z+z^2/2)^{n-k}}{(n-k)!}.$ You only get a $z^n$ term if you pick $z$ from $n-2k$ of the $n-k$ factors and $z^2/2$ from the other $k$ factors. Since there are $\binom{n-k}{k}$ ways to do this, the coefficient of $z^n$ in this term is $\binom{n-k}{k} \dfrac{1}{2^k (n-k)!} =\dfrac{1}{2^k k! (n-2k)!}$ .

Therefore $\displaystyle a_n = \sum_{k=0}^{ \lfloor n/2 \rfloor } \frac{n!}{2^k k! (n-2k)!}$ (which of course is the same as what Kaj wrote in his answer).

This is one of the first (if not the very first, I haven't got the book handy) examples in R.P. Stanley's Enumerative Combinatorics, Vol 1. There is no closed formula for the answer $a_n$ for a given value of$~n$, but better than to describe $a_n$ by a summation involving amongst others binomial coefficients, you have the power series expression $$ \sum_{n\in\Bbb N}a_n\frac{X^n}{n!}=\exp\Bigl(X+\frac{X^2}2\Bigr). $$ See also this answer. There are fast methods to find the coefficients of such power series expression. In the current case they lead to the recurrence $a_n=(n-1)a_{n-2}+a_{n-1}$ and initial values $a_0=a_1=1$ that will rapidly give you the initial part of the sequence $(a_n)_{n\in\Bbb N}=(1,1,2,4,10,26,76,232,764,\ldots)$ Note that this requires fewer arithmetic operations, and simpler ones with smaller numbers, to compute all of $(a_0,\ldots,a_n)$, than the summation expression $a_n=\sum_{k=0}^{\lfloor n/2 \rfloor}\binom{n}{2k}\frac{(2k)!}{k!2^k}$ requires to compute just$~a_n$.

The recurrence is also easily understood directly: if $n$ is a fixed point of an involution of $S_n$ then what remains is an involution of $S_{n-1}$, and otherwise $n$ is interchanged with one of the $n-1$ previous elements, and an involution of the remaining $n-2$ elements remains to be chosen.