%I M3541 #34 Oct 21 2023 19:20:42
%S 1,4,18,32,160,324,1456,2048,13122,25600,117128,209952,913952,2119936,
%T 9447840,13107200,86093440,172186884,774840976,1310720000,6964002864,
%U 13718968384,62761410632,88159684608,557885504000,835308258304,5083731656658,8988257288192,45753584909920,89261680665600,411782264189296,564050001920000
%N Generalized Euler phi function (for p=3).
%C For n >= 2, a(n) is the number of n X n circulant invertible matrices over GF(3). - Yuval Dekel (dekelyuval(AT)hotmail.com), Aug 22 2003
%D N. J. A. Sloane and Simon Plouffe, The Encyclopedia of Integer Sequences, Academic Press, 1995 (includes this sequence).
%H J. T. B. Beard Jr. and K. I. West, <a href="http://dx.doi.org/10.1090/S0025-5718-1974-0364196-5">Factorization tables for x^n-1 over GF(q)</a>, Math. Comp., 28 (1974), 1167-1168.
%t p = 3; numNormalp[n_] := Module[{r, i, pp}, pp = 1; Do[r = MultiplicativeOrder[p, d]; i = EulerPhi[d]/r; pp *= (1-1/p^r)^i, {d, Divisors[n]}]; Return[pp]]; numNormal[n_] := Module[{t, q, pp}, t=1; q=n; While[0==Mod[q, p], q /= p; t += 1]; pp = numNormalp[q]; pp *= p^n/n; Return[pp]]; a[n_] := If[n==1, 1, n*numNormal[n]]; Array[a, 40] (* _Jean-François Alcover_, Dec 10 2015, after _Joerg Arndt_ *)
%o (PARI)
%o p=3; /* global */
%o num_normal_p(n)=
%o {
%o my( r, i, pp );
%o pp = 1;
%o fordiv (n, d,
%o r = znorder(Mod(p,d));
%o i = eulerphi(d)/r;
%o pp *= (1 - 1/p^r)^i;
%o );
%o return( pp );
%o }
%o num_normal(n)=
%o {
%o my( t, q, pp );
%o t = 1; q = n;
%o while ( 0==(q%p), q/=p; t+=1; );
%o /* here: n==q*p^t */
%o pp = num_normal_p(q);
%o pp *= p^n/n;
%o return( pp );
%o }
%o a(n)=if ( n==1, 1, n * num_normal(n) );
%o v=vector(66,n,a(n))
%o /* _Joerg Arndt_, Jul 03 2011 */
%Y Cf. A003473 (p=2), A192037 (p=5).
%K nonn
%O 1,2
%A _N. J. A. Sloane_
%E Terms > 86093440 from _Joerg Arndt_, Jul 03 2011