Generating functions for tail length and rho-length I am trying to obtain generating functions for tail length and rho length of a random point in a random mapping. 
Let $\phi:\{1,2,\ldots,n\}\to \{1,2,\ldots,n\}$ be a random function. Consider the directed graph whose nodes are the elements $\{1,2,\ldots,n\}$ and whose edges are the ordered pairs $(x,\phi(x))$, for all $x\in \{1,2,\ldots,n\}$. We start from any $u_0$ and keep iterating $\phi$, i.e. we consider the sequence $u_1=\phi(u_0), u_2=\phi(u_1), \ldots$. In fact, starting from any $u_0$, the iteration structure of $\phi$ is described by a simple path that connects to a cycle. The length of this path (measured by the number of edges) is called the tail length of $u_0$ and is denoted by $\lambda(u_0)$. The length of the cycle (measured by the number of edges or nodes) is called the cycle length of $u_0$ and is denoted by $\mu(u_0)$. We also call rho-length of $u_0$ the quantity $\rho(u_0)=\lambda(u_0) +\mu(u_0)$
In the paper entitled "Random mapping statistics" and in the Theorem 3, the authors obtained the expectations for these parameters. How we can obtain generating functions for these parameters?
 A: Cycle size plus tail length.
We will show that with these parameters being cumulative we are indeed
justified in adding the generating functions for the two cases that we
considered separately above. This matches  what is being stated in the
cited paper.
Using the  same notation  as in  the companion answer  we get  for the
generating function
$$G(z, v) = \exp\left(\sum_{q\ge 1} \frac{T(v^q z, v)^q}{q}\right).$$
where  $T(z,v)$  is  the  generating  function from  the  tail  length
computation.
This yields for  the generating function $H(z)$ of  the expected cycle
size plus tail length same as in the other two answers
$$H(z) = \left.\frac{\partial}{\partial v} G(z,v)\right|_{v=1}.$$
This works out to
$$H(z) = \left. \exp\left(\sum_{q\ge 1} \frac{T(v^q z, v)^q}{q}\right)
\\ \times
\left(\sum_{q\ge 1} 
\frac{qT(v^q z)^{q-1} 
\left(qz v^{q-1}\frac{\partial}{\partial z} T(z,v) +
\frac{\partial}{\partial v} T(z, v)\right)}{q}\right)
\right|_{v=1}.$$
Put $$S(z) = 
\left.\frac{\partial}{\partial v} T(z,v)
\right|_{v=1}$$
and use the functional equation to obtain
$$\left.\frac{\partial}{\partial v} T(z,v)
\right|_{v=1}
= \left.z\exp T(vz, v)
\left(z \frac{\partial}{\partial z} T(z,v)+
\frac{\partial}{\partial v} T(z,v)
\right)\right|_{v=1}
\\ = T(z) (z T'(z) + S(z))
= z T(z) T'(z) + T(z) S(z)$$
so that
$$S(z) = \frac{z T(z) T'(z)}{1-T(z)}
= \frac{T(z)^2}{(1-T(z))^2}.$$
Entering this into the equation for $H(z)$ we find
$$H(z) = \exp\left(\sum_{q\ge 1} \frac{T(z)^q}{q}\right)
\left(\sum_{q\ge 1} T(z)^{q-1} 
\left(T'(z)q z + \frac{T(z)^2}{(1-T(z))^2}\right)\right)
\\ = 
\left(\frac{z T'(z)}{(1-T(z))^2} 
+ \frac{T(z)^2}{(1-T(z))^3}\right)
\exp\log\frac{1}{1-T(z)}
\\ = 
\left(\frac{T(z)+T(z)^2}{(1-T(z))^3}\right)
\exp\log\frac{1}{1-T(z)}
\\ = \frac{T(z)+T(z)^2}{(1-T(z))^4}.$$
This is the sum of the two generating functions as predicted above.
We get for the closed form
$$Q_n = 
n! \sum_{q=0}^{n-2} \frac{n^q}{q!} (n-q)^2
+ n! \frac{n^{n-1}}{(n-1)!}
\\ = n^n +
n! \sum_{q=0}^{n-2} \frac{n^q}{q!} (n-q)^2.$$
The sequence starts with
$$ 1, 12, 153, 2272, 39225, 776736, 17398969, 435538944, 
\\ 12058401393, 366021568000,\ldots$$
The Maple code for this was as follows where the reader is asked to
note certain  optimizations that have  been made in comparison  to the
first answer.

Q :=
proc(n)
    option remember;
    local ind, d, gf, pos, q, x, seen, traj;

    if n = 1 then return v fi;

    gf := 0;

    for ind from n^n to 2*n^n-1 do
        d := convert(ind, base, n);
        d := map(l->l+1, [seq(d[q], q=1..n)]);

        q := 0;
        for pos to n do
            x := pos; traj := [];

            do
                for seen to nops(traj) do
                    if traj[seen] = x then break fi;
                od;

                if seen <= nops(traj) then break fi;

                traj := [op(traj), x];
                x := d[x];
            od;

            q := q + nops(traj);
        od;

        gf := gf+v^q;
    od;


    gf;
end;

EX :=
proc(n)
    local T;

    T := solve(TT=z*exp(TT), TT);
    n!*coeftayl((T+T^2)/(1-T)^4, z=0, n);
end;

EX2 :=
proc(n)
    n! * residue(exp(w*n)*(1/w^(n-1)+1/w^n)*1/(1-w)^3, w=0);
end;

EX3 :=
proc(n)
    n^n + n! * add(n^q/q!*(n-q)^2, q=0..n-2);
end;

