The best way to find all the lattice points in a sphere

Question

Given a bunch of arbitrary vectors (stored in matrix A) and radius r, I would like to find all integer signed linear combinations of those vectors that fall into the sphere of radius r. I would then save the necessary coordinates in the matrix V. So, for example, if a linear combination

K=[0; 1; 0] 

the earth inside my sphere, that is, something like

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

and etc.

The vectors in A are undoubtedly the simplest basis for this lattice, and the largest vector will have a length of 1. Not sure if this limits the vectors in any useful way, but I suspect that this is possible. - They will not have the same directions as the less ideal base.

I have tried several approaches already, but none of them seem particularly satisfactory. I can't seem to find a good pattern to cross the grid.

My current approach involves starting in the middle (i.e. with a linear combination of all 0s) and alternately go through the necessary coordinates. It includes saving many additional vectors for tracking, so I can go through all the octants (in the 3D case) of the coordinates and find them one at a time. This implementation seems terribly complex and not very flexible (in particular, it does not seem to be easily generalizable for arbitrary measurement numbers, although this is not strictly necessary for the current purpose, it would be nice)

Is there a good way to find all the necessary points?

(* Ideally effective and elegant **. If REALLY necessary, it doesn’t matter to have a lot of extra points outside the sphere, but preferably not so much. I definitely need all the vectors inside the sphere - if that matters a lot, I’m most interested in a 3D case.

** I am sure my current implementation is not a single one.)

---

Related questions I found:

Finding all the points in a sphere of radius r around an arbitrary coordinate is actually a much more general case than what I'm looking for. I deal only with periodic lattices, and my sphere is always focused on 0, coinciding with one point on the lattice. But I do not have a list of points, but rather a matrix of vectors with which I can generate all points.

How to efficiently list all points of a sphere in an n-dimensional grid is a case for a completely regular hypercubic lattice and Manhattan distance. I am looking for completely arbitrary lattices and Euclidean distance (or, for efficiency, obviously, the square of this).

Answer 1

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

---

Step 1: Find out the maximum ranges. Thanks for ElKamina for the vague suggestion, as well as for answering another question on math.se from chappers: https://math.stackexchange.com/a/1230160/49989

 function I=findMaxComponents(A,r) //returns matrix I, takes matrix A, scalar r [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^TA iAT=pinv(A'); //finds the (pseudo-)inverse of A^T I=[]; //initializes I as an empty vector for i=1:vecs t=r*(iATA*U(:,i))/norm(iAT*U(:,i)) //put it all together as per above link I=[I,t(i)]; //take only the maximized component and store it in I. end I=[-I;I]; //I want to go from minimum to maximum value. endfunction 

In my question, I asked only for a general basis, i.e. for n dimensions, a set of n arbitrary, but linearly independent vectors. The above code, due to the use of pseudo-inversions, works for arbitrary-form matrices and, similarly, Scilab "A" returns conjugate transposition, and not just transposition of A, so it should work the same for complex matrices.

In the last step, I put the appropriate minimal components.

For one such A, as an example, this gives me the following in the Scilab 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. 

For each component, the above values ​​are the maximum ones that still land on the sphere, so I can safely reset the part after the decimal point to get the maximum integers. Therefore, for a given matrix A, I will 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 that I will land on all points inside the sphere. Further:

---

Step 2: Using the information above, find all the 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 it, it remains only to check whether this combination really lies inside or outside the sphere. All I have to do for this:

---

Step 3: Create Actual 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 fall within the scope of radius r.

In the above example, the output is quite long: initially 315 possible points for a sphere of radius 3, I get 163 remaining points.

First 4: (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 rest of the work is optimization. Presumably, some of these cycles could be done faster, and especially as the number of measurements increases, I have to generate a lot of points that I have to drop, so there may be a better way than to take (distorted at this coordinate frame) bounding the hyperfield of the n-sphere as a starting point.

Answer 2

Imagine only K as X.

The problem can be represented as:

(a11x1 + a12x2 ..) ^ 2 + (a21x1 + a22x2 ..) ^ 2 ... <r ^ 2

(x1, x2, ...) will not form a sphere.

Answer 3

Without taking into account any statements, I believe that 1) if the set of vectors does not have a maximum rank, then the number of solutions is infinite; 2) if the set has the maximum rank, then the image of the linear transformation generated by the vectors is a subspace (for example, a plane) of the target space that intersects the sphere in a sphere with a lower dimension; 3) it follows that you can reduce the problem to a linear transformation 1-1 (kxk-matrix on k-dimensional space); 4) since the matrix is ​​reversible, you can “pull” the sphere to an ellipsoid in the space containing the lattice points, and as a bonus you will get a good geometric description of the ellipsoid (main axis theorem); 5) your problem now becomes an exact definition of the lattice points inside the ellipsoid.

The last problem is related to the old problem (counting the lattice points inside the ellipse), which was examined by Gaus, which received a good approximation. Determining the lattice points inside an ellipse (oid) is probably not such a neat task, but it can probably be reduced one dimension at a time (the cross section of the ellipsoid and the plane is another ellipsoid).