4.改進LLL方法－遞迴LLL

原本LLL方法是針對整個Lattice進行LLL運算，而遞迴LLL一開始先對前兩個向量進行LLL運算，之後再加入第三個向量進行LLL運算...，每次都新增一個向量進行LLL運算直到整個Lattice都化簡為止。

範例出自https://thomas-plantard.github.io/pdf/talk10-scn.pdf

這個LLL副程式略有修改，參數m代表將第1~m列向量進行LLL運算。
(%i1)
LLL(v,m):=block
  (vstar:copymatrix(v),
   delta:3/4,
   mu:zeromatrix(m,m),
   for i:1 thru m do
     (for j:1 thru i do
        (mu[ i ][j]:v[ i ].vstar[j]/vstar[j].vstar[j],
         vstar[ i ]:v[ i ]-sum(mu[ i ][j]*vstar[j],j,1,i-1)
        )
     ),
   k:2,
   while k<=m do
     (for i:k-1 thru 1 step -1 do
        (q:round(mu[k][ i ]),
         if q#0 then
           (v:rowop(v,k,i,q),
            mu:rowop(mu,k,i,q)
           )
        ),
      if (vstar[k].vstar[k])>=(delta-mu[k][k-1]^2)*(vstar[k-1].vstar[k-1]) then
        (k:k+1
        )
      else
        (v:rowswap(v,k,k-1),
         vstar[1]:v[1],
         for i:1 thru m do
           (for j:1 thru i do
             (mu[ i ][j]:v[ i ].vstar[j]/vstar[j].vstar[j],
              vstar[ i ]:v[ i ]-sum(mu[ i ][j]*vstar[j],j,1,i-1)
             )
           ),
         if k>2 then
           (k:k-1
           )
        )
     ),
  return(v)
)$

(%i2)
B:matrix([1031,0,0,0,0],
         [354,1,0,0,0],
         [322,0,1,0,0],
         [916,0,0,1,0],
         [426,0,0,0,1]);
(%o2)　103135432291642601000001000001000001

(%i3)
for i:2 thru length(B) do
  (B: LLL(B,i),
   print(B)
  );
−3113322916426−3−32000001000001000001

−10−279164262−6−1100113−4000001000001

−1−4404262−30−3014540025−7000001

−10−2−232−2−210123−1301−3400−2115
(%o3)　done

另一個範例出自http://www1.uwindsor.ca/sac2012/system/files/2a_Plantard.pdf
(%i4)
B:matrix([86670401,0,0,0,0,0,0,0],
         [38009011,1,0,0,0,0,0,0],
         [10117311,0,1,0,0,0,0,0],
         [38269415,0,0,1,0,0,0,0],
         [45874978,0,0,0,1,0,0,0],
         [33538152,0,0,0,0,1,0,0],
         [61611560,0,0,0,0,0,1,0],
         [66174289,0,0,0,0,0,0,1]);
(%o4)　866704013800901110117311382694154587497833538152616115606617428901000000001000000001000000001000000001000000001000000001

(%i5)
for i:2 thru length(B) do
  (B: LLL(B,i),
   print(B)
  );
−3227−14111101173113826941545874979335381546161156066174289−316513018000000001000000001000000001000000001000000001000000001

−24183−920382664154587497833538152616115606617428915324244000000−21576343000000001000000001000000001000000001000000001

8−40010345874978335381526161156066174289−27−45126260000−66−47−180000042−38−23−53000000001000000001000000001000000001

−1724−3813335381526161156066174289−31−20−740000−5−6−45−1315000643−14−350001717−3624000000001000000001000000001

−7−14−28−546161156066174289−6−915−1411−1600−14−3145−6120011−11013−10−4004−15−6−14−121200−7−64−21213000000001000000001

591434766174289−4−1−48−2−9−100−840−4−11−7500−669−49701−7−124−53−2010−8−10−3−2−10−124−3−67−4000000001

−223−254−256−30−1−42311−1−8−3−643−41−81−54−1−10303−3−6−2−40−21−12−301113−3450−7−12−2115−31−7−24
(%o5)　done




執行結果第三列向量和pdf檔所寫的不一樣，我還以為我的LLL程式出了問題，我還另外用maple再寫一次來測試，結果是相同的。

maple檔下載遞迴LLL1.rar
−223−254−256−30−1−42311−1−8−3−643−41−81−54−1−10303−3−6−2−40−21−12−301113−3450−7−12−2115−31−7−24−223−254−256−3−3−1−42311−1−86−643−41−8124−1−10303−6−6−2−40−21−12−301113−3440−7−12−2113−31−7−24





＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
在 http://www1.uwindsor.ca/sac2012/system/files/2a_Plantard.pdf 還提到另外一種遞迴LLL方法，先將兩個向量為一組進行LLL運算，再以四個向量為一組進行LLL運算，再以八個向量為一組進行LLL運算…，直到全部的Lattice向量都化簡完為止。

這個LLL副程式略有修改，參數m1m2代表將第m1~m2列向量進行LLL運算。
(%i1)
LLL(v,m1,m2):=block
  (vstar:copymatrix(v),
   delta:3/4,
   m:length(v),
   mu:zeromatrix(m,m),
   for i:m1 thru m2 do
     (for j:m1 thru i do
        (mu[ i ][j]:v[ i ].vstar[j]/vstar[j].vstar[j],
         vstar[ i ]:v[ i ]-sum(mu[ i ][j]*vstar[j],j,1,i-1)
        )
     ),
   k:m1+1,
   while k<=m2 do
     (for i:k-1 thru 1 step -1 do
        (q:round(mu[k][ i ]),
         if q#0 then
           (v:rowop(v,k,i,q),
            mu:rowop(mu,k,i,q)
           )
        ),
      if (vstar[k].vstar[k])>=(delta-mu[k][k-1]^2)*(vstar[k-1].vstar[k-1]) then
        (k:k+1
        )
      else
        (v:rowswap(v,k,k-1),
         vstar[m1]:v[m1],
         for i:m1 thru m2 do
           (for j:m1 thru i do
             (mu[ i ][j]:v[ i ].vstar[j]/vstar[j].vstar[j],
              vstar[ i ]:v[ i ]-sum(mu[ i ][j]*vstar[j],j,1,i-1)
             )
           ),
         if k>m1+1 then
           (k:k-1
           )
        )
     ),
  return(v)
)$

(%i2)
B:matrix([86670401,0,0,0,0,0,0,0],
         [38009011,1,0,0,0,0,0,0],
         [10117311,0,1,0,0,0,0,0],
         [38269415,0,0,1,0,0,0,0],
         [45874978,0,0,0,1,0,0,0],
         [33538152,0,0,0,0,1,0,0],
         [61611560,0,0,0,0,0,1,0],
         [66174289,0,0,0,0,0,0,1]);
(%o2)　866704013800901110117311382694154587497833538152616115606617428901000000001000000001000000001000000001000000001000000001

