When the product of dice rolls yields a square

For $1\le i\le6,\;$ let $a_i$ be the number of dice which have the digit $i$ appearing.

The product of the rolls will be a perfect square when $a_2+a_6,\;$ $a_3+a_6,\;$ and $a_5$ are all even;

so we can consider two cases:

$\textbf{1)}$ When $a_2, a_3, a_6$ are all odd, we get the exponential generating function

$\;\;\;\displaystyle\underbrace{\big(1+x+\frac{x^2}{2!}+\cdots\big)^2}_{a_1, a_4}\underbrace{\big(x+\frac{x^3}{3!}+\frac{x^5}{5!}+\cdots\big)^3}_{a_2, a_3, a_6}\underbrace{\big(1+\frac{x^2}{2!}+\frac{x^4}{4!}+\cdots\big)}_{a_5}$

$\;\;\;\displaystyle=e^{2x}\left(\frac{e^x-e^{-x}}{2}\right)^3\left(\frac{e^x+e^{-x}}{2}\right)=\color{red}{\frac{1}{16}\big(e^{6x}-2e^{4x}-e^{-2x}+2\big)}$

$\textbf{2)}$ When $a_2, a_3, a_6$ are all even, we get the exponential generating function

$\;\;\;\displaystyle\underbrace{\big(1+x+\frac{x^2}{2!}+\cdots\big)^2}_{a_1,a_4}\underbrace{\big(1+\frac{x^2}{2!}+\frac{x^4}{4!}+\cdots\big)^4}_{a_2,a_3, a_5, a_6}$

$\;\;\;\displaystyle=e^{2x}\left(\frac{e^x+e^{-x}}{2}\right)^4=\color{red}{\frac{1}{16}\big(e^{6x}+4e^{4x}+6e^{2x}+e^{-2x}+4\big)}$

Adding the two cases gives the generating function

$\;\;\;\displaystyle g_e(x)=\frac{1}{16}\big[2e^{6x}+2e^{4x}+6e^{2x}+6\big]=\color{red}{\frac{1}{8}\big[e^{6x}+e^{4x}+3e^{2x}+3\big]}$

$\hspace{.3 in}\displaystyle=1+\frac{1}{8}\sum_{n=1}^{\infty}\left(6^n+4^n+3\cdot2^n\right)\frac{x^n}{n!},\;\;$ so there are

$\displaystyle \hspace{.5 in}\color{blue}{\frac{1}{8}\big(6^n+4^n+3\cdot2^n\big)}$ ways to roll $n$ dice and get a product which is a perfect square.

You can bypass the cumbersome transition matrix method in this problem, and get a quicker explicit solution using generating functions and some algebra tricks.

Your goal is to count the number of 2s, 3s, and 5s that occur in the prime factorization $ (2^A 3^B 5^C)$ of the product of N rolls of a die. The exponents of these primes increase in an additive fashion with each roll. We can model this by treating the exponents $(A,B,C)$ as position vectors on a three-dimensional lattice.

Step 1. Each outcome of each roll of the die has a prime factorization, and has the incremental effect of shifting your position on the lattice by a displacement vector. Below we tabulate the possible outcomes of a roll of a die, and the exponents that occur in the prime factorization of the number rolled).

1: (0,0,0) 2:(1,0,0) 3:(0,1,0) 4:(2,0,0) 5:(0,0,1) 6:(1,1,0)

Since we want to keep track only of the parity of these exponents we can regard (2,0,0) as equivalent to (0,0,0). Note that there are now two ways to roll this null displacement vector, and one way to roll the other displacement vectors.

Step 2. We now use algebraic expressions as a notational device for modeling a displacement on a lattice. To each of the outcomes listed above we use its prime factorization to construct an associated algebraic expression. (Admittedly it seems like we are chasing our tails, toggling back and forth between additive and multiplicative descriptions of numbers, but be patient. ) :)

$1\to 2^0 3^0 5^0 \to a^0 b^ 0 c^0 $ (null displacement on exponent lattice)

$ 2\to 2^1 3^0 5^0 \to a^1 b^0 c^0 = a$

$\ldots$

$ 4\to 2^2 3^0 5^0 \to a^2 \to a^0$ (null displacement) because in our case we care only about parity on the exponent lattice

$\ldots$

$ 6\to a*b$

