LLL Algorithm
[2001.08.27]LLLアルゴリズム
■実数をできるだけ高さの低い(分子や分母の小さい)有理数で近似するために、LLL(Lenstra, Lenstra and Lov\'asz)アルゴリズムを紹介する。
ここでは、R^n上の線形独立なn個のベクトルで張られたZ-加群を格子(Lattice)と呼ぶ。
つまり、b1,...,bnをR^nのある基底列ベクトルとするとき、
L = {Σi=1n xibi: xi \in Z}
は、b1,...,bnを基底とする格子である。
■The Gram-Schmidt-Process
R^nの内積< , >を使って、R^nの基底b1,...,bnから、直交基底b*1,...,b*nを求めることができる。
具体的には、各i=1,...,nに対して、
b*i = bi - Σj=1i-1 μi,jb*i
ただし、
μi,j = <b*i,bj>/<b*j,b*j>
とすれば良い。
■処理A
概要:指定されたμk,lについて、|μk,l|<= 1/2を満たすようにする。
入力:整数k,l; n次正方行列B,U
出力:n次正方行列B,U
処理: [x]は、xに最も近い整数とする。例えば、[(2m+1)/2]=mに固定する。
If |μk,l| > 1/2 then
r := [μk,l].
bk := bk - rbl.
For j=1 to l-1 do μk,j := μk,j - rμl,j.
μk,l := μk,l - r.
Endif
■処理B
概要:bk-1とbkを交換する。
入力:整数k; n次正方行列B,U; n次ベクトルA=(A1,...,An)
出力:n次正方行列B,U; n次ベクトルA
処理:
u := μk,k-1; v := Ak + u2Ak-1; μk,k-1 := uAk-1/v;
Ak := Ak-1Ak/v; Ak-1 := v;
ベクトルbk-1とbkを交換する。
For j=1 to k-2 do 要素μk-1,jとμk,jを交換する。
For i=k+1 to n do
t := μi,k-1,
μi,k-1 := μk,k-1μi,k-1 + μi,k - μi,kμk,k-1u.
μi,k := t - uμi,k.
Endif
■LLLアルゴリズム
概要:任意に与えられた基底Bから、簡約した(reduced)基底を求める。
入力:n次正方行列B; Bの列ベクトルb1,...,bnは基底を表す。
出力:n次正方行列B; Bの列ベクトルb1,...,bnは簡約した基底を表す。
処理:
Gram-Schmidt-Processにより、μk,k-1とb*iを計算する。
For i=1 to n do Ai := < b*i , b*i > .
k := 2.
Repeat
l=k-1に対して、処理Aを実行する。
If Ak < (3/4 - μi,k-12)Ak-1 then
処理Bを実行する。
If k > 2 then k := k - 1 .
Else
For l=k-2 to 1 do 処理Aを実行する。
k := k + 1 .
Endif
Until (k > n).
■LLLアルゴリズムの応用として、円周率πに近い４次以下の代数的数を求めてみよう。
つまり、|x1π4 + x2π3 + x3π2 + x4π + x5 |が小さくなる整数の組み(x1,x2,x3,x4,x5)を求めてみよう。
ここで、各整数xiは10^6程度とする。
B = ((1 0 0 0 [106π4])t
(0 1 0 0 [106π3])t
(0 0 1 0 [106π2])t
(0 0 0 1 [106π])t
(0 0 0 0 106)t)
とする。Bから簡約された基底B'は、
B' = ((2 5 0 -6 9)t
(6 -5 12 1 -551)t
(4 2 -5 12 14)t
(3 14 9 -1 -6)t
(16 5 -12 -1 0)t)
となる。B-1B'を計算すると、
B-1B' = ((2 5 0 -6 -331)t
(6 -5 12 1 -551)t
(4 2 -5 12 -440)t
(3 14 9 -1 -812)t
(16 5 -12 -1 -1692)t)
となる。
B-1B'の第１列(2 5 0 -6 -331)tに対応する４次方程式
2x4+5x3-6x-331 = 0
の実数解は、
β=-4.5280612244...
と
α=3.1415926291...
であり、
|α - π| <= 0.0000000245
となる。
[2021.09.26追記]
Pari/GPを使うと、以下のように、LLL-algorithmを簡単に実行できる。
[Pari/GPによる計算]
gp> A=[1,0,0,0,0;0,1,0,0,0;0,0,1,0,0;0,0,0,1,0;floor(Pi^4*10^6),floor(Pi^3*10^6),floor(Pi^2*10^6),floor(Pi*10^6),10^6]
%1 =
[ 1 0 0 0 0]
[ 0 1 0 0 0]
[ 0 0 1 0 0]
[ 0 0 0 1 0]
[97409091 31006276 9869604 3141592 1000000]
gp> B=qflll(A,1)
%2 =
[ -2 4 -6 -12 5]
[ -5 2 5 -3 19]
[ 0 -5 -12 7 9]
[ 6 12 -1 13 -7]
[331 -440 551 1152 -1143]
gp> polroots(-2*x^4-5*x^3+6*x+331)
time = 2 ms.
%3 = [-4.3508990352028642112674596118328069927 + 0.E-38*I, 3.1415926291137477851684864709196875781 + 0.E-38*I, -0.64534679695544178695051342954344029270 - 3.4192741655585576260837808960532422770*I, -0.64534679695544178695051342954344029270 + 3.4192741655585576260837808960532422770*I]~
gp> real(polroots(-2*x^4-5*x^3+6*x+331)[1])
time = 1 ms.
%4 = -4.3508990352028642112674596118328069927
gp> real(polroots(-2*x^4-5*x^3+6*x+331)[2])
time = 1 ms.
%5 = 3.1415926291137477851684864709196875781
[参考文献]
- [1]Nigel P. Smart, "The Algorithmic Resolution of Diophantine Equations", LMSST 41, Cambridge University Press, 1998, ISBN0-521-64633-2, p64-76.
- [2]Henri Cohen, "A Course in Computational Algebraic Number Theory", GTM 138, Springer-Verlag New York Inc., 1996, ISBN-387-55640-0, p79-105.
| Last Update: 2021.09.26 |
| H.Nakao |