from mpmath import mp
mp.dps = 80
prob = lambda k: mp.log(1 + 1/(k*(k + 2))) / mp.log(2)
mp.exp(mp.nsum(lambda n: (mp.zeta(2*n) - 1)/n
* mp.fsum((-1)**(j+1)/mp.mpf(j)
for j in range(1, 2*int(n))),
[1, mp.inf]) / mp.log(2)) # K_0
mp.nsum(lambda k: prob(k)/k, [1, mp.inf]) ** -1 # K_{-1}Each row was computed from (3) for $p=0$ and from (4) for $p\ne0$, twice: once at $150$ working decimal digits with $680$ grouped zeta terms, and once at $210$ working decimal digits with $920$ grouped zeta terms.
The two values agreed in their first $105$ significant digits, and what is stored is the union of the two decimal values. Agreement between two truncations is evidence and not a proof, which is what the rigour says. The rows $K_0,K_{-1},\ldots,K_{-10}$ were compared with the decimal expansions in the OEIS entries named in their comments and with the table on MathWorld's Khinchin harmonic mean page [6]. The two half-integer rows have no OEIS entry to compare with and were checked against direct partial sums of (2), with the positive tail bounded by $\sum_{k>N}k^{p-2}/\log 2$; at $N=200000$ the resulting intervals contain the stored values. The monotonicity in (8) was checked on the stored values.