zpp(g) = {sum(i=1, #g, my(f=factor(g[i])); sum(j=1, #f~, x^f[j, 1]*(y^f[j, 2]-1)))}
A316089(n) = {my(g=znstar(n).cyc); my(p=zpp(g), r=#g); prod(i=1, poldegree(p), my(u=Vec(r + polcoeff(p, i))); r!/prod(j=1, #u, u[j]!))}
zp(g) = {sum(i=1, #g, my(f=factor(g[i])); sum(j=1, #f~, x^f[j, 1]*y^f[j, 2]))}
aut(p, q)={my(s=0, d=0, m=1); forstep(i=poldegree(q), 1, -1, my(c=polcoeff(q, i)); if(c, s+=i*c*d + (i-1)*c*(d+c); m*=prod(i=1, c, p^i-1); d+=c)); s+=d*(d-1)/2; m*p^s}
A258615(n) = {my(p=zp(znstar(n).cyc)); prod(i=1, poldegree(p), aut(i, polcoeff(p, i)))}
A302257(n) = A258615(n)*A316089(n)/(#znstar(n)[2])!
A327790(n) = prod(i=1,#znstar(n)[2],eulerphi(znstar(n)[2][i]))
a(n) = A302257(n)/A327790(n)