|
%I #23 Apr 30 2019 08:22:30
%S 1,1,0,1,1,3,13,24,67,252,1795,4038,16812,61750,349806,3485026,
%T 10391070,49433135,240064988,1282012986,9167581934,131550811985,
%U 459677212302,2707382738558,14318807586215,94084166753923,601900541189696,5894253303715121
%N Number of primitive (aperiodic, or Lyndon) necklaces with n beads of unlabeled colors such that the numbers of beads per color are distinct.
%H Alois P. Heinz, <a href="/A261599/b261599.txt">Table of n, a(n) for n = 0..300</a>
%H F. Ruskey, <a href="http://combos.org/necklace">Necklaces, Lyndon words, De Bruijn sequences, etc.</a>
%H F. Ruskey, <a href="/A000011/a000011.pdf">Necklaces, Lyndon words, De Bruijn sequences, etc.</a> [Cached copy, with permission, pdf format only]
%H Eric Weisstein's World of Mathematics, <a href="http://mathworld.wolfram.com/Necklace.html">Necklace</a>
%H Wikipedia, <a href="https://en.wikipedia.org/wiki/Lyndon_word">Lyndon word</a>
%H Wikipedia, <a href="https://en.wikipedia.org/wiki/Necklace_(combinatorics)">Necklace (combinatorics)</a>
%H <a href="/index/Ne#necklaces">Index entries for sequences related to necklaces</a>
%F a(n) = (1/n) * Sum_{d | n} moebius(n/d) * A007837(d) for n>0. - _Andrew Howroyd_, Dec 21 2017
%e a(4) = 1: 0001.
%e a(5) = 3: 00001, 00011, 00101.
%e a(6) = 13: 000001, 000011, 000101, 000112, 000121, 000122, 001012, 001021, 001022, 001102, 001201, 001202, 010102.
%e a(7) = 24: 0000001, 0000011, 0000101, 0000111, 0000112, 0000121, 0000122, 0001001, 0001011, 0001012, 0001021, 0001022, 0001101, 0001102, 0001201, 0001202, 0010011, 0010012, 0010021, 0010022, 0010101, 0010102, 0010201, 0010202.
%p with(numtheory):
%p b:= proc(n, i, g, d, j) option remember; `if`(i*(i+1)/2<n or g>0
%p and g<d, 0, `if`(n=0, `if`(d=g, 1, 0), b(n, i-1, g, d, j)+
%p `if`(i>n, 0, binomial(n/j, i/j)*b(n-i, i-1, igcd(i, g), d, j))))
%p end:
%p a:= n-> `if`(n=0, 1, add(add((f-> `if`(f=0, 0, f*b(n$2, 0, d, j)))(
%p mobius(j)), j=divisors(d)), d=divisors(n))/n):
%p seq(a(n), n=0..30);
%t a[0] = 1; a[n_] := With[{P = Product[1 + x^k/k!, {k, 1, n}] + O[x]^(n+1) // Normal}, DivisorSum[n, MoebiusMu[n/#]*#!*Coefficient[P, x, #]&]/n];
%t Table[a[n], {n, 0, 30}] (* _Jean-François Alcover_, May 28 2018, after _Andrew Howroyd_ *)
%o (PARI) a(n)={if(n==0, 1, my(p=prod(k=1, n, (1+x^k/k!) + O(x*x^n))); sumdiv(n, d, moebius(n/d)*d!*polcoeff(p, d))/n)} \\ _Andrew Howroyd_, Dec 21 2017
%Y Cf. A007837, A261531, A261600.
%K nonn
%O 0,6
%A _Alois P. Heinz_, Aug 25 2015
|