Step 3. Now we are ready to introduce generating functions! To account for all six outcomes of a dice roll, set $F(a,b,c)= 2+ a + b + a*b + c$ which correctly double-counts the null outcome. The magic of generating functions as a book-keeping device is that for example the 3rd power of the multivariable polynomial $F(a,b,c)$ tells you the number of ways that three rolls of the die can add up to a given displacement vector on the exponent lattice. That is, the expansion of the third power $F(a,b,c)^3=8 + 12 a + 6 a^2 + a^3 + 12 b + 24 a b + 15 a^2 b + 3 a^3 b + 6 b^2 + 15 a b^2 + 12 a^2 b^2 + 3 a^3 b^2 + b^3 + 3 a b^3 + 3 a^2 b^3 + a^3 b^3 + 12 c + 12 a c + 3 a^2 c + 12 b c + 18 a b c + 6 a^2 b c + 3 b^2 c + 6 a b^2 c + 3 a^2 b^2 c + 6 c^2 + 3 a c^2 + 3 b c^2 + 3 a b c^2 + c^3$

tells you that in three rolls there are $6$ ways to get the algebraic expression $a^2 b c$ (which corresponds to the displacement vector $(2,1,1)$ on the exponent lattice, which in turn corresponds to increasing the cumulative product of numbers rolled by the factor $2^2 \times 3 \times 5$.)

Now, since you want to keep track of the general $n_{th}$ power of $F(a,b,c)$ you can introduce the geometric series $g(a,b,c;t)= \frac{1}{ 1- t F} = 1+ t F+ t^2 F^2 + t^3 F^3 \ldots$ whose expansion in powers of $t$ consists of the powers of $F$.

Step 4. Recall that you want to lump together all lattice positions that have the same parity mod 2. A sneaky trick for accomplishing that "topological identification" of congruent lattice points is to symmetrize the function $g$ to make it even in each of the variables $(a,b,c)$, so that only even powers of $a,b,c$ occur in its Taylor expansion. This symmetrized function is an eight-term expression $S(a,b,c;t)= \frac{1}{8}[g(a,b,c;t)+ g(-a,b,c;t)+ g(a,b,c) + g( a,-b, c) +\ldots]$

In principle we could expand $8S(a,b,c;t) $ out as a Taylor series: $8 + 16 t + 8 (4 + a^2 + b^2 + a^2 b^2 + c^2) t^2 + (64 + 48 a^2 + 48 b^2 + 96 a^2 b^2 + 48 c^2) t^3 +\ldots$

but actually all we care about is the coefficient-weighted total number of different terms that occur for each power of $t$.

Step 5. This latter can be found by the simple trick of replacing $a=1, b=1, c=1$ which is a massive simplification. So go back and do that from the very beginning, where $g$ is introduced, repeat the symmetrization process, insert the numerical values $a=1,b=1,c=1$ and collect terms in an expression that now depends only on $t$. We see that $ 8 S(a,b,c;t) =3 + 1/(1 - 6 t) + 1/(1 - 4 t) + 3/(1 - 2 t)$.

Step 6. Now at last it is clear (by expanding the last expression in its geometric series) why you get the answer $\frac{6^n + 4^n + 3* 2^n}{8}$.

Note: This answer can be regarded as supplement to the nice answer of @MathWonk. Here we put a strong focus on generating functions.

Intro: A typical representation of one roll of a six-sided die is given by \begin{align*} x^1+x^2+x^3+x^4+x^5+x^6 \end{align*} The exponents of $x$ represent the pips of the die, the coefficients the number of occurrences of the respective event. Since we want to count the number of squares in $n$ rolls, we also keep track of the prime factors $2,3$ and $5$ which occur in the numbers $1,\ldots,6$. We use the variables $a,b$ and $c$ to mark the number of these prime factors. We obtain the generating function \begin{align*} x+ax^2+bx^3+a^2x^4+cx^5+abx^6 \end{align*} The variable $a$ represents the occurrence of $2$, $b$ represents $3$ and $c$ the prime $5$. Since the number $6=2\cdot3$ we count the prime factor $2$ and $3$ by multiplying the term $x^6$ with $a b$. This is similarly done for all other faces of the die.

We also want to keep track of the number of rolls, so we introduce a variable $t$ and multiply each term with it. This way we can define as basic building block the generating function \begin{align*} A(a,b,c;t;x)=(x+ax^2+bx^3+a^2x^4+cx^5+abx^6)t \end{align*} A generating function representing $n\geq 1$ rolls is \begin{align*} A_n(a,b,c;t;x)&:= \left(A(a,b,c;t;x)\right)^n\\ &=(x^1+ax^2+bx^3+a^2x^4+cx^5+abx^6)^nt^n \end{align*}