(%i3)
d:length(B)$
for i:1 thru ceiling(log(d)/log(2)) do
  (for j:1 thru ceiling(d/2^i) do
     (m1: (j-1)*2^i+1,
      m2:j*2^i,
      if m2>d then
        (m2:d),
      B: LLL(B,m1,m2),
      print("第",m1,"到",m2,"列向量進行LLL運算"),
      print(B)
     )
  );
第1到2列向量進行LLL運算
−3227−14111101173113826941545874978335381526161156066174289−316513018000000001000000001000000001000000001000000001000000001

第3到4列向量進行LLL運算
−3227−141111391−812145874978335381526161156066174289−3165130180000000040363949000000−1067−1044000000001000000001000000001000000001

第5到6列向量進行LLL運算
−3227−141111391−81211348108946161156066174289−3165130180000000040363949000000−1067−1044000000002830−2009000000−38712748000000001000000001

第7到8列向量進行LLL運算
−3227−141111391−8121134810894329437−3165130180000000040363949000000−1067−1044000000002830−2009000000−3871274800000000224829000000−2093−27

第1到4列向量進行LLL運算
8−400−103134810894329437−27−45126−260000−66−47−180000042−38−2353000000002830−2009000000−3871274800000000224829000000−2093−27

第5到8列向量進行LLL運算
8−400−103−22−9332−1−27−45126−260000−66−47−180000042−38−235300000000−72597250000−192−131110000−365−63−93000048−23−213

第1到8列向量進行LLL運算
−2−32−4−2−535301−264−111836−3−1−401−15−41−81−13−336402−2−21−23−110113−4−501−37−5−2−1−5371−1−14
(%o4)　done




紅色數字都是執行結果和pdf檔所寫的不一樣的地方，只是再用maple測試，執行結果和我的相同，或許原作者所使用的LLL程式多了其他的設定(=075改為=099)

maple下載遞迴LLL2.rar
−2−32−4−2−535301−264−111836−3−1−401−15−41−81−13−336402−2−21−23−110113−4−501−37−5−2−1−5371−1−14−223−254−256−3−3−1−42311−1−86−643−41−8124−1−10303−6−6−2−40−21−12−301113−3440−7−12−2113−31−7−24



參考資料：
Recursive Lattice Reduction，Thomas Plantard, Willy Susilo
https://thomas-plantard.github.io/pdf/PlaSus10.pdf
簡報檔
https://thomas-plantard.github.io/pdf/talk10-scn.pdf

Lattice Reduction for Modular Knapsack，Thomas Plantard, Willy Susilo, Zhenfei Zhang
https://thomas-plantard.github.io/pdf/PlaSusZha12.pdf
簡報檔
http://www1.uwindsor.ca/sac2012/system/files/2a_Plantard.pdf

5.改進LLL方法－PotLLL

原本的LLL只考慮相鄰的兩列向量( vi、vi−1 )是否滿足Lovász condition( vi2(43−2ii−1)vi−12 )，違反的話則將vk−1vk交換。而PotLLL考慮是否滿足(Pot(B)Pot(klB))，違反的話則將第l列向量插入第k列向量的位置。


kl(i)=ili−1for　i<k or i>lfor　i=k,　for　k<i ≦ l.代表將第l列向量插入第k列向量。
B是論文用來代表lattice，我的maxima程式用v表示。
Pot(B):=ni=1vol(L(b1bi))2=ni=1bi2(n−i+1) 計算lattice的Potential值
以上式子就是將lattice B分成很多sub lattice
B1=L(b1)，由b1向量所組成的lattice，其體積的平方為b12
B_2=L(b_1,b_2) ，由 b_1,b_2 向量組成的lattice，其體積的平方為 \Vert\; b_1^* \Vert\;^2 \cdot \Vert\; b_2^* \Vert\;^2
\ldots
B_n=L(b_1,b_2,\ldots,b_n) ，由 b_1,b_2,\ldots,b_n 向量組成的lattice，其體積的平方為 \Vert\; b_1^* \Vert\;^2 \cdot \Vert\; b_2^* \Vert\;^2 \ldots \Vert\; b_n^* \Vert\;^2
將全部sub lattice體積的平方全部乘起來 \displaystyle \prod_{i=1}^n \Vert\; b_i^* \Vert\;^{2(n-i+1)} 就是Potential值


以 B=\left[ \matrix{9&2&7 \cr 8&6&1 \cr 3&2&6 } \right] 為例，計算 B 的Potential值
先計算GramSchmit正交化後的向量 b^*=\left[ \matrix{\displaystyle 9&2&7\cr \frac{253}{134}&\frac{311}{67}&-\frac{503}{134}\cr -\frac{8080}{5253}&\frac{9494}{5253}&\frac{7676}{5253}} \right] ，注意各向量彼此正交(垂直)。

各向量長度的平方為 \matrix{\displaystyle \Vert\; b_1^*\Vert\;^2=9^2+2^2+7^2=134　　　　　　　　　　 \cr \Vert\; b_2^*\Vert\;^2=\left( \frac{253}{134} \right)^2+\left( \frac{311}{67} \right)^2+\left( -\frac{503}{134} \right)^2=\frac{5253}{134}　 \cr \Vert\; b_3^*\Vert\;^2=\left( -\frac{8080}{5253} \right)^2+\left( \frac{9494}{5253} \right)^2+\left( \frac{7676}{5253} \right)^2=\frac{40804}{5253}}
B_1=\left[ \matrix{9&2&7} \right] ， b_1 向量長度就是 \Vert\; b_1^* \Vert\; ，平方為 \Vert\; b_1^* \Vert\;^2=134
B_2=\left[ \matrix{9&2&7 \cr 8&6&1} \right] ， b_1,b_2 所張成的平行四邊形面積等同於 b_1^*長、b_2^* 寬所張成的長方形面積，平方為 \displaystyle \Vert\; b_1^* \Vert\;^2 \cdot \Vert\; b_2^* \Vert\;^2=134 \cdot \frac{5253}{134}=5253
B_3=\left[ \matrix{9&2&7 \cr 8&6&1 \cr 3&2&6} \right] ， b_1,b_2,b_3 所張成的平行六面體體積等同於 b_1^*長、 b_2^*寬、 b_3^* 高所張成的長方體體積，平方為 \displaystyle \Vert\; b_1^* \Vert\;^2 \cdot \Vert\; b_2^* \Vert\;^2 \cdot \Vert\; b_3^* \Vert\;^2=134 \cdot \frac{5253}{134} \cdot \frac{40804}{5253}=40804
將全部sub lattice體積的平方全部乘起來 134 \cdot 5253 \cdot 40804=28722017208 就是Potential值


＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
以下程式計算Potential值

要先載入eigen.mac才能使用gramschmidt指令
(%i1)　load("eigen.mac");
(%o1)　C:/Program Files (x86)/Maxima-sbcl-5.37.1/share/maxima/5.37.1/share/matrix/eigen.mac

