今のところもっと幸せになる方法を見つけました。まだ改善の余地があるかもしれませんので、より良い方法がある場合、またはこのコードでエラーを見つけた場合は、ぜひ共有してください。ここに私が今持っているものがあります：（すべて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 の " A'" は、単に転置するのではなく、共役転置を返すAため、複素行列に対しても同様に機能するはずです。
最後のステップでは、対応する最小限のコンポーネントを配置します。
例としてそのような A の 1 つを使用すると、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 に、2 番目のコンポーネントで -3 から 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 つ: (各列は 1 つ)
- 0.2425356 0.2425356 1.2126781 - 0.9701425
- 2.4253563 - 2.6678919 - 2.4253563 - 2.4253563
1.6977494 0. 0.2425356 0.4850713
したがって、残りの作業は最適化です。おそらく、これらのループのいくつかはより高速に作成でき、特に次元数が増えるにつれて、破棄しなければならない非常に多くのポイントを生成する必要があるため、バウンディングを取得するよりも良い方法があるかもしれませnんn-1-開始点としての球。