The closed disk of radius has area . But it happens to contain only points with integer coordinates. Here is a picture of one quarter of this disk.

The radius is somewhat notable is that the discrepancy of 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 . The curve in red is , which is an experimentally found upper bound for the discrepancy in this range of .

On the scale up to , the upper bound is , and the radii bumping against this bound are and . The exponent begins to resemble the conjectural in the Gauss circle problem.

Finally, over the interval the upper bound comes out as . The exponent 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);
```