After the first read of the post, my impression was also, that the next case concerns $GL_n(\mathbb F_3)$, and
i gave up for a while. The case of finding the elements $A$ with $A^2=1$ in
$$
G = GL_3(\mathbb F_q)\ ,
$$
$q$ being and *odd* prime power, can be restated and solved theoretically. After doing this, the general case can also be solved by the same lines. (Odd characteristic.)

Let us consider such an $A\in G$, a $3\times 3$ matrix for the beginning.

Since $-A\in G$ we can chose now or later the class we want to count.
(The $\det:G\to \mathbb{F}_q^\times$ with the image $\pm 1$ separates the cases. This is because our $n=3$ is odd.)

The characteristic polynomial $P$ of $A$ has three non-zero roots.

The minimal polynomial is $X^2-1$ or a divizor, i.e. one among $X\pm1$, (since $A^2=1$).

In case of three equal eigenvalues, we can consider the case with $1,1,1$ being the three ones.
Then $P=(X-1)(X-1)(X-1)$. Subtracting from this $(X^2-1)(X-1)$ we see that $2(X-1)(X-1)$ also anihilates $A$.
But $2$ is not zero, since $q$ was odd. So $(X-1)(X-1)$ anihilates $A$. Using again $X^2-1=(X-1)(X+1)$, we see that
$2(X-1)$ anihilates $A$. So $X-1$ does it, so $A=1$. (This is too pedestrian, but shows where the odd characteristic is needed. In characteristic two, all standard Jordan blocks are idempotents. The solution in this case would be to also consider conjugations of such block compositions. In odd characteristic we have a diagonal shape and less problems.)

The case with $-1,-1,-1$ reduces to the above one, by passing from $A$ to $-A$, so we get $A=-1$.

Let us consier all other cases. Then the eigenvalues of $A$ are the two divisors $\pm1$ of $X^2-1$
and a further value $a$, say. The set of eigenvalues $1,-1,a$ is stable with respect to $x\to1/x$.
So $a$ is either $-1$ or $+1$.
In both cases we will restrict to the half of the $A$'s by imposing the condition, that the eigenvalues are
$$
-1,1,1\ .
$$
Since the minimal polynomial $X^2-1=(X+1)(X-1)$ splits in two distinct factors, odd characteristic again,
we have a diagonalizable matrix $A$, so there exists an invertible matrix $S$ so that we can write:
$$
A = S^{-1}DS\ ,
$$
$D$ being the diagonal matrix with entries $-1,1,1$.
We let $S$ run in the space of all invertible matrices. When do we have the same writing
$$
T^{-1} D T = A = S^{-1} D S\ ,
$$
for $S,T$ invertible?
This means that $ST^{-1}$ commutes with $D$, so it is of the shape / is in the $G$--subgroup of matrices

```
* 0 0
0 * *
0 * *
```

which has
$$ (q-1)\cdot(q^2-1)(q^2-q)=F(q,1)\cdot F(q,2) $$
elements, where $F(q,n)$ is the product over $k$ with $0 \le k < n$ from the factors $q^n-q^k$.
$F$ is here inspired from the word factorial.

We may use the convention, that an empty product is one, so $F(q,0)=1$.

So $S,T$ are generating by conjugation as above the same $A$, iff they are in the right classes w.r.t. the
above subgroup. (The right classes to be considered are maybe the left sided ones.)

So we expect a number of $A$'a equal to
$$
2 + 2\frac {F(q,3)}{F(q,1)F(q,2)}
=2
+2p^2(p²+p+1)\ .
$$
And indeed:

```
def F(q,n):
return prod( [ q^n - q^k for k in [ 0 .. n-1 ] ] ) # with no checks
for p in [ 3,5 ]:
print "%s -> %s" % ( p, 2 + 2*F(p,3) / F(p,1) / F(p,2) )
print "... or simply"
for p in [ 3,5 ]:
print "%s -> %s" % ( p, 2 + 2 * p^2 * (p^2+p+1) )
```

And we get:

```
3 -> 236
5 -> 1552
... or simply
3 -> 236
5 -> 1552
```

For a general $n$ we have to implement over the field with (odd) $q$ elements the following lines,
some $q$--combinatorial / factorial formula:
$$
\sum_{0\le s, t\le n\ ,\ s+t=n}
\frac{F(q,n)}{F(q,s)F(q,t)}
=2+
\sum_{0< s, t< n\ ,\ s+t=n}
\frac{F(q,n)}{F(q,s)F(q,t)}
\ .
$$
For four, $n=4$, just an example, the terms in the sum correspond to diagonal matrices $D$ with diagonal
$(-1,-1,-1,-1)$, and $(-1,-1,-1,1)$ and $(-1,-1,1,1)$ and $(-1,1,1,1)$ and $(1,1,1,1)$, and the corresponding
cosets have to be taken with respect to diagonal block matrices for the splittings $4=4+0=3+1=2+2=1+3=0+4$, correspondingly

```
**** *000 **00 ***0 ****
**** 0*** **00 ***0 ****
**** 0*** 00** ***0 ****
**** 0*** 00** 000* ****
```

and it is also possible to construct the matrices $A$ after constructing coset representatives.