Let's consider approximation of $\prod_{p=1}^n N(p,a)$, $n\to \infty$, where the number of fixed necklaces of length n composed of $a$ types of beads $N(n,a)$: $\frac {(a-1)^{n+1}} {(a-3) \cdot n!} \prod_{p=1}^n \frac {1-a^p} {1-a}$, $a > 3$.
Please see Error term in formula for products of necklaces at MO.
My question is how to derive the such approximation and what the error estimation is?
I am interested in that due to $\frac {1-a^p} {1-a}$ looks like a generation functions of Mahonian triangles (A008302, OEIS) and I would like to understand why so?