Best way to find all points of lattice in sphere

Question

Given a bunch of arbitrary vectors (stored in a matrix A) and a radius r, I'd like to find all integer-valued linear combinations of those vectors which land inside a sphere of radius r. The necessary coordinates I would then store in a Matrix V. So, for instance, if the linear combination

K=[0; 1; 0]

lands inside my sphere, i.e. something like

if norm(A*K) <= r then
   V(:,1)=K
end

etc.

The vectors in A are sure to be the simplest possible basis for the given lattice and the largest vector will have length 1. Not sure if that restricts the vectors in any useful way but I suspect it might. - They won't have as similar directions as a less ideal basis would have.

I tried a few approaches already but none of them seem particularly satisfying. I can't seem to find a nice pattern to traverse the lattice.

My current approach involves starting in the middle (i.e. with the linear combination of all 0s) and go through the necessary coordinates one by one. It involves storing a bunch of extra vectors to keep track of, so I can go through all the octants (in the 3D case) of the coordinates and find them one by one. This implementation seems awfully complex and not very flexible (in particular it doesn't seem to be easily generalizable to arbitrary numbers of dimension - although that isn't strictly necessary for the current purpose, it'd be a nice-to-have)

Is there a nice* way to find all the required points?

(*Ideally both efficient and elegant**. If REALLY necessary, it wouldn't matter THAT much to have a few extra points outside the sphere but preferably not that many more. I definitely do need all the vectors inside the sphere. - if it makes a large difference, I'm most interested in the 3D case.

**I'm pretty sure my current implementation is neither.)

---

Similar questions I found:

Find all points in sphere of radius r around arbitrary coordinate - this is actually a much more general case than what I'm looking for. I am only dealing with periodic lattices and my sphere is always centered at 0, coinciding with one point on the lattice. But I don't have a list of points but rather a matrix of vectors with which I can generate all the points.

How to efficiently enumerate all points of sphere in n-dimensional grid - the case for a completely regular hypercubic lattice and the Manhattan-distance. I'm looking for completely arbitary lattices and euclidean distance (or, for efficiency purposes, obviously the square of that).

Answer 1

Offhand, without proving any assertions, I think that 1) if the set of vectors is not of maximal rank then the number of solutions is infinite; 2) if the set is of maximal rank, then the image of the linear transformation generated by the vectors is a subspace (e.g., plane) of the target space, which intersects the sphere in a lower-dimensional sphere; 3) it follows that you can reduce the problem to a 1-1 linear transformation (kxk matrix on a k-dimensional space); 4) since the matrix is invertible, you can "pull back" the sphere to an ellipsoid in the space containing the lattice points, and as a bonus you get a nice geometric description of the ellipsoid (principal axis theorem); 5) your problem now becomes exactly one of determining the lattice points inside the ellipsoid.

The latter problem is related to an old problem (counting the lattice points inside an ellipse) which was considered by Gauss, who derived a good approximation. Determining the lattice points inside an ellipse(oid) is probably not such a tidy problem, but it probably can be reduced one dimension at a time (the cross-section of an ellipsoid and a plane is another ellipsoid).

Answer 2

I found a method that makes me a lot happier for now. There may still be possible improvements, so if you have a better method, or find an error in this code, definitely please share. Though here is what I have for now: (all written in SciLab)

---

