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;