%I #14 Mar 29 2022 15:10:57
%S 1,1,4,20,139,1192,12318,148318,2041754,31616757,544005172,
%T 10296204096,212589150300,4755177958104,114545293676588,
%U 2956316416222300,81386676426000157,2380590235918735576,73729207700492304684,2410324868012471929670,82944575892433740648996
%N a(0) = 1; a(n) = Sum_{k=1..n} binomial(n,k) * sigma_0(k) * a(n-k).
%H Seiichi Manyama, <a href="/A340903/b340903.txt">Table of n, a(n) for n = 0..413</a>
%F E.g.f.: 1 / (1 - Sum_{i>=1} Sum_{j>=1} x^(i*j) / (i*j)!).
%F E.g.f.: 1 / (1 - Sum_{k>=1} sigma_0(k) * x^k/k!). - _Seiichi Manyama_, Mar 29 2022
%t a[0] = 1; a[n_] := a[n] = Sum[Binomial[n, k] DivisorSigma[0, k] a[n - k], {k, 1, n}]; Table[a[n], {n, 0, 20}]
%t nmax = 20; CoefficientList[Series[1/(1 - Sum[Sum[x^(i j)/(i j)!, {j, 1, nmax}], {i, 1, nmax}]), {x, 0, nmax}], x] Range[0, nmax]!
%o (PARI) my(N=40, x='x+O('x^N)); Vec(serlaplace(1/(1-sum(k=1, N, numdiv(k)*x^k/k!)))) \\ _Seiichi Manyama_, Mar 29 2022
%o (PARI) a(n) = if(n==0, 1, sum(k=1, n, numdiv(k)*binomial(n, k)*a(n-k))); \\ _Seiichi Manyama_, Mar 29 2022
%Y Cf. A000005, A000670, A129921, A295739, A340904.
%K nonn
%O 0,3
%A _Ilya Gutkovskiy_, Jan 26 2021