|
|
10#
大 中
小 發表於 2015-5-20 22:56 只看該作者
3.改進LLL方法－modified LLL或MLLL
原本的LLL需要線性獨立的向量，若Lattice中有線性相依的向量，那在計算Gram Schmidt正交化時會發生錯誤，Pohst在1987年發表了modified
LLL演算法，就算有線性相依的向量也能處理。
maxima程式碼我參考Computation with finitely presented groups的第374頁的虛擬碼，只是我寫的maxima程式碼一直有錯誤，而我找不出錯誤在哪裡，所以我將我目前所寫的程式碼公布出來，看有沒有網友能幫我除錯。
115.6.19補充，利用claude修正程式碼
| maxima的go無法跳出while do迴圈 | maxima和c/c++的round執行結果不一致 |
虛擬碼的Goto 99雖然maxima也有對應的go指令
但實際使用時go無法從while do跳出來(出現do loop: 'go' not within 'block': 99錯誤訊息)
MLLL():=block(
x:1,
99,
while x<5 do
(x:x+1,
if x=3 then
(go(99))
else
(print("x=",x))
)
)$
MLLL();
\(x=2\)
do loop: 'go' not within 'block': 99
#0: MLLL()
-- an error. To debug this try: debugmode(true); |
maxima的\(round(-2.5)=-2\)
c/c++的\(round(-2.5)=-3\)
在第二個範例\(\left[\matrix{4&-1\cr5&4\cr-2&-4}\right]\)
以maxima的round執行結果\(\left[\matrix{-1&1\cr2&1\cr0&0}\right]\)
和書上執行結果不一致\(\left[\matrix{-1&1\cr1&2\cr0&0}\right]\) |
當時我使用startagain變數取代Goto 99，但程式沒寫好導致無法得到正確結果，改用claude修正程式碼如下。
Example 7.2. If \(n=1\) in MLLL, then MLLL performs the version of the Euclidean algorithm in which remainders under division are chosen to have minimum absolute value. For example, suppose that \(c_1=36,c_2=84,\) and \(c_3=100\). Table 8.7.1 summarizes the changes in the values of \(b_1,b_2,\) and \(b_3\) during the execution of MLLL\((c_1,c_2,c_3;b_1,b_2,b_3)\). The fourth column indicates the "places" in the procedure at which the changes occur.
Table 8.7.1
| \(b_1\) | \(b_2\) | \(b_3\) | Place |
\(\matrix{36\cr36\cr12\cr12\cr12\cr12\cr4\cr4}\)
|
\(\matrix{84\cr12\cr36\cr0\cr100\cr4\cr12\cr0}\)
|
\(\matrix{100\cr100\cr100\cr100\cr0\cr0\cr0\cr0}\)
|
\(\matrix{ \cr2\cr5\cr2\cr4\cr2\cr5\cr4}\)
|
Example 7.3. Now let us consider an example in \(\mathbb{Z}^2\). Let \(c_1=(4,-1),c_2=(5,4),\) and \(c_3=(-2,-4)\). Table 8.7.2 shows the changes in the values of \(b_1,b_2,\) and \(b_3\) in the execution of MLLL.
Table 8.7.2
| \(b_1\) | \(b_2\) | \(b_3\) | Place |
\(\matrix{(4, -1)\cr(4, -1)\cr(4, -1)\cr(4, -1)\cr(-1, 1)\cr(-1, 1)\cr(-1, 1)\cr(-1, 1)\cr(-1, 1)\cr(-1, 1)}\)
|
\(\matrix{(5, 4)\cr(1, 5)\cr(1, 5)\cr(-1, 1)\cr(4, -1)\cr(1, 2)\cr(1, 2)\cr(-1, 1)\cr(0, 0)\cr(1, 2)}\)
|
\(\matrix{(-2, -4)\cr(-2, -4)\cr(-1, 1)\cr(1, 5)\cr(1, 5)\cr(1, 5)\cr(-1, 1)\cr(1, 2)\cr(1, 2)\cr(0, 0)}\)
|
\(\matrix{ \cr2\cr2\cr5\cr5\cr2\cr2\cr5\cr2\cr4}\)
|
Example 7.4. Our last example is in \(\mathbb{Z}^3\). Suppose the input to MLLL consists of the rows of the matrix \(\begin{bmatrix} 48 & -124 & 292 \\ 171 & -142 & 141 \\ -291 & 254 & -277 \end{bmatrix}\).
Then the computation proceeds as shown in Table 8.7.3.
Table 8.7.3
| \(b_1\) | \(b_2\) | \(b_3\) | Place |
\(\matrix{(48, -124, 292)\cr(48, -124, 292)\cr(123, -18, -151)\cr(123, -18, -151)\cr(123, -18, -151)\cr(123, -18, -151)\cr(51, -30, 5)\cr(51, -30, 5)\cr(51, -30, 5)\cr(51, -30, 5)\cr(51, -30, 5)\cr(-12, 20, -40)\cr(-12, 20, -40)\cr(-12, 20, -40)\cr(-12, 20, -40)\cr(-12, 20, -40)\cr(18, -8, -6)\cr(18, -8, -6)\cr(18, -8, -6)\cr(18, -8, -6)\cr(18, -8, -6)\cr(18, -8, -6)\cr(18, -8, -6)}\)
|
\(\matrix{(171, -142, 141)\cr(123, -18, -151)\cr(48, -124, 292)\cr(171, -142, 141)\cr(171, -142, 141)\cr(51, -30, 5)\cr(123, -18, -151)\cr(21, 42, -161)\cr(21, 42, -161)\cr(192, -100, -20)\cr(-12, 20, -40)\cr(51, -30, 5)\cr(39, -10, -35)\cr(39, -10, -35)\cr(-18, 52, -126)\cr(18, -8, -6)\cr(-12, 20, -40)\cr(39, -10, -35)\cr(3, 6, -23)\cr(3, 6, -23)\cr(-18, 8, 6)\cr(0, 0, 0)\cr(3, 6, -23)}\)
|
\(\matrix{(-291, 254, -277)\cr(-291, 254, -277)\cr(-291, 254, -277)\cr(-291, 254, -277)\cr(51, -30, 5)\cr(171, -142, 141)\cr(171, -142, 141)\cr(171, -142, 141)\cr(192, -100, -20)\cr(21, 42, -161)\cr(21, 42, -161)\cr(21, 42, -161)\cr(21, 42, -161)\cr(-18, 52, -126)\cr(39, -10, -35)\cr(39, -10, -35)\cr(39, -10, -35)\cr(-12, 20, -40)\cr(-12, 20, -40)\cr(-18, 8, 6)\cr(3, 6, -23)\cr(3, 6, -23)\cr(0, 0, 0)}\)
|
\(\matrix{ \cr2\cr5\cr2\cr2\cr5\cr5\cr2\cr2\cr5\cr2\cr5\cr2\cr2\cr5\cr2\cr5\cr5\cr2\cr2\cr5\cr2\cr4}\)
|
不使用maxima內建指令round(-5/2)=-2，書上需要round(-5/2)=-3，所以另外寫ROUND指令
(%i1) ROUND(x):=if x>=0 then floor(x+1/2) else ceiling(x-1/2)$
(%i2)
ADJUST_MU(m,p):=block
(if abs(mu[m,p])>1/2 then
(r:ROUND(mu[m,p]),/*r:round(mu[m,p]),*/
b[m]:b[m]-r*b[p],
mu[m,p]:mu[m,p]-r,
for j:1 thru p-1 do
(mu[m,j]:mu[m,j]-r*mu[p,j])
)
)$
(%i3)
MLLL(c):=block
([bstar,s,h,k,mu,b,ZeroVector,H,i,t,m,nu,B,C,restart],
h:length(c[1]),
k:length(c),
mu:zeromatrix(k,k),
b:zeromatrix(k,k),
ZeroVector:create_list(0,i,1,h),
bstar:zeromatrix(k,h),
B:create_list(0,i,1,k),
b:c,
s:k,
i:1,
/* 修正3: 外層while不再用startagain=false限制，改用restart局部旗標 */
while i<=s do
(if b[ i ]=ZeroVector then
(if i<s then ([b[ i ],b[s]]:[b[s],b[ i ]], print(" P1",b)),
s:s-1
)
else /*b[ i ]不為零向量*/
(bstar[ i ]:b[ i ],
for j:1 thru i-1 do (mu[i,j]: (b[ i ].bstar[j])/B[j], bstar[ i ]:bstar[ i ]-mu[i,j]*bstar[j]),
B[ i ]:bstar[ i ].bstar[ i ],
if i=1 then
(i:2)
else /*i>1*/
(t:i, m:i,
restart:false, /* 修正3: 每次進入i>1區塊重設restart */
while m<=t and restart=false do /* 修正1+2: P4後立即退出內層while */
(ADJUST_MU(m,m-1), print(" P2",b),
nu:mu[m,m-1], C:B[m]+nu^2*B[m-1],
if C>=3/4*B[m-1] then
(for p:m-2 thru 1 step -1 do (ADJUST_MU(m,p), print("P 3",b)),
m:m+1
)
else
(if b[m]=ZeroVector then /* 修正1: P4與P5互斥 */
(if m<s then
([b[m],b[s]]:[b[s],b[m]], print("P 4",b)),
s:s-1, i:m,
restart:true /* 修正2: 設旗標讓內層while退出 */
)
else /*b[m]不為零向量，才執行Lovász交換(P5)*/
(if C#0 then
(mu[m,m-1]:nu*B[m-1]/C, B[m]:B[m-1]*B[m]/C,
for j:m+1 thru t do
(temp:matrix([1,mu[m,m-1]],[0,1]).matrix([0,1],[1,-nu]).matrix([mu[j,m-1]],[mu[j,m]]),
mu[j,m-1]:temp[1][1], mu[j,m]:temp[2][1]
)
),
B[m-1]:C,
[b[m-1],b[m]]:[b[m],b[m-1]], print("P 5",b),
if B[m-1]=0 then (t:m-1),
for j:1 thru m-2 do ([mu[m-1,j],mu[m,j]]:[mu[m,j],mu[m-1,j]]),
bstar[m-1]:b[m-1],
for j:1 thru m-2 do (bstar[m-1]:bstar[m-1]-mu[m-1,j]*bstar[j]),
if m<=t then
(bstar[m]:b[m],
for j:1 thru m-1 do (bstar[m]:bstar[m]-mu[m,j]*bstar[j])
),
if m>2 then (m:m-1)
)
)/*C<3/4B[m-1]*/
),/*While m<=t and restart=false*/
/* 修正3: restart=true時不遞增i，讓外層while從i=m重新開始(模擬Goto 99) */
if restart=false then (i:i+1)
)/* i>1 */
)/* b[ i ]#0 */
),
return(b)
)$
(%i5)
c:matrix([36],
[84],
[100])$
MLLL(c);
\(P2\left[\matrix{36\cr12\cr100}\right]\)
\(P5\left[\matrix{12\cr36\cr100}\right]\)
\(P2\left[\matrix{12\cr0\cr100}\right]\)
\(P4\left[\matrix{12\cr100\cr0}\right]\)
\(P2\left[\matrix{12\cr4\cr0}\right]\)
\(P5\left[\matrix{4\cr12\cr0}\right]\)
\(P2\left[\matrix{4\cr0\cr0}\right]\)
(%o5) \(\left[\matrix{4\cr0\cr0}\right]\)
(%i7)
c:matrix([4,-1],
[5,4],
[-2,-4])$
MLLL(c);
\(P2\) \(\left[\matrix{4&-1\cr 1&5\cr -2&-4}\right]\)
\(P2\) \(\left[\matrix{4&-1\cr 1&5\cr -1&1}\right]\)
\(P5\) \(\left[\matrix{4&-1\cr -1&1\cr 1&5}\right]\)
\(P2\) \(\left[\matrix{4&-1\cr -1&1\cr 1&5}\right]\)
\(P5\) \(\left[\matrix{-1&1\cr 4&-1\cr 1&5}\right]\)
\(P2\) \(\left[\matrix{-1&1\cr 1&2\cr 1&5}\right]\)
\(P2\) \(\left[\matrix{-1&1\cr 1&2\cr -1&1}\right]\)
\(P5\) \(\left[\matrix{-1&1\cr -1&1\cr 1&2}\right]\)
\(P2\) \(\left[\matrix{-1&1\cr 0&0\cr 1&2}\right]\)
\(P4\) \(\left[\matrix{-1&1\cr 1&2\cr 0&0}\right]\)
\(P2\) \(\left[\matrix{-1&1\cr 1&2\cr 0&0}\right]\)
(%o7) \(\left[\matrix{-1&1\cr 1&2\cr 0&0}\right]\)
(%i9)
c:matrix([48,-124,292],
[171,-142,141],
[-291,254,-277])$
MLLL(c);
\(P2 \begin{bmatrix} 48 & -124 & 292 \\ 123 & -18 & -151 \\ -291 & 254 & -277 \end{bmatrix}\)
\(P5 \begin{bmatrix} 123 & -18 & -151 \\ 48 & -124 & 292 \\ -291 & 254 & -277 \end{bmatrix}\)
\(P2 \begin{bmatrix} 123 & -18 & -151 \\ 171 & -142 & 141 \\ -291 & 254 & -277 \end{bmatrix}\)
\(P2 \begin{bmatrix} 123 & -18 & -151 \\ 171 & -142 & 141 \\ 51 & -30 & 5 \end{bmatrix}\)
\(P5 \begin{bmatrix} 123 & -18 & -151 \\ 51 & -30 & 5 \\ 171 & -142 & 141 \end{bmatrix}\)
\(P2 \begin{bmatrix} 123 & -18 & -151 \\ 51 & -30 & 5 \\ 171 & -142 & 141 \end{bmatrix}\)
\(P5 \begin{bmatrix} 51 & -30 & 5 \\ 123 & -18 & -151 \\ 171 & -142 & 141 \end{bmatrix}\)
\(P2 \begin{bmatrix} 51 & -30 & 5 \\ 21 & 42 & -161 \\ 171 & -142 & 141 \end{bmatrix}\)
\(P2 \begin{bmatrix} 51 & -30 & 5 \\ 21 & 42 & -161 \\ 192 & -100 & -20 \end{bmatrix}\)
\(P5 \begin{bmatrix} 51 & -30 & 5 \\ 192 & -100 & -20 \\ 21 & 42 & -161 \end{bmatrix}\)
\(P2 \begin{bmatrix} 51 & -30 & 5 \\ -12 & 20 & -40 \\ 21 & 42 & -161 \end{bmatrix}\)
\(P5 \begin{bmatrix} -12 & 20 & -40 \\ 51 & -30 & 5 \\ 21 & 42 & -161 \end{bmatrix}\)
\(P2 \begin{bmatrix} -12 & 20 & -40 \\ 39 & -10 & -35 \\ 21 & 42 & -161 \end{bmatrix}\)
\(P2 \begin{bmatrix} -12 & 20 & -40 \\ 39 & -10 & -35 \\ -18 & 52 & -126 \end{bmatrix}\)
\(P5 \begin{bmatrix} -12 & 20 & -40 \\ -18 & 52 & -126 \\ 39 & -10 & -35 \end{bmatrix}\)
\(P2 \begin{bmatrix} -12 & 20 & -40 \\ 18 & -8 & -6 \\ 39 & -10 & -35 \end{bmatrix}\)
\(P5 \begin{bmatrix} 18 & -8 & -6 \\ -12 & 20 & -40 \\ 39 & -10 & -35 \end{bmatrix}\)
\(P2 \begin{bmatrix} 18 & -8 & -6 \\ -12 & 20 & -40 \\ 39 & -10 & -35 \end{bmatrix}\)
\(P2 \begin{bmatrix} 18 & -8 & -6 \\ -12 & 20 & -40 \\ 39 & -10 & -35 \end{bmatrix}\)
\(P5 \begin{bmatrix} 18 & -8 & -6 \\ 39 & -10 & -35 \\ -12 & 20 & -40 \end{bmatrix}\)
\(P2 \begin{bmatrix} 18 & -8 & -6 \\ 3 & 6 & -23 \\ -12 & 20 & -40 \end{bmatrix}\)
\(P2 \begin{bmatrix} 18 & -8 & -6 \\ 3 & 6 & -23 \\ -18 & 8 & 6 \end{bmatrix}\)
\(P5 \begin{bmatrix} 18 & -8 & -6 \\ -18 & 8 & 6 \\ 3 & 6 & -23 \end{bmatrix}\)
\(P2 \begin{bmatrix} 18 & -8 & -6 \\ 0 & 0 & 0 \\ 3 & 6 & -23 \end{bmatrix}\)
\(P4 \begin{bmatrix} 18 & -8 & -6 \\ 3 & 6 & -23 \\ 0 & 0 & 0 \end{bmatrix}\)
\(P2 \begin{bmatrix} 18 & -8 & -6 \\ 3 & 6 & -23 \\ 0 & 0 & 0 \end{bmatrix}\)
(%o9) \(\left[\matrix{18&-8&-6\cr 3&6&-23\cr 0&0&0}\right]\)
－－－－－－－－－－－－－－－－－－－－
另外 http://www.numbertheory.org/calc/krm_calc.html有calc的數學軟體
在 http://www.numbertheory.org/calc/下載calc_win32.exe，輸入以下指令(紅色文字)
CALC
A NUMBER THEORY CALCULATOR
K.R.MATTHEWS, 16th January 2015
Type exit to quit, help for information:
> mlll()
Do you wish to use an existing matrix from a file? (Y/N)
Enter y or n : n
enter the matrix of integers (first row non-zero) :
Enter the number of rows and number of columns
3 2
Enter row 1:
4 -1
Enter row 2:
5 4
Enter row 3:
-2 -4
The matrix entered is
4 -1
5 4
-2 -4
VERBOSE? (Y/N)
Enter y or n : n
enter the parameters m1 and n1 (normally 1 and 1) : 1 1
L =
0 0 0
1 0 0
1 0 0
D[0] = 1, D[1] = 2, D[2] = 9, D[3] = 0,
The corresponding transformation matrix is
-1 1 1
-2 3 3
4 -6 -7
The corresponding reduced basis is
-1 1
1 2
>
程式輸出\( (-1,1),(1,2) \)結果是正確的但對除錯沒有太大的幫助。
參考資料
MLLL原始的論文
http://www.researchgate.net/prof ... 73c47f645000000.pdf
有MLLL演算法虛擬碼
http://www.numbertheory.org/pdfs/mlll.pdf
－－－－－－－－－－－－－－－－－－－－
115.7.15新增
在MLLL運算過程中，再加上一個轉換矩陣\(H\)，\(H\)是從單位矩陣\(I\)開始，隨著MLLL化簡,同步記錄了所有列運算。
操作1：大小約化(Size Reduction) | 操作2：兩列交換(Swap) | 操作3：遇到零向量 |
將某一列減去另一列的整數倍(\(b_i\leftarrow b_i-r\cdot b_j\))。
將\(H\)對應的某一列減去另一列的整數倍(\(H_i\leftarrow H_i-r\cdot H_j\)) | 當不滿足Lovász條件時，交換相鄰的兩列。
將\(H\)交換對應相鄰的兩列。 | 將零向量推至矩陣底部，\(H\)對應的列也推至矩陣底部。 |
當MLLL執行結束後所得到的化簡矩陣\(R\)和轉換矩陣\(H\)，符合\(HM=R\)且\(|\;det(H)|\;=1\)。
而這個轉換矩陣\(H\)可以用來解線性方程組\(Ax=b\)。關鍵在於將「線性方程組」的問題，轉換為「尋找一組向量之間的整數線性關係（Integer Linear Relation）」。
方法 | 範例 |
問題敘述 |
| 假設\(A\)是一個\(m\times n\)矩陣，\(b\)則是一個\(m \times 1\)的行向量，要找出\(n\times 1\)整數解\(x\)，滿足\(Ax=b\)。 | \(A=\left[\matrix{-8&5&7&-7&3&-7&4&9&-6\cr
1&-2&0&-10&-4&3&8&5&2\cr
-7&3&6&5&1&2&5&0&-6\cr
-9&-3&4&9&-2&6&1&-10&-9\cr
-2&1&-5&-4&3&7&-8&-8&-5\cr
-1&1&-8&4&-8&-1&-9&8&6}\right]\)，\(b=\left[\matrix{3\cr-1\cr-1\cr-7\cr9\cr8}\right]\)
求整數解\(x\)滿足\(Ax=b\)。 |
步驟1：將原問題轉換成求整數線性關係，計算\([A|b]^T\) |
\(A\)矩陣由\(n\)個長度為\(m\)的行向量\(A_1,A_2,\ldots,A_n\)組成，\(b\)則是一個\(m \times 1\)的行向量。
方程組\(Ax=b\)可以展開寫成行向量的線性組合形式\(x_1 A_1+x_2 A_2+\ldots+x_n A_n=b\)
移項後可得\(x_1 A_1+x_2 A_2+\ldots+x_n A_n-1 \cdot b = 0\)，
求解\(Ax=b\)其實等同於去尋找一組係數\((x_1,x_2,\ldots,x_n,-1)\)，使得這\(n+1\)個向量的線性組合結果為零向量。
將矩陣\(A\)和矩陣\(b\)結合在一起，新矩陣大小為\(m\times (n+1)\)，其中前\(n\)行是\(A\)的行向量，最後一行是\(b\)。
MLLL是以列運算進行Lattice化簡，所以要將新矩陣轉置，就是在對\(A_1,A_2,\ldots,A_n\)和\(b\)進行線性組合。 | \([A|b]^T=\left[\matrix{-8&1&-7&-9&-2&-1\cr
5&-2&3&-3&1&1\cr
7&0&6&4&-5&-8\cr
-7&-10&5&9&-4&4\cr
3&-4&1&-2&3&-8\cr
-7&3&2&6&7&-1\cr
4&8&5&1&-8&-9\cr
9&5&0&-10&-8&8\cr
-6&2&-6&-9&-5&6\cr
3&-1&-1&-7&9&8}\right]\) |
步驟2：將\([A|b]^T\)利用MLLL化簡，得到轉換矩陣\(H\) |
| MLLL執行完後產生化簡矩陣\(R\)和轉換矩陣\(H\)，滿足\(H\cdot [A|b]^T=R\)，其中\(H\)記錄了所有列變換步驟。 | \([R,H]\) :MLLL_H\(([A|b]^T)\)
\(R=\left[\matrix{0&0&-1&0&0&0\cr
1&0&0&0&0&0\cr
0&0&0&-1&0&0\cr
0&1&0&0&0&0\cr
0&0&0&0&0&1\cr
0&0&0&0&1&0\cr
0&0&0&0&0&0\cr
0&0&0&0&0&0\cr
0&0&0&0&0&0\cr
0&0&0&0&0&0}\right]\) |
\(H=\left[\matrix{-37818&-19676&-85649&14204&43543&-24334&48460&0&0&0\cr
-189557&-98622&-429305&71196&218253&-121971&242900&0&-1&0\cr
214507&111603&485811&-80567&-246980&138025&-274871&0&1&0\cr
-438519&-228151&-993149&164704&504904&-282166&561922&0&-2&0\cr
440143&228996&996827&-165314&-506774&283211&-564003&0&2&0\cr
622969&324116&1410888&-233982&-717277&400851&-798278&0&3&0\cr
3216146&1673283&7283868&-1207958&-3703021&2069439&-4121200&1&17&0\cr
-1489972&-775197&-3374459&559621&1715531&-958726&1909263&0&-7&0\cr
-31640&-16463&-71657&11884&36430&-20359&40544&0&0&1\cr
-86456&-44981&-195803&32472&99544&-55630&110785&0&0&0}\right]\) |
步驟3：求得整數解\(x\) |
如果\(Ax=b\)存在整數解，在化簡矩陣\(R\)中，必然會出現某一列全為0(因為\(x_1A_1+x_2A_2+\ldots+x_nA_n-b=0\))。
從\(R\)矩陣最後一列開始，符合第\(i\)列全為\(0\)，\(H\)矩陣第\(i\)列第\(n+1\)行為\(\pm1\)
(1)該值為\(+1\)，則整數解\(x\)為\(H\)矩陣第\(i\)行前\(n\)個數字加上負號。
(2)該值為\(-1\)，則整數解\(x\)為\(H\)矩陣第\(i\)行前\(n\)個數字。 | \(R\)矩陣第9行全為\(0\)，\(H\)矩陣第9行第\(10\)行為\(+1\)
\(H[9]=\left[\matrix{-31640&-16463&-71657&11884&36430&-20359&40544&0&0&1}\right]\)
整數解\(x\)為\(H\)矩陣第\(9\)行前\(9\)個數字加上負號
\(x=\left[\matrix{31640&16463&71657&-11884&-36430&20359&-40544&0&0}\right]\) |
參考資料
Computational Algebraic Number Theory
作者：Michael Pohst
https://books.google.com.tw/book ... 20reduction&f=false
因為程式關係計算出來的轉換矩陣\(H\)和整數解\(x\)和書上不同，但都是符合\(Ax=b\)的解。
\(H^t=\begin{pmatrix}
-37818 & -19676 & -85649 & 14204 & 43543 & -24334 & 48460 & 0 & 0 & 0 \\
-386767 & -201227 & -875938 & 145266 & 445318 & -248864 & 495604 & 1 & 0 & 0 \\
-1221509 & -635526 & -2766436 & 458787 & 1406428 & -785976 & 1565244 & 3 & 0 & 0 \\
-1219885 & -634681 & -2762758 & 458177 & 1404558 & -784931 & 1563163 & 3 & 0 & 0 \\
828011 & 430797 & 1875254 & -310993 & -953360 & 532781 & -1061015 & -2 & 0 & 0 \\
-803061 & -417816 & -1818748 & 301622 & 924633 & -516727 & 1029044 & 2 & 0 & 0 \\
-86456 & -44981 & -195803 & 32472 & 99544 & -55630 & 110785 & 0 & 0 & 0 \\
2816502 & 1465368 & 6378727 & -1057851 & -3242880 & 1812269 & -3609071 & -7 & 0 & 0 \\
9467074 & 4925523 & 21440740 & -3555742 & -10900253 & 6091559 & -12131128 & -23 & 1 & 0 \\
-11297648 & -5877935 & -25586565 & 4243288 & 13007950 & -7269435 & 14476828 & 28 & 0 & 1
\end{pmatrix}\)，\(x=\begin{pmatrix}11297648 \\5877935 \\25586565 \\-4243288 \\-13007950 \\7269435 \\-14476828 \\-28 \\0\end{pmatrix}\)
(%i1) ROUND(x):=if x>=0 then floor(x+1/2) else ceiling(x-1/2)$
(%i2)
ADJUST_MU(m,p):=block
(if abs(mu[m,p])>1/2 then
(/*r:round(mu[m,p]),*/
r:ROUND(mu[m,p]),
b[m]:b[m]-r*b[p],
H[m]:H[m]-r*H[p], /* H 同步更新 */
mu[m,p]:mu[m,p]-r,
for j:1 thru p-1 do
(mu[m,j]:mu[m,j]-r*mu[p,j])
)
)$
與MLLL完全相同的演算法，新增一個轉換矩陣H
(%i3)
_MLLL_H(c):=block
([bstar,s,h,k,mu,b,ZeroVector,H,i,t,m,nu,B,C,restart],
h:length(c[1]),
k:length(c),
mu:zeromatrix(k,k),
b:zeromatrix(k,k),
ZeroVector:create_list(0,i,1,h),
bstar:zeromatrix(k,h),
B:create_list(0,i,1,k),
b:c,
H:ident(k),/* H 初始化為單位矩陣 */
s:k,
i:1,
/* 修正3: 外層while不再用startagain=false限制，改用restart局部旗標 */
while i<=s do
(if b[ i ]=ZeroVector then
(if i<s then
([b[ i ],b[s]]:[b[s],b[ i ]],
[H[ i ],H[s]]:[H[s],H[ i ]] /* H 同步交換 */
/*,print(" P1",b)*/
),
s:s-1
)
else /*b[ i ]不為零向量*/
(bstar[ i ]:b[ i ],
for j:1 thru i-1 do (mu[i,j]: (b[ i ].bstar[j])/B[j], bstar[ i ]:bstar[ i ]-mu[i,j]*bstar[j]),
B[ i ]:bstar[ i ].bstar[ i ],
if i=1 then
(i:2)
else /*i>1*/
(t:i, m:i,
restart:false, /* 修正3: 每次進入i>1區塊重設restart */
while m<=t and restart=false do /* 修正1+2: P4後立即退出內層while */
(ADJUST_MU(m,m-1), /*print(" P2",b),*/
nu:mu[m,m-1], C:B[m]+nu^2*B[m-1],
if C>=3/4*B[m-1] then
(for p:m-2 thru 1 step -1 do (ADJUST_MU(m,p)/*,print("P 3",b)*/),
m:m+1
)
else
(if b[m]=ZeroVector then /* 修正1: P4與P5互斥 */
(if m<s then
([b[m],b[s]]:[b[s],b[m]],
[H[m],H[s]]:[H[s],H[m]] /* H 同步交換 */
/*print("P 4",b)*/
),
s:s-1, i:m,
restart:true /* 修正2: 設旗標讓內層while退出 */
)
else /*b[m]不為零向量，才執行Lovász交換(P5)*/
(if C#0 then
(mu[m,m-1]:nu*B[m-1]/C, B[m]:B[m-1]*B[m]/C,
for j:m+1 thru t do
(temp:matrix([1,mu[m,m-1]],[0,1]).matrix([0,1],[1,-nu]).matrix([mu[j,m-1]],[mu[j,m]]),
mu[j,m-1]:temp[1][1], mu[j,m]:temp[2][1]
)
),
B[m-1]:C,
[b[m-1],b[m]]:[b[m],b[m-1]], /*print("P 5",b),*/
[H[m-1],H[m]]:[H[m],H[m-1]], /* H 同步交換(P5也只是換位置) */
if B[m-1]=0 then (t:m-1),
for j:1 thru m-2 do ([mu[m-1,j],mu[m,j]]:[mu[m,j],mu[m-1,j]]),
bstar[m-1]:b[m-1],
for j:1 thru m-2 do (bstar[m-1]:bstar[m-1]-mu[m-1,j]*bstar[j]),
if m<=t then
(bstar[m]:b[m],
for j:1 thru m-1 do (bstar[m]:bstar[m]-mu[m,j]*bstar[j])
),
if m>2 then (m:m-1)
)
)/*C<3/4B[m-1]*/
),/*While m<=t and restart=false*/
/* 修正3: restart=true時不遞增i，讓外層while從i=m重新開始(模擬Goto 99) */
if restart=false then (i:i+1)
)/* i>1 */
)/* b[ i ]#0 */
),
return([b,H])
)$
A.x=b的A矩陣
(%i4)
A:matrix([-8,5,7,-7,3,-7,4,9,-6],
[1,-2,0,-10,-4,3,8,5,2],
[-7,3,6,5,1,2,5,0,-6],
[-9,-3,4,9,-2,6,1,-10,-9],
[-2,1,-5,-4,3,7,-8,-8,-5],
[-1,1,-8,4,-8,-1,-9,8,6]);
(%o4) \(\left[\matrix{-8&5&7&-7&3&-7&4&9&-6\cr
1&-2&0&-10&-4&3&8&5&2\cr
-7&3&6&5&1&2&5&0&-6\cr
-9&-3&4&9&-2&6&1&-10&-9\cr
-2&1&-5&-4&3&7&-8&-8&-5\cr
-1&1&-8&4&-8&-1&-9&8&6}\right]\)
A.x=b的b矩陣
(%i5) b:matrix([3],[-1],[-1],[-7],[9],[8]);
(%o5) \(\left[\matrix{3\cr-1\cr-1\cr-7\cr9\cr8}\right]\)
計算[A|b]^T
(%i6) M:transpose(addcol(A,b));
(%o6) \(\left[\matrix{-8&1&-7&-9&-2&-1\cr
5&-2&3&-3&1&1\cr
7&0&6&4&-5&-8\cr
-7&-10&5&9&-4&4\cr
3&-4&1&-2&3&-8\cr
-7&3&2&6&7&-1\cr
4&8&5&1&-8&-9\cr
9&5&0&-10&-8&8\cr
-6&2&-6&-9&-5&6\cr
3&-1&-1&-7&9&8}\right]\)
執行MLLL_H()副程式
(%i7) [R,H]:MLLL_H(M);
(%o7) \([\left[\matrix{0&0&-1&0&0&0\cr
1&0&0&0&0&0\cr
0&0&0&-1&0&0\cr
0&1&0&0&0&0\cr
0&0&0&0&0&1\cr
0&0&0&0&1&0\cr
0&0&0&0&0&0\cr
0&0&0&0&0&0\cr
0&0&0&0&0&0\cr
0&0&0&0&0&0}\right]\),\(\left[\matrix{-37818&-19676&-85649&14204&43543&-24334&48460&0&0&0\cr
-189557&-98622&-429305&71196&218253&-121971&242900&0&-1&0\cr
214507&111603&485811&-80567&-246980&138025&-274871&0&1&0\cr
-438519&-228151&-993149&164704&504904&-282166&561922&0&-2&0\cr
440143&228996&996827&-165314&-506774&283211&-564003&0&2&0\cr
622969&324116&1410888&-233982&-717277&400851&-798278&0&3&0\cr
3216146&1673283&7283868&-1207958&-3703021&2069439&-4121200&1&17&0\cr
-1489972&-775197&-3374459&559621&1715531&-958726&1909263&0&-7&0\cr
-31640&-16463&-71657&11884&36430&-20359&40544&0&0&1\cr
-86456&-44981&-195803&32472&99544&-55630&110785&0&0&0}\right] ]\)
找整數解x
(%i10)
n:length(R)$
ZeroVector:create_list(0,i,1,length(R[1]))$
for i:n thru 1 step -1 do
(if R[ i ]=ZeroVector and H[i,n]=1 then
(print("在R[",i,"]全為0，H[",i,",",n,"]=+1，整數解x=",x:-rest(H[ i ],-1))),/*x再加上負號*/
if R[ i ]=ZeroVector and H[i,n]=-1 then
(print("在R[",i,"]全為0，H",i,",",n,"]=-1，整數解x=",x:rest(H[ i ],-1)))
)$
在\(R[9]\)全為\(0\)，\(H[9,10]=+1\)，整數解\(x=[31640,16463,71657,-11884,-36430,20359,-40544,0,0]\)
驗證整數解x是否符合A.x=b
(%i11) is(A.x=b);
(%o11) true
－－－－－－－－－－－－－－－－－－－－
115.7.15新增
Example 8.1. Let B be Matrix 8.8.1. This matrix was obtained by forming a product UV, where U was an 18-by-14 matrix, V was a 14-by-18 matrix, and the entries of U and V were chosen randomly from the set \({-2,-1,0,1,2}\). The rank of B is 14.
Matrix 8.8.1
由\(18\times14\)和\(14\times 18\)矩陣相乘，可確定矩陣的秩為14。
Matrix 8.8.2
因為當年電腦運算速度較慢，要先將矩陣轉換成Hermite normal form後再執行MLLL加快速度，但Hermite normal form不在本文討論範圍。
Matrix 8.8.3
執行MLLL結果。
Matrix 8.8.1矩陣
(%i12)
c:matrix([13,-1,14,-6,-10,-11,-1,2,-1,-4,-9,-1,2,8,4,3,-6,-8],
[5,12,15,2,-6,-6,5,-4,-10,-10,-3,7,1,-7,9,4,5,-3],
[-10,-3,5,-2,-5,15,1,-11,-8,-2,-13,-9,9,1,5,5,-4,-16],
[-4,3,1,6,-11,11,3,-8,-8,6,-5,6,15,-2,-1,-1,-6,6],
[-2,-3,13,-3,2,-9,-2,-4,-11,-9,2,-4,-2,1,-7,-12,2,8],
[-4,-3,-10,9,10,5,-6,-2,4,-1,0,-5,-10,-2,3,-1,9,8],
[0,1,7,12,-9,-3,7,5,-4,-15,1,2,-2,-10,-5,2,2,3],
[8,-4,-18,-7,7,-1,-4,-2,21,6,10,3,-7,-3,6,13,8,-19],
[7,-3,3,-6,-2,-8,-6,-1,9,-4,0,-5,-12,7,11,-4,8,-10],
[12,10,-5,-4,21,-18,1,3,9,10,5,-9,-21,2,17,2,11,11],
[6,-12,-2,-5,11,1,-4,-11,-1,6,8,-1,-7,-9,7,-14,8,-4],
[5,9,3,-2,16,-13,7,5,0,1,11,-5,-16,-7,10,-4,9,6],
[2,2,-2,-3,6,0,8,-1,3,3,-5,-9,-7,-7,10,9,6,-12],
[0,2,-5,5,0,3,-7,-5,3,-8,4,17,1,-13,4,7,6,-13],
[14,-1,-5,2,2,0,-1,-1,7,1,2,1,-2,3,8,14,-1,-1],
[-6,-3,7,0,-19,4,-4,-3,-6,1,-13,3,13,9,-9,-10,-6,4],
[5,-10,1,1,7,4,-6,-1,-5,-6,12,4,2,3,-3,-2,-10,6],
[-2,9,-3,-18,-3,-12,6,-5,8,11,-3,-4,-4,5,4,0,10,-6]);
(%o12) \(\left[\matrix{13&-1&14&-6&-10&-11&-1&2&-1&-4&-9&-1&2&8&4&3&-6&-8\cr
5&12&15&2&-6&-6&5&-4&-10&-10&-3&7&1&-7&9&4&5&-3\cr
-10&-3&5&-2&-5&15&1&-11&-8&-2&-13&-9&9&1&5&5&-4&-16\cr
-4&3&1&6&-11&11&3&-8&-8&6&-5&6&15&-2&-1&-1&-6&6\cr
-2&-3&13&-3&2&-9&-2&-4&-11&-9&2&-4&-2&1&-7&-12&2&8\cr
-4&-3&-10&9&10&5&-6&-2&4&-1&0&-5&-10&-2&3&-1&9&8\cr
0&1&7&12&-9&-3&7&5&-4&-15&1&2&-2&-10&-5&2&2&3\cr
8&-4&-18&-7&7&-1&-4&-2&21&6&10&3&-7&-3&6&13&8&-19\cr
7&-3&3&-6&-2&-8&-6&-1&9&-4&0&-5&-12&7&11&-4&8&-10\cr
12&10&-5&-4&21&-18&1&3&9&10&5&-9&-21&2&17&2&11&11\cr
6&-12&-2&-5&11&1&-4&-11&-1&6&8&-1&-7&-9&7&-14&8&-4\cr
5&9&3&-2&16&-13&7&5&0&1&11&-5&-16&-7&10&-4&9&6\cr
2&2&-2&-3&6&0&8&-1&3&3&-5&-9&-7&-7&10&9&6&-12\cr
0&2&-5&5&0&3&-7&-5&3&-8&4&17&1&-13&4&7&6&-13\cr
14&-1&-5&2&2&0&-1&-1&7&1&2&1&-2&3&8&14&-1&-1\cr
-6&-3&7&0&-19&4&-4&-3&-6&1&-13&3&13&9&-9&-10&-6&4\cr
5&-10&1&1&7&4&-6&-1&-5&-6&12&4&2&3&-3&-2&-10&6\cr
-2&9&-3&-18&-3&-12&6&-5&8&11&-3&-4&-4&5&4&0&10&-6}\right]\)
經MLLL_H(c)化簡可立即得到Matrix 8.8.3結果，不需要再先化簡Hermite normal form再MLLL來加速了
(%i13) MLLL_H(c);
(%o13)
\([\left[\matrix{0&-1&0&0&-1&-2&0&1&1&0&0&1&0&-2&-1&-1&-1&-2\cr
0&1&1&0&-3&-2&0&-1&1&1&-2&-1&1&1&1&0&-1&-1\cr
-2&1&0&0&2&1&2&0&-2&2&-1&-2&0&-1&0&-1&0&2\cr
-2&2&-1&1&-2&1&-1&0&1&1&-1&1&1&1&-1&-1&2&1\cr
-1&-2&-3&2&0&0&-1&0&2&-1&-1&-1&-1&0&-1&2&-1&1\cr
0&-2&-1&0&0&2&0&1&0&-1&-1&0&-2&1&3&0&-2&0\cr
-1&-2&0&2&-1&2&0&1&-1&-2&-1&-1&0&1&-1&0&-2&2\cr
-1&3&0&0&-1&0&1&2&0&0&0&3&-1&-1&2&-1&1&0\cr
-2&1&1&0&0&1&0&-2&-2&-2&0&0&1&0&0&1&0&0\cr
1&-2&-3&2&0&0&-1&-1&2&2&1&1&1&-2&-2&0&0&1\cr
0&1&0&2&1&-2&-2&0&1&1&0&-1&0&1&2&0&-3&1\cr
1&-1&-2&1&2&2&-2&0&1&0&-2&0&0&0&-1&3&2&-1\cr
-1&1&0&1&0&-1&-1&-1&-1&-1&0&2&-1&-1&0&-2&2&3\cr
0&0&0&-2&-2&-1&1&0&-1&1&-1&2&2&1&-2&0&-2&2\cr
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\cr
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\cr
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\cr
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0}\right],\)
\(\left[\matrix{645790&-2095817&-1739381&2461627&2183923&810421&-1832214&-849278&1462227&-2347825&-1596116&2374546&3067484&1567022&1317684&0&0&0\cr
-5349&17360&14408&-20390&-18090&-6713&15177&7035&-12112&19448&13221&-19669&-25409&-12980&-10915&0&0&0\cr
8612409&-27950310&-23196799&32828842&29125316&10807972&-24434845&-11326173&19500606&-31311157&-21286183&31667513&40908695&20898175&17572950&0&0&1\cr
184332&-598221&-496482&702637&623370&231324&-522980&-242414&417372&-670154&-455589&677781&875570&447284&376114&0&0&0\cr
1312928&-4260918&-3536263&5004631&4440043&1647633&-3724998&-1726632&2972793&-4773265&-3244997&4827590&6236374&3185847&2678929&0&0&0\cr
1134188&-3680846&-3054844&4323311&3835585&1423327&-3217885&-1491573&2568084&-4123443&-2803230&4170372&5387368&2752133&2314226&0&0&0\cr
1422220&-4615609&-3830632&5421231&4809645&1784787&-4035078&-1870362&3220257&-5170606&-3515120&5229453&6755508&3451046&2901931&0&0&0\cr
2665909&-8651820&-7180405&10161935&9015535&3345530&-7563633&-3505937&6036275&-9692146&-6588987&9802453&12662997&6468882&5439582&0&0&0\cr
931364&-3022607&-2508552&3550182&3149675&1168797&-2642437&-1224837&2108838&-3386056&-2301934&3424593&4423955&2259974&1900377&0&0&0\cr
28150667&-91358862&-75821457&107304918&95199509&35327122&-79868154&-37020922&63740017&-102344184&-69576384&103508974&133714865&68308133&57439244&0&1&2\cr
-2150276&6978408&5791591&-8196441&-7271775&-2698447&6100696&2827828&-4868755&7817518&5314563&-7906489&-10213755&-5217688&-4387473&0&0&0\cr
-26812140&87014872&72216249&-102202714&-90672901&-33647364&76070531&35260628&-60709266&97477856&66268121&-98587262&-127356904&-65060174&-54708086&0&-1&-2\cr
27395592&-88908379&-73787728&104426720&92646010&34379556&-77725882&-36027925&62030344&-99599046&-67710164&100732593&130128284&66475931&55898574&0&1&2\cr
28961421&-93990050&-78005159&110395361&97941310&36344562&-82168403&-38087147&65575767&-105291755&-71580224&106490091&137565931&70275447&59093528&0&1&2\cr
-55387255&179751216&149180921&-211125545&-187307801&-69507139&157142915&72839741&-125410337&201365157&136893558&-203656914&-263087884&-134398237&-113013382&1&-2&-4\cr
-57950352&188069374&156084409&-220895579&-195975647&-72723647&164414852&76210469&-131213819&210683520&143228437&-213081331&-275262525&-140617643&-118243184&0&-2&-4\cr
8097896&-26280533&-21811001&30867617&27385343&10162294&-22975085&-10649537&18335622&-29440600&-20014527&29775667&38464772&19649699&16523126&0&0&1\cr
4332703&-14061155&-11669774&16515432&14652274&5437240&-12292606&-5697936&9810305&-15751919&-10708587&15931193&20580220&10513389&8840545&0&0&0}\right] ]\)
|
