with(combinat); paux := proc(n, k, b) option remember; if n = 0 and k = 0 then return 1 fi; if n <= 0 or k <= 0 or n < k then return 0 fi; if b = 0 then return 0 fi; if n = k then return 1 fi; if k = 1 then if n <= b then return 1 else return 0; fi; fi; paux(n-1, k-1, b) + paux(n-k, k, b-1) end; p := (n, k, a, b) -> paux(n-(a-1)*k, k, b-a+1); F := (k, b) -> 1/k!*add(B[q]*z^q/q!, q=1..b)^k; EX := proc(n, k, a, b, r) option remember; local curF, m, m1, res, Bsubs; if n = 1 then return p(r, k, a, b); fi; Bsubs := [seq(B[q]=add(w^l, l=max(a-q, 0)..b-q), q=1..b)]; res := 0; for m1 from 0 to k-1 do curF := coeftayl(subs(Bsubs, F(k-m1, b)), z=0, n-1); for m from 0 to r do res := res + (n-1)!*p(m, m1, a, b)* coeff(curF, w, r-m); od; od; res; end; T := (n,k,r) -> EX(n, k, 1, n+r-1, r); OEISDATA := (mx, r) -> seq(seq(T(n, k, r), k=1..n+r-1), n=1..mx); FCF := proc(alpha, beta, k, a, b) option remember; local res, q; if alpha < 1 or beta < 0 then return 0 fi; if k = 1 then if max(a-alpha, 0) <= beta and beta <= b-alpha then return 1/alpha! fi; return 0; fi; res := 0; for q to a-1 do res := res + add(`if`(alpha-q>0 and beta-qq>=0, FCF(alpha-q, beta-qq, k-1, a, b), 0), qq=a-q..b-q)/k/q!; od; for q from a to b do res := res + add(`if`(alpha-q>0 and beta-qq>=0, FCF(alpha-q, beta-qq, k-1, a, b), 0), qq=0..b-q)/k/q!; od; res; end; EX2 := (n, k, a, b, r) -> `if`(n=1, p(r,k,a,b), (n-1)!*add(add(p(m, m1, a, b)*FCF(n-1, r-m, k-m1, a, b), m1=0..k-1), m=0..r)); ST2 := (n, k) -> EX2(n, k, 1, n, 1); EX1N := (n, k, r) -> add(add(p(m,m1,1,r)*stirling2(n-1,k-m1) *binomial(k+r-m-m1-1, k-m1-1), m1=0..k-1), m=0..r); T2 := (n,k,r) -> EX2(n, k, 1, n+r-1, r); OEISDATA2 := (mx, r) -> seq(seq(T2(n, k, r), k=1..n+r-1), n=1..mx); VERIF := proc(mx, r) local lA, lB; lA := [OEISDATA(mx, r)]; lB := [OEISDATA2(mx, r)]; {seq(lA[p]-lB[p], p=1..nops(lA))}; end;