{-# LANGUAGE NoImplicitPrelude,RebindableSyntax #-}
module Algebra.Linear.Integral
  (
    divModUpperTriangular
  , divModLowerTriangular
  , hermite
  ) where


import           Algebra.Linear
import           Auxiliary     (minimumAmong,adorn)

import qualified Algebra.Ring
import           NumericPrelude

import           Control.Arrow (first,second,(***))
import           Data.Function (on)
import           Data.List     (partition)


divModUpperTriangular :: Vector Integer -> Matrix Integer -> (Vector Integer,Vector Integer)
divModUpperTriangular x u = reverse *** reverse $ divModLowerTriangular (reverse x) (reverse . map reverse $ u)

divModLowerTriangular :: Vector Integer -> Matrix Integer -> (Vector Integer,Vector Integer)
divModLowerTriangular x l = go l x [] [] where
  go []        []        q r = (q,r)
  go (lᵢ : l') (xᵢ : x') q r = go l' x' (q ++ [qᵢ]) (r ++ [rᵢ]) where
    (lFirst,lᵢᵢ : _) = splitAt (length q) lᵢ
    y = xᵢ - innerProduct q lFirst
    (qᵢ,rᵢ)
      | lᵢᵢ == 0  = error "Algebra.Linear.divModLowerTriangular: zero on diagonal"
      | otherwise = y `divMod` lᵢᵢ

-- Compute the Hermite normal form of the matrix,
-- together with the unimodular basis transformation matrix.
-- In particular, if (h,u) = hermite l, then u l = h, u is unimodular,
-- and h is upper triangular.
hermite :: Matrix Integer -> (Matrix Integer,Matrix Integer)
hermite []          = ([],[])
hermite xs@([] : _) = (xs,identity (length xs))
hermite vs          = case minimumAmong (compare `on` (abs . head . fst)) nonZero of
  Nothing                    -> first (map (0 :)) $ hermite (map tail vs)
  Just ((v,i),[])            -> case hermite (map (tail . fst) startZero) of
    (h,u) -> ((positive v :) . map (0 :) $ h,changeSign v `matrixProduct` rowSwap n (1,i) `matrixProduct` shift u)
  Just ((v@(v₀ : _),i),rest) -> let
    (reduced,translates) = unzip . flip map rest $ \ (x@(x₀ : _),j) -> let
      (c,_) = x₀ `divMod` v₀ in ((zipWith (\ vᵢ xᵢ -> xᵢ - c * vᵢ) v x,j),(j,c))
    (h,u) = hermite $ v : map fst reduced ++ map fst startZero
      in (h,
          u
            `matrixProduct` permute (i : map snd reduced ++ map snd startZero)
            `matrixProduct` affine n i (map (second negate) translates)
          )
 where
  (startZero,nonZero) = partition ((==) 0 . head . fst) . flip zip [1 ..] $ vs
  positive v = map (* signum (head v)) v
  changeSign v = matrixFromFunction n n f where
    f i j
      | i == j && i == 1 = signum (head v)
      | i == j           = 1
      | otherwise        = 0
  shift u = (1 : replicate (n - 1) 0) : map (0 :) u
  n = length vs


-- adjoint :: Matrix Integer -> Matrix Integer
-- adjoint m = strip . fst . hermite . adorn $ m where
--   n = length m
--   strip = map (drop n)

-- Alternative implementation of Hermite normal form
{-
hermite' :: Matrix Integer -> Matrix Integer
hermite' []          = []
hermite' xs@([] : _) = xs
hermite' vs          = case minimumAmong (compare `on` (abs . head)) nonZero of
  Nothing                -> map (0 :) $ hermite' (map tail vs)
  Just (v,[])            -> (positive v :) . map (0 :) $ hermite' (map tail startZero)
  Just (v@(v₀ : _),rest) -> hermite' $ v : reduced ++ startZero where
    reduced = map (\ x@(x₀ : _) -> let (c,_) = x₀ `divMod` v₀ in zipWith (\ vᵢ xᵢ -> xᵢ - c * vᵢ) v x) rest
 where
  (startZero,nonZero) = partition ((==) 0 . head) vs
  positive v = map (* signum (head v)) v
-}

-- Examples

la :: Matrix Integer
la = 
  [
    [ 6, 2,-4,-4, 2,-2,-2]
  , [ 2, 6,-4,-4, 2,-2,-2]                                                                                                                                                                
  , [-4,-4, 8, 4,-4, 0, 4]                                                                                                                                                                 
  , [-4,-4, 4, 8,-4, 4, 0]                                                                                                                                                                 
  , [ 2, 2,-4,-4, 6,-2,-2]                                                                                                                                                                
  , [-2,-2, 0, 4,-2, 6,-2]                                                                                                                                                                
  , [-2,-2, 4, 0,-2,-2, 6]
  ]

l :: Matrix Integer
l = 
  [
    [2,1,1,1,1,1,1]                                                                                                                                                                    
  , [1,2,1,1,1,1,1]                                                                                                                                                                    
  , [1,1,2,0,1,1,0]                                                                                                                                                                    
  , [1,1,0,2,1,0,1]                                                                                                                                                                    
  , [1,1,1,1,2,1,1]                                                                                                                                                                   
  , [1,1,1,0,1,2,1]                                                                                                                                                                    
  , [1,1,0,1,1,1,2]
  ]