|
%I #33 May 09 2023 11:41:02
%S 1,8,21,84,206,717,1810,5462,13859,38497,96113,253206,620480,1566292,
%T 3770933,9212041,21768608,51795427,120279052,279849177,639379257,
%U 1459282932,3283758256,7369471795,16351101855,36147590987,79162129897,172646751524,373527250619,804631686843,1721283389932,3666041417241
%N Number of commuting 4-tuples of elements from S_n, divided by n!.
%C Euler transform of A001001.
%H Alois P. Heinz, <a href="/A226313/b226313.txt">Table of n, a(n) for n = 1..9000</a>
%H Tad White, <a href="http://arxiv.org/abs/1304.2830">Counting Free Abelian Actions</a>, arXiv preprint arXiv:1304.2830, 2013
%F a(n) ~ exp(2^(7/4) * Pi^(3/2) * Zeta(3)^(1/4) * n^(3/4) / (3^(3/2) * 5^(1/4)) - sqrt(5*Zeta(3)*n) / (2^(3/2)*Pi) + (sqrt(Pi) * 5^(1/4) / (2^(15/4) * 3^(3/2) * Zeta(3)^(1/4)) - sqrt(3) * 5^(5/4) * Zeta(3)^(3/4) / (2^(15/4) * Pi^(7/2))) * n^(1/4) - 25*Zeta(3) / (16*Pi^6) + (5 - 2*Zeta(3)) / (192*Pi^2)) * Pi^(1/4) * Zeta(3)^(1/8) / (2^(13/8) * 3^(1/4) * 5^(1/8) * n^(5/8)). - _Vaclav Kotesovec_, Mar 26 2018
%p with(numtheory):
%p b:= proc(n) option remember; add(d*sigma(d), d=divisors(n)) end:
%p a:= proc(n) option remember; `if`(n=0, 1, add(add(d*
%p b(d), d=divisors(j)) *a(n-j), j=1..n)/n)
%p end:
%p seq(a(n), n=1..40); # _Alois P. Heinz_, Mar 06 2015
%t b[n_] := b[n] = DivisorSum[n, #*DivisorSigma[1, #]&];
%t a[n_] := a[n] = If[n == 0, 1, Sum[DivisorSum[j, #*b[#]&]*a[n-j], {j, 1, n}] /n];
%t Array[a, 40] (* _Jean-François Alcover_, Mar 27 2017, after _Alois P. Heinz_ *)
%t nmax = 40; Rest[CoefficientList[Series[Exp[Sum[Sum[Sum[d*DivisorSigma[1, d], {d, Divisors[k]}] * x^(j*k) / j, {k, 1, Floor[nmax/j] + 1}], {j, 1, nmax}]], {x, 0, nmax}], x]] (* _Vaclav Kotesovec_, Mar 31 2018 *)
%o (PARI) EulerT(v)={Vec(exp(x*Ser(dirmul(v, vector(#v, n, 1/n))))-1, -#v)}
%o seq(n)={my(v=vector(n,i,1)); for(k=1, 2, v=dirmul(v, vector(n,i,i^k))); EulerT(v)} \\ _Andrew Howroyd_, May 09 2023
%Y Column k=4 of A362826.
%Y Cf. A001001, A061256, A301777.
%K nonn
%O 1,2
%A _N. J. A. Sloane_, Jun 08 2013
|
