An unexpected application of non-trivial combinatorics PROBLEM STATEMENT

Given two finite sets $A$ and $B$, each containing $s \in \mathbb N$ elements, how many pairs of functions $f \colon A \rightarrow B$ and $g \colon B \rightarrow A$ are there, such that their composition $h = g \circ f$ has no fixed points?

THE ORIGIN OF THE PROBLEM

It is highly improbable that two people could happen to like each other by chance.

Which, under the following assumptions:


*

*We are considering a group of $2n$ people with $n$ men and $n$ women;

*Each person likes exactly one person of the opposite sex with equal probability;


reduces to

Given two finite sets $A$ and $B$, each containing $n \in \mathbb N$ elements, and two arbitrary functions $f \colon A \rightarrow B$ and $g \colon B \rightarrow A$, what is the expected number of pairs of elements $(a \in A, b \in B)$ such that $f(a) = b$ and $g(b) = a$?

WHAT I HAVE DONE
Here I am using $E$ for expectation, $P$ for probability and $N$ for number.
$$E(\text{number of pairs}) = \sum\limits_{k=0}^n k \cdot P(\text{there are exactly }k\text{ pairs})$$
$$P(\text{there are exactly }k\text{ pairs}) = \frac {N(\text{situations with }k\text{ pairs})} {N(\text{total situations})}$$
$$N(\text{total situations}) = (n^n)^2,$$
because there are $n^n$ functions from $A$ to $B$ and the same number in the reverse direction.
$$ \begin{aligned} 
N & (\text{situations with }k\text{ pairs}) = \\
 & = N(\text{pairs of functions }(f \colon A \rightarrow B, g \colon B \rightarrow A) \\
 & \quad \text{ whose composition has exactly }k\text{ fixed points}) = \\
 & = \sum [\text{over selections of }k\text{ elements from both }A\text{ and }B] \\
 & \quad N(\text{pairs of functions whose composition is identity only on the selection}) = \\
 & = {\binom {n} {k}}^2 \cdot k! \cdot N(\text{pairs of functions on two sets containing } \\
 & \quad s = n - k\text{ elements whose composition has no fixed points})
\end{aligned} $$
This is how I reduced the original problem to the problem in question. If the solution to that problem was $N(s)$, the answer to the original problem would be
$$E(\text{number of pairs}) = \frac {1} {(n^n)^2} \sum\limits_{k=0}^n k \cdot {\binom {n} {k}}^2 \cdot k! \cdot N(n - k)$$
WHAT I TRIED TO DO
Approach: fix $g$, count $f$s.
$$ \begin{multline}
N(s) = \sum [\text{over all functions }g \colon B \rightarrow A] \\
N(\text{functions }f \colon A \rightarrow B\text{ such that the composition }h = g \circ f \\
\text{ has no fixed points})
\end{multline} $$
Problem: all such $f$s have to send every $a$ to anywhere but $g^{-1}(a)$. Therefore, there are $\prod\limits_{a \in A} (n - \left|g^{-1}(a)\right|)$ such functions. The sets $\left\{g^{-1}(a)\right\}_{a \in A}$ form a partition of $B$. The obvious next idea is to classify all $g$s based on that partition, but the formulas start to get really scary.
UPD: I wrote a simulation that actually counted directly the expectation for $n=1..5$ (the $O(n\cdot(n^n)^2)$ approach). The expected value is always exactly 1. I have no idea how to interpret that, given the crazy formula for expectation.
I would mostly appreciate an idea for an approach, not a complete solution.
 A: To begin with, the expected number $E(H)$  of hits (first choice matchings) is $1$, and more is true: One  has