計算Potential值副程式
(%i2)
Pot(v):=block
(vstar:gramschmidt(v),
n:length(v[1]),
product((vstar[ i ].vstar[ i ])^(n-i+1),i,1,n)
)$

計算v的Potential值
(%i3)
v:matrix([9,2,7],
              [8,6,1],
              [3,2,6]);
Pot(v);
(%o3)　 \left[ \matrix{9&2&7 \cr 8&6&1 \cr 3&2&6} \right]
(%o4)　28722017208



＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
以下程式為PotLLL執行結果

GramSchmit正交化
(%i1)
GramSchmit(v):=block
([m:length(v)],
vstar:copymatrix(v),
mu:zeromatrix(m,m),
for i:1 thru m do
  (for j:1 thru i do
     (mu[ i ][j]:v[ i ].vstar[j]/vstar[j].vstar[j],
      vstar[ i ]:v[ i ]-sum(mu[ i ][j]*vstar[j],j,1,i-1)
     )
  )
)$

檢查是否違反Size condition
(%i2)
SizeReduce(i,j):=block
([q:round(mu[ i ][j])],
if q#0 then
   (print("round(μ[",i,"][",j,"])=",round(mu[ i ][j]),"違反Size condition"),
    print("v=",v,"  第",i,"列-",q,"*第",j,"列  v=",rowop(v,i,j,q)),
    v:rowop(v,i,j,q)
   )
)$

第l列移動到第k列
(%i3)
MoveRow(v,k,l) := block
([n:length(v)],
head(v, n) :=if n > 0 then head(submatrix(length(v), v), n - 1) else v,
tail(v, n) :=if n > 0 then tail(submatrix(1, v), n - 1) else v,
v:addrow(head(v, n - k+1), v[l], tail(v, k-1)),
if k<l then v:submatrix(l+1,v) else v:submatrix(l,v)
)$

PotLLL主程式
(%i4)
PotLLL(v):=block(
[l:1,n:length(v),delta:3/4],
while l<=n do
  (if l>1 then
     (print("第",l,"列進行SizeReduce,μ=",mu)
     ),
   for i:1 thru l-1 do
     (SizeReduce(l,i)
     ),
   GramSchmit(v),
   P:1,
   Pmin:1,
   k:1,
   for j:l-1 thru 1 step -1 do
     (P: P*(vstar[l].vstar[l]+sum(mu[l,i]^2*vstar[ i ].vstar[ i ],i,j,l-1))/(vstar[j].vstar[j]),
      if P< Pmin then
        (k:j,
         Pmin: P
        )
     ),
   print("第",l,"列移到第",k,"列有最小的potential比值",Pmin),
   if delta> Pmin then
     (print("potential比值",Pmin,"比",delta,"小"),
      tempv:v,
      v:MoveRow(v,k,l),
      print(tempv,"第",l,"列移到第",k,"列",v),
      GramSchmit(v),
      print("l值",l,"=>",k),
      l:k
     )
   else
     (print("potential比值",Pmin,"比",delta,"大,l值",l,"=>",l+1),
      l:l+1
     )
  ),
return(v)
)$

(%i6)
v:matrix([9,2,7],
              [8,6,1],
              [3,2,6]);
PotLLL(v);
(%o5)　 \left[ \matrix{9&2&7 \cr 8&6&1 \cr 3&2&6} \right]
第1列移到第1列有最小的potential比值1
potential比值1比 \displaystyle \frac{3}{4} 大,l值1=>2

第2列進行SizeReduce,μ= \left[ \matrix{\displaystyle 1&0&0 \cr \frac{91}{134}&1&0 \cr \frac{73}{134}&-\frac{1015}{5253}&1} \right]
round(μ[2][1])=1違反Size condition
v=\left[ \matrix{9&2&7 \cr 8&6&1 \cr 3&2&6} \right]   第2列-1*第1列   \left[ \matrix{9&2&7 \cr -1&4&-6 \cr 3&2&6} \right]

第2列移到第1列有最小的potential比值 \displaystyle \frac{53}{134}
potential比值 \displaystyle \frac{53}{134} 比 \displaystyle \frac{3}{4} 小
\left[ \matrix{9&2&7\cr -1&4&-6\cr 3&2&6} \right] 第2列移到第1列 \left[ \matrix{-1&4&-6 \cr 9&2&7 \cr 3&2&6} \right]
l值2=>1

第1列移到第1列有最小的potential比值1
potential比值1比 \displaystyle \frac{3}{4} 大,l值1=>2

第2列進行SizeReduce, \mu=\left[ \matrix{\displaystyle 1&0&0\cr -\frac{43}{53}&1&0\cr -\frac{31}{53}&\frac{2536}{5253}&1} \right]
round(μ[2][1])=-1違反Size condition
v=\left[ \matrix{-1&4&-6\cr 9&2&7\cr 3&2&6} \right]   第2列 - -1* 第1列   v=\left[ \matrix{-1&4&-6\cr 8&6&1\cr 3&2&6} \right]

第2列移到第1列有最小的potential比值1
potential比值1比 \displaystyle \frac{3}{4} 大,l值2=>3

第3列進行SizeReduce, \mu=\left[ \matrix{\displaystyle 1&0&0\cr \frac{10}{53}&1&0\cr -\frac{31}{53}&\frac{2536}{5253}&1} \right]
round(μ[3][1])=-1違反Size condition
v=\left[ \matrix{-1&4&-6\cr 8&6&1\cr 3&2&6} \right]   第3列 - -1* 第1列   v=\left[ \matrix{-1&4&-6\cr 8&6&1\cr 2&6&0} \right]

第2列移到第1列有最小的potential比值1
potential比值1比 \displaystyle \frac{3}{4} 大,l值2=>3

第3列進行SizeReduce, \mu=\left[ \matrix{\displaystyle 1&0&0\cr -\frac{9}{20}&1&0\cr \frac{13}{10}&-\frac{186}{409}&1} \right]
round(μ[3][1])=1違反Size condition
v=\left[ \matrix{2&6&0\cr -3&-2&-6\cr 8&6&1} \right]   第3列 -1* 第1列   v=\left[ \matrix{2&6&0\cr -3&-2&-6\cr 6&0&1} \right]

第3列移到第1列有最小的potential比值 \displaystyle \frac{65440}{278409}
potential 比值 \displaystyle \frac{65440}{278409} 比 \displaystyle \frac{3}{4} 小
\left[ \matrix{-1&4&-6 \cr 8&6&1 \cr 2&6&0} \right] 第3列移到第1列  \left[ \matrix{2&6&0 \cr -1&4&-6 \cr 8&6&1} \right]
l值3=>1

