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.

Unusually few lattice points in the disk
Unusually few lattice points in the 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}.

Radius from 1 to 100
Radius up to 100

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.

Radius up to 1000
Radius up to 1000

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.

Radius from 1000 to 3000
Radius from 1000 to 3000

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))];
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;
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);

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s