$$E(H\,|\, f)=1$$ 
for every individual $f:\>A\to B$.
Proof. Assume that an $f$ is given. We have a hit at $b\in B$ when $g(b)\in f^{-1}(\{b\})$. The probability  that this is the case is equal to
$$P(b\,|\, f):={\#f^{-1}(\{b\})\over n}\ .$$
Therefore $E(H\,|\, f)$ becomes
$$E(H\,|\, f)=\sum_{k=1}^n P(b_k\,|\,f)={1\over n}\sum_{k=1}^n \#f^{-1}(\{b_k\})=1\ .$$
When $n=5$ there are  $5^{10}=9\,765\,625$ cases in all. Mathematica has found the following probability distribution of $H$:
$$(0.283124, 0.469117, 0.213893, 0.0323789, 0.00147456, 0.000012288)\ .$$
These values are to be interpreted as follows: With $28.3\%$ probability we have no hits, with $46.9\%$ probability exactly one hit, and so on. The last value is $=5!/5^{10}$ and comes from bijective $f$'s with $g=f^{-1}$.
(One might have conjectured that $H$ is approximatively Poisson distributed with mean $1$, but this seems not  to be the case.)
A: Here's a not-carefully checked attempt to solve a simpler, related question: What is the probability that there is at least one fixed point. If it's wrong, I hope it at least helps. Perhaps you can then compute the probability of at least two, at least three, etc., and then find the expected value you’re looking for.
For $J\subseteq N=\{1,\dots,n\}$, let $P_J$ be the probability that for a random pair of functions $f,g:N\to N$, $J$ is contained in the set of points fixed by $g\circ f$.
If $|J|=1$, WLOG $J=\{1\}$, and $P_J=\dfrac{1}{n}$. ($g$ must take $f(1)$ to $1$).
If $|J|=2$, WLOG $J=\{1,2\}$. For the set of fixed points to include $J$, it must be the case that $f(1)\ne f(2)$ (probability $\dfrac{n-1}{n}$) and $g(f(i))=i$ for $i=1,2$ (probability $\dfrac{1}{n^2}$. So $P_J=\dfrac{n-1}{n^3}$.
If $J=\{1,\dots,k\}$, it must be the case that $f(1),f(2),\dots,f(k)$ are distinct (probability $\dfrac{n!}{(n-k)!\,n^k}$) and $g(f(i))=i$ for $i=1\dots k$ (probability $\dfrac{1}{n^k}$. In this case, $P_J=\dfrac{n!}{(n-k)!\,n^{2k}}$.
By inclusion-exclusion, the probability of at least one fixed point is
$${n\choose1}P_{\{1\}}-{n\choose2}P_{\{1,2\}}+{n\choose3}P_{\{1,2,3\}}\cdots \pm {n\choose n}P_{\{1,\dots,n\}} = \sum_{k=1}^n(-1)^{k+1}{n\choose k} \dfrac{n!}{(n-k)!\,n^{2k}}\textrm.$$
Mathematica evaluates the sum as $1 - n^{-2 n} U(-n, 1, n^2)$, where $U(a,b,z)$ is the confluent hypergeometric function $\dfrac{1}{\Gamma(a)}\int_0^\infty e^{-zt}t^{a-1}(1+t)^{b-a-1}\,dt$. Mathematica doesn’t find its limit as $n\to\infty$, but it sure looks like it could be $1-\dfrac{1}{e}$, which suggests there might be a much easier formulation of this as a derangement problem.

A: As Steve points out, this can be reformulated as a more standard derangement problem. We can produce all pairs of functions $f \circ g$ as follows. Choose $f \circ g = h$ a derangement of $N =  \{ 1, 2 , \dots, n\}$. Then let $f$ be any bijection you like. This forces $g$ completely. Moreover, all pairs $f, g$ arise like this as clearly both $f$ and $g$ must be bijections.
So the total number of such pairs is just 
$$n! \left[ \frac{n!}{e} \right]$$
where the brackets denote the `nearest integer' function.
Edit: This is wrong, thanks to Steve for pointing out my misread of question. Let's try again...
Edit 2: The below is also false unfortunately. In particular I don't consider the case that $h|K$ is not a bijection. Nor do I count the ways that $h$ can be decomposed. The idea might be salvageable, but will probably be somewhat complicated.
Fix some $h = f \circ g$ and denote its image by $K = \{ 1, \dots, k\}$.
If we restrict the domain of $f \circ g|_K:K \to K$ to coincide with its image, then it is a bijection. By the above analysis, there are 
$$k! \left[ \frac{k!}{e} \right]$$
such maps. Of course, we can set the rest of $h$ in any way we like. That is, for each $x \in N \setminus K$ we can choose $h(x) \in K$ arbitrarily. Thus, there are 
$$k! \left[ \frac{k!}{e} \right] \cdot (n-k)^k $$
pairs. We're not done though, since we began by making a choice of $Im h$. Let's sum over every possible such choice
$$ \sum_{k=1}^n {n \choose k} k! \left[ \frac{k!}{e} \right] \cdot (n-k)^k $$