第1列移到第1列有最小的potential比值1
potential比值1比 \displaystyle \frac{3}{4} 大,l值1=>2
第2列進行SizeReduce, \mu=\left[ \matrix{\displaystyle 1&0&0 \cr \frac{11}{20}&1&0 \cr \frac{13}{10}&-\frac{186}{409}&1} \right]
round(μ[2][1])=1違反Size condition
v=\left[ \matrix{2&6&0 \cr -1&4&-6 \cr 8&6&1} \right]   第2列-1*第1列   v=\left[ \matrix{2&6&0 \cr -3 & -2 &-6 \cr 8&6&1} \right]

第2列移到第1列有最小的potential比值1
potential比值1比 \displaystyle \frac{3}{4} 大,l值2=>3
第3列進行SizeReduce, \mu=\left[ \matrix{\displaystyle 1&0&0 \cr -\frac{9}{20}&1&0 \cr \frac{13}{10}&-\frac{186}{409}&1} \right]
round(μ[3][1])=1違反Size condition
v=\left[ \matrix{2&6&0 \cr -3&-2&-6 \cr 8&6&1} \right]   第3列-1*第1列   v=\left[ \matrix{2&6&0 \cr -3&-2&-6 \cr 6&0&1} \right]
第3列移到第1列有最小的potential比值 \displaystyle \frac{6179}{8180}
potential比值 \displaystyle \frac{6179}{8180} 比 \displaystyle \frac{3}{4} 大,l值3=>4
(%o6)　 \left[ \matrix{2&6&0 \cr -3&-2&-6 \cr 6&0&1} \right]


參考資料
PotLLL: A Polynomial Time Version of LLL With Deep Insertions
https://arxiv.org/abs/1307.7534
PDF檔https://arxiv.org/pdf/1307.7534v1.pdf

6-1.改進LLL方法－Jacobi方法

之前的Gaussian/Lagrange方法針對兩個向量進行化簡
https://math.pro/db/viewthread.php?tid=1876&page=1#pid10232
化簡後的向量符合(1)  \displaystyle \Bigg\vert\; \frac{v_1 \cdot v_2}{v_1 \cdot v_1} \Bigg\vert\;\le \frac{1}{2} 　(2) \Vert\; v_1 \Vert\; \le \Vert\; v_2 \Vert\; 兩個條件
擴展到更多向量就是將任兩個向量進行Gaussian/Lagrange化簡，任兩個向量滿足
(1)  \displaystyle \Bigg\vert\; \frac{v_i \cdot v_j}{v_i \cdot v_i} \Bigg\vert\;\le \frac{1}{2} 　(2) \Vert\; v_i \Vert\; \le \Vert\; v_j \Vert\; ， 1 \le i<j\le n 兩個條件



(%i1)
Jacobi(v):=block
([i,j,n:length(v),t,End],
do
  (End:true,
   for i:1 thru n-1 do
     (for j:i+1 thru n do
        (if v[ i ].v[ i ]>v[j].v[j] then
           (print("第",i,"列長度",sqrt(v[ i ].v[ i ]),"比第",j,"列長度",sqrt(v[j].v[j]),"長"),
            print(v,"第",i,"列和第",j,"列交換",v:rowswap(v,i,j)),
            End:false
           ),
        t:round(v[ i ].v[j]/v[ i ].v[ i ]),
        if t#0 then
           (print("t=round(v",i,".v",j,"/v",i,".v",i,")=round(",v[ i ].v[j],"/",v[ i ].v[ i ],")=",t,"不等於0"),
            print(v,"第",j,"列-",t,"*第",i,"列",v:rowop(v,j,i,t)),
            End:false
           )
        )
     ),
   if End=true then
     (print("化簡後v=",v,
             ",",vi.vj/vi.vi,"=",genmatrix(lambda([i,j],if i<j then v[ i ].v[j]/v[ i ].v[ i ] else "*"),n,n),
             ",向量長度=",genmatrix(lambda([i,j],sqrt(v[ i ].v[ i ])),n,1)),
      print("所有向量符合(1)|vi.vj/vi.vi|≦1/2，(2)∥vi∥≦∥vj∥，1≦i<j≦n兩個條件，程式結束"),
      return(v)
     ),
   print("------------------------")
  )
)$

(%i2)
v:matrix([1,1,1],
             [-1,0,2],
             [3,5,6]);
(%o2)　 \left[ \matrix{1&1&1\cr -1&0&2\cr 3&5&6} \right]

(%i3)　Jacobi(v);
\displaystyle t=round(v1.v3/v1.v1)=round(14/3)=5不等於0
\left[ \matrix{1&1&1\cr -1&0&2\cr 3&5&6} \right] 第3列-5*第1列 \left[ \matrix{1&1&1\cr -1&0&2\cr -2&0&1} \right]
\displaystyle t=round(v2.v3/v2.v2)=round(4/5)=1不等於0
\left[ \matrix{1&1&1\cr -1&0&2\cr -2&0&1} \right] 第3列-1*第2列 \left[ \matrix{1&1&1\cr -1&0&2\cr -1&0&-1} \right]
------------------------
第1列長度 \sqrt{3}比第3列長度 \sqrt{2}長
\left[ \matrix{1&1&1\cr -1&0&2\cr -1&0&-1} \right] 第1列和第3列交換 \left[ \matrix{-1&0&-1\cr -1&0&2\cr 1&1&1} \right]
\displaystyle t=round(v1.v3/v1.v1)=round(-2/2)=-1不等於0
\left[ \matrix{-1&0&-1\cr -1&0&2\cr 1&1&1} \right] 第3列--1*第1列 \left[ \matrix{-1&0&-1\cr -1&0&2\cr 0&1&0} \right]
第2列長度 \sqrt{5}比第3列長度1長
\left[ \matrix{-1&0&-1\cr -1&0&2\cr 0&1&0} \right] 第2列和第3列交換 \left[ \matrix{-1&0&-1\cr 0&1&0\cr -1&0&2} \right]
------------------------
第1列長度\sqrt{2}比第2列長度1長
\left[ \matrix{-1&0&-1\cr 0&1&0\cr -1&0&2} \right] 第1列和第2列交換 \left[ \matrix{0&1&0\cr -1&0&-1\cr -1&0&2} \right]
------------------------
化簡後 v=\left[ \matrix{0&1&0\cr -1&0&-1\cr -1&0&2} \right] , \frac{vi.vj}{vi^2}=\left[ \matrix{*&0&0\cr *&*&-\frac{1}{2} \cr *&*&*} \right] ,向量長度 =\left[ \matrix{1 \cr \sqrt{2} \cr \sqrt{5}} \right]
所有向量符合(1)|vi.vj/vi.vi|≦1/2，(2)∥vi∥≦∥vj∥，1≦i<j≦n兩個條件，程式結束

(%o3)　 \left[ \matrix{0&1&0\cr -1&0&-1\cr -1&0&2} \right]

(%i4)
v:matrix([1,1,7,2],
         [9,8,4,6],
         [1,8,5,7],
         [2,3,1,1]);
(%o4)　 \left[ \matrix{1&1&7&2\cr 9&8&4&6\cr 1&8&5&7\cr 2&3&1&1} \right]

