%I #6 Aug 13 2012 23:20:08
%S 1,1,4,25,200,1888,20158,237357,3032188,41554144,605964370,9345693140,
%T 151727166822,2583300560490,45984983349166,853637181574329,
%U 16489023127843088,330789284905928356,6880312907650893934,148151276593976715612,3297947033016039111690
%N G.f. satisfies: A(x) = 1 + x*A(x)^3*A(x*A(x)^3).
%F a(n) = coefficient of x^n in (1+x*A(x))^(3*n+1)/(3*n+1) where A(x) = Sum_{n=0} a(n)*x^n.
%F Recurrence:
%F Let A(x)^m = Sum_{n>=0} a(n,m)*x^n with a(0,m)=1, then
%F a(n,m) = Sum_{k=0..n} m*binomial(3*n+m,k)/(3*n+m) * a(n-k,k).
%F G.f. A(x) = F(x,1) where F(x,n) satisfies: F(x,n) = F(x,n-1)*(1 + x*F(x,n)^2*F(x,n+1)) for n>0 with F(x,0)=1.
%e G.f.: A(x) = 1 + x + 4*x^2 + 25*x^3 + 200*x^4 + 1888*x^5 + 20158*x^6 +...
%e G.f. A(x) is the unique solution to variable A in the infinite system of simultaneous equations:
%e A = 1 + x*A^2*B;
%e B = A*(1 + x*B^2*C);
%e C = B*(1 + x*C^2*D);
%e D = C*(1 + x*D^2*E);
%e E = D*(1 + x*E^2*F); ...
%e where B(x) = A(x)*A(x*A(x)^3), C(x) = A(x)*B(x*A(x)^3), D(x) = A(x)*C(x*A(x)^3), ...
%e Expansions of a few of the functions described above begin:
%e B(x) = 1 + 2*x + 12*x^2 + 100*x^3 + 998*x^4 + 11258*x^5 + 139398*x^6 +...
%e C(x) = 1 + 3*x + 24*x^2 + 253*x^3 + 3090*x^4 + 41646*x^5 + 604636*x^6 +...
%e D(x) = 1 + 4*x + 40*x^2 + 512*x^3 + 7452*x^4 + 118016*x^5 + 1990284*x^6 +...
%e E(x) = 1 + 5*x + 60*x^2 + 905*x^3 + 15340*x^4 + 280400*x^5 + 5417498*x^6 +...
%e ALTERNATE GENERATING METHOD.
%e Suppose functions A=A(x), B=B(x), C=C(x), etc., satisfy:
%e A = 1 + x*A^3*B,
%e B = 1 + x*(A*B)^3*C,
%e C = 1 + x*(A*B*C)^3*D,
%e D = 1 + x*(A*B*C*D)^3*E, etc.,
%e then B(x) = A(x*A(x)^3), C(x) = B(x*A(x)^3), D(x) = C(x*A(x)^3), etc.,
%e where A(x) = 1 + x*A(x)^3*A(x*A(x)^3) is the g.f. of this sequence.
%e Expansions of a few of the functions described above begin:
%e B(x) = 1 + x + 7*x^2 + 64*x^3 + 681*x^4 + 8058*x^5 + 103570*x^6 +...
%e C(x) = 1 + x + 10*x^2 + 121*x^3 + 1630*x^4 + 23678*x^5 + 364984*x^6 +...
%e D(x) = 1 + x + 13*x^2 + 196*x^3 + 3209*x^4 + 55660*x^5 + 1010248*x^6 +...
%e E(x) = 1 + x + 16*x^2 + 289*x^3 + 5580*x^4 + 112860*x^5 + 2367358*x^6 +...
%o (PARI) {a(n)=local(A=1+x); for(i=0, n, A=1+x*A^3*subst(A, x, x*A^3+x*O(x^n))); polcoeff(A, n)}
%o (PARI) /* a(n) = [x^n] (1+x*A(x))^(3*n+1)/(3*n+1): */
%o {a(n)=local(A=1+x); for(i=0, n, A=sum(m=0,n,polcoeff((1+x*A+x*O(x^m))^(3*m+1)/(3*m+1),m)*x^m)+x*O(x^n));polcoeff(A,n)}
%o (PARI) {a(n, m=1)=if(n==0, 1, if(m==0, 0^n, sum(k=0, n, m*binomial(3*n+m, k)/(3*n+m)*a(n-k, k))))}
%o for(n=0,31,print1(a(n),", "))
%Y Cf. A088717, A215506, A215507.
%K nonn
%O 0,3
%A _Paul D. Hanna_, Aug 13 2012