Gauss par tête et queue

Voici une des 257 implémentations de l'algorithme de Gauss... Petite originalité : on utilise des listes dynamiques et on ne les parcourt que par «tête et queue».

C'est assez efficace...Quelques microsecondes pour le calcul du déterminant d'une matrice de taille 1000.

Haskell
det matrice taille 1000 : 592 microsec
 
warming up
estimating clock resolution...
mean is 1.373373 us (640001 iterations)
found 3269 outliers among 639999 samples (0.5%)
  2863 (0.4%) high severe
estimating cost of a clock call...
mean is 33.83300 ns (9 iterations)
found 1 outliers among 9 samples (11.1%)
  1 (11.1%) high severe
 
benchmarking detMat/1000
collecting 100 samples, 1 iterations each, in estimated 5.147409 s
mean: 25.38584 ms, lb 25.30688 ms, ub 25.58190 ms, ci 0.950
std dev: 592.3872 us, lb 196.0864 us, ub 1.124442 ms, ci 0.950
found 8 outliers among 100 samples (8.0%)
  5 (5.0%) high mild
  3 (3.0%) high severe
variance introduced by outliers: 17.070%
variance is moderately inflated by outliers

Voici le code avec des noms les plus explicites possibles:

Haskell
-- déterminant
 
rotatePivotDet :: (Eq a, Num a) => Matrice a -> Matrice a
rotatePivotDet (Mat (ligne:lignes))
  | ((head ligne) /= 0) || (all (== 0) (head (List.transpose (ligne:lignes)))) = 
    Mat (ligne:lignes)
  | otherwise       = 
        rotatePivotDet (Mat (lignes ++  [map (* ((-1) ^ (length lignes))) ligne]))
 
determinant :: (Eq a, Fractional a) => Matrice a -> a
determinant (Mat []) = error "Matrice vide"
determinant (Mat [[a]]) = a
determinant m = 
  if pivot == 0 then 
    0
  else
    pivot * (determinant (Mat lignes'))
  where
    Mat ((pivot:l):lignes) = rotatePivotDet m    
    lignes' = map transVec lignes
    transVec (t:q)
      | t == 0     = q
      | otherwise  = fmap (/ pivot) (tail $ zipWith (-) (map (* pivot) (t:q)) (map (* t) (pivot:l)))
 
 
-- inverse
 
rotatePivot :: (Eq a, Num a) => [[a]] -> [[ a]]
rotatePivot (ligne:lignes)
  | ((head ligne) /= 0) || (all (== 0) (head (List.transpose (ligne:lignes)))) = 
    ligne:lignes
  | otherwise       = 
       rotatePivot  (lignes ++ [ligne])
                      
                      
                      
triangle :: (Eq a, Num a) => [[a]] -> [[a]]
triangle [] =  []
triangle m  = (pivot:l):(triangle lignes')
  where
    (pivot:l):lignes = rotatePivot m    
    lignes' = map transVec lignes
    transVec (t:q)
      | t == 0     = q
      | otherwise  = tail $ zipWith (-) (map (* pivot) (t:q)) (map (* t) (pivot:l))
 
 
completeLigne :: (Num a) =>  Int -> [a] -> [a]
completeLigne t l = [0 | x <- [0..(t - (length l) - 1)]] ++ l 
 
complete :: (Num a) =>  [[a]] -> [[a]]
complete m = map (completeLigne (2 * (length m))) m 
 
retourne_complet :: [[a]] -> [[a]]
retourne_complet m = reverse (map reverse m)
 
unitet :: (Num a) => Int ->  [[a]]
unitet n = [[if i == j then 1 else 0 | j <- [0..n-1]] | i <- [0..n-1]]
 
 
collet :: (Num a) => [[a]] -> [[a]] -> [[a]]
collet m n  =  zipWith (++) m n
 
tab :: (Num a) => [[a]] -> [[a]]
tab m = collet m (unitet (length m))
         
carre_g :: [[a]] -> [[a]]
carre_g  m = map (take long) m 
  where long = length m
        
carre_d :: [[a]] -> [[a]]
carre_d m = map (drop long) m 
  where long = length m
 
nm :: (Num a, Eq a) => [[a]] -> [[a]]
nm a = collet (retourne_complet m1) (reverse m2) 
  where (m1,m2) = (carre_g (m a),carre_d (m a))
          where m a = complete $ triangle $ tab a
                
nnm :: (Fractional a, Eq a) => [[a]] -> [[a]]
nnm a = map (\ligne -> (map (\x -> x / (head ligne))) ligne) (triangle (nm a))
 
inverse :: (Fractional a, Eq a) => Matrice a -> Matrice a
inverse (Mat a) = 
  if (determinant (Mat a) == 0) then  
    error "Matrice non inversible"
  else
    Mat (reverse $ carre_d $ complete (nnm a))

courtesy of webmatter.de