Flajolet and Sedgewick generating function for Hertzsprung Problem The Hertzsprung Problem goes as follows: In how many can we place exactly $n$ non-attacking kings on a $n \times n$ chessboard such that there is exactly $1$ king in each row and column where $n \in \mathbb{N}$.
My main question is: How did Flajolet and Sedgewick get the below generating function for the Hertzsprung problem?
$$\sum_{n=0}^\infty n! x^n \frac{(1-x)^n}{(1+x)^n}$$
Flajolet and Sedgewick discuss this generating function briefly and somewhat vaguely on page $373$ in Enumerative Combinatorics. They give a sketch of a sketch. However I'm totally not sure how they got the generating function. I have tried to extract the same generating function for a long time using restricted permutation, I mean the Hertzsprung problem reduces to finding number of permutations of $[n]$ such that no two adjacent digits in the permutation are consecutive.
Moreover using the restricted permutations argument we get a nice closed form as follows
$$n!+\sum_{k=1}^n {(-1)^k}(n-k)!\sum_{i=1}^k 2^{i} \binom{k-1}{i-1}\binom{n-k+1}{i}$$
by simply using Principle of inclusion-exclusion and stars and bars argument.
However I don't know how Flajolet and Sedgewick got the generating function
$$\sum_{n=0}^\infty n! x^n \frac{(1-x)^n}{(1+x)^n}$$
Any idea how to approach the problem for finding the generating function?  I have already searched all the references in (oeis.org/A002464) but no reference gives a proof on how they got the generating function. All the references just show or give an approach on how to get the closed form of the double sum which is very easy to get.
Your help would be highly appreciated.
Thanks.

Here is my guess at what Flajolet and Sedgewick had in mind.
We want to find a generating function for the number of permutations $\sigma$ of $[1..n]$ such that $|\sigma_{i+1} - \sigma_i| \ne 1$.  We start by solving an associated problem:  in preparation to apply inclusion / exclusion, we say that a permutation $\sigma$ of $[1..n]$ has "property $i$" if $|\sigma_{i+1} - \sigma_i| = 1$ for $1 \le i < n$, and we define $s_{j,n}$ to be the number of permutations with $j$ of the properties. We define the bivariate generating function of $s_{j,n}$ by
$$f(u,z) = \sum_{n \ge 0} \sum_{j \ge 0} s_{n,j} u^j z^n$$
Once we have $f(u,z)$, it's easy to find the generating function for the number of permutations which have none of the properties. By inclusion / exclusion, the GF is simply $f(-1,z)$.
Consider an arbitrary permutation of $[1..n]$.  We can break the permutation into a sequence of chunks where each chunk is either a single integer or consists of two or more consecutive integers in ascending or descending order.  For example, here is a permutation of $[1..12]$ with seven chunks: four single integers, one chunk of length $2$ and two chunks of length $3$.
$$9\; 2\; \overleftarrow{\boxed{12\; 11\; 10}}\; 8\; \overrightarrow{\boxed{5\; 6\; 7}} \; 1\; \overrightarrow{\boxed{3\; 4}}$$
How many permutations consist of $m$ chunks?  If the lengths of the chunks are $t_1, t_2, t_3, \dots, t_m$, then we must have
$$t_1+t_2+t_3+\dots+t_m=n$$
with $t_i \ge 1$ for all $i$.  Taking into account that each chunk of length greater than $1$ may be either ascending or descending and the chunks can be ordered in $m!$ ways, a GF for the number of permutations of $[1..n]$ having $m$ chunks is
$$m!\; (z + 2z^2 + 2z^3 + 2z^4 + \dots)^m$$
Any permutation may be broken into $m$ chunks for exactly one value of $m$, so a GF for the number of permutations of $[1..n]$ is
$$\sum_{m \ge 0} m!\; (z + 2z^2 + 2z^3 + 2z^4 + \dots)^m$$
Observing that a chunk of length $j$ contributes $j-1$ to $s_{j,n}$, we can modify the previous GF to produce the bivariate GF
$$f(u,z) = \sum_{m \ge 0} m!\; (z + 2u z^2 + 2 u^2 z^3 + 2u^3 z^4 + \dots)^m$$
Applying the formula for the sum of an infinite geometric series and simplifying,
$$f(u,z) = \sum_{m \ge 0} m!\; \left( z + \frac{2 u z^2}{1 - uz} \right)^m$$
We are now in a position to apply  inclusion / exclusion to find the GF for the number of permutations $\sigma$ of $[1..n]$ with no instances of $|\sigma_{i+1} - \sigma_i| = 1$:
$$f(-1,z) = \sum_{m \ge 0} m!\; \left( z\; \frac{1-z}{1+z} \right)^m $$
