module Math.Lattices.LLL (
lll,
lllDelta,
Basis(..)
) where
import Prelude hiding ((*>))
import Data.Array
import Data.Ratio
import Math.Algebra.LinearAlgebra hiding ((!))
import Math.Lattices.Internal
import Math.LinearAlgebra.GramSchmidt
type Basis = Array Int [Rational]
type GSO = Array (Int, Int) Rational
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