module Math.Lattices.LLL (
lll,
lllDelta,
closeVector,
Basis(..)
) where
import Data.Array
import Data.Ratio
import Math.Algebra.LinearAlgebra hiding ((!))
import Math.LinearAlgebra.GramSchmidt
type Basis = Array Int [Rational]
type GSO = Array (Int, Int) Rational
norm2 v = v <.> v
rnd x = floor $ x + 1%2
lll :: [[Rational]] -> Basis
lll basis = lllDelta basis $ 3%4
lllDelta :: [[Rational]] -> Rational -> Basis
lllDelta basis delta = lllLoop b' delta bb' mu_arr 1 n
where
n = length basis 1
(b, mu) = gramSchmidtOrthogonalization basis
bb = map norm2 b
b' = listArray (0, n) basis
bb' = listArray (0, n) bb
mu_arr = array ( ((0, 0), (n, n) ) ) [ ( (i,j), m ) | i <- [0..n],
j <- [0..n],
let m = (basis !! i <.> (b !! j)) / (norm2 $ b !! j) ]
sizeReduction :: Int -> Basis -> GSO -> (Basis, GSO)
sizeReduction k b mu = sizeReduction' indices k b mu
where
indices = reverse $ [0..k1]
sizeReduction' (l:ls) k b mu = sizeReduction' ls k b' mu'
where
(b', mu') = sizeReduction'' k l b mu
sizeReduction' [] _ b mu = (b, mu)
sizeReduction'' k l b mu = (b', mu'')
where
r = toRational $ round $ mu ! (k, l)
b_k = b ! k
b_l = b ! l
b_k' = b_k <-> (r *> b_l)
b' = b // [ (k, b_k') ]
mu' = mu // [ update | j <- [0..l1],
update <- [ ( (k, j), mu ! (k,j) (r * mu ! (l,j)) ) ] ]
mu'' = mu' // [ ( (k, l), mu' ! (k, l) r) ]
lovaszCondition :: Array Int Rational -> Int -> Rational -> GSO -> Bool
lovaszCondition bb k delta mu = (bb ! k) >= (delta m^2)*(bb ! (k1))
where
m = mu ! (k, k1)
swapBaseVectors b bb mu_ k n = (b', bb', mu'')
where
b' = b // [ (k 1, b ! k), (k, b ! (k1)) ]
m = mu_ ! (k, k1)
bb_k1 = bb ! (k1)
bb_k = bb ! k
btmp = bb_k + m^2*bb_k1
bb' = bb // [ (k, bb_k1*bb_k/btmp), (k1, btmp) ]
mu = mu_ // [ ( (k, k1), m*bb_k1/btmp ) ]
mu' = mu // [ update | j <- [0..k2],
update <- [ ( (k1, j), mu!(k,j) ), ( (k, j), mu!(k1, j)) ] ]
mu'' = mu' // [ update | i <- [k+1..n],
update <- [ ( (i, k1), update_i_k1 i), ( (i, k), update_i_k i) ] ]
where
update_i_k1 i = (mu' ! (k, k1)) * (mu' ! (i, k1)) + (mu' ! (i, k)) m*(mu' ! (i, k)) * (mu' ! (k, k1))
update_i_k i = (mu' ! (i, k1)) m * (mu' ! (i, k))
lllLoop :: Basis -> Rational -> Array Int Rational -> GSO -> Int -> Int -> Basis
lllLoop b delta bb mu k n | k > n = b
| isLovasz = lllLoop b' delta bb mu' (k+1) n
| otherwise = lllLoop b'' delta bb' mu'' nextk n
where
(b', mu') = sizeReduction k b mu
isLovasz = lovaszCondition bb k delta mu'
(b'', bb', mu'') = swapBaseVectors b' bb mu' k n
nextk = max 1 $ k 1
closeVector :: [[Rational]] -> [Ratio Integer] -> [Rational]
closeVector basis x = foldl1 (<+>) $ babaiNP (reverse $ [0..d]) basis' b' x
where
d = length basis 1
b' = listArray (0, d) $ gramSchmidtBasis basis
basis' = listArray (0, d) basis
projectTo v b = (v <.> b) / (norm2 b)
vsum zero = foldl (<+>) zero
babaiNP [] _ _ _ = []
babaiNP (i:is) b b' w = y_i : recurse
where
l_i = projectTo w $ b' ! i
delta = toRational $ rnd $ l_i
y_i = delta *> b ! i
w_i1 = w <-> (l_i delta) *> (b' ! i) <-> y_i
recurse = babaiNP is b b' w_i1