(%i5)　Jacobi(v);
\displaystyle t=round(v1.v2/v1.v1)=round(57/55)=1不等於0
\left[ \matrix{1&1&7&2\cr 9&8&4&6\cr 1&8&5&7\cr 2&3&1&1}\right] 第2列-1*第1列 \left[ \matrix{1&1&7&2\cr 8&7&-3&4\cr 1&8&5&7\cr 2&3&1&1} \right]
\displaystyle t=round(v1.v3/v1.v1)=round(58/55)=1不等於0
\left[ \matrix{1&1&7&2\cr 8&7&-3&4\cr 1&8&5&7\cr 2&3&1&1}\right] 第3列-1*第1列 \left[ \matrix{1&1&7&2\cr 8&7&-3&4\cr 0&7&-2&5\cr 2&3&1&1} \right]
第1列長度 \sqrt{55}比第4列長度 \sqrt{15}長
\left[ \matrix{1&1&7&2\cr 8&7&-3&4\cr 0&7&-2&5\cr 2&3&1&1}\right] 第1列和第4列交換 \left[ \matrix{2&3&1&1\cr 8&7&-3&4\cr 0&7&-2&5\cr 1&1&7&2} \right]
\displaystyle t=round(v1.v4/v1.v1)=round(14/15)=1不等於0
\left[ \matrix{2&3&1&1\cr 8&7&-3&4\cr 0&7&-2&5\cr 1&1&7&2}\right] 第4列-1*第1列 \left[ \matrix{2&3&1&1\cr 8&7&-3&4\cr 0&7&-2&5\cr -1&-2&6&1} \right]
第2列長度 \sqrt{138}比第3列長度 \sqrt{78}長
\left[ \matrix{2&3&1&1\cr 8&7&-3&4\cr 0&7&-2&5\cr -1&-2&6&1}\right] 第2列和第3列交換 \left[ \matrix{2&3&1&1\cr 0&7&-2&5\cr 8&7&-3&4\cr -1&-2&6&1} \right]
\displaystyle t=round(v2.v3/v2.v2)=round(75/78)=1不等於0
\left[ \matrix{2&3&1&1\cr 0&7&-2&5\cr 8&7&-3&4\cr -1&-2&6&1}\right] 第3列-1*第2列 \left[ \matrix{2&3&1&1\cr 0&7&-2&5\cr 8&0&-1&-1\cr -1&-2&6&1} \right]
第2列長度 \sqrt{78}比第4列長度 \sqrt{42}長
\left[ \matrix{2&3&1&1\cr 0&7&-2&5\cr 8&0&-1&-1\cr -1&-2&6&1}\right] 第2列和第4列交換 \left[ \matrix{2&3&1&1\cr -1&-2&6&1\cr 8&0&-1&-1\cr 0&7&-2&5} \right]
------------------------
\displaystyle t=round(v1.v3/v1.v1)=round(14/15)=1不等於0
\left[ \matrix{2&3&1&1\cr -1&-2&6&1\cr 8&0&-1&-1\cr 0&7&-2&5}\right] 第3列-1*第1列 \left[ \matrix{2&3&1&1\cr -1&-2&6&1\cr 6&-3&-2&-2\cr 0&7&-2&5} \right]
\displaystyle t=round(v1.v4/v1.v1)=round(24/15)=2不等於0
\left[ \matrix{2&3&1&1\cr -1&-2&6&1\cr 6&-3&-2&-2\cr 0&7&-2&5}\right] 第4列-2*第1列 \left[ \matrix{2&3&1&1\cr -1&-2&6&1\cr 6&-3&-2&-2\cr -4&1&-4&3} \right]
第3列長度 \sqrt{53}比第4列長度 \sqrt{42}長
\left[ \matrix{2&3&1&1\cr -1&-2&6&1\cr 6&-3&-2&-2\cr -4&1&-4&3}\right] 第3列和第4列交換 \left[ \matrix{2&3&1&1\cr -1&-2&6&1\cr -4&1&-4&3\cr 6&-3&-2&-2} \right]
\displaystyle t=round(v3.v4/v3.v3)=round(-25/42)=-1不等於0
\left[ \matrix{2&3&1&1\cr -1&-2&6&1\cr -4&1&-4&3\cr 6&-3&-2&-2}\right] 第4列--1*第3列 \left[ \matrix{2&3&1&1\cr -1&-2&6&1\cr -4&1&-4&3\cr 2&-2&-6&1} \right]
------------------------
\displaystyle t=round(v2.v4/v2.v2)=round(-33/42)=-1不等於0
\left[ \matrix{2&3&1&1\cr -1&-2&6&1\cr -4&1&-4&3\cr 2&-2&-6&1}\right] 第4列--1*第2列 \left[ \matrix{2&3&1&1\cr -1&-2&6&1\cr -4&1&-4&3\cr 1&-4&0&2} \right]
第3列長度 \sqrt{42}比第4列長度 \sqrt{21}長
\left[ \matrix{2&3&1&1\cr -1&-2&6&1\cr -4&1&-4&3\cr 1&-4&0&2}\right] 第3列和第4列交換 \left[ \matrix{2&3&1&1\cr -1&-2&6&1\cr 1&-4&0&2\cr -4&1&-4&3} \right]
------------------------
\displaystyle t=round(v1.v3/v1.v1)=round(-8/15)=-1不等於0
\left[ \matrix{2&3&1&1\cr -1&-2&6&1\cr 1&-4&0&2\cr -4&1&-4&3}\right] 第3列--1*第1列 \left[ \matrix{2&3&1&1\cr -1&-2&6&1\cr 3&-1&1&3\cr -4&1&-4&3} \right]
第2列長度 \sqrt{42}比第3列長度2\sqrt{5}長
\left[ \matrix{2&3&1&1\cr -1&-2&6&1\cr 3&-1&1&3\cr -4&1&-4&3}\right] 第2列和第3列交換 \left[ \matrix{2&3&1&1\cr 3&-1&1&3\cr -1&-2&6&1\cr -4&1&-4&3} \right]
------------------------
化簡後 v=\left[ \matrix{2&3&1&1\cr 3&-1&1&3\cr -1&-2&6&1\cr -4&1&-4&3} \right] , \frac{vi.vj}{vi^2}=\left[ \matrix{ *&\frac{7}{15}&-\frac{1}{15}&-\frac{2}{5} \cr *&*&\frac{2}{5}&-\frac{2}{5} \cr *&*&*&-\frac{19}{42} \cr *&*&*&*} \right] ,向量長度 =\left[ \matrix{\sqrt{15} \cr 2\sqrt{5} \cr \sqrt{42} \cr \sqrt{42}} \right]
所有向量符合(1)|vi.vj/vi.vi|≦1/2，(2)∥vi∥≦∥vj∥，1≦i<j≦n兩個條件，程式結束
(%o5)　 \left[ \matrix{2&3&1&1\cr 3&-1&1&3\cr -1&-2&6&1\cr -4&1&-4&3} \right]

