### general linear group finite field

I am trying to compute the number of 3 x 3 matrices A over the field of p elements (p prime) such that A^2 = I.
I wrote the following code which gives me the answers for p=2 and p=3, but it is taking for ever to compute when for p=5. Is there a better way to do this job which can speed up the computations? Many thanks for help.

def num_invol(n):
G = GL(n,GF(3));
sum = 0
for A in G:
if A^2 == G.one():
sum = sum +1
return sum

for p in (2..11):
if p.is_prime():
print (p, num_invol(p))