Only the start of an answer, but too long for a comment: $H_k$ is proportional to the $z^{4k-2}$ coefficient of the $\wp$ function associated to an elliptic curve of $j$-invariant $1728$. The Weierstrass differential equation for $\wp$ gives various ways to compute these coefficients, and thus the $H_k$, without recourse to high-precision floating-point arithmetic. For example, writing $\wp$ as the inverse function of an elliptic integral yields the following **gp** code that recovers all the tabulated numbers; increase $N$ to get more terms: N = 20; x = 1 / serreverse(intformal(1/sqrt(1-r^4/16+O(r^(4*N+1)))))^2; H = vector(N, n, (4*n)! * polcoeff(x,4*n-2) / (4*n-1)) This connection might yield at least some of the results that David observed experimentally. For a start, most of the divisibility by high powers of small primes $p \equiv 3 \bmod 4$ is explained by the factor $k!/(k-1)$ together with observation that in the power-series expansion of $\int (1-(r^4/16))^{-1/2} dr = \sum_m a_m r^m$ the valuation of $a_m$ is not as negative as the usual $-v_p(m)$ — this must be a manifestation of the supersingularity of the curve at such primes; and the primes in the denominator of $H_k$ should be amenable to a similar analysis. **EDIT** Here's an alternative algorithm for producing the $H_k$, via an explicit recursion starting from the differential equation ${\wp'}^2 = \wp^3 - \frac14\wp$. It's convenient to differentiate again and cancel the common factor of $\wp'$, obtaining $\wp'' = 6 \wp^2 - \frac18$. Then if we let the $z^{4k-2}$ coefficient of $\wp$ be $h_k = (k-1) H_k / k!$, and equate coefficients in the formula for $\wp''$, we find $$ h_k = \frac6{(k+1)(k-6)} \sum_{0<j<k} h_j h_{k-j} $$ for each $k>4$. In **gp**: N = 20; v = vector(N); v[1] = 1/80 \\ = polcoeff(x,2) for(n=2, N, v[n] = 6*sum(m=1,n-1,v[m]*v[n-m]) / ((4*n+1)*(4*n-6))) for(n=1, N, v[n] *= (4*n)!/(4*n-1)) v (Yes, this agrees with the previous calculation through $H_{80}$.) This should make it easier to account for the denominators of $H_k$.