cp's OEIS Frontend

A226313 Number of commuting 4-tuples of elements from S_n, divided by n!.

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