|
%I M0114 N0044 #27 Jan 29 2022 01:02:47
%S 2,1,2,3,4,8,8,18,18,38,28,142,72,234,360,669,520,2606,1608,7338,8856,
%T 19370,16768,94308,67556,216200,277512,815310,662368,4499852,2311468,
%U 8465496,13045076,31592762,40937592,159769394,103197488,401912086
%N Number of equivalence classes of binary sequences of primitive period n.
%C The number of equivalence classes of primitive sequences of period p, taking values in a set with b elements, is given by: N'(p) = sum_{d|p} mobius(p/d)*N(d) where N denotes the number of equivalence classes in the set of all sequences with period p, taking b values (see A002729). - Pab Ter (pabrlos2(AT)yahoo.com), Oct 22 2005
%D N. J. A. Sloane, A Handbook of Integer Sequences, Academic Press, 1973 (includes this sequence).
%D N. J. A. Sloane and Simon Plouffe, The Encyclopedia of Integer Sequences, Academic Press, 1995 (includes this sequence).
%D R. C. Titsworth, Equivalence classes of periodic sequences, Illinois J. Math., 8 (1964), 266-270.
%H Vincenzo Librandi, <a href="/A002730/b002730.txt">Table of n, a(n) for n = 1..50</a>
%H <a href="/index/Lu#Lyndon">Index entries for sequences related to Lyndon words</a>
%F Reference gives formula.
%p with(numtheory): E:=proc(k,L) if(L=1) then RETURN(1) else RETURN(order(k,L)) fi end; M:=proc(k,L) local s,EkL: EkL:=E(k,L): if(k>1) then s:=(k^EkL-1)/(k-1): RETURN(L*EkL/igcd(L,s)) else RETURN(L*EkL/igcd(L,EkL)) fi end; C:=proc(k,t,p) local u: RETURN(add(M(k,p/igcd(p,u*(k-1)+t))^(-1),u=0..p-1)) :end; N:=proc(p) options remember: local s,t,k: if(p=1) then RETURN(2) fi: s:=0: for t from 0 to p-1 do for k from 1 to p-1 do if igcd(p,k)=1 then s:=s+2^C(k,t,p) fi od od: RETURN(s/(p*phi(p))):end; Nprimitive:=proc(p) options remember: local d: RETURN(add(mobius(p/d)*N(d),d=divisors(p))): end; seq(Nprimitive(p),p=1..51); (Pab Ter)
%t max = 38; m[k_, n_] := (s = 1; Do[ If[ Mod[s, n] == 0, Return[e], s = s + k^e ] , {e, 1, max}]); c[k_, t_, n_] := Sum[ m[k, n/GCD[n, u*(k-1) + t]]^(-1), {u, 0, n-1}]; (* b = A002729 *) b[n_] := b[n] = (s = 0; Do[ If[ GCD[n, k] == 1 , s = s + 2^c[k, t, n]] , {k, 1, n-1}, {t, 0, n-1}]; s / (n*EulerPhi[n]) ); b[0] = 1; b[1] = 2; a[n_] := Sum[ MoebiusMu[n/d]*b[d], {d, Divisors[n]}]; Table[a[n], {n, 1, max}] (* _Jean-François Alcover_, Dec 06 2011, after Maple *)
%Y Cf. A002729.
%K nonn,easy,nice
%O 1,1
%A _N. J. A. Sloane_
%E More terms from Pab Ter (pabrlos2(AT)yahoo.com), Oct 22 2005
|
