I believe the recurrence should be $f(n)=n(f(n-1)+1)$ with $f(2)=2.$ You can absorb the $(-1)^{j+k}$ in the addition of the products, so each term just takes one multiply. This gives $0,2,9,40,205,1236$, which is OEIS A038156 and matches your author's expression except that the upper limit of the sum should be $n-1$, not $n$. Another formula is $f(n)=\lfloor n!(e-1) \rfloor -1$
We can verify your author's formula from the recurrence by induction. The base case $f(2)=2=2!\frac 1{1!}$ works. Now assume we have $f(n)=n!\sum_{k=1}^{n-1}\frac 1{k!}$ We have $$\begin {align} f(n+1) &=(n+1)\left(f(n)+1\right)\\
&=(n+1)\left(n!\sum_{k=1}^{n-1}\frac 1{k!}+1\right)\\
&=(n+1)!\left(\sum_{k=1}^{n-1}\frac 1{k!}+ \frac 1{n!}\right)\\
&==(n+1)!\left(\sum_{k=1}^{n}\frac 1{k!}\right) \end {align}$$
OK, so this is a weird proof, but there's really only two parts of this proof that actually matter here. First:
[...] right multiplication by an elementary matrix is a column operation, and effect of column operations on the determinant is well known.
What this means is that when you multiply $A$ by $E$ on the right, you are doing a column operation on $A$, whether that be switching two columns, doing a column addition, or multiplying a column by a scalar. All of these operations will multiply the determinant by some known number, $k$ (i.e. switching two columns multiplies determinant by $k=-1$, column addition multiplies determinant by $k=1$, and multiplying a column by the scalar $c$ multiplies determinant by $k=c$). Thus, for any matrix $A$, the determinant of $AE$ will be $k$ times the original determinant of $A$. We can write this in an equation as follows:
$$\det AE=k(\det A)$$
Second:
[...] for a column operation the corresponding elementary matrix can be obtained from the identity matrix I by this column operation. So, its determinant is 1 (determinant of I) times the effect of the column operation.
Now, this is really confusing at first, but it can be understood in terms of our $\det AE=k(\det A)$ above. See, this equation works for any matrix $A$, which means we could also substitute the identity matrix $I$ for $A$ into this equation. Therefore, we get:
$$\det IE=k(\det I)\rightarrow \det E=k$$
(Note that this uses the fact that $IE=E$ and $\det I=1$.)
Therefore, we now know the value of $k$ is $\det E$. Thus, we can substitute that back into our original equation $\det AE=k(\det A)$ to get:
$$\det AE=(\det E)(\det A)$$
Since scalar multiplication is commutative, we can switch the right side around to get the final lemma:
$$\det AE=(\det A)(\det E)$$
Best Answer
Usually you define the determinant via
$$\det(A)=\sum_{\sigma\in S_n}\bigg(sgn(\sigma)\prod_{i=1}^nA_{i,\sigma(i)}\bigg)$$
So the inner product
$$\prod_{i=1}^n A_{i, \sigma(i)}$$
does exactly $n$ multiplications. Multiplying by $sgn(\sigma)$ is $O(1)$ so all in all
$$O\bigg(sgn(\sigma)\prod A_{i, \sigma(i)}\bigg)=O(n)$$
for a fixed $\sigma$. I also assume that calculating $sgn(\sigma)$ is $O(1)$. Now $S_n$ has exactly $n!$ elements so
$$O\bigg(\sum_{\sigma\in S_n}sgn(\sigma)\prod A_{i, \sigma(i)}\bigg)=O\bigg(\sum_{\sigma\in S_n}n\bigg)=O\bigg(n\cdot\sum_{\sigma\in S_n}1\bigg)=O(n\cdot n!)$$
which finally gives $O(n\cdot n!)$, not $O(n!)$.
Although I didn't prove that it is not $O(n!)$ (see: the big-O notation ) it actually is not $O(n!)$ because $\det(A)$ does precisely $n!\cdot(n+1)$ arithmetic operations (I ignore the fact that you have to generate the set of all permutations $S_n$ because it can be done in $O(n!)$) by the former argumentation. The $n!$ comes from summation over $S_n$, the $n$ from inner $\prod$ product and $+1$ from multiplying by $sgn(\sigma)$. Or in other words $\det(A)$ is $\Theta(n\cdot n!)$ - the big-Theta notation.
Also note the problem with the big-O notation: constant function is $O(n!)$. Does it mean it is slow? No. Whoever gave you the exercise most likely meant some other growth describing notation, e.g. the big-Theta.