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));;