Where $\pi_1$, $\pi_2$, $\pi_3$ are elements of the symmetric group $S_n$, I'm interested in computing the following sum: $$ f(\pi_1,\pi_2,\pi_3) = \sum_{\widetilde{\pi}_1\in S_n}\sum_{\widetilde{\pi}_2\in S_n}\sum_{\widetilde{\pi}_3\in S_n}\delta(\widetilde{\pi}_1\pi_1\widetilde{\pi}^{-1}_1\widetilde{\pi}_2\pi_2\widetilde{\pi}^{-1}_2\widetilde{\pi}_3\pi_3\widetilde{\pi}^{-1}_3,\mathrm{id}). $$ Here $\mathrm{id}$ is the identity element of $S_n$ and $\delta$ is the Kronecker delta. Note that $f$ is manifestly invariant under conjugation of $\pi_1$, $\pi_2$, or $\pi_3$, so it is a function only of the cycle structures of $\pi_1$, $\pi_2$, and $\pi_3$. In particular, let $a^{(1)}_k$ be the number of cycles of size $k$ in the cycle decomposition of $\pi_1$. (Define $a^{(2)}_k$ and $a^{(3)}_k$ similarly.) I'm specifically after a formula for $f(\pi_1,\pi_2,\pi_3)$ in terms of the numbers $a^{(i)}_k$.
Something I've thought of: you can rewrite $f$ in terms of irreducible characters as $$ f(\pi_1,\pi_2,\pi_3)=n!^2\sum_q\frac{1}{d_q}\chi_q(\pi_1)\chi_q(\pi_2)\chi_q(\pi_3), $$ where $q$ labels irreducible characters $\chi_q$, and $d_q=\chi_q(\mathrm{id})$ is the dimension of the $q$-th representation. From here you can write $\chi_q$ in terms of the $a_k$ using character polynomials. This doesn't really satisfy what I'm trying to do, however, for two reasons: (1) it seems difficult to work out the character polynomials for large $n$, and (2) I'd really like a formula that works for all $n$.
To be clear about what I'm looking for, consider the generalization of $f$ to two arguments. You can get the following formula for the two-argument case $$ f(\pi_1,\pi_2) = \sum_{\widetilde{\pi}_1\in S_n}\sum_{\widetilde{\pi}_2\in S_n}\delta(\widetilde{\pi}_1\pi_1\widetilde{\pi}^{-1}_1\widetilde{\pi}_2\pi_2\widetilde{\pi}^{-1}_2,\mathrm{id}) = n!\prod_k a^{(1)}_k!k^{a^{(1)}_k}\delta_{a^{(1)}_k,a^{(2)}_k} . $$ I want a formula like this, but for the three-argument function $f(\pi_1,\pi_2,\pi_3)$.
In the end I'd be interested in the generalizations of $f$ to more than three arguments, but those can be reduced to the 3-argument case. For example $$ f(\pi_1,\pi_2,\pi_3,\pi_4)=\frac{1}{n!^2}\sum_{\pi\in S_n}f(\pi_1,\pi_2,\pi)f(\pi^{-1},\pi_3,\pi_4). $$
The equation $\pi_1\pi_2=e$ is equivalent to $\pi_1=\pi_2^{-1}$, which is nice in this context because inversion plays nicely with conjugacy classes (indeed, in $S_n$ it fixes conjugacy classes). The equation $\pi_1\pi_2\pi_3=e$ on the other hand is equivalent to $\pi_1\pi_2=\pi_3^{-1}$, which is not nice for us because multiplication does not play nice with conjugacy classes: it is difficult to describe (in general) how often a permutation of cycle type $\lambda$ is produced by multiplying permutations of cycle types $\mu$ and $\nu$.
Let $C(\pi)$ be the centralizer of $\pi$. Your first formula is
$$ f(\pi_1,\pi_2)=\begin{cases} |G||C(\pi)| & \pi_1\sim\pi_2 \\ 0 & \pi_1\not\sim\pi_2 \end{cases} $$
which is not hard to see by counting: pick anything for $\bar{\pi}_1$, then the valid $\bar{\pi}_2$ for $\bar{\pi}_1\pi_1\bar{\pi}_1^{-1}=\bar{\pi}_2\pi_2^{-1}\bar{\pi}_2^{-1}$ are in a coset of $C(\pi_2^{-1})$ (which is conjugate to $C(\pi_2)$ since $\pi_2^{-1}\sim\pi_2$). The explicit formula for the size of a centralizer is classical, indeed $C(\pi)$ is a direct product of wreath products $C_k\wr S_{c_k(\pi)}$, where $C_k$ is cyclic (generated by a $k$-cycle) and $c_k(\pi)$ is the number of $k$-cycles in $\pi$.
Let $K(\pi)$ be the conjugacy class of $\pi$, so $|C(\pi)||K(\pi)|=|S_n|$ by orbit-stabilizer. The map $G\to K(\pi_1)$ given by $\bar{\pi}_1\mapsto \bar{\pi_1}\pi_1\bar{\pi}_1^{-1}$ is a $|C(\pi_1)|$-to-$1$ map. Therefore we may rewrite
$$ f(\pi_1,\pi_2,\pi_3)=|C(\pi_1)||C(\pi_2)||C(\pi_3)|g(\pi_1,\pi_2,\pi_3), $$
$$ g(\pi_1,\pi_2,\pi_3):=\#\{(\sigma_1,\sigma_2,\sigma_3)\in K(\pi_1)\times K(\pi_2)\times K(\pi_3)\mid \sigma_1\sigma_2=\sigma_3\} $$
or more simply (note $K(\pi_3^{-1})=K(\pi_3)$),
$$ g(\pi_1,\pi_2,\pi_3)=\#\{(\alpha,\beta)\in K(\pi_1)\times K(\pi_2)\mid \alpha\beta\in K(\pi_3)\} . $$
Abusing notation, we may view $K(\pi)$ as the sum of permutations from $\pi$\s conjugacy class as an element of the group ring center $Z(\mathbb{Z}[G])$. Or $K_\lambda$, if we use cycle types $\lambda$. Indeed, these sums form a nice integral basis. Then $K_\mu K_\nu = \sum C_{\mu\nu}^{\lambda} K_\lambda$ for some structure constants $C_{\mu\nu}^{\lambda}$, called connection coefficients. Summing all the coefficients of elements of $K_{\lambda}$ on both sides of this equation yields $g(\mu,\nu)=C_{\mu\nu}^{\lambda}|K_{\lambda}|$.
Therefore it suffices to figure out these coefficients $C_{\mu\nu}^{\lambda}$, which is discussed in this MO thread.