Generating functions for tail length and rho-length

Let me point out that the cited paper is a landmark brilliant event containing profound insights that surpass the MSE question answer format.

As an example we do the calculation of the expected cycle length $\mu(u_0)$ and the tail length $\lambda(u_0)$ using the labelled tree function. We refer the reader to the paper as regards the asymptotics.

We will provide a closed form of the exponential generating function of the quantities that are involved.

The species of labelled trees has the specification $$\def\textsc#1{\dosc#1\csod} \def\dosc#1#2\csod{{\rm #1{\small #2}}}\mathcal{T} = \mathcal{Z} \times \textsc{SET}(\mathcal{T})$$ which gives the functional equation $$T(z) = z \exp T(z).$$

We have that $$T(z) = \sum_{n\ge 1} n^{n-1} \frac{z^n}{n!}.$$

Cycle length.

To compute the expected cycle length note that these random mappings are sets of cycles of trees having combinatorial specification

$$\textsc{SET} \left(\sum_{q\ge 1} \textsc{CYC}_{=q}(\mathcal{T}(\mathcal{Z}))\right).$$

Now a tree on $n$ nodes that goes into a cycle of size $q$ contributes $nq$ to the total count so we get the species $$\textsc{SET} \left(\sum_{q\ge 1} \textsc{CYC}_{=q}(\mathcal{T}(\mathcal{V}^q\mathcal{Z}))\right).$$

It follows that the bivariate generating function of mappings by node count and cycle size is $$G(z, v) = \exp\left(\sum_{q\ge 1} \frac{T(v^q z)^q}{q}\right).$$

