The problem is:
For $f:X\rightarrow X$, where $X$ has $n$ distinct elements, how many different $f$ exist that meet $f(f(f(x)))=x$ for all $x\in X$?
I can calculate for small number of $n$'s by the following way (for example, $n=8$):
case-1) $f$ is an identity function ($f(x)=x$) for all $x$.
case-2) $f$ is an identity function for 5 elements in $X$ and for the rest of 3, it is cyclic, i.e. $f(a)=b,f(b)=c,f(c)=a$.
case-3) $f$ is an identity function for 2 elements in $X$ and for the rest of 6, it can be divided into two groups of each having 3, and each groups are cyclic.
As there are only two different ways of being cyclic (1,2,3 or 1,3,2), the number of functions for $n=8$ can be calculated as:
$$1+\binom85\cdot2+\binom82\binom63\frac12\cdot2^2=1233$$
I could also generalize it for any $n$ as:
$$n=3m+p\quad (p=0,1,2)$$ $$a_n=\sum_{k=0}^m\left(\binom{n}{n-3k}\cdot2^k\cdot\frac1{k!}\prod_{i=0}^{k-1}\binom{n-3k}{3}\right)$$
And then I found that $a_n$ are $1,1,3,9,21,81,351,1233,\dots$ that is the sequence in this link https://oeis.org/A001470/internal. In the link, it shows a simple recurrence equation of
$$a_n=a_{n-1}+(n-1)(n-2)a_{n-3}$$
But I can't figure out how one can reach to that recurrence equation.
So my question is "how one can reach to the above recurrence equation"?