Calcul matriciel

À l'horizon de ce module : une structure en CAML autour de la méthode du pivot de Gauss (ou algorithme 方程 ) tout en fonctionnel...Pas une seule boucle pour calculer des l-réduites échelonnées, des inverses de matrices, des résolutions de systèmes dans des matrices à coefficients dans des anneaux quelconques, tout en récursion terminale. On terminera par un travail sur les matrices à coefficients dans Z/nZ avec comme fil rouge le chiffrement de Hill. En amphi, quelques exercices sur les bases, les changements de bases.

Le poly et ses sources TEX ainsi que le diaporama et ses sources TEX (version 6 nov 2012). Le diaporama d'algèbre linéaire et ses sources TEX (version 18 déc 2012). Les exercices du TD d'algèbre linéaire.

Calcul matriciel avec des listes et des fonctions récursives terminales

Voici le résultat de nos trois semaines de TP : il reste quelques petites fonctions à créer pendant les vacances...

OCaml
include List ;;
open List ;;
 
 
type 'a anneau =
  {zero      : 'a;
   un        : 'a;
   som       : 'a -> 'a -> 'a;
   prod      : 'a -> 'a -> 'a;
   sous      : 'a -> 'a -> 'a;
   div       : 'a -> 'a -> 'a;
   to_string : 'a -> string;
   egal      : 'a -> 'a -> bool;
   ordre     : 'a -> 'a -> bool;
  };;
 
 
exception Non_simplifiable;;
 
let entier =
  {zero      = 0;
   un        = 1;
   som       = ( + );
   prod      = ( * );
   sous      = ( - );
   div       = (fun x y -> if (x mod y = 0) 
                   then x/y else raise Non_simplifiable);
   to_string = string_of_int;
   egal      = (=);
   ordre     = (>=);
  };;
 
 
let reel =
  {zero      = 0.;
   un        = 1.;
   som       = ( +. );
   prod      = ( *. );
   sous      = ( -. );
   to_string = string_of_float;
   div       = ( /. );
   egal      = (=);
   ordre     = (>=);
  };;
 
let xor = fun (x:bool) y ->
  match x,y with
  |x,y when x = y -> false
  |_ -> true;;
 
let booleen =
  {zero      = false;
   un        = true;
   som       = ( xor );
   prod      = ( && );
   sous      = ( xor );
   div       = (fun x y -> if y then x else raise Non_simplifiable);
   egal      = (=);
   ordre     = (fun x y ->  if x then x else y);
   to_string = (fun x -> if x then "vrai" else "faux");
  };;
 
 
 
(*              joli affichage            *)
 
let rec printv a v =
match v with
  |[] -> print_string " | ";
  |t::q ->
    Printf.printf "%7s" (" "^(a.to_string t)^" ");
    (printv a q);;
 
let rec printm a m = 
match m with 
  |[] -> print_newline();
  |t::q ->
    print_string " |";
    printv a t;
    print_newline();
    printm a q;
;; 
 
 
(* pour créer un vecteur ligne de longueur longueur rempli de motif
    [x;x;x] = x::[x;x] *)
 
(*version non terminale *)
let rec make_vec0 = fun longueur motif ->
  match longueur with
  |0 -> []
  |_ -> motif::(make_vec0 (longueur - 1) motif);;
 
 
(*version terminale*)
let make_vec = fun longueur motif ->
  let rec aux = fun long acc ->
      match long with
      |0 -> acc
      |_ -> aux (long-1) (motif::acc)
  in aux longueur [];;
 
 
(* pour créer une matrice de taille row*col remplie de x *)
let make_mat row col x =
  make_vec row (make_vec col x);;
 
 
 
(* matrice remplie de 1  *)
let attila a (r,c) =
  make_mat r c a.un;; 
 
(* matrice remplie de neutres *)
let zeros a (r,c) =
  make_mat r c a.zero;; 
 
 
 
(* Variations autour de la longueur d'une liste *)
 
(* long t::liste = 1 + long liste *)
let rec longueur0 = fun liste ->
  match liste with
  |[] -> 0
  |t::q -> 1 + longueur0 q;;
 
(* la même en rec terminale *)
let longueur1 = fun liste ->
  let rec long_aux = fun l acc ->
    match l with
    |[] -> acc
    |t::q -> long_aux q (1 + acc)
  in long_aux liste 0;;
 
(* la même en utilisant le pliage gauche *)
let compte = fun somme_partielle y ->
  1+somme_partielle;;
 
let longueur liste = fold_left compte 0 liste;;
 
(* sinon, on utilise le List.length de caml ;-) *)
 
 
 
(* variations autour de nth *)
 
let rec neme = fun liste no ->
  match no with
  |0 -> hd liste
  |_ -> neme (tl liste) (no -1);;
 
(* sinon, on utilise nth de caml...*)
 
  
(* nb de lignes *)
let n_rows mat =
  length mat;;
 
 
