When the product of dice rolls yields a square Succinct Question:
Suppose you roll a fair six-sided die $n$ times.
What is the probability that the product of the rolls is a square?
Context:
I used this as one question in a course for elementary school teachers when $n=2$, and thought the generalization might be a good follow-up question for secondary school math teachers. But I encountered quite a bit of difficulty in tackling it, and I am wondering if there is a neater solution than what I have already seen, and to what deeper phemonena it connects.
Known:
Since the six sides of a die are $1, 2, 3, 2^2, 5,$ and $2\cdot3$, the product of the rolls is always of the form $2^{A}3^{B}5^{C}$, and the question is now transformed into the probability that $A, B, C$ are all even. The actual "probability" component is mostly for ease of phrasing; its only contribution is a $6^n$ in the denominator, and my true question is of a more combinatorial nature: namely,

In how many ways can the product of $n$ rolls yield a square?

One approach that I have seen involves first creating an $8 \times 8$ matrix corresponding to the eight cases around the parity of $A, B, C$; one can then take the dot product of each roll with this matrix, and hope to spot a pattern. In this way, one may discover the formula:
$$\frac{6^n + 4^n + 3\cdot2^n}{8}$$
and the "probability" version is simply this formula with another $6^n$ multiplied in the denominator.
As for proving this: Some guesswork around linear combinations of the numerator yields a formula for each of the eight cases concerning $A,B,C$ parity, and one can then prove all eight of them by induction. And so I "know" the answer in the sense that I have all eight of the formulae (and the particular one listed above is correct) but they were not found in a particularly organized fashion.

My Actual Question:
What is a systematic way to deduce the formula, given above, for the number of ways the product of $n$ rolls yields a square, and to what deeper phenomena does this connect?

 A: 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*}

A: Using diagonalization on Wolfram Alpha, I was able to confirm your result.  I used the matrix 
$$M=\left(\begin{array}{cccccccc} 
2&1&1&1&1&0&0&0\\
1&2&1&0&1&1&0&0\\ 
1&1&2&0&1&0&1&0 \\ 
1&0&0&2&0&1&1&1 \\ 
1&1&1&0&2&0&0&1\\
0&1&0&1&0&2&1&1\\ 
0&0&1&1&0&1&2&1\\ 
0&0&0&1&1&1&1&2
 \end{array} \right)$$
which gives one step transitions between square root parts of products (no square root part , $\sqrt{2}$, $\sqrt{3}$, $\sqrt{5}$, $\sqrt{6}$, $\sqrt{10}$, $\sqrt{15}$, $\sqrt{30}$ respectively).
For the number you require, you want the first entry of $M^n\ \vec{b}$ where $\vec{b}=\left(\begin{array}{c}1\\0\\0\\0\\0\\0\\0\\0\end{array}\right)$
As mentioned above, using diagonalization with Wolfram Alpha (the diagonalizing matrix was actually quite nice--integer entries; the inverse had denominators of eighths), I was able to confirm your result.  
I don't know if this method is an improvement over your description.  I hope it helps!
A: There's a slick approach to this based on bijections, though it loses a lot of the generality of the generating function methods.
Let $S=\{1,2,3,6\}$ and $T=\{4,5\}$.  We will divide roll sequences into classes based on whether the sequence contains elements from $S, T,$ or both.
Class 1: Sequences consisting only of rolls in $T$.  Swapping the first die roll between $4$ and $5$ gives a bijection between squares and non-squares, so exactly half the sequences in this class give a square.
Class 2: Sequences consisting only of rolls in $S$.  Now we can divide the sequences into groups of $4$ that share the same last $n-1$ rolls.  Each group has one square product, so exactly $1/4$ the sequences in this class give a square.
Class 3: Sequences containing both a roll in $S$ and a roll in $T$.  To each sequence we assign a "type", consisting of (1): The location of the first roll in $S$, (2): the location of the first roll in $T$, and (3): The remaining $n-2$ rolls.  Once the type is fixed, there's $8$ choices for the remaining roll, and exactly one of them gives a square.  
So the number of square sequences is 
$$\frac{1}{2} |\textrm{Class } 1| + \frac{1}{4} |\textrm{Class } 2| + \frac{1}{8} |\textrm{Class } 3| = \frac{1}{2} 2^n + \frac{1}{4} 4^n + \frac{1}{8} (6^n-2^n-4^n)$$
A: 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.
