%% syntaxe : M i j any --> depose any dans M en a_ij
/put_ij {
5 dict begin
   /a exch def
   /j exch def
   /i exch def
   /M exch def
   /L M i get_Li def
   L j a put
   M i L put_Li
end
} def

%% syntaxe : M i j get_ij --> le coeff c_ij
/get_ij {
   3 1 roll   %% j M i
   get_Li     %% j L_i
   exch get
} def

%% syntaxe : M i L put_Li --> remplace dans M la ligne Li par L
/put_Li {
   put
} def

%% syntaxe : M i get_Li --> la ligne Li de M
/get_Li {
   get
} def

%% syntaxe : M i --> depose sur la pile la colonne d'indice i de M
/get_Ci {
4 dict begin
   /j exch def
   /M exch def
   M dimmatrix pop /m exch def
   /i 0 def
   [
   m {
      M i j get_ij
      /i i 1 add store
   } repeat
   ]
end
} def

%% syntaxe : M i j exch_l --> echange les lignes i et j
/exch_L {
3 dict begin
   /j exch def
   /i exch def
   /M exch def
   M i get_Li
   M i
   M j get_Li
   put_Li
   M j 3 -1 roll
   put_Li
end
} def

%% syntaxe : B M i j elimine_coeff --> remplace la ligne j par Lj - (a_ji/a_ii) Li
%% et applique la meme transformation au vecteur colonne B
/elimine_coeff {
4 dict begin
   /j exch def
   /i exch def
   /M exch def
   /B exch def
   %% le coeff
   M j i get_ij
   M i i get_ij
   div neg
   /k exch def
   
   B j
   B j get 
   B i get k mul add
   put

   M j
   M i get_Li {k mul} apply
   M j get_Li 
   {add} fuzapply
   put_Li
end
} def

%% syntaxe : B M i echelonne --> echelonne avec le coeff a_ii les lignes Li+1 .. Ln
/echelonne {
4 dict begin
   /i exch def
   /M exch def
   /B exch def
   /j i 1 add def
   M length 1 sub i sub {
      B M i j elimine_coeff
      /j j 1 add store
   } repeat
end
} def

%% syntaxe : B A pivot_gauss --> A est echelonne
%% ATTENTION : il est suppose que A s'echelonne "bien" (determinant non nul)
/pivot_gauss {
4 dict begin
   /A exch def
   /B exch def
   /i 0 def
   
   A length 1 sub {
      A i i get_ij 0 ne               %% si a_ii <> 0
         {B A i echelonne}            %% on echelonne
         {                            %% sinon on recherche j tq a_ji <> 0
            /j i 1 add def
            {
            A j i get_ij 
            0 ne
               {                      %% on a trouve
                  A i j exch_L        %% on echange les lignes Li et Lj
                  B i get             %% idem pour B
                  B i
                  B j get put
                  B j 3 -1 roll put
                  B A i echelonne     %% puis on echelonne
                  exit                %% et on sort de la boucle
               }
               {                      %% on a pas trouve
                  /j j 1 add store    %% on voit la suite
               }
            ifelse
            } loop
         }
      ifelse
      /i i 1 add store                %% colonne suivante
   } repeat
end
} def

%% B A 0 echelonne
%% B A 1 echelonne

%% B A 0 1 elimine_coeff
%% B A 0 2 elimine_coeff

%% syntaxe : B A solve_trig --> X ou A est une matrice carre
%% triangulaire superieure et X est l'unique vecteur solution de AX = B
/solve_triang {
6 dict begin
   /A exch def
   /B exch def
   /N A length def
   /sols N array def
   /i N 1 sub def
   
   N {
   %% la solution x_i
   /j i 1 add def
   
   sols i
   B i get
   N 1 sub i sub {
      A i j get_ij
      sols j get
      mul sub
      /j j 1 add store
   } repeat
   A i dup get_ij div
   put
   /i i 1 sub store
   } repeat
   sols
end
} def

%% syntaxe : B A solve_syst --> X ou A est une matrice carre
%% de determinant non nul
/solve_syst {
2 dict begin
   /A exch def
   /B exch def
   B A pivot_gauss
   B A solve_triang
end
} def

%% syntaxe : M view_square_matrix --> depose les coeffs de la matrice carree M sur la pile
/view_square_matrix {
2 dict begin
   /M exch def
   /i 0 def
   (matrix)
   M length {
      M i get 
      aload pop ()
      /i i 1 add store
   } repeat
   (fin matrix)
} def

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%