Step 1: Figure out the maximal ranges as defined by a bounding n-parallelotope aligned with the axes of the lattice vectors. Thanks for ElKamina's vague suggestion as well as this reply to another of my questions over on math.se by chappers: https://math.stackexchange.com/a/1230160/49989

    function I=findMaxComponents(A,r) //given a matrix A of lattice basis vectors
                                      //and a sphere radius r,
                                      //find the corners of the bounding parallelotope
                                      //built from the lattice, and store it in I.
        [dims,vecs]=size(A); //figure out how many vectors there are in A (and, unnecessarily, how long they are)
        U=eye(vecs,vecs); //builds matching unit matrix
        iATA=pinv(A'*A); //finds the (pseudo-)inverse of A^T A
        iAT=pinv(A'); //finds the (pseudo-)inverse of A^T
        I=[]; //initializes I as an empty vector
        for i=1:vecs                           //for each lattice vector,
            t=r*(iATA*U(:,i))/norm(iAT*U(:,i)) //find the maximum component such that
                                               //it fits in the bounding n-parallelotope
                                               //of a (n-1)-sphere of radius r
            I=[I,t(i)];                        //and append it to I
        end
        I=[-I;I]; //also append the minima (by symmetry, the negative maxima)
    endfunction

In my question I only asked for a general basis, i.e, for n dimensions, a set of n arbitrary but linearly independent vectors. The above code, by virtue of using the pseudo-inverse, works for matrices of arbitrary shapes and, similarly, Scilab's "A'" returns the conjugate transpose rather than just the transpose of A so it equally should work for complex matrices.

In the last step I put the corresponding minimal components.

For one such A as an example, this gives me the following in Scilab's console:

     A  =
     
        0.9701425  - 0.2425356    0.
        0.2425356    0.4850713    0.7276069
        0.2425356    0.7276069  - 0.2425356
    
    r=3;

    I=findMaxComponents(A,r)
    
     I  =
     
      - 2.9494438  - 3.4186986  - 4.0826424
        2.9494438    3.4186986    4.0826424
    
     I=int(I)
     
     I  =
     
      - 2.  - 3.  - 4.
        2.    3.    4.

The values found by findMaxComponents are the largest possible coefficients of each lattice vector such that a linear combination with that coefficient exists which still land on the sphere. Since I'm looking for the largest such combinations with integer coefficients, I can safely drop the part after the decimal point to get the maximal plausible integer ranges. So for the given matrix A, I'll have to go from -2 to 2 in the first component, from -3 to 3 in the second and from -4 to 4 in the third and I'm sure to land on all the points inside the sphere (plus superfluous extra points, but importantly definitely every valid point inside) Next up:

---

Step 2: using the above information, generate all the candidate combinations.

    function K=findAllCombinations(I) //takes a matrix of the form produced by
                                      //findMaxComponents() and returns a matrix
                                      //which lists all the integer linear combinations
                                      //in the respective ranges.
        v=I(1,:); //starting from the minimal vector
        K=[];
        next=1; //keeps track of what component to advance next
        changed=%F; //keeps track of whether to add the vector to the output
        
        while or(v~=I(2,:)) //as long as not all components of v match all components of the maximum vector
            if v <= I(2,:) then //if each current component is smaller than each largest possible component
                if ~changed then
                    K=[K;v]; //store the vector and
                end
                v(next)=v(next)+1; //advance the component by 1
                next=1; //also reset next to 1
                changed=%F;
            else
                v(1:next)=I(1,1:next); //reset all components smaller than or equal to the current one and
                next=next+1; //advance the next larger component next time
                changed=%T;
            end
        end
        K=[K;I(2,:)]'; //while loop ends a single iteration early so add the maximal vector too
                       //also transpose K to fit better with the other functions
    endfunction

So now that I have that, all that remains is to check whether a given combination actually does lie inside or outside the sphere. All I gotta do for that is:

---

Step 3: Filter the combinations to find the actually valid lattice points

    function points=generatePoints(A,K,r)
        possiblePoints=A*K; //explicitly generates all the possible points
        points=[];
        for i=possiblePoints
            if i'*i<=r*r then //filter those that are too far from the origin
                points=[points i];
            end
        end
    endfunction

And I get all the combinations that actually do fit inside the sphere of radius r.

For the above example, the output is rather long: Of originally 315 possible points for a sphere of radius 3 I get 163 remaining points.

The first 4 are: (each column is one)

  - 0.2425356    0.2425356    1.2126781  - 0.9701425
  - 2.4253563  - 2.6678919  - 2.4253563  - 2.4253563
    1.6977494    0.           0.2425356    0.4850713

so the remainder of the work is optimization. Presumably some of those loops could be made faster and especially as the number of dimensions goes up, I have to generate an awful lot of points which I have to discard, so maybe there is a better way than taking the bounding n-parallelotope of the n-1-sphere as a starting point.