(* nb de colonnes *)
let n_cols = fun mat ->
  length (hd mat);;
 
(* extraction de la ligne no i *)
let extract_line = fun mat i ->
  nth mat i;;
 
(* VARIATIONS AUTOUR DE L'EXTRACTION DE LA COLONNE NO j*)
(* explore le peigne et on extrait le j-eme
   de chaque tête *)
let rec extract_col1 = fun mat j ->
  match mat with
  |[]          -> []
  |tete::queue -> (nth tete j)::(extract_col1 queue j);; 
 
(* la même en version terminale *)
let extract_col2 = fun mat j ->
  let rec aux = fun m acc ->
    match m with
    |[]   -> rev acc
    |t::q -> aux q ((nth t j)::acc)
  in aux mat [];;
 
(* Pour le fun on peut  utiliser fold_left : on colle
   l'élément j d'une ligne aux autres par pliage
   On commence par coller le j-eme d'une liste
   à une liste initiale *)
let jth_col = fun j liste_ini liste_parcourue  ->
  (nth liste_parcourue j)::liste_ini;;
 
let extract_col = fun mat j ->
  rev (fold_left (jth_col j) [] mat);;
 
(* extraction de la 1ère colonne par pliage : on fabrique une fonction qui
   extrait les têtes de listes et les accumule*)
let extrait_tete = fun  accumulateur_de_tetes ligne_source  ->
  (hd ligne_source)::accumulateur_de_tetes;;
 
(* on extrait les têtes par pliage *)
let tete_col mat =
  rev (fold_left extrait_tete [] mat);;
 
(* queue_col : la même chose que extract_col mais
   je prends la queue au lieu de la tête avec pattern-matching (mais comme
   dirait M. Faucou, les vrai(e)s informaticien(ne)s font du fold_left) *)
let queue_col2 = fun mat ->
  let rec aux = fun m acc ->
    match m with
    |[]   -> rev acc
    |t::q -> aux q ((tl t)::acc)
  in aux mat [];;
 
(* Version fold_left : étape 1 -> Extraction de la matrice privée de sa 1ère colonne
    avec le  principe du pliage *)
let extrait_queue = fun accumulateur_de_queues ligne_source ->
  (tl ligne_source)::accumulateur_de_queues;;
(* Étape 2 -> On extrait les queues par pliage *)
let queue_col mat =
  rev (fold_left extrait_queue [] mat);;
 
 
(* Option tout en un : on met hd ou tl en option *)
let extrait_extremite = fun extremite accumulateur_d_extremites ligne_source ->
  (extremite ligne_source)::accumulateur_d_extremites;;
(* Étape 2 -> On extrait les queues par pliage *)
let extrem_col = fun extrem mat -> 
  rev (fold_left (extrait_extremite extrem) [] mat);;
 
 
(* Transposée d'une matrice : on extrait la 1ère colonne et on la colle
   à la fonction appliquée aux colonnes restantes.
   Attention ! Le cas terminal est une liste remplie
   de lignes vides != liste vide *)
(* Version 1 : non terminale *)
let rec transpose_n_t = fun mat ->
  match mat with
  |[]::_ -> [] 
  |_     -> tete_col(mat)::(transpose_n_t (queue_col mat));;
  
(* Version 2 : terminale *)
let transpose = fun mat ->
  let rec transp_aux = fun mat acc -> 
    match mat with 
    |[]::_ -> rev acc 
    |_     -> transp_aux (queue_col mat) ((tete_col mat)::acc)
  in transp_aux mat [];;
 
(* Opération terme à terme sur deux lignes
   [x1;x2;x3;...] [y1;y2;y3;...] -> [x1 op y1; x2 op y2; x3 op y3;...]
   On utilise rev_map2 pour rester tail-recursive yeah...*)
let vecs_op = fun op v1 v2 ->
  rev (rev_map2 op v1 v2);;
 
(* opération terme à terme sur deux matrices
   [ligne1; ligne2;...] [ligne'1; ligne'2; ...] -> [vecs_op op ligne1 ligne'1 ; vecs_op op ligne2 ligne'2 ; ....] *)
let mats_op = fun op m1 m2 ->
  vecs_op (vecs_op op) m1 m2;;
 
 
(* opération sur les coefs d'une ligne
   [op v1 ; op v2 ; ...] *)
let map_ligne = fun op ligne ->
  rev (rev_map (fun x -> op x) ligne);;
 
(* opération sur les coefs d'une matrice *)
let map_mat = fun  op mat ->
  rev (rev_map (fun ligne -> (map_ligne op  ligne)) mat);;
 
 
(* som de tous les termes d'une ligne (v1 + v2 +...) *)
let som_vec = fun a vec  ->
  fold_left a.som a.zero vec;;
 
let prod_vec = fun a vec  ->
  fold_left a.prod a.un vec;;
 
 
(* Multiplication d'une ligne par une matrice :
   plm ligne1 est une fonction qui à une ligne l2 associe
   la somme du produit terme à terme de ligne1 avec l2.
   Un map de cette fonction sur la transposée de m2 fait
   ce qu'on veut en fin de fonction.
   (version 1) -> 
*)
let ligne_prod_mat = fun a m2 ligne1 ->
  let f = fun s l c -> a.som s (a.prod l c) in  
  let plm = fun l1 l2 -> fold_left2 f a.zero l1 l2 in
  rev (rev_map (plm ligne1) (transpose m2));;
 
(* ou bien (version 2) avec du filtrage par motif avec la transposée -> *)
let ligne_prod_mat_pm = fun a m2 ligne1 ->
  let tr_m2 = transpose m2 in
  let rec prod_aux = fun tr_m acc -> 
    match tr_m with
    |[]                  -> rev acc
    |tete_mat::queue_mat -> prod_aux queue_mat ( (som_vec a (vecs_op a.prod ligne1 tete_mat))::acc )
  in prod_aux tr_m2 [];;
 
(* ou bien (version 3) avec du filtrage par motif sans la transposée -> *)
let ligne_prod_mat_pm_bis = fun a m2 ligne1 ->
  let rec prod_aux = fun m acc -> 
    match m with
    |[]::_ -> rev acc
    |_     -> prod_aux (queue_col m) ( (som_vec a (vecs_op a.prod ligne1 (tete_col m)))::acc )
  in prod_aux m2 [];;
 
(* Multiplication de deux  matrices :
   on map la fonction partielle ligne_prod_mat a m2
   sur m1 (qui est la liste de ses lignes) *)
let prod_mat a m1 m2 =
  rev (rev_map (ligne_prod_mat a m2) m1);;
  
(* (aij) telle que aij = f(i,j) *)
 
(* Avec Ocaml 4 (oct 2012) :
  let mat_fonc f r c =
  mapi (fun i x -> mapi (fun j z -> f i j) x) (attila entier (r,c));;*)
 
(* Avec Ocaml 3 :  *)
let mat_fonc  = fun f nb_lignes nb_cols ->
  let rec aux_ligne = fun f i j acc ->
    match j with
    |j when (j = nb_cols)  -> rev acc
    |_ -> aux_ligne f i (j+1)  ((f i j)::acc) (* construit la ligne i*)
  in
  let rec aux_col = fun f i j acc_de_lignes ->
    match i with
    |i when (i = nb_lignes) -> rev acc_de_lignes
    |_ -> aux_col f (i+1) j ((aux_ligne f i 0 [])::acc_de_lignes)
  in
  aux_col f 0 0 [];;
 
(* matrice identité de taille n *)
let unite a = fun n ->
  mat_fonc (fun i j  -> if i=j then a.un else a.zero) n n;;
 
 
(* puissance naïve *)
 
let rec pow1 = fun a mat n ->
  match n with
  |0 -> unite a  (n_rows mat)
  |_ ->  prod_mat a mat (pow1 a mat (n-1));;
 
(* puissance dicho *)
 
let rec pow = fun a mat n ->
  match n with
  |0 -> unite a  (n_rows mat)
  |1 -> mat
  |n when (n mod 2 = 0) -> pow a (prod_mat a mat mat)  (n/2)
  |_ ->  prod_mat a mat (pow a mat (n-1));;
 
 
(* trace par mapping : lent car on parcourt beaucoup de chemin inutilement *)
let trace1 = fun a mat ->
  let diag = mats_op (a.prod) mat (unite a (length mat)) in
  fold_left
    (fun s l -> a.som s (hd l))
    a.zero
    (rev_map (fun x -> filter (fun y -> y > 0) x) diag)
;;
 
(* par pliage avec la matrice ((i,j)) pour bricoler un parcours de tableau t(i,j) : lent *)
let trace3 = fun a mat ->
  let n = length mat in
  let mat_coeff = mat_fonc (fun i j -> (i,j)) n n in
  let coeff_diag = fold_left2 (fun s (i,j) x -> if (i=j) then a.som s x else s) a.zero in 
   fold_left2 (fun s x y -> a.som s (coeff_diag x y)) (a.zero) mat_coeff mat;;
 
 
(* on parcourt la matrice en diagonale : n appels à queue_col. On choisit celle-ci : *)
let trace = fun a mat ->
  let rec trace_aux matrice acc =
    match matrice with
    |[]                 -> acc
    |ligne::queue_ligne -> trace_aux (queue_col queue_ligne) ((a.som) (hd ligne) acc)
  in trace_aux mat a.zero;;
 
(* matrice élémentaire : on utilise mat_fonc sur la fonction qui renvoie 1 si ligne=i et colonne = j et 0 sinon *)
let elem a = fun i j n ->
  mat_fonc (fun no_ligne no_col -> if (no_ligne,no_col) = (i,j)  then a.un else a.zero) n n;;
 
(* kelem : idem mais avec k au li de 1 :*)
let kelem a = fun i j k n ->
  mat_fonc (fun no_ligne no_col -> if (no_ligne,no_col) = (i,j)  then k else a.zero) n n;;
(* dans ce cas elem a i j n = kelem a i j a.un n *)
 
 
(* dilatation : * Diag(1 1 1 k 1 1 1...)   Li <- k Li *)
let dilat1 a i k  = fun mat ->
  let n = length mat in
  let d = mat_fonc (fun no_ligne no_col ->
    if (no_ligne , no_col) = (i , i)
    then k
    else if no_ligne = no_col
         then a.un
         else a.zero) n n
  in
  prod_mat a d mat;;
 
 
(* Plus rapide :  on cr&eacute;e une matrice ad hoc en copiant mat sauf &agrave; la
   ligne i qui est remplac&eacute;e par k*ligne i de mat *)
let dilat a i k = fun mat ->
  let li = rev (rev_map (fun x -> a.prod k x) (nth mat i)) (* li <- k*Li *)
  in
  let rec cree_mat = fun m panier no_l-> (* on prend les lignes de m qu'on place dans notre panier sauf Li *)
    match m with
    |[]   -> rev panier (* les premiers &eacute;taient les derniers *)
    |ligne_vue::reste_matrice ->
             if no_l = i  (* on colle la i-eme ligne *)
             then cree_mat reste_matrice (       li::panier) (no_l+1) (* qui sera k*Li *)
             else cree_mat reste_matrice (ligne_vue::panier) (no_l+1) (* sinon on colle la ligne courante de m *)
  in
  cree_mat mat [] 0;; (* on commence &agrave; la ligne 0 avec un panier vide *)
(* on parcourt n lignes et n &eacute;l&eacute;ments dans une ligne : on &eacute;vite prod_mat et ses n^3 op&eacute;rations  *)
 
  
(* transvection : (Id +k E(ij))*Mat  Li <- Li+k.Lj *)
let transvec0 a i j k = fun mat ->
  let n = length mat in
  let t = mats_op (a.som) (unite a n) (kelem a i j k n) in
  prod_mat a t mat;;
(* Pb : unit&eacute;, kelem et a.som entra&icirc;nent 3 parcours de tableaux *)
 
(*  En fait, on peut se d&eacute;brouiller avec un seul appel &agrave;  mat_fonc *)
let transvec1 a i j k = fun mat ->
  if i = j 
  then dilat a i (a.som a.un k) mat (* si i=j alors Li <- Li+k*Li = (1+k)*Li  *)
  else
    let n = length mat
    in
    let t = mat_fonc
      (fun r l -> if (r=l)
                  then a.un
                  else if (r,l) = (i,j)
                       then k
                       else a.zero
      ) n n
    in
    prod_mat a t mat;;
(* PB : plus qu'un seul mat_fonc mais tjs prod_mat en n^3 *)
 
 
(* On peut alors fabriquer une matrice ad hoc avec li = Li + k*Lj *)
let transvec a i j k = fun mat ->
  let klj = rev (rev_map (fun grain -> a.prod k grain) (nth mat j)) (* lj <- k*Lj et rev (rev_map) = map mais terminal *)
  in
  let li = vecs_op (a.som) (nth mat i) klj        (* li <- Li + k*Lj *)
  in
  let rec cree_mat = fun matrice hotte no_l-> (* on copie matrice dans hotte sauf &agrave; la ligne i qui re&ccedil;oit li *)
    match matrice with
    |[]                       -> rev hotte  (* la premi&egrave;re ligne est tout au fond de la hotte : on retourne celle-ci *)
    |ligne_lue::reste_matrice -> if no_l = i
             then cree_mat reste_matrice (       li::hotte) (no_l+1)  (* on ins&egrave;re li au bon endroit (no_l = i) *)
             else cree_mat reste_matrice (ligne_lue::hotte) (no_l+1)  (* sinon on ins&egrave;re la ligne de mat inchang&eacute;e *)
  in
  cree_mat mat [] 0;; (* on commence &agrave; la ligne 0 avec une hotte vide *)
 
 
(* &eacute;change de lignes : cf poly -> produit par quatre matrices  *)
let swap0 a i j = fun mat ->
  let s = mats_op a.som in
  let n = length mat in
  let p = s (unite a n) (s (kelem a i i (a.sous a.zero a.un) n) (s (kelem a j j (a.sous a.zero a.un) n) (s (elem a i j n) (elem a j i n)))) in
  prod_mat a p mat;;
(* ici encore, trop gourmand -> on cherche &agrave; exprimer la matrice d'&eacute;change avec un seul mat_fonc*)
 
(* Plus &eacute;conomique : directement avec un seul mat_fonc *)
let swap1 a i j = fun mat ->
  let n = length mat in
  let s = mat_fonc
    (fun r c -> if (r,c) = (i,j) or (r,c) = (j,i)
                then a.un
                else if (r != i && r != j) && r=c
                     then a.un
                     else a.zero
    ) n n
  in
  prod_mat a s mat
;;
 
(* encore plus &eacute;conomique : matrice ad hoc *)
let swap i j  = fun mat ->
  let rec remplit_panier = fun panier1 panier2 no_l->
    match panier1 with
    |[]                           -> rev panier2 (* on retourne le panier car les premiers sont au fond *)
    |grappe_vue::reste_du_panier1 -> if no_l = i
             then remplit_panier reste_du_panier1 ((nth mat j)::panier2) (no_l+1) (* on met lj dans Li*)
             else if no_l = j
                  then remplit_panier reste_du_panier1 ((nth mat i)::panier2) (no_l+1) (* on met li dans Lj *)
                  else remplit_panier reste_du_panier1 ( grappe_vue::panier2) (no_l+1) (* on met la grappe1 sinon *)
  in
  remplit_panier mat [] 0 (* le panier2 est vide, le panier1 est mat, on commence &agrave; la grappe no 0 *)
;; 
 
(* &eacute;limination ligne i *)
let elim no_ligne = fun mat ->
  let rec cree_mat = fun mat_a_cloner panier no_ligne_courant ->
    match mat_a_cloner with
    |[]                            -> rev panier (* les premi&egrave;res grappes sont au fond du panier *)
    |ligne_vue::reste_mat_a_cloner -> if no_ligne_courant = no_ligne
             then cree_mat reste_mat_a_cloner panier (*inchang&eacute;*) (no_ligne_courant+1)
             else cree_mat reste_mat_a_cloner (ligne_vue::panier) (no_ligne_courant+1)
  in
  cree_mat mat [] 0
;;
 
 
 
(*
 
  *********************** G A U S S     J O R D A N **********************************
  
*)
 
let concat_mat = fun mat1 mat2 ->
  let collage_lignes = fun l1 l2 ->  rev_append (rev l1) l2  (* = append l1 l2 mais terminal *)
  in
  rev (rev_map2 collage_lignes mat1 mat2);; (* = map2 collage_lignes mat1 mat2 mais terminal *)
  (* j'applique (map in technical english :-) ) la fonction collage_lignes aux lignes des matrices 1 et 2
      maporum lignosa comme on dit &agrave; Poudlard... *)
 
(* teste si une matrice est carr&eacute;e : renvoie vrai si mat est carr&eacute;e  *)
let est_carree = fun mat -> 
  n_rows mat = n_cols mat;;
(* pour agr&eacute;menter le debugger *)
exception Matrice_non_carree;;
 
(* colle l'identit&eacute; pour le tableau de recherche de l'inverse par Gauss-Jordan
   Uniquement  pour les matrices carr&eacute;es car c'est
   pour obtenir l'inverse : pas la peine de se fatiguer sinon *)
let mat_to_tab a mat =
  if est_carree mat then
    let n = n_rows mat in
    let u = unite a n in
    concat_mat mat u
  else raise Matrice_non_carree;;
 
 
(* extrait le demi-tableau de droite *)
(* on commence par extraire la demi-ligne de droite
   en extrayant la demi_ligne de gauche de la ligne retourn&eacute;e
   if you see what I mean *)
let bout_ligne n = fun ligne ->
  let engil = rev ligne (* on retourne la ligne *)
  in 
  let rec fabrique_ligne = fun ligne_aux panier no_col ->
    if no_col = n then panier
    else fabrique_ligne (tl ligne_aux) ((hd ligne_aux)::panier) (no_col + 1)
  in
  fabrique_ligne engil [] 0;; (* panier vide et on commence colonne no 0 *)
 
(* on applique &agrave; tout le tableau : maporum lignosa *)
let tab_to_mat = fun mat -> 
  rev (rev_map (bout_ligne (n_rows mat)) mat);; (* rev (rev_map ) = map mais terminal *)
    
 
(*
                     TRIANGULATION DE GAUSS
*)
 
let permut_circi = fun mat i ->
  rev_append (rev            (* on colle de fa&ccedil;on terminale ...*)
                (elim i mat) (*... &agrave; la matrice priv&eacute;e de sa ligne i ...*)
  ) [nth mat i];;            (*... la ligne i tout en bas *)
 
 
(* cherche un pivot sur la colonne no pos*)
let trouve_pivot = fun a mat pos -> 
  let rec cherche = fun cpt matrix -> (* r&eacute;cursion sur le nombre de permutations *)
    if cpt = (length matrix) then matrix (* on a tout permut&eacute; et tjs pas de pivot -> on abandonne *)
    else 
      if a.egal (nth (nth matrix pos) pos)  (a.zero) (* si m[pos][pos] est nul...*)
      then
        let mat_perm = permut_circi matrix pos  (* on permute *)
        in
        cherche (cpt + 1) mat_perm (* et on tente notre chance avec la nouvelle ligne *)  
      else matrix (* si m[pos][pos] != 0 alors on a notre candidat pivot *)
  in
  cherche pos mat(* on cherche &agrave; partir de la ligne no pos dans la mat de d&eacute;part *)
;; 
 
(*  Li <- kj*Li - ki*Lj *)
let comb = fun a i ki j kj mat ->
  transvec a i j (a.sous a.zero kj) (dilat a i ki mat);;
 
(* on fait le vide sous mat[pos][pos] *)
let vide_sous_pivot = fun a mat pos ->
  let p = trouve_pivot a mat pos (* on met un pivot en position pos pos *)
  in
  let pivot = nth (nth p pos) pos (* on appelle ce pivot pivot *)
  in
  if a.egal pivot a.zero then mat (* s'il est nul, on ne pourra pas faire grand chose avec *)
  else
    let rec fait_vide = fun i m -> (* on annule mat[i][pos] *)
      if i = (n_rows m) then m (* si on est arriv&eacute; en bas de la matrice, c'est fini *)
      else  (* sinon on combine : Li <- lpos,pos * Li - li,pos * Lpos : *)
        fait_vide (i+1) (comb a i pivot pos (nth (nth m i) pos)  m) 
    in
    fait_vide (pos+1) p (* on commence sous la ligne o&ugrave; se trouve le pivot *)
;;
 
 
 
 
 
(* on peut enfin trigonaliser....*)
let tri_gauss a = fun mat -> 
  let rec tri_aux = fun i m  -> (* on applique vide_sous_pivot de haut en bas de gauche &agrave; droite*) 
    if i = (n_rows m) then m (* on est tout en bas : c'est fini *)
    else
      let v = vide_sous_pivot a m i (* v est la matrice rang&eacute;e jusqu'&agrave; la colonne i *)
      in
      tri_aux (i+1) v (* et on va voir la colonne d'apr&egrave;s *)
  in tri_aux 0 mat ;;
 
(* on adapte vide_sous_pivot pour vider aussi au-dessus *)
let vide_autour_pivot a = fun mat pos ->
  let p = trouve_pivot a mat pos in
  let pivot = nth (nth p pos) pos in
  if pivot = a.zero then mat
  else
    let rec vide i m  =
      if i = (n_rows m) then m
      else
        if i = pos then vide (i+1) m
        else 
          vide (i+1) (comb a i pivot pos (nth (nth m i) pos)  m) 
    in vide 0 p;;
 
(* tri_gauss devient l&eacute;*)
let le a = fun mat ->
  let rec le_aux = fun i m ->
    if i = (n_rows m) then m
    else
      let v = vide_autour_pivot a m i in
      le_aux (i+1) v 
  in le_aux 0  mat;;
 
 
(* simplifie si possible une ligne par le premier élément non nul
   recontré; (normalement le candidat pivot) *)  
let simp_vec a = fun vec ->
  let rec cherche = fun v -> 
  match v with
  |[] -> vec
  |t::q -> if a.egal t a.zero then cherche q
    else map_ligne (fun x -> (a.div x t)) vec
  in cherche vec;;
 
 
(* lr&eacute; : on simplifie les lignes de la l&eacute; *)
let lre a = fun mat -> 
  rev(rev_map (fun ligne -> simp_vec a ligne) (le a mat));;
 
 
 
 
(* extrait la diagonale d'une matrice carr&eacute;e*)
 
let diag mat =
  let rec aux m acc =
  match m with
  |[] -> rev acc
  |_  -> aux (queue_col (tl m)) ((hd (hd  m))::acc)
  in aux mat [];;
 
 
 
(* inverse d'une matrice *)
exception Matrice_non_inversible;;
 
let inverse a mat =
  if est_carree mat then
    let r = le a (mat_to_tab a mat) in
    let d = diag r in
    if mem a.zero d then
      raise Matrice_non_inversible
    else
      let rec reduit m di i =
        match di with
        |[] -> tab_to_mat m
        |t::q -> reduit (dilat a i (a.div (a.un) t) m) q (i+1)
      in reduit r d 0
  else raise Matrice_non_inversible;;
 
 
 
(*
            D&Eacute;TERMINANT
  *)
 
 
 
(* permutation circulaire des lignes i &agrave; n  *)
let compte_permut_circi a mat i =
  let parite = if ((n_rows mat) -i) mod 2 = 0 then (a.sous a.zero a.un) else a.un in
  ( (rev_append (rev (elim i mat)) [nth mat i]) , parite);;
 
 
let nb_permut a mat i =
  snd (compte_permut_circi a mat i);;
 
 
(* pour choisir le pivot max *)
let max = fun a x y ->
  if (a.ordre) x y then x else y;;
 
let abs = fun a x ->
  let oppx = a.sous a.zero x in
  max a x oppx;;
  
let max_abs = fun a x y ->
  if a.ordre (abs a x) (abs a y) then x else y;;
 
exception Max_indefini;;
 
let rec max_abs_vec = fun a vec ->
  match vec with
  |[] -> raise Max_indefini
  |t::[] -> t
  |t::q -> max_abs a t (max_abs_vec a q);;
 
 
 
 
 
 
let compte_trouve_pivot a mat pos =
  let rec trouve i m nb =
  match m with 
  |tete::queue ->
    if i = (length m) then (m,nb)
    else 
      if a.egal (nth (nth m pos) pos)  (a.zero)
      then
        let c = compte_permut_circi a m pos in
        trouve (i+1) (fst c) (a.prod nb (snd c)) 
      else (m,nb)
  |_ -> (m,nb)
  in trouve 0 mat a.un;;
 
 
(*let compte_trouve_pivot_max a mat pos =
  let maxi = max_abs_vec a (coupe pos (extract_col mat pos)) in
  let rec trouve i m nb =
  match m with 
  |tete::queue ->
    if i = (length m) then (m,nb)
    else 
      if not (a.egal (nth (nth m pos) pos)  maxi)
      then
        let c = compte_permut_circi a m pos in
        trouve (i+1) (fst c) (a.prod nb (snd c)) 
      else (m,nb)
  |_ -> (m,nb)
  in trouve 0 mat a.un;;
*)
 
let compte_vide_sous_pivot a mat pos =
  let c = compte_trouve_pivot a mat pos in
  let p = fst c in
  let pivot = nth (nth p pos) pos in
  let k = snd c in
  if pivot = a.zero then (mat,a.un)
  else
    let rec vide i m ct =
      if i = (n_rows m) then (m,ct)
      else
        vide (i+1) (comb a i pivot pos (nth (nth m i) pos)  m) (a.prod ct pivot)
    in vide (pos+1) p k;;
 
 
let compte_tri_gauss a mat =
  let rec aux i m k =
    if i = (n_rows m) then (m,k)
    else
      let v = compte_vide_sous_pivot a m i in
      aux (i+1) (fst v) (a.prod k (snd v))
  in aux 0 mat a.un;;
 
  
 
  
(* D&eacute;terminant de la matrice par Gauss
   Attention aux permutations et aux mult de lignes...*)
 
let det a mat =
  let g = compte_tri_gauss a mat in
    a.div (prod_vec a (diag (fst g))) (snd g);;
 
 
 
 
(*  Syst&egrave;me  AX = B *)
 
exception Dim_erreur;;
 
  
let mat_float_of_int = fun m -> 
  map (fun x -> map (fun z -> float_of_int z) x) m;;
 
exception Pas_de_solution_unique;;
 
let systeme1 a m v =
  let r = lre a (concat_mat m v) in
  let d = diag r in
  if not (est_carree m) or (mem a.zero d) then
    raise  Pas_de_solution_unique
  else
    extract_col r (length m)
;;
 
  
exception Pas_de_solution;;
 
(* extraction d'une sous-liste Li..Lk *)
 
let sous_liste i k l =
  let rec local n acc =
    if n > k then rev acc
    else local (n+1) ((nth l n)::acc)
  in local i [];;
 
let sous_matrice m i1 i2 j1 j2 =
  sous_liste i1 i2 (rev (rev_map  (sous_liste j1 j2) m));;
 
  
(* AX = B dans tous les cas de figure *)
 
let systeme2 a m v =
  let rr = lre a (concat_mat m v) in
  let rrr = sous_matrice rr 0 ((n_rows rr) -1)0 ((n_cols rr) -2) in
  let d = diag rrr in
  let rang = length (filter (fun x -> (a.egal x a.un)) d) in
  let dim = min (n_rows m) (n_cols m) in
  let r = 
  if rang = dim then
   sous_matrice rr 0 (rang-1) 0 ((n_cols rr) -1) 
  else
    rr
  in
    let vx = extract_col r (n_cols m) in
    if (rang = dim) && (not (mem a.zero d)) then
      sous_matrice r 0 ((n_rows r) -1) (n_cols m) (n_cols m) 
    else
      let bout = sous_liste rang ((n_rows r) -1) vx in 
      if not (a.egal (prod_vec a bout) a.zero) then
        raise Pas_de_solution
      else
        sous_matrice r 0 (rang -1) rang ((n_cols r) - 1) 
;;
      
(*
 
  Matrice inverse dans Z/nZ - HILL
 
*)
 
 
 
(* un peu d'arithm&eacute;tique *)
 
let rem a b =
  let m = a mod b in
  if m >= 0 then m else m + b;;
 
type classe  =
  {representant : int ; modulo : int};;
 
let (%) a n = {representant = (rem a n); modulo = n};;
 
exception Classes_incompatibles;;
 
let (+%) a b =
  {modulo =
      if a.modulo = b.modulo
      then a.modulo 
      else raise Classes_incompatibles ;
   representant =
      rem ( (a.representant) + (b.representant)) (a.modulo)
  };; 
 
let classe_canon n a = {representant = (rem a n); modulo = n};;
 
let est_dans_classe n c =
   c.representant mod c.modulo = rem n c.modulo;; 
 
 
let loi_induite loi =
  fun a  b -> 
  {modulo =
      if a.modulo= b.modulo
      then a.modulo 
      else raise Classes_incompatibles ;
   representant =
      rem (loi (a.representant) (b.representant)) (a.modulo)
  };; 
 
let ( +% ) = loi_induite  ( + ) ;;
let ( *% ) = loi_induite  ( * ) ;;
let ( -% ) = loi_induite  ( - ) ;;
 
let egal_mod x y = 
  let xr = x.representant in
  let xm = x.modulo in
  let yr = y.representant in
  let ym = y.modulo in
  (xr) % (xm)  = (yr) % (ym);;
 
 
let sup_mod x y = 
  let xr = x.representant in
  let xm = x.modulo in
  let yr = y.representant in
  let ym = y.modulo in
  (xr) % (xm)  >= (yr) % (ym);;
 
(*
(9%5) *% (7%5);;
*)
 
(* Euclide &eacute;tendu : coeff de B&eacute;zout *)
 
 
let rec euclide2 a b =
  let rec euc(u,v,r,u',v',r') =
    if r' = 0 then
      [|u;v;r|]
    else
      let q = r/r' in
      euc(u',v',r',u - q * u', v - q * v', r - q * r')
  in euc(1,0,a,0,1,b);;
 
 
 
exception Non_inversible ;;
 
 
let inv_mod a =
  let b = euclide2 a.representant a.modulo in
  if b.(2) = 1 then
     (b.(0))%(a.modulo)
  else
    raise Non_inversible;;
 
let string_of_q c =
(string_of_int (rem c.representant c.modulo))^"%"^(string_of_int c.modulo);;
 
let quotient n  =
  {zero      = 0%n ;
   un        = 1%n ;
   som       = ( +% );
   prod      = ( *% );
   sous      = (-%);
   egal      = egal_mod ;
   div       = (fun x y -> ( *% ) x (inv_mod y));
   ordre     = sup_mod;
   to_string = string_of_q;
  };;
 
 
 
(*
        Rationnels
 *)
 
 
let rec pgcd_pos a b =
  if a >= b then
    if b = 0 then a 
    else pgcd_pos  b (a mod b)
  else pgcd_pos b a;;
 
let pgcd a b = pgcd_pos (abs entier  a) (abs entier b);;
 
 
type rationnel = {
num : int;
den : int;
};;
 
  
let frac a b =
  let d = pgcd a b in
  if (a < 0 && b < 0) or ( a >= 0 && b < 0) then
    {num = -a/d ; den = -b/d}
  else 
    {num = a/d ; den = b/d}
;;
 
let (//) a b = frac a b ;;
 
  
let r_plus r1 r2 = 
(r1.num * r2.den + r2.num * r1.den) // (r1.den * r2.den);; 
r_plus (2 // 3) (5 // 6);;
 
let r_mult r1 r2 =
(r1.num * r2.num) // (r1.den * r2.den);;
 
let r_sup r1 r2 = 
(r1.num * r2.den) >= (r2.num * r1.den);;
 
let r_egal r1 r2 = 
(r1.num // r1.den) = (r2.num //r2.den);;
 
 
 
let r_inv r = 
if r.num = 0 then raise Non_simplifiable
else r.den // r.num;;
 
 
let r_to_string r =
  let s = (r.num // r.den) in
  if s.num = 0 then "0"
  else if s.den = 1 then (string_of_int s.num)
  else (string_of_int s.num)^"/"^(string_of_int s.den);;
 
  
let r_opp r =
(- r.num) // (r.den);;
 
let rat = {
   zero = 0//1;
   un   = 1//1;
   som  = r_plus;
   prod = r_mult;
   sous  = (fun x y -> r_plus  x (r_opp y));
   div = (fun x y -> r_mult x (r_inv y));
   to_string = r_to_string;
   egal = r_egal;
   ordre = r_sup;
  };;
 
(* syst&egrave;mes et inverse *)
 
let rat_of_int x = x // 1;;
 
let mat_rat_of_int = fun m -> 
  map (fun x -> map (fun z -> rat_of_int z) x) m;;
 
 
 
tri_gauss reel (attila reel (50,50));;

courtesy of webmatter.de