In fact these are only introductory notes, giving some background knowledge. We can instead start with

Main part: Let $A_n(a,b,c;t;x)$ be a generating function of a six-sided die representing $n$ rolls, which keeps track of the prime factors $2,3$ and $5$ of the pips and the number of rolls. It is given for $n\geq 1$ by \begin{align*} A_n(a,b,c;t;x)&=(x^1+ax^2+bx^3+a^2x^4+cx^5+abx^6)^nt^n\\ &=t^n\sum_{{i_1+i_2+\ldots+i_6=n}\atop{i_j\geq 0,1\leq j \leq 6}}\binom{n}{i_1,i_2,\ldots,i_6} x^{i_1+2i_2+\ldots+6i_6}a^{i_2+2i_4+i_6}b^{i_3}c^{i_5}\tag{1} \end{align*} with $\binom{n}{i_1,i_2,\ldots,i_6}=\frac{n!}{i_1!i_2!\cdots i_6!}$ the multinomial coefficients.

Since we want to count the rolls giving square numbers we are looking for a generating function $B_n(a,b,c;t;x)$, which is based upon $A_n(a,b,c;t;x)$ but additionally fulfills, that the exponents of $a,b$ and $c$ are even. In fact, this was the rationale for introducing these variables.

In order to obtain even exponents of $a,b$ and $c$ we need according to the representation in (1)

\begin{align*} i_2+2i_4+i_6&\equiv 0(2)\\ i_3&\equiv 0(2)\tag{2}\\ i_5&\equiv 0(2) \end{align*}

Now recall, that each function $f(x)$ can be represented as sum of an even and odd function via \begin{align*} f(x)&=f_e(x)+f_o(x)\\ &=\frac{f(x)+f(-x)}{2}+\frac{f(x)-f(-x)}{2} \end{align*} The even part $G_e(x)$ of a generating function $G(x)=\sum_{n=0}^{\infty}g_nx^n$ contains even powers of $x$ only, since \begin{align*} G_e(x)=\frac{G(x)+G(-x)}{2}=\sum_{n=0}^{\infty}g_{2n}x^{2n} \end{align*}

We need according to (2) an even generating function in the variables $a,b$ and $c$ which leads to \begin{align*} B_n(a,b,c;t;x)=\frac{1}{8}&\left(A_n(a,b,c;t;x)+A_n(-a,b,c;t;x)\right.\\ &+A_n(a,-b,c;t;x)+A_n(-a,-b,c;t;x)\\ &+A_n(a,b,-c;t;x)+A_n(-a,b,-c;t;x)\tag{3}\\ &\left.+A_n(a,-b,-c;t;x)+A_n(-a,-b,-c;t;x)\right)\\ \end{align*}

Note, we need the variables $a,b$ and $c$ for the derivation of the appropriate generating function $B_n(a,b,c;t;x)$. We don't need the variables to count the number of occurrences of squares. We also don't need to differentiate the pips, so we also don't need $x$ any longer.

We simply need the variable $t$ which counts the number of rolls and we want to add up all terms for a specific $t^n$. This way we count all occurrences of squares in $n$ rolls.

We obtain: The generating function $C_n(t)$ representing all occurrences of squares when rolling a die $n$ times is $(n\geq 1):$ \begin{align*} C_n(t)&=B_n(1,1,1;t;1)\\ &=\frac{1}{8}\left(A_n(1,1,1;t;1)+A_n(-1,1,1;t;1)\right.\\ &\qquad+A_n(1,-1,1;t;1)+A_n(-1,-1,1;t;1)\\ &\qquad+A_n(1,1,-1;t;1)+A_n(-1,1,-1;t;1)\\ &\qquad\left.+A_n(1,-1,-1;t;1)+A_n(-1,-1,-1;t;1)\right)\\ &=\frac{1}{8}\left((6t)^n+(2t)^n+(2t)^n+(2t)^n+(4t)^n+0+0+0\right)\\ &=\frac{1}{8}\left(6^n+4^n+3\cdot2^n\right)t^n \end{align*}

Note, it's convenient to use the coefficient of operator $[x^n]$ to denote the coefficient of $x^n$ of a series.

We conclude, the number of occurrences of squares when rolling a die $n$ times and multiplying the resulting pips is the coefficient of $t^n$ of the generating function $C_n(t)$ \begin{align*} [t^n]C_n(t)=\frac{1}{8}\left(6^n+4^n+3\cdot2^n\right)\qquad\qquad n\geq 1 \end{align*}