module LinAlg where

import ListNumSyntax
import Data.Number.IReal
import Data.Number.IReal.FoldB
import Data.Ratio
import Data.List
import Data.Bits

type Vector a = [a]
type Matrix a = [Vector a]

-- Ad hoc choice of precision restriction
prec0 x = prec 200 x

-- Multiplication  ------------------------------------------

mulM :: (VarPrec a, Num a) => Matrix a -> Matrix a -> Matrix a
mulM = map . apply 

apply :: (Num a, VarPrec a) => Matrix a -> Vector a -> Vector a
apply ass bs = map (`dot` bs) (transpose ass)

-- Multiplication with scalar --------------------------------------------------

smul :: Num a => a -> Vector a -> Vector a
smul x = map (x *)

smulM :: Num a => a -> Matrix a -> Matrix a
smulM x = map (smul x)

-- Scalar product --------------------------------------------------------------

dot :: (Num a, VarPrec a) => Vector a -> Vector a -> a
dot xs ys = prec0 (bsum (xs * ys))

-- Norms of vector ------------------------------------------------------------

norm2, norm1, normInf :: (Floating a, Ord a, VarPrec a) => Vector a -> a
norm2 xs = sqrt (dot xs xs)

norm1 = bsum . map abs 

normInf = maximum . map abs

-- L_1 and L_infinity norms of matrix -----------------------------------------
normM1, normMInf ::  (Ord a, Floating a, VarPrec a) => Matrix a -> a
normM1 = maximum . map norm1

normMInf = normM1 . transpose

-- Condition number of matrix (in L_1 norm) -----------------------------------

condM :: (Floating a, Ord a, VarPrec a) => Matrix a -> a
condM ass = normM1 ass * normM1 (inverse ass)

-- LQ factorization ------------------------------------------------------------

-- LQ factorization of quadratic matrix.
-- We choose LQ rather than QR since it fits better with
-- head/tail access pattern of lists. Returns L and Q, where
-- L is on triangular form, i.e. only non-zero elements present, and
-- Q is represented as a list of normal vectors of Householder reflectors.

lq :: (Floating a, Ord a, VarPrec a) => Matrix a -> (Matrix a, [Vector a])
lq ass@((x:xs):xss@(_:_)) = (bs : lss, vs : vss)
    where h:hs = map head ass
          d = norm2 (h:hs)
          vs = (if h >= 0 then h + d else h - d): hs
          v2 = 2*d*(d+abs h)  -- d1 = ||vs||^2
          -- Multiply ass by Householder reflector matrix 
          -- generated by vector vs, i.e. compute  A - 2*A*v*v^t
          bs : bss =  prec0 (zipWith f ass vs)
             where f as x = as - apply ass vs * repeat (2*x/v2)
          (lss,vss) = lq (map tail bss)
lq ass = (ass, [])

-- Solve quadratic, non-singular system Ax = b, using LQ factorization.
solve :: (Floating a, Ord a, VarPrec a) => Matrix a -> Vector a -> Vector a
solve ass = qIter vss . subst lss
   where (lss,vss) = lq ass
         
         -- Solve triangular system (matrix is lower triangular)
         subst :: Fractional a => [[a]] -> [a] -> [a]
         subst [] _ = []
         subst ((x:xs):xss) (b:bs) = y : subst xss (bs - xs * repeat y)
               where y = b/x

         -- Multiply by sequence of Householder matrices
         qIter :: (Fractional a, VarPrec a) => [Vector a] -> Vector a -> Vector a
         qIter [] bs = bs
         qIter (vs : vss) (b : bs) = hh vs (b : qIter vss bs)
         
         -- Vector rs multiplied by Householder reflector matrix generated by vs
         hh :: (Fractional a, VarPrec a) => Vector a -> Vector a -> Vector a
         hh vs rs = rs - vs * repeat x
            where x = 2 * dot rs vs/dot vs vs

inverse :: (Floating a, Ord a, VarPrec a) => Matrix a -> Matrix a
inverse ass = map (solve ass) bss
  where bss = map (\k -> replicate k 0 ++ 1 : replicate (n-1-k) 0) [0..n-1]  
        n = length ass      

-- Notorious ill-conditioned example: Hilbert matrix 
hilbert n = [[fromRational (1 % (i+j-1)) | i <- [1..n] ] | j <- [1..n] ]

vandermonde xs = take n (iterate (*xs) (replicate n 1))
  where n = length xs

 
{-
Examples:
> condM (hilbert 15) ?? 20
0.15391915629553122413e22
(0.82 secs, 283246752 bytes)

>  let ass = hilbert 25 :: Matrix IReal
> mapM_ (? 100) $ solve ass (map sqrt [1..25])
705.1533284963290956608362822423912734910227236293299614122851194081918876326343926875640023119369034821
-439673.9873388597705868990725370476651285605876080956807569044940939165590176737169191270235259441503432119
68204344.9731886365740080624537609996623689444810205308753901423565119434568104831725210699973202379088129882
-4664448582.3034346633655587322862893452286509473730850571024608318992483060505363965128396539508414113288932734
177397301669.4470302993210422259749351478977835134529306395373729499456573463841842107447990766834809423011040738
-4254091706828.1686190716920574829426672777980678548199026748623525425695787052836566266674697317136802722120565833
69545239200122.9144040744433369547417865224879339613940156554671219472545236203917004334667205391996866968931595183
-816847000729018.7737530264497363783287163634866491132650857885691558585008629694744116239295210128489892702899000614
7154341228591462.3012615296054988779018448134073649486823700047095232218238925733189512085896861380925718869840128164
-48009667586828894.5767332565013430439910972214304035550971986293103555764342546810272410148080254252305995120781866206
251844392677531274.0709127325911434838350769774189068254774415143641327544366238651541772916524464117349357235144334991
-1048143734335817624.0336204855464617426412722915751934481683873315855051574469784131062190560275498658327486102529792226
3498210324096944636.4524065506620210607308978049694297040490003525118886804565590254308130573820526159628678426741651309
-9431171401929169455.3424006600477669947185012388808800099904639824329048176939584968420634515763762886683395660091124001
20625644400616390081.3023145571266820730894216095277381736236614083175100283703764598631503021932601076440164701624585386
-36637408250801684992.2327515887414617772781813348542150962395077849459653041754367086153854944131054328104108083193973575
52765488024278928905.4679536245314145289284377844729508688894134162341828418605873919241149475390504553871953987733953659
-61295595789243212173.9536012643337759381753284186801798055604026908049933788464084009039536683226486562275743371985047603
56896779843800653456.9314422266376319337514520994749231451564802652135775892050259875530397079988163136075067115610230454
-41573863983516368926.7055317586515189490950549358051162588877518503235490221753772940736406787440880748304465115113504359
23365653654180951042.9851530663181890194238095089360096426483401824061668322486605839625160258109763064486458688933759486
-9740723523526602431.8575162294625153120767034736388134987019551900158949982443457241958170904379968964746678508149489593
2835293607475155284.1978846389721895478771840783053407794046273623991370837498430393230149403586055906293606569603644531
-514097728749287006.1850108351735689777337751827160575523829387979908417439866364790805580176949298837083538483210154166
43696874386454381.7499373427978656981157281081737498936218190682725431997586473877133365258605530660231535107025476866
(1.06 secs, 291352880 bytes)

-}