6-2.改進LLL方法－Jacobi方法

原本的Gauss/Lagrange方法	本論文的Gauss/Lagrange方法
1   do\{\; 2　if(v_1的長度>v_2的長度) 3　　 \{\; v_1和v_2兩個向量交換 \}\; 4　t=\lfloor\; v_1.v_2/v_1.v_1 \rceil\; 5  if(t!=0) 6    \{\; v_2=v_2−tv_1; \}\; 7　 \}\; 8  while(t!=0);	1  G=V^TV 2  U=I_2 3   if G(1,1)<G(2,2) 4    swap G(:,1) and G(:,2) 5    swap G(1,: ) and G(2,: ) 6    swap Z(:,1) and G(:,2) 7  end 8 9  while G(1,1)>G(2,2) 10   q=\lfloor\; G(1,2)/G(2,2) \rceil\; 11   G(:,2)=G(:,2)-q \times G(:,1) 12   G(2,: )=G(2,: )-q \times G(1,: ) 13   U(:,2)=U(:,2)-q \times U(:,1) 14 end 15 16 return U

註：說明論文所用的符號
1.G=V^TV代表gram矩陣，其中G(i,i)=v_i.v_i,G(i,j)=v_i.v_j，對角線元素為每個向量長度的平方，其他元素為兩向量的內積。
2.B代表lattice，為文章一致改用V表示lattice。
3.Z代表unimodular矩陣，但一般常用U表示unimodular矩陣，其中det(U)=+1或-1。

Jacobi方法以列計算

Jacobi方法以行計算

1  G=VV^T 2  U=I_n 3 4  while not all pairs (v_i,v_j) satisfy \Vert\; v_i \Vert\; \le \Vert\; v_j \Vert\; and \displaystyle |\; v_i v_j |\; \le \frac{\Vert\; v_i \Vert\;^2}{2} 5    for i=1 to n-1 6      for j=i+1 to n 7        q=G(i,j)/G(i,i) 8        if |\; q |\;>1/2 9          G(:,j)=G(:,j)-\lfloor\; q \rceil\; \times G(:,i) 10         G(j,: )=G(j,: )-\lfloor\; q \rceil\; \times G(i,: ) 11         U(j,: )=U(j,: )-\lfloor\; q \rceil\; \times U(i,: ) 12       end 13       if G(i,i)>G(j,j) 14         swap G(:,i) and G(:,j) 15         swap G(i,: ) and G(j,: ) 16         swap U(i,: ) and U(j,: ) 17       end 18     end 19   end 20 end 21 22 return UV

1  G=V^TV 2  U=I_n 3 4  while not all pairs (v_i,v_j) satisfy \Vert\; v_i \Vert\; \le \Vert\; v_j \Vert\; and \displaystyle |\; v_i v_j |\; \le \frac{\Vert\; v_i \Vert\;^2}{2} 5    for i=1 to n-1 6      for j=i+1 to n 7        q=G(i,j)/G(i,i) 8        if |\; q |\;>1/2 9          G(:,j)=G(:,j)-\lfloor\; q \rceil\; \times G(:,i) 10         G(j,: )=G(j,: )-\lfloor\; q \rceil\; \times G(i,: ) 11         U(:,j)=U(:,j)-\lfloor\; q \rceil\; \times U(:,i) 12       end 13       if G(i,i)>G(j,j) 14         swap G(:,i) and G(:,j) 15         swap G(i,: ) and G(j,: ) 16         swap U(:,i) and U(:,j) 17       end 18     end 19   end 20 end 21 22 return VU

註：
論文演算法以行向量運算，本文章卻以列向量運算，所以將兩個演算法並列出來，其中第1,11,16,22行有差異。
論文演算法第7行q=G(i,j)/G(j,j)會導致lattice尚未化簡完程式就提早結束，更正為q=G(i,j)/G(i,i)。

參考資料：
http://www.cas.mcmaster.ca/~qiao/
Filip Jeremic and Sanzheng Qiao. A Parallel Jacobi-type Lattice Basis Reduction Algorithm. International Journal of Numerical Analysis and Modeling, Series B. Vol. 5, No. 1-2, 2014, 1-12.
http://www.cas.mcmaster.ca/~qiao/publications/JQ14.pdf
Filip Jeremic and Sanzheng Qiao. A GPU Implementation of a Jacobi Method for Lattice Basis Reduction. Proceedings of International Workshop on Data-Intensive Scientific Discovery (DISD), 2013, Shanghai University, Shanghai, China, August 1-4, 2013.
http://www.cas.mcmaster.ca/~qiao/publications/JQ13.pdf



以列向量運算的Jacobi方法
(%i1)
JacobibyRow(V):=block
([End,G,n:length(V),U,i,j],
G:V.transpose(V),
U:ident(n),
do
(End:true,
for i:1 thru n-1 do
  (for j:i+1 thru n do
     (q:G[i,j]/G[i,i],
      if abs(q)>1/2 then
        (print("G(",i,",",j,")/G(",i,",",i,")=",q,">",1/2,"  ,  round(",q,")=",round(q)),
         print("G=",G,"  第",j,"行-",round(q),"*第",i,"行  G=",G:columnop(G,j,i,round(q))),
         print("G=",G,"  第",j,"列-",round(q),"*第",i,"列  G=",G:rowop(G,j,i,round(q))),
         print("U=",U,"  第",j,"列-",round(q),"*第",i,"列  U=",U:rowop(U,j,i,round(q))),
         End:false
        ),
      if G[i,i]>G[j,j] then
        (print("G(",i,",",i,")=",G[i,i],">G(",j,",",j,")=",G[j,j]),
         print("G=",G,"  第",i,"行和第",j,"行交換  G=",G:columnswap(G,i,j)),
         print("G=",G,"  第",i,"列和第",j,"列交換  G=",G:rowswap(G,i,j)),
         print("U=",U,"  第",i,"列和第",j,"列交換  U=",U:rowswap(U,i,j)),
         End:false
        )
     )
    ),
   if End=true then
     (print("G=",G,",",vi.vj/vi.vi,"=",genmatrix(lambda([i,j],if i<j then G[i,j]/G[i,i] else "*"),n,n),
             ",向量長度=",genmatrix(lambda([i,j],sqrt(G[i,i])),n,1)),
      print("所有向量符合(1)|vi.vj/vi.vi|≦1/2，(2)∥vi∥≦∥vj∥，1≦i<j≦n兩個條件，程式結束"),
      print("化簡後向量V=U.V=",U,V,"=",U.V),
      return(U.V)
     ),
    print("------------------------")
  )
)$

列向量的lattice
(%i2)
V:matrix([1,1,1],
             [-1,0,2],
             [3,5,6]);

(%o2)　 \left[ \matrix{1&1&1\cr -1&0&2\cr 3&5&6} \right]