This yields for the generating function $H(z)$ of the expected cycle size $$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)^q}{q}\right) \left(\sum_{q\ge 1} \frac{qT(v^q z)^{q-1} T'(v^q z)q v^{q-1} z}{q}\right) \right|_{v=1} \\ = \exp\left(\sum_{q\ge 1} \frac{T(z)^q}{q}\right) \left(\sum_{q\ge 1} T(z)^{q-1} T'(z)q z\right) \\ = \frac{z T'(z)}{(1-T(z))^2} \exp\log\frac{1}{1-T(z)} \\ = \frac{z T'(z)}{(1-T(z))^3}.$$

We get from the functional equation $$T'(z) = \exp(T(z)) + z \exp(T(z)) T'(z)$$ so that $$T'(z) = \frac{\exp(T(z))}{1-z\exp(T(z))} = \frac{T(z)/z}{1-zT(z)/z} = \frac{1}{z} \frac{T(z)}{1-T(z)}$$ and therefore $$H(z) = \frac{T(z)}{(1-T(z))^4}.$$

Extracting coefficients via Lagrange inversion we have $$Q_n = n! [z^n] H(z) = n! \frac{1}{2\pi i} \int_{|z|=\epsilon} \frac{1}{z^{n+1}} \frac{T(z)}{(1-T(z))^4} dz.$$

Put $T(z)=w$ so that $z=w/\exp(w) = w\exp(-w)$ and $dz = \exp(-w) - w\exp(-w)$ to get $$n! \frac{1}{2\pi i} \int_{|w|=\gamma} \frac{\exp(w(n+1))}{w^{n+1}} \frac{w}{(1-w)^4} (\exp(-w) - w\exp(-w)) \; dw \\ = n! \frac{1}{2\pi i} \int_{|w|=\gamma} \frac{\exp(wn)}{w^{n+1}} \frac{w}{(1-w)^4} (1 - w) \; dw \\ = n! \frac{1}{2\pi i} \int_{|w|=\gamma} \frac{\exp(wn)}{w^{n}} \frac{1}{(1-w)^3} \; dw.$$

This gives the closed form $$Q_n = \frac{1}{2} n! \sum_{q=0}^{n-1} \frac{n^q}{q!} (n+1-q) (n-q).$$

The average in question is $$\frac{Q_n}{n^n \times n}.$$

The sequence $Q_n$ starts with $$1, 10, 117, 1648, 27425, 528336, 11581885, \\ 284878336, 7772592897, 233010784000,\ldots$$

The Maple code for verification by total enumeration of this sequence was as follows.

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

    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
            seen := {}; x := pos; traj := [];

            while not(x in seen) do
                traj := [op(traj), x];
                seen := seen union {x};

                x := d[x];
            od;

            cycinit := 1;
            while traj[cycinit] <> x do
                cycinit := cycinit + 1;
            od;

            q := q + nops(traj)-cycinit+1;
        od;

        gf := gf+v^q;
    od;


    gf;
end;

EX :=
proc(n)
    local T;

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

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

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

Addendum Mon Sep 22 2019. A better way to extract the coefficients from $H(z)$ is not to compute the derivative of $T(z)$ and use

$$Q_n = n! [z^n] H(z) = n! \frac{1}{2\pi i} \int_{|z|=\epsilon} \frac{1}{z^{n+1}} \frac{z T'(z)}{(1-T(z))^3} \; dz \\ = n! \frac{1}{2\pi i} \int_{|z|=\epsilon} \frac{1}{z^{n}} \frac{T'(z)}{(1-T(z))^3} \; dz.$$

Put $T(z)=w$ so that $z=w/\exp(w) = w\exp(-w)$ and $dw = T'(z) \; dz$ to get

$$n! \frac{1}{2\pi i} \int_{|w|=\gamma} \frac{\exp(nw)}{w^{n}} \frac{1}{(1-w)^3} \; dw.$$

We then continue as before. The standard Lagrange inversion does not help here, we get

$$Q_n = (n-1)! [z^{n-1}] H'(z) \\ = (n-1)! \frac{1}{2\pi i} \int_{|z|=\epsilon} \frac{1}{z^{n}} \left( \frac{1}{(1-T(z))^4} + \frac{4T(z)}{(1-T(z))^5} \right) T'(z) \; dz \\ = (n-1)! \frac{1}{2\pi i} \int_{|w|=\gamma} \frac{\exp(nw)}{w^{n}} \left( \frac{1}{(1-w)^4} + \frac{4w}{(1-w)^5} \right) \; dw.$$

This yields

$$(n-1)! \sum_{q=0}^{n-1} \frac{n^q}{q!} {n-1-q+3\choose 3} + 4(n-1)! \sum_{q=0}^{n-2} \frac{n^q}{q!} {n-2-q+4\choose 4} \\ = (n-1)! \sum_{q=0}^{n-1} \frac{n^q}{q!} {n-q+2\choose 3} + 4(n-1)! \sum_{q=0}^{n-2} \frac{n^q}{q!} {n-q+2\choose 4}.$$

With

$${n-q+2\choose 3} + 4 {n-q+2\choose 4} = {n-q+2\choose 3} + (n-q-1) {n-q+2\choose 3}$$

we find

$$(n-1)! \sum_{q=0}^{n-1} \frac{n^q}{q!} (n-q) {n-q+2\choose 3}.$$

Tail length.

This requires a modification to the tree function in order to count the total tail length for all nodes in the tree.

We obtain the functional equation $$T(z, v) = z\sum_{q\ge 0} \frac{T(vz,v)^q}{q!} = z\exp(T(vz, v)).$$

This is because when we attach a set of trees to the new root all nodes in the tree have their tail length incremented by one.

We then have the a simple species for average tail lengths: $$\textsc{SET} \left(\textsc{CYC}(\mathcal{T})\right).$$

It follows that the bivariate generating function of mappings by node count and tail length is $$G(z, v) = \exp\log\frac{1}{1-T(z,v)} = \frac{1}{1-T(z,v)}.$$

This yields for the generating function $H(z)$ of the expected tail length (same as before) $$H(z) = \left.\frac{\partial}{\partial v} G(z,v)\right|_{v=1}.$$

This works out to $$H(z) = \left.\frac{1}{(1-T(z,v))^2} \frac{\partial}{\partial v} T(z,v) \right|_{v=1} = \frac{1}{(1-T(z,v))^2} \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} \\ = \frac{T(z)}{(1-T(z))^2} (z T'(z) + H(z) (1-T(z))^2) \\ = \frac{z T(z) T'(z)}{(1-T(z))^2} + H(z) T(z).$$

This yields $$H(z) = \frac{z T(z) T'(z)}{(1-T(z))^3} = \frac{T(z)^2}{(1-T(z))^4}.$$

The Lagrange inversion computation can be carried out as before with an extra $w$ in the integrand and it produces $$Q_n = \frac{1}{2} n! \sum_{q=0}^{n-2} \frac{n^q}{q!} (n-q) (n-1-q).$$

The sequence $Q_n$ starts with $$0, 2, 36, 624, 11800, 248400, 5817084, 150660608, \\ 4285808496, 133010784000,\ldots$$

The Maple code for verification by total enumeration of this sequence was as follows.

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

    if n = 1 then return 1 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
            seen := {}; x := pos; traj := [];

            while not(x in seen) do
                traj := [op(traj), x];
                seen := seen union {x};

                x := d[x];
            od;

            cycinit := 1;
            while traj[cycinit] <> x do
                cycinit := cycinit + 1;
            od;

            q := q + cycinit-1;
        od;

        gf := gf+v^q;
    od;


    gf;
end;

EX :=
proc(n)
    local T;

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

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

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

The labelled tree function recently appeared at this MSE link and at this MSE link II which is closely related to the present subject.

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;