Recall the cycle index of the cyclic group
$$Z(C_n) = \frac{1}{n} \sum_{d|n} \varphi(d) a_d^{n/d}.$$
Supposing that we have $n = b+w$ beads with $b$ black and $w$ white
and there is a run of black beads of length at least $q$ (maximum run
length at least $q$) we get from the Polya Enumeration Theorem the
closed form
$$Q_{b,w} = [B^b W^w] Z(C_n; B+W)
- [B^b W^w] \sum_{m=1}^{\lfloor n/2 \rfloor}
Z\left(C_m; \sum_{p=1}^{q-1} B^p \sum_{p=1}^w W^p\right).$$
Here we have the boundary case that $Q_{0,w}$ is one if $q=0$ and zero
otherwise and that $Q_{b,0}$ is one if $b\ge q$ and zero otherwise.
Note that by doing the substitution we find
$$[B^b W^w] Z(C_n; B+W)
= \frac{1}{n} \sum_{d|(b,w)}
\varphi(d) {n/d\choose b/d}.$$
When we ask about binary necklaces on $n$ beads with maximum run length
for black beads being at least $q$ we find the quantity
$$P_n = \sum_{k=0}^n Q_{k, n-k}.$$
These formulas have an implicit extra parameter $q.$ We present an
implementation in Maple that also includes an enumeration routine that
was used to verify these statistics. We can compute some example
sequences. Here is the count for necklaces with the same number of
black and white beads with the maximum run length at least zero (admit
all, this is $Q_{b,b,0}$):
$$1, 2, 4, 10, 26, 80, 246, 810, 2704, \ldots$$
which points us to OEIS A003239 where these
data are confirmed. As an additional sanity check we can compute the
count of all binary necklaces (admit all, which is $P_{n,0}$):
$$1, 2, 3, 4, 6, 8, 14, 20, 36, 60, 108, 188, 352, 632,
\\ 1182, 2192, 4116, 7712, \ldots$$
which points to OEIS A000031 for
confirmation. We obtain one minus this sequence when we compute
$P_{n,1}$ which is correct as well (omit only the white monochrome
necklace). The sequence with maximum run length at least three is
given by
$$0, 0, 1, 2, 3, 5, 9, 17, 31, 60, 113, 220, 419, 813,
\\ 1565, 3033, 5855, \ldots$$
Maximum run length at least four is
$$0, 0, 0, 1, 2, 3, 5, 9, 17, 33, 63, 124, 241, 475,
\\ 930, 1831, 3593, \ldots$$
To conclude we present the sequence where we have four more black
beads than white ones and the maximum run length for black is at least
four (this is $Q_{w+4,w,4}$):
$$1, 3, 10, 34, 116, 411, 1464, 5292, 19246,
\\ 70533, 259766, 961423, \ldots$$
We have posted these values here to encourage the reader to present
their own implementation and verify them, or perhaps simplify the
formula from the introduction. The Maple code was as follows (PET
first, followed by ENUM):
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;