Counting integer points inside a sphere of radius R and dimension D

Question

I am trying to write an efficient algorithm that counts the number of points inside a Sphere of Radius R and Dimension D. The sphere is always at the origin. Suppose we have a sphere of dimension 2 (circle) with radius 5.

My strategy is to generate all possible points within the first quadrant, so for the above example we know that (1,2) is in the circle, so must all + / - combinations of that point which is simply dimension squared. So for each point found in a single quadrant of an n-dimensional sphere we add 2 ^ dimension to the total count.

I'm not sure if there is a much more efficient solution to this problem but this is what I have so far in terms of implementation.

int count_lattice_points(const double radius, const int dimension) {
int R = static_cast<int>(radius);

int count  = 0;

std::vector<int> points;
std::vector<int> point;

for(int i = 0; i <= R; i++)
    points.push_back(i);

do {
    for(int i = 0; i < dimension - 1; i++)
        point.push_back(points.at(i));

    if(isPointWithinSphere(point, radius)) count += std::pow(2,dimension);
    point.clear();

}while(std::next_permutation(points.begin(), points.end()));

return count + 3;
}

What can I fix or improve in this situation ?

Answer 1

For 2D case this is Gauss's circle problem. One possible formula:

N(r) = 1 + 4 * r + 4 * Sum[i=1..r]{Floor(Sqrt(r^2-i^2))}

(central point + four quadrants, 4*r for points at the axis, others for in-quadrant region).

Note that there is no known simple closed math expression for 2D case.

In general your idea with quadrants, octants etc is right, but checking all the points is too expensive.

One might find the number of ways to compose all squares from 0 to r^2 from 1..D integer squares (extension of (4) formula).

Note that combinatorics would help to make calculation faster. For example, it is enough to find the number of ways to make X^2 from D natural squares, and multiply by 2^D (different sign combinations); find the number of ways to make X^2 from D-1 natural squares, and multiply by D*2^(D-1) (different sign combinations + D places for zero addend) etc

Example for D=2, R=3

addends: 0,1,4,9
possible sum     compositions    number of variants        
0               0+0             1
1               0+1,1+0         2*2=4
2               1+1             4      
4               0+4,4+0         2*2=4
5               1+4,4+1         2*4=8  
8               4+4             4
9               0+9,9+0         2*2=4
-------------------------------------
                                29

Answer 2

I presented my algorithm for 2D here (with some source code and an ugly but handy illustration): https://stackoverflow.com/a/42373448/5298879

It's around 3.4x faster than MBo's counting points between the origin and the edge of the circle in one of the quarters.

You just imagine an inscribed square and count only one-eighth of what's outside that square inside that circle.

public static int gaussCircleProblem(int radius) {
    int allPoints=0; //holds the sum of points
    double y=0; //will hold the precise y coordinate of a point on the circle edge for a given x coordinate.
    long inscribedSquare=(long) Math.sqrt(radius*radius/2); //the length of the side of an inscribed square in the upper right quarter of the circle
    int x=(int)inscribedSquare; //will hold x coordinate - starts on the edge of the inscribed square
    while(x<=radius){
        allPoints+=(long) y; //returns floor of y, which is initially 0
        x++; //because we need to start behind the inscribed square and move outwards from there
        y=Math.sqrt(radius*radius-x*x); // Pythagorean equation - returns how many points there are vertically between the X axis and the edge of the circle for given x
    }
    allPoints*=8; //because we were counting points in the right half of the upper right corner of that circle, so we had just one-eightth
    allPoints+=(4*inscribedSquare*inscribedSquare); //how many points there are in the inscribed square
    allPoints+=(4*radius+1); //the loop and the inscribed square calculations did not touch the points on the axis and in the center
    return allPoints;
}

Answer 3

An approach similar to that described by MBo, including source code, can be found at https://monsiterdex.wordpress.com/2013/04/05/integer-lattice-in-n-dimensional-sphere-count-of-points-with-integer-coordinates-using-parallel-programming-part-i/.

The approach consists in finding partitions of the radius, and then for each partition in the sphere compute the number of ways it can be represented in the sphere by both permuting coordinates and flipping the signs of nonzero coordinates.