with(combinat);

CAT := n -> 1/(n+1)*binomial(2*n,n);

CAT_PARS :=
proc(n)
option remember;
local recurse, res;

    res := [];

    recurse :=
    proc(sofar, oc, cc)
        if oc = 0 and cc = 0 then
            res := [op(res), sofar];
            return;
        fi;

        if oc > 0 then
            recurse([op(sofar), `(`], oc-1, cc+1);
        fi;
        if cc > 0 then
            recurse([op(sofar), `)`], oc, cc-1);
        fi;
    end;

    recurse([], n, 0);
    convert(res, `set`);
end;

BALCOUNT :=
proc(parlist)
local posA, posB, intervA, intervB, len, subseq, count,
    res, n;

    subseq := [];

    n := nops(parlist);

    for posA to n do
        for posB from posA+1 to n do
            intervA := parlist[posA..posB];
            len := posB - posA + 1;

            if type(len, even) and
            intervA in CAT_PARS(len/2) then
                subseq := [op(subseq), [posA, posB]];
            fi;
        od;
    od;

    res := [];

    subseq :=
    sort(subseq,
         (a, b) -> a[2]-a[1] < b[2]-b[1]);
    count := nops(subseq);

    for posA to count do
        intervA := subseq[posA];
        for posB from posA+1 to count do
            intervB := subseq[posB];

            if intervB[1] <= intervA[1] and
            intervB[2] >= intervA[2] then
                break;
            fi;
        od;

        if posB = count + 1 then
            res := [op(res), intervA];
        fi;
    od;

    add((iv[2]-iv[1]+1)/2, iv in res);
end;

GFP :=
proc(n)
option remember;
local gf, idx, d, parlist;

    gf := 0;

    for idx from 2^n to 2*2^n-1 do
        d := convert(idx, base, 2);

        parlist :=
        subs([0=`(`, 1=`)`], d[1..n]);

        gf := gf + u^BALCOUNT(parlist);
    od;

    gf;
end;

GFA :=
proc(n)
local C;

    C := (1-sqrt(1-4*u*z^2))/(2*u*z^2);

    coeftayl(C/(1-z*C)^2, z=0, n);
end;

GFB :=
n -> add(u^k*add(binomial(2*k,q)*binomial(n-2*k+1, k-q)
         *(n+2*q+1-4*k), q=0..k), k=0..floor(n/2));

GFC :=
n -> add(u^k*(n+1-2*k)^2/(n+1)*binomial(n+1,k),
         k=0..floor(n/2));