module Main where

import List
import Monad

main = run cgs >> putStrLn "\n" >> run mgs where
    run f = sequence $ liftM (run' f) ks
    ks = [1e-1,1e-5,1e-10]

    run' f k = mprint $ res where
        q = f [[1,k,0,0,0],[1,0,k,0,0],[1,0,0,k,0],[1,0,0,0,k]]
        res = (transpose q ?* q) ?- ident (length q)

    --print matrix row by row.
    mprint xs = sequence (liftM print $ transpose xs) >> putStrLn "\n"

--------------------------------------------------------------------------------------- vector utils
infix  7 .& -- inner product
infixl 5 .- -- vector subtraction
infixl 6 .* -- vector scaling

(.&) xs ys = foldl1 (+) $ zipWith (*) xs ys
(.*) xs y = map (*y) xs
(.-) = zipWith (-)

norm xs = sqrt $ xs .& xs
unit xs = xs .* 1/norm xs

----------------------------------------------------------------------------------------matrix utils
infixl 6 ?* -- matrix multiplication
infixl 5 ?- -- matrix subtraction

(?*) xs ys = map m [0..length(ys)-1] where
    xs' = transpose xs
    m i = map (yi .&) xs' where yi = ys !! i
(?-) = zipWith (.-)

ident n = [ map (delta i) [1..n] | i <- [1..n] ] where
    delta i j = if(i==j) then 1 else 0

----------------------------------------------------------classical gram schmidt... functional style
cgs as = map q [0..length(as)-1] where
    v i = foldl (.-) ai $ map (\j -> q j .* (ai .& q j)) [0..i-1]
        where ai = as !! i
    q i = unit $ v i

-----------------------------------------------------------modified gram schmidt... functional style
mgs as = map q [1..length(as)] where
    v i 1 = as !! (i-1)
    v j i = v j (i-1) .- q (i-1) .* (v j (i-1) .& q (i-1))
    q i = unit $ v i i