# Lattice points in a disk

The closed disk of radius ${72}$ has area ${\pi \cdot 72^2\approx 16286}$. But it happens to contain only ${16241}$ points with integer coordinates. Here is a picture of one quarter of this disk.

The radius ${72}$ is somewhat notable is that the discrepancy of ${45}$ between the area and the number of integer points is unusually large. Here is the plot of the absolute value of the difference |area-points| as a function of integer radius ${n}$. The curve in red is ${y = 1.858 r^{0.745}}$, which is an experimentally found upper bound for the discrepancy in this range of ${n}$.

On the scale up to ${n=1000}$, the upper bound is ${4.902 r^{0.548}}$, and the radii bumping against this bound are ${449}$ and ${893}$. The exponent ${0.548}$ begins to resemble the conjectural ${0.5+\epsilon}$ in the Gauss circle problem.

Finally, over the interval ${1000\le n\le 3000}$ the upper bound comes out as ${6.607n^{0.517}}$. The exponent ${0.517}$ looks good.

This little numerical experiment in Matlab involved using the convex hull function convhull on log-log scale. The function identifies the vertices of the convex hull, which is a polygon. I pick the side of the polygon lying over the midpoint of the range; this yields a linear upper bound for the set of points. On the normal scale, this line becomes a power function. Matlab code is given below; it’s divided into three logical steps.

##### Find the difference between area and the number of lattice points
a = 1000;
b = 3000;
R = a:b;
E = [];
for n = a:b
[X,Y] = meshgrid(1:n, 1:n);
pts = 4*n + 1 + 4*nnz(X.^2+Y.^2<=n^2);
E = [E, max(1,abs(pts - pi*n^2))];
end
##### Pick a suitable side of log-log convex hull
ix = convhull(log(R), log(E));
k = numel(ix);
while (R(ix(k))<(a+b)/2)
k = k-1;
end
##### Plot the result and output the parameters of the upper bound
R1 = R(ix(k)); E1 = E(ix(k));
R2 = R(ix(k+1)); E2 = E(ix(k+1));
b = log(E1/E2)/log(R1/R2);
a = E1/R1^b;
plot(R, E, '.');
hold on
plot(R, a*R.^b , 'r-');
axis tight
hold off
fprintf('a = %.3f, b = %.3f\n', a, b);