How many hash collision will exist where the digest is a pre-image? This question follows on from the question Circular hash collision where digest is pre-image: Could hash(a) = b, hash(b)=c, then hash(c)=a?

The question was answered but the issue of the number of collisions was also raised. I have moved my answer to this query as a separate question. The other question begins:

For any SHA-family of algorithms (including Keccak) except for SHA1, would it be theoretically possible for there to exist a pre-image a whose hash digest b is the pre-image to the hash digest c which is the pre-image to a thus completing a circle (i.e. circular hash or collision of pre-image a with the digest a)? 

${{hash(A)\to B}}$, ${{hash(B)\to C}}$ , ${{hash(C)\to A}}$

The answer to this is yes. My question is, how many such collisions would typically exist?
The expected number of 3-cycles in a random function on a large domain is $1/3$. More generally, for a random function on a set of size $n$, the expected number of cycles of length $\ell$ is asymptotically equivalent to $1/\ell$ as $n \to \infty$.

This result, along with a general approach to the combinatorical properties of random functions, can be found in the paper Random Mapping Statistics by Flajolet and Odlyzko (likely, the result has been derived elsewhere using different methods too).

Let me sketch how this result is obtained. Recall that the functional graph of a mapping is a set of components, each of which consists of a cycle to which trees are connected.

Let $t(z)$ be the exponential generating function for the number of (rooted) trees. It can be shown that (this is known as Cayley's formula)

$$t(z) = \sum_{n = 1}^\infty \frac{n^{n - 1}}{n!} z^n.$$

Since components are cycles of trees, the generating function for the number of components is

$$\log\left(\frac{1}{1-t(z)}\right).$$

Also, the generating function for the number of functional graphs would then be the exponential of the above. In order to count the number of cycles of length $\ell$, one can construct a bivariate generating function for the number of functional graphs in which the second variable "marks" $\ell$-cycles. That is, instead of the above generating function for the number of components, we define

$$\log\left(\frac{1}{1-t(z)}\right) + (u- 1)\frac{t(z)^\ell}{\ell}.$$

That is, the term (in the Taylor series in the variable $t(z)$) corresponding to cycles of length $\ell$ is multiplied by the marker variable $u$. The desired bivariate generating function at the level of the functional graph is then

$$\xi(z, u) = \exp\left(\log\left(\frac{1}{1-t(z)}\right) + (u- 1)\frac{t(z)^\ell}{\ell}\right).$$

One can see that the generating function for the sum of number of cycles of length $\ell$ over all functional graphs of fixed size can then be computed as

$$\Xi(z) = \frac{\partial\xi(z, u)}{\partial u}\Bigg|_{u = 1} = \frac{z^\ell}{\ell(1 - t(z))}.$$

It remains to determine the Maclaurin series coefficients of the above, at least for terms of very high order. This is possible by singularity analysis. We use the fact (by Proposition 1 in the paper of Flajolet and Odlyzko) that (as $z \to 1/e$)

$$1 - t(z) \sim \sqrt{2(1 - ez)}.$$

Hence,

$$\Xi(z) \sim \frac{z^\ell}{\ell\sqrt{2(1 - ez)}}.$$

The singularity analysis technique (Theorem 1 in Flajolet and Odlyzko) then shows that (with $\Xi_n / n!$ the coefficient of order $n$)

$$\Xi_n \overset{\star_1}{\sim} e^n\,\frac{n!}{\sqrt{2}\ell}\,\frac{\sqrt{n}}{n\Gamma(1/2)} \sim \frac{e^n n!}{\ell\sqrt{2\pi n}} \overset{\star_2}{\sim} \frac{n^n}{\ell},$$

where $\star_1$ is the actual singularity analysis and $\star_2$ follows from Stirling's approximation $n! \sim \sqrt{2\pi n}\, n^n/e^n$. In order to obtain the expected value of the number of $\ell$-cycles of a random function on a set of size $n$, it suffices to divide by the number of such functions ($n^n$): $\Xi_n \sim 1/\ell$.