This is a front-end for the Online Encyclopedia of Integer Sequences, made by Christian Perfect. The idea is to provide OEIS entries in non-ancient HTML, and then to think about how they're presented visually. The source code is on GitHub.
G.f. = x + x^2 + 2*x^3 + 5*x^4 + 12*x^5 + 33*x^6 + 90*x^7 + 261*x^8 + ... a(4)=5 with the following series-reduced planted trees: (oooo), (oo(oo)), (o(ooo)), (o(o(oo))), ((oo)(oo)). - _Michael Somos_, Jul 25 2003
Method 1: a := [1,1]; for n from 3 to 30 do L := series( mul( (1-x^k)^(-a[k]),k=1..n-1)/(1-x^n)^b, x,n+1); t1 := coeff(L,x,n); R := series( 1+2*add(a[k]*x^k,k=1..n-1)+2*b*x^n, x, n+1); t2 := coeff(R,x,n); t3 := solve(t1-t2,b); a := [op(a),t3]; od: A000669 := n-> a[n]; Method 2, more efficient: with(numtheory): M := 1001; a := array(0..M); p := array(0..M); a[1] := 1; a[2] := 1; a[3] := 2; p[1] := 1; p[2] := 3; p[3] := 7; Method 2, cont.: for m from 4 to M do t1 := divisors(m); t3 := 0; for d in t1 minus {m} do t3 := t3+d*a[d]; od: t4 := p[m-1]+2*add(p[k]*a[m-k],k=1..m-2)+t3; a[m] := t4/m; p[m] := t3+t4; od: # A000669 := n-> a[n]; A058757 := n->p[n]; # Method 3: b:= proc(n, i) option remember; `if`(n=0, 1, `if`(i<1, 0, add(binomial(a(i)+j-1, j)* b(n-i*j, i-1), j=0..n/i))) end: a:= n-> `if`(n<2, n, b(n, n-1)): seq(a(n), n=1..40); # Alois P. Heinz, Jan 28 2016
b[n_, i_] := b[n, i] = If[n==0, 1, If[i<1, 0, Sum[Binomial[a[i]+j-1, j]* b[n-i*j, i-1], {j, 0, n/i}]]];
a[n_] := If[n<2, n, b[n, n-1]];
Array[a, 40] (* Jean-François Alcover, Jan 08 2021, after Alois P. Heinz *)
{a(n) = my(A, X); if( n<2, n>0, X = x + x * O(x^n); A = 1 / (1 - X); for(k=2, n, A /= (1 - X^k)^polcoeff(A, k)); polcoeff(A, n)/2)}; /* Michael Somos, Jul 25 2003 */
from collections import Counter def A000669_list(n): list = [1] + [0] * (n - 1) for i in range(1, n): for p in Partitions(i + 1, min_length=2): m = Counter(p) list[i] += prod(binomial(list[s - 1] + m[s] - 1, m[s]) for s in m) return list print(A000669_list(20)) # M. Eren Kesim, Jun 21 2021
E.g.f.: A(x) = x + x^2/2! + 4*x^3/3! + 26*x^4/4! + 236*x^5/5! + 2752*x^6/6! + ...
where exp(A(x)) = 1 - x + 2*A(x), and thus
Series_Reversion(A(x)) = x - x^2/2! - x^3/3! - x^4/4! - x^5/5! - x^6/6! + ...
O.g.f.: G(x) = x + x^2 + 4*x^3 + 26*x^4 + 236*x^5 + 2752*x^6 + 39208*x^7 + ...
where
G(x) = x/2 + x/(2*(2-x)) + x/(2*(2-x)*(2-2*x)) + x/(2*(2-x)*(2-2*x)*(2-3*x)) + x/(2*(2-x)*(2-2*x)*(2-3*x)*(2-4*x)) + x/(2*(2-x)*(2-2*x)*(2-3*x)*(2-4*x)*(2-5*x)) + ...
From _Gus Wiseman_, Dec 28 2019: (Start)
A rooted tree is series-reduced if it has no unary branchings, so every non-leaf node covers at least two other nodes. The a(4) = 26 series-reduced rooted trees with 4 labeled leaves are the following. Each bracket (...) corresponds to a non-leaf node.
(1234) ((12)34) ((123)4)
(1(23)4) (1(234))
(12(34)) ((124)3)
(1(24)3) ((134)2)
((13)24) (((12)3)4)
((14)23) ((1(23))4)
((12)(34))
(1((23)4))
(1(2(34)))
(((12)4)3)
((1(24))3)
(1((24)3))
(((13)2)4)
((13)(24))
(((13)4)2)
((1(34))2)
(((14)2)3)
((14)(23))
(((14)3)2)
(End)
M:=499; a:=array(0..500); a[0]:=0; a[1]:=1; a[2]:=1; for n from 0 to 2 do lprint(n,a[n]); od: for n from 2 to M do a[n+1]:=(n+2)*a[n]+2*add(binomial(n,k)*a[k]*a[n-k+1],k=2..n-1); lprint(n+1,a[n+1]); od: Order := 50; t1 := solve(series((exp(A)-2*A-1),A)=-x,A); A000311 := n-> n!*coeff(t1,x,n); # second Maple program: b:= proc(n, i) option remember; `if`(n=0, 1, `if`(i<1, 0, add(combinat[multinomial](n, n-i*j, i$j)/j!* a(i)^j*b(n-i*j, i-1), j=0..n/i))) end: a:= n-> `if`(n<2, n, b(n, n-1)): seq(a(n), n=0..40); # Alois P. Heinz, Jan 28 2016 # faster program: b:= proc(n, i) option remember; `if`(i=0 and n=0, 1, `if`(i<=0 or i>n, 0, i*b(n-1, i) + (n+i-1)*b(n-1, i-1))) end: a:= n -> `if`(n<2, n, add(b(n-1, i), i=0..n-1)): seq(a(n), n=0..40); # Peter Luschny, Feb 15 2021
nn = 19; CoefficientList[ InverseSeries[ Series[1+2a-E^a, {a, 0, nn}], x], x]*Range[0, nn]! (* Jean-François Alcover, Jul 21 2011 *)
a[ n_] := If[ n < 1, 0, n! SeriesCoefficient[ InverseSeries[ Series[ 1 + 2 x - Exp[x], {x, 0, n}]], n]]; (* Michael Somos, Jun 04 2012 *)
a[n_] := (If[n < 2,n,(column = ConstantArray[0, n - 1]; column[[1]] = 1; For[j = 3, j <= n, j++, column = column * Flatten[{Range[j - 2], ConstantArray[0, (n - j) + 1]}] + Drop[Prepend[column, 0], -1] * Flatten[{Range[j - 1, 2*j - 3], ConstantArray[0, n - j]}];]; Sum[column[[i]], {i, n - 1}] )]); Table[a[n], {n, 0, 20}] (* Peter Regner, Oct 05 2012, after a formula by Felsenstein (1978) *)
multinomial[n_, k_List] := n!/Times @@ (k!); b[n_, i_] := b[n, i] = If[n == 0, 1, If[i<1, 0, Sum[multinomial[n, Join[{n-i*j}, Array[i&,j]]]/j!*a[i]^j *b[n-i*j, i-1], {j, 0, n/i}]]]; a[n_] := If[n<2, n, b[n, n-1]]; Table[ a[n], {n, 0, 40}] (* Jean-François Alcover, Feb 07 2016, after Alois P. Heinz *)
sps[{}]:={{}};sps[set:{i_,_}]:=Join@@Function[s,Prepend[#,s]&/@sps[Complement[set,s]]]/@Cases[Subsets[set],{i,_}];
mtot[m_]:=Prepend[Join@@Table[Tuples[mtot/@p],{p,Select[sps[m],1Gus Wiseman, Dec 28 2019 *)
(* Lengthy but easy to follow *)
lead[, n /; n < 3] := 0
lead[h_, n_] := Module[{p, i},
p = Position[h, {_}];
Sum[MapAt[{#, n} &, h, p[[i]]], {i, Length[p]}]
]
follow[h_, n_] := Module[{r, i},
r = Replace[Position[h, {_}], {a__} -> {a, -1}, 1];
Sum[Insert[h, n, r[[i]]], {i, Length[r]}]
]
marry[, n /; n < 3] := 0
marry[h_, n_] := Module[{p, i},
p = Position[h, _Integer];
Sum[MapAt[{#, n} &, h, p[[i]]], {i, Length[p]}]
]
extend[a_ + b_, n_] := extend[a, n] + extend[b, n]
extend[a_, n_] := lead[a, n] + follow[a, n] + marry[a, n]
hierarchies[1] := hierarchies[1] = extend[hier[{}], 1]
hierarchies[n_] := hierarchies[n] = extend[hierarchies[n - 1], n] (* Daniel Geisler, Aug 22 2022 *)
a(n):=if n=1 then 1 else sum((n+k-1)!*sum(1/(k-j)!*sum((2^i*(-1)^(i)*stirling2(n+j-i-1,j-i))/((n+j-i-1)!*i!),i,0,j),j,1,k),k,1,n-1); /* Vladimir Kruchinin, Jan 28 2012 */
{a(n) = local(A); if( n<0, 0, for( i=1, n, A = Pol(exp(A + x * O(x^i)) - A + x - 1)); n! * polcoeff(A, n))}; /* Michael Somos, Jan 15 2004 */
{a(n) = my(A); if( n<0, 0, A = O(x); for( i=1, n, A = intformal( 1 / (1 + x - 2*A))); n! * polcoeff(A, n))}; /* Michael Somos, Oct 25 2014 */
{a(n) = n! * polcoeff(serreverse(1+2*x - exp(x +x^2*O(x^n))), n)}
for(n=0, 30, print1(a(n), ", ")) \\ Paul D. Hanna, Oct 27 2014
\p100 \\ set precision
{A=Vec(sum(n=0, 600, 1.*x/prod(k=0, n, 2 - k*x + O(x^31))))}
for(n=0, 25, print1(if(n<1,0,round(A[n])),", ")) \\ Paul D. Hanna, Oct 27 2014
from functools import lru_cache from math import comb @lru_cache(maxsize=None) def A000311(n): return n if n <= 1 else -(n-1)*A000311(n-1)+comb(n,m:=n+1>>1)*(0 if n&1 else A000311(m)**2) + (sum(comb(n,i)*A000311(i)*A000311(n-i) for i in range(1,m))<<1) # Chai Wah Wu, Nov 10 2022
From _Andrew Howroyd_, Nov 26 2020: (Start) In the following examples of series-parallel networks, elements in series are juxtaposed and elements in parallel are separated by '|'. The unit element is denoted by 'o'. a(1) = 1: (o). a(2) = 2: (oo), (o|o). a(3) = 5: (ooo), (o(o|o)), ((o|o)o), (o|o|o), (o|oo). a(4) = 15: (oooo), (oo(o|o)), (o(o|o)o), ((o|o)oo), ((o|o)(o|o)), (o(o|oo)), (o(o|o|o)), ((o|oo)o), ((o|o|o)o), (o|o|o|o), (o|o|oo), (oo|oo), (o|ooo), (o|o(o|o)), (o|(o|o)o). (End)
terms = 25; A[] = 1; Do[A[x] = Exp[Sum[(1/k)*(A[x^k] + 1/A[x^k] - 2 + x^k), {k, 1, terms + 1}]] + O[x]^(terms + 1) // Normal, terms + 1]; CoefficientList[A[x], x] // Rest (* Jean-François Alcover, Jun 29 2011, updated Jan 12 2018 *)
EulerT(v)={Vec(exp(x*Ser(dirmul(v,vector(#v,n,1/n))))-1, -#v)}
seq(n)={my(p=x+O(x^2)); for(n=2, n, p=x*Ser(EulerT(Vec(p^2/(1+p)+x, 1-n)))); Vec(p)} \\ Andrew Howroyd, Nov 27 2020
a(2) = 2 since given two 1-ohm resistors, a series circuit yields 2 ohms, while a parallel circuit yields 1/2 ohms.
r:= proc(n) option remember; `if`(n=1, {1}, {seq(seq(seq(
[f+g, 1/(1/f+1/g)][], g in r(n-i)), f in r(i)), i=1..n/2)})
end:
a:= n-> nops(r(n)):
seq(a(n), n=1..15); # Alois P. Heinz, Apr 02 2015
r[n_] := r[n] = If[n == 1, {1}, Union @ Flatten @ {Table[ Table[ Table[ {f+g, 1/(1/f+1/g)}, {g, r[n-i]}], {f, r[i]}], {i, 1, n/2}]}]; a[n_] := Length[r[n]]; Table[a[n], {n, 1, 15}] (* Jean-François Alcover, May 28 2015, after Alois P. Heinz *)
\\ not efficient; just to show the method
N=10;
L=vector(N); L[1]=[1];
{ for (n=2, N,
my( T = Set( [] ) );
for (k=1, n\2,
for (j=1, #L[k],
my( r1 = L[k][j] );
for (i=1, #L[n-k],
my( r2 = L[n-k][i] );
T = setunion(T, Set([r1+r2, r1*r2/(r1+r2) ]) );
);
);
);
T = vecsort(Vec(T), , 8);
L[n] = T;
); }
for(n=1, N, print1(#L[n], ", ") );
\\ Joerg Arndt, Mar 07 2015
D^3(1) = (12*x^2+56*x+52)/(x-1)^6. Evaluated at x = 0 this gives a(4) = 52. a(3) = 8: The 8 possible increasing plane trees on 3 vertices with vertices of outdegree k >= 1 coming in 2 colors, B or W, are ....................................................... .1B..1B..1W..1W.....1B.......1W........1B........1W.... .|...|...|...|...../.\....../..\....../..\....../..\... .2B..2W..2B..2W...2...3....2....3....3....2....3....2.. .|...|...|...|......................................... .3...3...3...3......................................... G.f. = x + 2*x^2 + 8*x^3 + 52*x^4 + 472*x^5 + 5504*x^6 + 78416*x^7 + ...
read transforms; t1 := 2*ln(1+x)-x; t2 := series(t1,x,10); t3 := seriestoseries(t2,'revogf'); t4 := SERIESTOLISTMULT(%); # N denotes all series-parallel networks, S = series networks, P = parallel networks; spec := [ N, N=Union(Z,S,P),S=Set(Union(Z,P),card>=2), P=Set(Union(Z,S), card>=2)}, labeled ]: A006351 := n->combstruct[count](spec,size=n); A006351 := n -> add(combinat[eulerian2](n-1,k)*2^(n-k-1),k=0..n-1): seq(A006351(n), n=1..18); # Peter Luschny, Nov 16 2012
max = 18; f[x_] := 2*Log[1+x]-x; Rest[ CoefficientList[ InverseSeries[ Series[ f[x], {x, 0, max}], x], x]]*Range[max]! (* Jean-François Alcover, Nov 25 2011 *)
a(n):=if n=1 then 1 else ((n-1)!*sum(binomial(n+k-1,n-1)* sum((-1)^(j)*binomial(k,j)*sum((binomial(j,l)*(j-l)!*2^(j-l)*(-1)^l* stirling1(n-l+j-1,j-l))/(n-l+j-1)!,l,0,j),j,1,k),k,1,n-1)); /* Vladimir Kruchinin, Jan 24 2012 */
x='x+O('x^66); Vec(serlaplace(serreverse( 2*log(1+x) - 1*x ))) \\ Joerg Arndt, May 01 2013
# uses[eulerian2 from A201637] def A006351(n): return add(A201637(n-1, k)*2^(n-k-1) for k in (0..n-1)) [A006351(n) for n in (1..18)] # Peter Luschny, Nov 16 2012
In the following examples of series-parallel networks, elements in series are juxtaposed and elements in parallel are separated by '|'. The unit element is denoted by 'o'. a(1) = 1: (o). a(2) = 2: (oo), (o|o). a(3) = 4: (ooo), (o(o|o)), (o|o|o), (o|oo). a(4) = 11: (oooo), (oo(o|o)), (o(o|o)o), ((o|o)(o|o)), (o(o|oo)), (o(o|o|o)), (o|o|o|o), (o|o|oo), (oo|oo), (o|ooo), (o|o(o|o)).
\\ here B(n) gives A003430 as a power series. EulerT(v)={Vec(exp(x*Ser(dirmul(v,vector(#v,n,1/n))))-1, -#v)} B(n)={my(p=x+O(x^2)); for(n=2, n, p=x*Ser(EulerT(Vec(p^2/(1+p)+x)))); p} seq(n)={my(q=subst(B((n+1)\2), x, x^2), s=x^2+q^2/(1+q), p=x+O(x^2), t=p); for(n=1, n\2, t=x + q*(1 + p); p=x + x*Ser(EulerT(Vec(t+(s-subst(t,x,x^2))/2))) - t); Vec(p+t-x+B(n))/2}
A(x) = x + 1/2*x^3 + 7/24*x^4 + 3/4*x^5 + 169/180*x^6 + ... For n=4 there are two unlabeled networks: ..o.....o--o ./.\.../....\ o...o o------o .\./ ..o which can be labeled in 3 (resp. 4) ways, for a total of 7.
Q := x; for d from 1 to 30 do Q := Q+c*x^(d+1)/(d+1)!; t1 := coeff(series(2*Q - (exp(Q)-1+log(1+x)), x, d+2), x, d+1); t2 := solve(t1,c); Q := subs(c=t2,Q); Q := series(Q,x,d+2); od: A058379 := n->coeff(Q,x,n)*n!; # method 1 Order := 50; t1 := solve(series((exp(A)-2*A-1),A)=-log(1+x),A); A058379 := n-> n!*coeff(t1,x,n); # method 2
CoefficientList[InverseSeries[Series[-1+E^(1+2*a-E^a), {a, 0, 20}], x], x]*Range[0, 20]! (* Jean-François Alcover, Jul 21 2011 *)
CoefficientList[Series[(-1 + Log[1+x] - 2*ProductLog[-Sqrt[1+x]/(2*Sqrt[E])])/2,{x,0,15}],x] * Range[0,15]! (* Vaclav Kotesovec, Jan 08 2014 *)
a(n):=sum((sum((m+k-1)!*sum(sum(((-1)^i*2^i*stirling2(m+j-i-1, j-i))/(i!*(m+j-i-1)!),i,0,j)/(k-j)!,j,0,k),k,0,m-1)) *stirling1(n,m), m,1,n); /* Vladimir Kruchinin, Feb 17 2012 */
Vec(serlaplace(-1/2 + log(1+x)/2 - lambertw(-exp(-.5)*sqrt(1+x)/2))) \/ 1 \\ Charles R Greathouse IV, Jun 16 2021
T(10,8) = 3 because the partitions of 10 with 8 parts are 31111111 and 22111111. The partition 31111111 corresponds to 2 graphs and the partition 22111111 corresponds to only one. T(n,m) = 1, if and only if m>=n-1. Because A000669(1)=A000669(2)=1, the partitions of n with all parts <=2 correspond to summands = 1. If there is only a summand (or partition), the total is equal to 1. It is clear that for m>=n-1 there is only one partition of n with exactly m parts. Triangle begins: 1, 1, 1, 2, 1, 1, 5, 3, 1, 1, 12, 7, 3, 1, 1, 33, 20, 8, 3, 1, 1, 90, 55, 22, 8, 3, 1, 1,
with(gfun): f := series((ln(1+x)-x^2/(1+x)), x, 21): egf := seriestoseries(f, 'revogf'): seriestolist(egf, 'Laplace');
lim = 19; Join[{1}, Drop[ CoefficientList[ InverseSeries[ Series[x + 2*(1 - Cosh[x]) , {x, 0, lim}], y] + InverseSeries[ Series[-Log[1 - x] - x^2/(1 - x),{x, 0, lim}], y], y], 2]*Range[2, lim]!] (* Jean-François Alcover, Sep 21 2011, after g.f. *)
m = 17; Rest[CoefficientList[InverseSeries[Series[Log[1+x]-x^2/(1+x), {x, 0, m}], x], x]*Table[k!,{k, 0, m}]](* Jean-François Alcover, Apr 18 2011 *)
h(n,k):=if n=k then 0 else (-1)^(n-k)*binomial(n-k-1,k-1); a(n):=if n=1 then 1 else -sum((k!/n!*stirling1(n,k)+sum(binomial(k,j)*sum((j)!/(i)!*stirling1(i,j)*h(n-i,k-j),i,j,n-k+j),j,1,k-1)+h(n,k))*a(k),k,1,n-1); /* Vladimir Kruchinin, Sep 08 2010 */
x='x+O('x^55);
s=-log(1-x)-x^2/(1-x);
A048174=Vec(serlaplace(serreverse(s)));
t=x+2*(1-cosh(x));
A058349=Vec(serlaplace(serreverse(t)));
A048172=A048174+A058349; A048172[1]-=1;
A048172 /* Joerg Arndt, Feb 04 2011 */
Comments