|
%I #42 Dec 12 2015 08:08:45
%S 1,1,1,2,1,6,1,8,3,10,1,36,1,14,15,128,1,648,1,800,21,22,1,13824,5,26,
%T 81,6272,1,972000,1,32768,33,34,35,3359232,1,38,39,20480000,1,
%U 96018048,1,247808,30375,46,1,4586471424,7,500000,51,1384448,1,204073344,55
%N a(n) is the denominator of c(n), where c(n) is calculated from Product_{i>=1}(1-c(i)*x^i) = exp(-(x^2)/(1-x))*(1-x).
%C c(1) = 1 and for n>1, c(n) satisfies Sum_{d|n} (1/d)*c(n/d)^d = 1 + 1/n.
%C c(p) = 1 for prime p and a(p) = 1 accordingly.
%p c := proc (n) option remember; 1+1/n-add(procname(n/d)^d/d, d = `minus`(numtheory:-divisors(n), {1})) end proc: c(1) := 1: a := denom(map(c, [`$`(1 .. 100)]));
%t nmax = 100; Remove[c]; Subscript[c, 1] = 1; Do[Subscript[c, k] = Subscript[c, k] /. (Flatten[Solve[SeriesCoefficient[E^(-x^2/(1 - x))*(1 - x), {x, 0, k}] == Coefficient[Expand[Product[1 - Subscript[c, i]*x^i, {i, 1, k}]], x^k], Subscript[c, k]]]), {k, 2, nmax}]; Table[Subscript[c, n], {n, 1, nmax}] // Denominator (* _Vaclav Kotesovec_, Dec 12 2015 *)
%o (PARI) lista(nn) = {vc = vector(nn); vc[1] = 1; for (n=2, nn, vc[n] = 1+1/n - sumdiv(n, d, if (d==1, 0, (vc[n/d]^d)/d)); print1(denominator(vc[n]), ", "););} \\ _Michel Marcus_, Nov 27 2015
%Y Cf. A259027 (numerators).
%K nonn,frac
%O 1,4
%A _Gevorg Hmayakyan_, Nov 26 2015
|