(%i3)　JacobibyRow(V);
\displaystyle G( 1,3)/G(1,1)=\frac{14}{3}>\frac{1}{2}  ,   \displaystyle round(\frac{14}{3})=5
G=\left[ \matrix{3&1&14\cr 1&5&9\cr 14&9&70} \right] 第3行-5*第1行   G=\left[ \matrix{3&1&-1\cr 1&5&4\cr 14&9&0} \right]
G=\left[ \matrix{3&1&-1\cr 1&5&4\cr 14&9&0} \right] 第3列-5*第1列   G=\left[ \matrix{3&1&-1\cr 1&5&4\cr -1&4&5} \right]
U=\left[ \matrix{1&0&0\cr 0&1&0\cr 0&0&1} \right] 第3列-5*第1列   U=\left[ \matrix{1&0&0\cr 0&1&0\cr -5&0&1} \right]
\displaystyle G( 2,3)/G(2,2)=\frac{4}{5}>\frac{1}{2}  ,   \displaystyle round(\frac{4}{5})=1
G=\left[ \matrix{3&1&-1\cr 1&5&4\cr -1&4&5} \right] 第3行-1*第2行   G=\left[ \matrix{3&1&-2\cr 1&5&-1\cr -1&4&1} \right]
G=\left[ \matrix{3&1&-2\cr 1&5&-1\cr -1&4&1} \right] 第3列-1*第2列   G=\left[ \matrix{3&1&-2\cr 1&5&-1\cr -2&-1&2} \right]
U=\left[ \matrix{1&0&0\cr 0&1&0\cr -5&0&1} \right] 第3列-1*第2列   U=\left[ \matrix{1&0&0\cr 0&1&0\cr -5&-1&1} \right]
G( 2,2)=5>G(3,3)=2
G=\left[ \matrix{3&1&-2\cr 1&5&-1\cr -2&-1&2} \right] 第2行和第3行交換   G=\left[ \matrix{3&-2&1\cr 1&-1&5\cr -2&2&-1} \right]
G=\left[ \matrix{3&-2&1\cr 1&-1&5\cr -2&2&-1} \right] 第2列和第3列交換   G=\left[ \matrix{3&-2&1\cr -2&2&-1\cr 1&-1&5} \right]
U=\left[ \matrix{1&0&0\cr 0&1&0\cr -5&-1&1} \right] 第2列和第3列交換   U=\left[ \matrix{1&0&0\cr -5&-1&1\cr 0&1&0} \right]
------------------------
\displaystyle G( 1,2)/G(1,1)=-\frac{2}{3}>\frac{1}{2}  ,   \displaystyle round(\frac{-2}{3})=-1
G=\left[ \matrix{3&-2&1\cr -2&2&-1\cr 1&-1&5} \right] 第2行--1*第1行   G=\left[ \matrix{3&1&1\cr -2&0&-1\cr 1&0&5} \right]
G=\left[ \matrix{3&1&1\cr -2&0&-1\cr 1&0&5} \right] 第2列--1*第1列   G=\left[ \matrix{3&1&1\cr 1&1&0\cr 1&0&5} \right]
U=\left[ \matrix{1&0&0\cr -5&-1&1\cr 0&1&0} \right] 第2列--1*第1列   U=\left[ \matrix{1&0&0\cr -4&-1&1\cr 0&1&0} \right]
G( 1,1)=3>G(2,2)=1
G=\left[ \matrix{3&1&1\cr 1&1&0\cr 1&0&5} \right] 第1行和第2行交換   G=\left[ \matrix{1&3&1\cr 1&1&0\cr 0&1&5} \right]
G=\left[ \matrix{1&3&1\cr 1&1&0\cr 0&1&5} \right] 第1列和第2列交換   G=\left[ \matrix{1&1&0\cr 1&3&1\cr 0&1&5} \right]
U=\left[ \matrix{1&0&0\cr -4&-1&1\cr 0&1&0} \right] 第1列和第2列交換   U=\left[ \matrix{ -4&-1&1\cr 1&0&0\cr 0&1&0} \right]
------------------------
\displaystyle G( 1,2)/G(1,1)=1>\frac{1}{2}  ,   \displaystyle round(1)=1
G=\left[ \matrix{1&1&0\cr 1&3&1\cr 0&1&5} \right] 第2行-1*第1行   G=\left[ \matrix{1&0&0\cr 1&2&1\cr 0&1&5} \right]
G=\left[ \matrix{1&0&0\cr 1&2&1\cr 0&1&5} \right] 第2列-1*第1列   G=\left[ \matrix{1&0&0\cr 0&2&1\cr 0&1&5} \right]
U=\left[ \matrix{-4&-1&1\cr 1&0&0\cr 0&1&0} \right] 第2列-1*第1列   U=\left[ \matrix{-4&-1&1\cr 5&1&-1\cr 0&1&0} \right]
------------------------
G=\left[ \matrix{1&0&0\cr 0&2&1\cr 0&1&5} \right], \displaystyle \frac{vi.vj}{vi^2}=\left[\matrix{\displaystyle *&0&0\cr *&*&\frac{1}{2}\cr *&*&*}\right],向量長度=\left[ \matrix{1\cr \sqrt{2}\cr \sqrt{5}} \right]
所有向量符合(1)|vi.vj/vi.vi|≦1/2，(2)∥vi∥≦∥vj∥，1≦i<j≦n兩個條件，程式結束
化簡後向量 V=U.V=\left[ \matrix{-4&-1&1\cr 5&1&-1\cr 0&1&0} \right] \left[\matrix{1&1&1\cr -1&0&2\cr 3&5&6}\right]=\left[ \matrix{0&1&0\cr 1&0&1\cr -1&0&2} \right]
(%o3)　 \left[ \matrix{0&1&0\cr 1&0&1\cr -1&0&2} \right]



－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－



以行向量運算的Jacobi方法
(%i1)
JacobibyColumn(V):=block
([End,G,n:length(V),U,i,j],
G:transpose(V).V,
U:ident(n),
do
(End:true,
for i:1 thru n-1 do
  (for j:i+1 thru n do
     (q:G[i,j]/G[i,i],
      if abs(q)>1/2 then
        (print("G(",i,",",j,")/G(",i,",",i,")=",q,">",1/2,"  ,  round(",q,")=",round(q)),
         print("G=",G,"  第",j,"行-",round(q),"*第",i,"行  G=",G:columnop(G,j,i,round(q))),
         print("G=",G,"  第",j,"列-",round(q),"*第",i,"列  G=",G:rowop(G,j,i,round(q))),
         print("U=",U,"  第",j,"行-",round(q),"*第",i,"行  U=",U:columnop(U,j,i,round(q))),
         End:false
        ),
      if G[i,i]>G[j,j] then
        (print("G(",i,",",i,")=",G[i,i],">G(",j,",",j,")=",G[j,j]),
         print("G=",G,"  第",i,"行和第",j,"行交換  G=",G:columnswap(G,i,j)),
         print("G=",G,"  第",i,"列和第",j,"列交換  G=",G:rowswap(G,i,j)),
         print("U=",U,"  第",i,"行和第",j,"行交換  U=",U:columnswap(U,i,j)),
         End:false
        )
     )
    ),
   if End=true then
     (print("G=",G,",",vi.vj/vi.vi,"=",genmatrix(lambda([i,j],if i<j then G[i,j]/G[i,i] else "*"),n,n),
             ",向量長度=",genmatrix(lambda([i,j],sqrt(G[i,i])),n,1)),
      print("所有向量符合(1)|vi.vj/vi.vi|≦1/2，(2)∥vi∥≦∥vj∥，1≦i<j≦n兩個條件，程式結束"),
      print("化簡後向量V=V.U=",V,U,"=",V.U),
      return(V.U)
     ),
    print("------------------------")
  )
)$

