给定一组任意向量（存储在矩阵a中）和半径r，我想找到这些向量的所有整数值线性组合，这些向量落在半径r的球体内。然后我将把必要的坐标存储在矩阵V中。因此，例如，如果线性组合
K=[0; 1; 0]
降落在我的球体内，即类似的东西
if norm(A*K) <= r then
V(:,1)=K
end
等等。
A中的向量肯定是给定晶格的最简单可能的基，并且最大的向量的长度为1。不确定这是否以任何有用的方式限制了向量，但我怀疑可能。-它们不会像不太理想的基础那样有相似的方向。
我已经尝试了几种方法，但似乎没有一种特别令人满意。我似乎找不到一个好的模式来遍历晶格。
我目前的方法包括从中间开始（即所有0的线性组合），并逐个通过必要的坐标。它包括存储一堆额外的向量来跟踪，这样我就可以遍历坐标的所有八分之一（在3D的情况下），并逐一找到它们。这种实现似乎非常复杂，也不太灵活（尤其是它似乎不容易推广到任意数量的维度——尽管这对当前的目的来说不是绝对必要的，但拥有会很好）
有没有一种很好的方法可以找到所有需要的分数？
（*理想情况下，既高效又优雅**。如果真的有必要，在球体外多几个点也没什么关系，但最好不要多。我确实需要球体内的所有矢量。-如果这有很大的不同，我最感兴趣的是3D情况。
**我很确定我目前的实现两者都不是。）
我发现了类似的问题：
找到任意坐标周围半径为r的球体中的所有点——这实际上是一个比我所寻找的更普遍的情况。我只处理周期晶格，我的球体总是以0为中心，与晶格上的一点重合。但我没有点的列表，而是一个向量矩阵，我可以用它生成所有的点。
如何有效地枚举n维网格中球面的所有点——完全正则超三次格和曼哈顿距离的情况。我在寻找完全任意的格和欧氏距离（或者，为了效率的目的，显然是它的平方）。
在不证明任何断言的情况下，我认为1）如果向量集不是最大秩的，那么解的数量是无限的；2） 如果该集合具有最大秩，则由向量生成的线性变换的图像是目标空间的子空间（例如，平面），该子空间在低维球体中与球体相交；3） 可以将问题简化为1-1线性变换（k维空间上的kxk矩阵）；4） 由于矩阵是可逆的，你可以将球体"拉回"到包含格点的空间中的椭球体，作为奖励，你可以得到椭球体的良好几何描述（主轴定理）；5） 现在，您的问题变成了确定椭球体内部晶格点的问题。
后一个问题与高斯考虑的一个老问题（计算椭圆内的格点）有关，高斯导出了一个很好的近似值。确定椭圆（oid）内的晶格点可能不是一个整洁的问题，但它可能可以一次减少一个维度（椭球和平面的横截面是另一个椭球）。
我找到了一个让我现在更快乐的方法。可能还有可能的改进，所以如果你有更好的方法，或者在这段代码中发现错误，请分享。尽管这是我现在所拥有的：（全部写在SciLab中）
步骤1：计算出由与晶格向量的轴对齐的边界n-等位多面体定义的最大范围。感谢ElKamina含糊其辞的建议，以及chappers对我在math.se上的另一个问题的回答：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
在我的问题中，我只要求一个一般的基础，即对于n维，一组n个任意但线性独立的向量。通过使用伪逆，上述代码适用于任意形状的矩阵，类似地，Scilab的"逆"；CCD_ 1"；返回共轭转置，而不仅仅是A的转置，因此它同样适用于复杂矩阵。
在最后一步中，我放置了相应的最小组件。
以一个这样的A为例，这在Scilab的控制台中给了我以下内容：
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.
由findMaxComponents找到的值是每个晶格矢量的最大可能系数，使得与该系数的线性组合仍然存在于球体上。由于我正在寻找具有整数系数的最大此类组合，因此我可以安全地删除小数点后的部分，以获得最大的合理整数范围。因此，对于给定的矩阵A，我必须从第一个分量的-2到2，从第二个分量的-3到3，从第三个分量的-4到4，我确信会落在球体内的所有点上（加上多余的额外点，但重要的是，肯定是里面的每个有效点）下一步：
步骤2：使用上述信息，生成所有候选组合。
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
现在我已经知道了，剩下的就是检查给定的组合是否真的位于球体内部或外部。我所要做的就是：
步骤3：过滤组合以找到实际有效的格点
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
我得到了所有的组合，它们确实适合半径r的范围。
对于上面的例子，输出相当长：对于半径为3的球体，在最初的315个可能点中，我得到了163个剩余点。
前4列为：（每列为一列）
- 0.2425356 0.2425356 1.2126781 - 0.9701425
- 2.4253563 - 2.6678919 - 2.4253563 - 2.4253563
1.6977494 0. 0.2425356 0.4850713
所以剩下的工作就是优化。据推测，其中一些循环可以更快地进行，尤其是随着维度数量的增加，我必须生成大量必须丢弃的点，所以也许有一种比以n-1-球体的边界n-平行四边形作为起点更好的方法。
让我们只将K表示为X。
问题可以表示为：
（a11x1+a12x2.）^2+（a21x1+a22x2.）^ 2…<r^2
（x1，x2，…）将不会形成球体。
这可以通过维度递归来完成——选择一个晶格超平面方向并索引所有与r半径球相交的超平面。每个这样的超平面的球交点本身就是一个低维的球。重复以下是Octave中的调用函数代码：
function lat_points(lat_bas_mx,rr)
% **globals for hyperplane lattice point recursive function**
clear global; % this seems necessary/important between runs of this function
global MLB;
global NN_hat;
global NN_len;
global INP; % matrix of interior points, each point(vector) a column vector
global ctr; % integer counter, for keeping track of lattice point vectors added
% in the pre-allocated INP matrix; will finish iteration with actual # of points found
ctr = 0; % counts number of ball-interior lattice points found
MLB = lat_bas_mx;
ndim = size(MLB)(1);
% **create hyperplane normal vectors for recursion step**
% given full-rank lattice basis matrix MLB (each vector in lattice basis a column),
% form set of normal vectors between successive, nested lattice hyperplanes;
% store them as columnar unit normal vectors in NN_hat matrix and their lengths in NN_len vector
NN_hat = [];
for jj=1:ndim-1
tmp_mx = MLB(:,jj+1:ndim);
tmp_mx = [NN_hat(:,1:jj-1),tmp_mx];
NN_hat(:,jj) = null(tmp_mx'); % null space of transpose = orthogonal to columns
tmp_len = norm(NN_hat(:,jj));
NN_hat(:,jj) = NN_hat(:,jj)/tmp_len;
NN_len(jj) = dot(MLB(:,jj),NN_hat(:,jj));
if (NN_len(jj)<0) % NN_hat(:,jj) and MLB(:,jj) must have positive dot product
% for cutting hyperplane indexing to work correctly
NN_hat(:,jj) = -NN_hat(:,jj);
NN_len(jj) = -NN_len(jj);
endif
endfor
NN_len(ndim) = norm(MLB(:,ndim));
NN_hat(:,ndim) = MLB(:,ndim)/NN_len(ndim); % the lowest recursion level normal
% is just the last lattice basis vector
% **estimate number of interior lattice points, and pre-allocate memory for INP**
vol_ppl = prod(NN_len); % the volume of the ndim dimensional lattice paralellepiped
% is just the product of the NN_len's (they amount to the nested altitudes
% of hyperplane "paralellepipeds")
vol_bll = exp( (ndim/2)*log(pi) + ndim*log(rr) - gammaln(ndim/2+1) ); % volume of ndim ball, radius rr
est_num_pts = ceil(vol_bll/vol_ppl); % estimated number of lattice points in the ball
err_fac = 1.1; % error factor for memory pre-allocation--assume max of err_fac*est_num_pts columns required in INP
INP = zeros(ndim,ceil(err_fac*est_num_pts));
% **call the (recursive) function**
% for output, global variable INP (matrix of interior points)
% stores each valid lattice point (as a column vector)
clp = zeros(ndim,1); % confirmed lattice point (start at origin)
bpt = zeros(ndim,1); % point at center of ball (initially, at origin)
rd = 1; % initial recursion depth must always be 1
hyp_fun(clp,bpt,rr,ndim,rd);
printf("%i lattice points foundn",ctr);
INP = INP(:,1:ctr); % trim excess zeros from pre-allocation (if any)
endfunction
关于NN_len（jj）*NN_hat（：；平行六面体；由晶格基中的矢量MLB形成。格基平行六面体的体积仅为prod（NN_len）——为了快速估计内部格点的数量，将半径为rr的ndim球的体积除以prod（NN_len）。这是递归函数代码：
function hyp_fun(clp,bpt,rr,ndim,rd)
%{
clp = the lattice point we're entering this lattice hyperplane with
bpt = location of center of ball in this hyperplane
rr = radius of ball
rd = recrusion depth--from 1 to ndim
%}
global MLB;
global NN_hat;
global NN_len;
global INP;
global ctr;
% hyperplane intersection detection step
nml_hat = NN_hat(:,rd);
nh_comp = dot(clp-bpt,nml_hat);
ix_hi = floor((rr-nh_comp)/NN_len(rd));
ix_lo = ceil((-rr-nh_comp)/NN_len(rd));
if (ix_hi<ix_lo)
return % no hyperplane intersections detected w/ ball;
% get out of this recursion level
endif
hp_ix = [ix_lo:ix_hi]; % indices are created wrt the received reference point
hp_ln = length(hp_ix);
% loop through detected hyperplanes (updated)
if (rd<ndim)
bpt_new_mx = bpt*ones(1,hp_ln) + NN_len(rd)*nml_hat*hp_ix; % an ndim by length(hp_ix) matrix
clp_new_mx = clp*ones(1,hp_ln) + MLB(:,rd)*hp_ix; % an ndim by length(hp_ix) matrix
dd_vec = nh_comp + NN_len(rd)*hp_ix; % a length(hp_ix) row vector
rr_new_vec = sqrt(rr^2-dd_vec.^2);
for jj=1:hp_ln
hyp_fun(clp_new_mx(:,jj),bpt_new_mx(:,jj),rr_new_vec(jj),ndim,rd+1);
endfor
else % rd=ndim--so at deepest level of recursion; record the points on the given 1-dim
% "lattice line" that are inside the ball
INP(:,ctr+1:ctr+hp_ln) = clp + MLB(:,rd)*hp_ix;
ctr += hp_ln;
return
endif
endfunction
这里面有一些八分之一/匹配的东西，但大多数应该很容易理解；M（：，jj）引用矩阵M的列jj；tic的意思是转置；[AB]连接矩阵A和B；A=[]声明了一个空矩阵。
更新/更好地优化了原始答案：
- "矢量化的"；代码中的递归函数，以避免大多数"；对于"；循环（这些循环使它慢了大约10倍；不过现在代码有点难以理解）
- 为内部点的INP矩阵预先分配的内存（这将其加速了另一个数量级；在此之前，Octave必须为每次调用最内部递归级别调整INP矩阵的大小——对于真正会减慢速度的大型矩阵/数组）
因为这个例程是项目的一部分，所以我也用Python对它进行了编码。从非正式测试来看，Python版本比这个（Octave）版本快2-3倍。
作为参考，这里是这个答案最初发布的旧的、慢得多的代码：
% (OLD slower code, using for loops, and constantly resizing
% the INP matrix) loop through detected hyperplanes
if (rd<ndim)
for jj=1:length(hp_ix)
bpt_new = bpt + hp_ix(jj)*NN_len(rd)*nml_hat;
clp_new = clp + hp_ix(jj)*MLB(:,rd);
dd = nh_comp + hp_ix(jj)*NN_len(rd);
rr_new = sqrt(rr^2-dd^2);
hyp_fun(clp_new,bpt_new,rr_new,ndim,rd+1);
endfor
else % rd=ndim--so at deepest level of recursion; record the points on the given 1-dim
% "lattice line" that are inside the ball
for jj=1:length(hp_ix)
clp_new = clp + hp_ix(jj)*MLB(:,rd);
INP = [INP clp_new];
endfor
return
endif