-
import Data.List (genericIndex)
a000081 = genericIndex a000081_list
a000081_list = 0 : 1 : f 1 [1,0] where
f x ys = y : f (x + 1) (y : ys) where
y = sum (zipWith (*) (map h [1..x]) ys) `div` x
h = sum . map (\d -> d * a000081 d) . a027750_row
-- Reinhard Zumkeller, Jun 17 2013
-
N := 30; P := PowerSeriesRing(Rationals(),N+1); f := func< A | x*&*[Exp(Evaluate(A,x^k)/k) : k in [1..N]]>; G := x; for i in [1..N] do G := f(G); end for; G000081 := G; A000081 := [0] cat Eltseq(G); // Geoff Bailey (geoff(AT)maths.usyd.edu.au), Nov 30 2009
-
N := 30: a := [1,1]; for n from 3 to N do x*mul( (1-x^i)^(-a[i]), i=1..n-1); series(%,x,n+1); b := coeff(%,x,n); a := [op(a),b]; od: a; A000081 := proc(n) if n=0 then 1 else a[n]; fi; end; G000081 := series(add(a[i]*x^i,i=1..N),x,N+2); # also used in A000055
spec := [ T, {T=Prod(Z,Set(T))} ]; A000081 := n-> combstruct[count](spec,size=n); [seq(combstruct[count](spec,size=n), n=0..40)];
# a much more efficient method for computing the result with Maple. It uses two procedures:
a := proc(n) local k; a(n) := add(k*a(k)*s(n-1,k), k=1..n-1)/(n-1) end proc:
a(0) := 0: a(1) := 1: s := proc(n,k) local j; s(n,k) := add(a(n+1-j*k), j=1..iquo(n,k)); # Joe Riel (joer(AT)san.rr.com), Jun 23 2008
# even more efficient, uses the Euler transform:
with(numtheory): a:= proc(n) option remember; local d, j; `if`(n<=1, n, (add(add(d*a(d), d=divisors(j)) *a(n-j), j=1..n-1))/ (n-1)) end:
seq(a(n), n=0..50); # Alois P. Heinz, Sep 06 2008
-
s[ n_, k_ ] := s[ n, k ]=a[ n+1-k ]+If[ n<2k, 0, s[ n-k, k ] ]; a[ 1 ]=1; a[ n_ ] := a[ n ]=Sum[ a[ i ]s[ n-1, i ]i, {i, 1, n-1} ]/(n-1); Table[ a[ i ], {i, 1, 30} ] (* Robert A. Russell *)
a[n_] := a[n] = If[n <= 1, n, Sum[Sum[d*a[d], {d, Divisors[j]}]*a[n-j], {j, 1, n-1}]/(n-1)]; Table[a[n], {n, 0, 30}] (* Jean-François Alcover, Feb 17 2014, after Alois P. Heinz *)
a[n_] := a[n] = If[n <= 1, n, Sum[a[n - j] DivisorSum[j, # a[#] &], {j, n - 1}]/(n - 1)]; Table[a[n], {n, 0, 30}] (* Jan Mangaldan, May 07 2014, after Alois P. Heinz *)
(* first do *) << NumericalDifferentialEquationAnalysis`; (* then *)
ButcherTreeCount[30] (* v8 onward Robert G. Wilson v, Sep 16 2014 *)
a[n:0|1] := n; a[n_] := a[n] = Sum[m a[m] a[n-k*m], {m, n-1}, {k, (n-1)/m}]/(n-1); Table[a[n], {n, 0, 30}] (* Vladimir Reshetnikov, Nov 06 2015 *)
terms = 31; A[] = 0; Do[A[x] = x*Exp[Sum[A[x^k]/k, {k, 1, j}]] + O[x]^j // Normal, {j, 1, terms}]; CoefficientList[A[x], x] (* Jean-François Alcover, Jan 11 2018 *)
-
g(m):= block([si,v],s:0,v:divisors(m), for si in v do (s:s+r(m/si)/si),s);
r(n):=if n=1 then 1 else sum(Co(n-1,k)/k!,k,1,n-1);
Co(n,k):=if k=1 then g(n) else sum(g(i+1)*Co(n-i-1,k-1),i,0,n-k);
makelist(r(n),n,1,12); /*Vladimir Kruchinin, Jun 15 2012 */
-
{a(n) = local(A = x); if( n<1, 0, for( k=1, n-1, A /= (1 - x^k + x * O(x^n))^polcoeff(A, k)); polcoeff(A, n))}; /* Michael Somos, Dec 16 2002 */
-
{a(n) = local(A, A1, an, i); if( n<1, 0, an = Vec(A = A1 = 1 + O(x^n)); for( m=2, n, i=m\2; an[m] = sum( k=1, i, an[k] * an[m-k]) + polcoeff( if( m%2, A *= (A1 - x^i)^-an[i], A), m-1)); an[n])}; /* Michael Somos, Sep 05 2003 */
-
N=66; A=vector(N+1, j, 1);
for (n=1, N, A[n+1] = 1/n * sum(k=1,n, sumdiv(k,d, d*A[d]) * A[n-k+1] ) );
concat([0], A) \\ Joerg Arndt, Apr 17 2014
-
from functools import lru_cache
from sympy import divisors
@lru_cache(maxsize=None)
def divisor_tuple(n): # cached unordered tuple of divisors
return tuple(divisors(n,generator=True))
@lru_cache(maxsize=None)
def A000081(n): return n if n <= 1 else sum(sum(d*A000081(d) for d in divisor_tuple(k))*A000081(n-k) for k in range(1,n))//(n-1) # Chai Wah Wu, Jan 14 2022
-
@CachedFunction
def a(n):
if n < 2: return n
return add(add(d*a(d) for d in divisors(j))*a(n-j) for j in (1..n-1))/(n-1)
[a(n) for n in range(31)] # Peter Luschny, Jul 18 2014 after Alois P. Heinz
-
[0]+[RootedTrees(n).cardinality() for n in range(1,31)] # Freddy Barrera, Apr 07 2019
Comments