行向量的lattice
(%i2)
V:matrix([1,-1,3],
             [1,0,5],
             [1,2,6]);
(%o2)　 \left[ \matrix{1&-1&3\cr 1&0&5\cr 1&2&6} \right]

(%i3)　JacobibyColumn(V);
\displaystyle G(1,3)/G(1,1)=\frac{14}{3}>\frac{1}{2}   ,   \displaystyle round(14/3)=5
G=\left[ \matrix{3&1&14\cr 1&5&9\cr 14&9&70} \right]   第3行-5*第1行   G=\left[ \matrix{3&1&-1\cr 1&5&4\cr 14&9&0} \right]
G=\left[ \matrix{3&1&-1\cr 1&5&4\cr 14&9&0} \right]   第3列-5*第1列   G=\left[ \matrix{3&1&-1\cr 1&5&4\cr -1&4&5} \right]
U=\left[ \matrix{1&0&0\cr 0&1&0\cr 0&0&1} \right]   第3行-5*第1行   U=\left[ \matrix{1&0&-5\cr 0&1&0\cr 0&0&1} \right]
\displaystyle G(2,3)/G(2,2)=\frac{4}{5}>\frac{1}{2}   ,   \displaystyle round(4/5)=1
G=\left[ \matrix{3&1&-1\cr 1&5&4\cr -1&4&5} \right]   第3行-1*第2行   G=\left[ \matrix{3&1&-2\cr 1&5&-1\cr -1&4&1} \right]
G=\left[ \matrix{3&1&-2\cr 1&5&-1\cr -1&4&1} \right]   第3列-1*第2列   G=\left[ \matrix{3&1&-2\cr 1&5&-1\cr -2&-1&2} \right]
U=\left[ \matrix{1&0&-5\cr 0&1&0\cr 0&0&1} \right]   第3行-1*第2行   U=\left[ \matrix{1&0&-5\cr 0&1&-1\cr 0&0&1} \right]
G(2,2)=5>G(3,3)=2
G=\left[ \matrix{3&1&-2\cr 1&5&-1\cr -2&-1&2} \right]   第2行和第3行交換   G=\left[ \matrix{3&-2&1\cr 1&-1&5\cr -2&2&-1} \right]
G=\left[ \matrix{3&-2&1\cr 1&-1&5\cr -2&2&-1} \right]   第2列和第3列交換   G=\left[ \matrix{3&-2&1\cr -2&2&-1\cr 1&-1&5} \right]
U=\left[ \matrix{1&0&-5\cr 0&1&-1\cr 0&0&1} \right]   第2行和第3行交換   U=\left[ \matrix{ 1&-5&0\cr 0&-1&1\cr 0&1&0} \right]
------------------------
\displaystyle G(1,2)/G(1,1)=-\frac{2}{3}>\frac{1}{2}   ,   \displaystyle round(-\frac{2}{3})=-1
G=\left[ \matrix{3&-2&1\cr -2&2&-1\cr 1&-1&5} \right]   第2行--1*第1行   G=\left[ \matrix{3&1&1\cr -2&0&-1\cr 1&0&5} \right]
G=\left[ \matrix{3&1&1\cr -2&0&-1\cr 1&0&5} \right]   第2列--1*第1列   G=\left[ \matrix{3&1&1\cr 1&1&0\cr 1&0&5} \right]
U=\left[ \matrix{1&-5&0\cr 0&-1&1\cr 0&1&0} \right]   第2行--1*第1行   U=\left[ \matrix{1&-4&0\cr 0&-1&1\cr 0&1&0} \right]
G(1,1)=3>G(2,2)=1
G=\left[ \matrix{3&1&1\cr 1&1&0\cr 1&0&5} \right]   第1行和第2行交換   G=\left[ \matrix{1&3&1\cr 1&1&0\cr 0&1&5} \right]
G=\left[ \matrix{1&3&1\cr 1&1&0\cr 0&1&5} \right]   第1列和第2列交換   G=\left[ \matrix{ 1&1&0\cr 1&3&1\cr 0&1&5} \right]
U=\left[ \matrix{1&-4&0\cr 0&-1&1\cr 0&1&0} \right]   第1行和第2行交換   U=\left[ \matrix{ -4&1&0\cr -1&0&1\cr 1&0&0} \right]
------------------------
\displaystyle G(1,2)/G(1,1)=1>\frac{1}{2}   ,   round(1)=1
G=\left[ \matrix{1&1&0\cr 1&3&1\cr 0&1&5} \right]   第2行-1*第1行   G=\left[ \matrix{1&0&0\cr 1&2&1\cr 0&1&5} \right]
G=\left[ \matrix{1&0&0\cr 1&2&1\cr 0&1&5} \right]   第2列-1*第1列   G=\left[ \matrix{1&0&0\cr 0&2&1\cr 0&1&5} \right]
U=\left[ \matrix{-4&1&0\cr -1&0&1\cr 1&0&0} \right]   第2行-1*第1行   U=\left[ \matrix{-4&5&0\cr -1&1&1\cr 1&-1&0} \right]
------------------------
G=\left[ \matrix{1&0&0\cr 0&2&1\cr 0&1&5}\right],\frac{vi.vj}{vi^2}=\left[ \matrix{*&0&0\cr *&*&\frac{1}{2}\cr *&*&*}\right] ,向量長度 =\left[ \matrix{1\cr \sqrt{2}\cr \sqrt{5}} \right]
所有向量符合(1)|vi.vj/vi.vi|≦1/2，(2)∥vi∥≦∥vj∥，1≦i<j≦n兩個條件，程式結束
化簡後向量 V=V.U=\left[ \matrix{1&-1&3\cr 1&0&5\cr 1&2&6} \right] \left[ \matrix{-4&5&0\cr -1&1&1\cr 1&-1&0} \right]=\left[ \matrix{0&1&-1\cr 1&0&0\cr 0&1&2} \right]
(%o3)　 \left[ \matrix{0&1&-1\cr 1&0&0\cr 0&1&2} \right]

[ 本帖最後由 bugmens 於 2017-12-31 13:54 編輯 ]