with(numtheory); with(combinat); pet_varinto_cind := proc(poly, ind) local subs1, subs2, polyvars, indvars, v, pot, res, k; res := ind; polyvars := indets(poly); indvars := indets(ind); for v in indvars do pot := op(1, v); subs1 := [seq(polyvars[k]=polyvars[k]^pot, k=1..nops(polyvars))]; subs2 := [v=subs(subs1, poly)]; res := subs(subs2, res); od; res; end; pet_cycleind_cyclic := proc(n) option remember; local d; 1/n*add(phi(d)*a[d]^(n/d), d in divisors(n)); end; pet_neckl_simple := proc(b,w) local n, d; n := b+w; 1/n*add(phi(d)*binomial(n/d, b/d), d in divisors(gcd(b,w))); end; Q := proc(b,w,q) option remember; local res, rep, p, m; if b=0 then if q=0 then return 1 else return 0 fi; fi; if w=0 then if b >= q then return 1 else return 0 fi; fi; rep := add(B^p, p=1..q-1)*add(W^p, p=1..w); for m to floor((b+w)/2) do res := res + pet_varinto_cind(rep, pet_cycleind_cyclic(m)); od; pet_neckl_simple(b, w) - coeff(coeff(expand(res), B, b), W, w); end; P := proc(n,q) option remember; local k; add(Q(k,n-k,q), k=0..n); end; QX := proc(b,w,q) option remember; local n, wpos, idx, maxrun, len, rot, orbit, orbits, src; if b=0 then if q=0 then return 1 else return 0 fi; fi; if w=0 then if b >= q then return 1 else return 0 fi; fi; n := b+w; orbits := table(); for wpos in choose(b+w,w) do maxrun := 0; for idx from 2 to w do len := wpos[idx] - wpos[idx-1] - 1; if len > maxrun then maxrun := len; fi; od; len := wpos[1]+n - wpos[w] - 1; if len > maxrun then maxrun := len; fi; if maxrun >= q then orbit := []; src := table([seq(`B`, idx=1..n)]); for idx to w do src[wpos[idx]] := `W`; od; for rot from 0 to n-1 do orbit := [op(orbit), [seq(src[1+((idx+rot) mod n)], idx = 0..n-1)]]; od; orbits[sort(orbit)[1]] := 1; fi; od; numelems(orbits); end; PX := proc(n,q) option remember; local k; add(QX(k,n-k,q), k=0..n); end; TRIANG := proc(mx, q) local n, b; seq(seq(Q(b, n-b, q), b=0..n), n=1..mx); end;