(GAP)
p := 3;
a := function(n)
local q, k, cnt, x;
q:=p^n; k:=GF(p, n); cnt:=0;
for x in k do
if Trace(k, GF(p), x)=0*Z(p) and Order(x)=q1 then
cnt := cnt+1;
fi;
od;
return cnt/n;
end;
for n in [1..16] do Print (a(n), ", "); od;
(Sage) # much more efficient
p=3; # choose characteristic
for n in xrange(1, 66):
F=GF(p^n, 'x');
g = F.multiplicative_generator(); # generator
vt = vector(ZZ, p); # stats: trace
m = p^n  1; # size of multiplicative group
## Compute all irreducible polynomials via Lyndon words:
for w in LyndonWords(p, n): # digits of Lyndon words range form 1, .., p
e = sum( (w[j]1) * p^j for j in xrange(0, n) )
if gcd(m, e) == 1: # primitive elements only
f = g^e;
t = f.trace().lift(); # trace (over ZZ)
vt[t] += 1;
print vt[0]; # choose index 0, 1, .., p1 for different traces
# Joerg Arndt, Oct 03 2012
