CLASS MATRICE AVEC TEMPLATE

/* LAMOME Julien 19/10/07 14:12
; voir notice ci dessous */

/*
Utilisation de la struct ci dessous.
Déclaration :
            - MATRIX<type> nom_matrice(n-ligne,m-colonne) crée une matrice de n lignes sur m colonne d'éléments 'type' de valeurs nul
            - MATRIX<type> nom_matrice(n-ligne/colonne) comme la présédente, mais une matrice carré.
            - MATRIX<type> nom_matrice crée une matrice vide
Acces données :
      lecture/écriture
           donnée de la matrice : nom_matrice(ligne,colonne)
           données de la matrice : nom_matrice[colonne][ligne] ATTENTION, la convention est inverse
      lecture
           sous matrice : nom_matrice(l1,l2,c1,c2) renvoie la matrice allant de la ligne l1 à l2 et de la colonne c1 à c2
           nombre de ligne :  nom_matrice.size_M().n
           nombre de colonnes : nom_matrice.size_M().m
           donnée de la matrice : nom_matrice(i,'r'|'c') renvoi la ligne(r) ou colonne(c) n° i
Fonctions membres :
          Id() : transforme la matrice en matrice identité, si elle est carré, sinon, ne fait rien
          inverse() : calcul l'inverse d'une matrice carré, quitte le programme sinon
          transpose() : prend la transposer de la matrice.
          affiche() : affiche la matrice, soit avec <stdio.h> soit avec <iostream> include précédent
          det() : calcul le déterminant de la matrice si carrée, renvoi -1 sinon.
          size() : retourne le nombre de lignes d'une matrice carré, et le nombre d'éléments d'une non carré
          size_M() : retourne la taille de la matrice dans la structure de type M_taille.
          Tr() : calcul la trace de la matrice (somme diagonal), renvoi 0 si non carré.
          clear() : rend la matrice vide.
          genre() : retourne le type de matrice : IS_VECTOR, si une des dimmension est de 1; IS_SCALAR, si les 2 dim sont à 1; 
                    VIDE si les dim sont nul; IS_MATRICE ; et MATRICE_DIM_ERR si les dimensions sont indéterminés
          push_back(type) : rempli une matrice par la fin, la rend carré. !! UTILISATION DECONSEILLER !!
          fct(T fonc(T)) : applique la fonction fonc a chaque membre de la matrice
Supporte :
         - ostream<<nom_matrice
         - ostream<<nom_matrice.size_M()
         - la plupart des opération algébrique de base : +,-,*,=
         - +,* par "type", et par "vector<type>"
         - les fonctions inverse(matrice), transpose(matrice).
         - la "procedure" TLU, qui à partir de la premiere matrice : A, retourne les matrice triangulaires basse et haute, L, et U, tel que A=L*U
Compilation tester avec succes avec gcc, jusqu'à la version 3.3.1 compatible cygwin */

#ifndef _MATRIX
#define _MATRIX

#include <vector>
#include <cmath>
#ifdef _STDIO_H_
#define __ENVOI(x) printf(x);
#define __ENVOI2(x,y) printf(x,y);
#else
#if defined(_CPP_IOSTREAM) || defined(__IOSTREAM__)   
#define __ENVOI(x) std::cout<<x;
#else
#define __ENVOI(x) throw (x);
#endif
#endif
struct M_taille
       {
       int n;
       int m;
       };
#define IS_MATRICE 2
#define IS_VECTOR 1
#define IS_SCALAR 0
#define VIDE 3
#define MATRICE_DIM_ERR -1
template <class T>struct MATRIX
      {
      // variable de la matrice : ses données + sa taille 
private:
      int n,m;
      std::vector<T> M;
      // constructeur, destructeur, operateur, et fonctions
public:
      MATRIX() {n=0;m=0;std::vector<T> M;};
      MATRIX(int ,int =0); //constructeur
      MATRIX(std::vector<T>);  //transformation d'un vecteur en matrice (1,n)
      MATRIX(T);// Attention si T=int !!!!! mettre "MATRIX operator=(T)" ne marche pas
      ~MATRIX(){M.clear();M.reserve(0);};
      int size() const;
      M_taille size_M() const;
      void clear();
      T& operator()(int ,int );
      const T operator()(int ,int )const;
      MATRIX operator()(int ,int ,int ,int );
      std::vector<T> operator()(int,char);
//      std::vector<T> operator[](int);
      T* operator[](int);
      const T* operator[](int)const;
      MATRIX operator=(MATRIX );
      MATRIX operator+=(MATRIX );
      MATRIX operator-=(MATRIX );
      MATRIX operator*=(MATRIX );
      MATRIX operator-() const ;
      MATRIX operator+() const {return *this;};
      MATRIX<T> fct(T t(T)) const;// pour appliquer n'importe quelle fonction
      T Tr() const;
      MATRIX& push_back(T);
      int genre() const;
      T det() const;
      MATRIX& Id();
      MATRIX& inverse();
      MATRIX& transpose();
      void affiche() const;
      void resize(int,int);
// Fonctions amies
      template <class _T>friend MATRIX<_T> operator+(const MATRIX<_T>&,_T);
      template <class _T>friend MATRIX<_T> operator*(const MATRIX<_T>&,const _T);
      template <class _T>friend MATRIX<_T> operator*(MATRIX<_T>&,MATRIX<_T>&);
      template <class _T>friend MATRIX<_T> operator+(MATRIX<_T>&,MATRIX<_T>&);
      template <class _T>friend MATRIX<_T> operator&(MATRIX<_T>&,MATRIX<_T>&);
      template <class _T>friend MATRIX<_T> operator-(MATRIX<_T>&,MATRIX<_T>&);
      template <class _T>friend bool       operator==(MATRIX<_T>,MATRIX<_T>);
      template <class _T>friend std::vector<_T> operator*(MATRIX<_T>,std::vector<_T>);
      template <class _T>friend int TLU(const MATRIX<_T>&,MATRIX<_T>&,MATRIX<_T>&);
//      template <class _T>friend std::ostream& operator<<(std::ostream&, MATRIX<_T>);
//      template <class _T>friend MATRIX<_T>& operator>>(std::istream&, MATRIX<_T>&);
      template <class _T,typename T1>friend MATRIX<_T>& operator>>(T1&, MATRIX<_T>&); 
      template <class T1,typename _T>friend T1& operator<<(T1&, MATRIX<_T>);
      };
/*------------fin de la struct MATRIX-----------------------------------------*/
/*------------------opération binaire sur les matrices------------------------*/
template <class T>MATRIX<T> operator+(const MATRIX<T>& ,const T );
template <class T>MATRIX<T> operator+(const T ,const MATRIX<T>& );
template <class T>MATRIX<T> operator*(const MATRIX<T>& ,const T );
template <class T>MATRIX<T> operator*(const T ,const MATRIX<T>& );
template <class T>MATRIX<T> operator+(const MATRIX<T>& ,const MATRIX<T>& );
template <class T>MATRIX<T> operator+(const MATRIX<T>& ,const MATRIX<T>& );
template <class T>MATRIX<T> operator-(const MATRIX<T>& ,const MATRIX<T>& );
template <class T>MATRIX<T> operator-(const MATRIX<T>& ,const MATRIX<T>& );
template <class T>MATRIX<T> operator&(MATRIX<T>& ,MATRIX<T>& ); //multiplication terme à terme
template <class T>MATRIX<T> operator&(const MATRIX<T>& ,const MATRIX<T>& );
template <class T>MATRIX<T> operator*(MATRIX<T>& ,MATRIX<T>& );
template <class T>MATRIX<T> operator*(const MATRIX<T>& ,const MATRIX<T>& );
template<class T> bool      operator==(const MATRIX<T>,const MATRIX<T>);
template<class T> bool      operator!=(MATRIX<T>,MATRIX<T>);
/*------------------opération binaire  matrice vecteur------------------------*/
template <class T>std::vector<T> operator*(MATRIX<T> ,std::vector<T> );
template <class T>std::vector<T> operator*(std::vector<T> ,MATRIX<T> );
/*-----------fonction supplémentaire utile au calculs de matrices-------------*/
template<class T>std::vector<T> operator+(std::vector<T> ,std::vector<T>);
template<class T>std::vector<T> operator*(T ,std::vector<T>);
template<class T>std::vector<T> operator*(std::vector<T>,T );
template<class T>MATRIX<T> transpose(MATRIX<T> );
template<class T>MATRIX<T> inverse(MATRIX<T> );
template<class T>  int     TLU(const MATRIX<T>& A,MATRIX<T>& L,MATRIX<T>& U);
#if defined(_CPP_IOSTREAM) || defined(__IOSTREAM__)   
//template<class T>std::ostream& operator<<(std::ostream& s, MATRIX<T> A);
//std::ostream& operator<<(std::ostream& s, M_taille A);
#endif   // IOSTREAM
/*--------------------- fonctions de flux pour matrices ----------------------*/
template <class _T,typename T1> MATRIX<_T>& operator>>(T1&, MATRIX<_T>&); 
template <class T1,typename _T> T1& operator<<(T1&, MATRIX<_T>);



typedef MATRIX<double> MATRIX2;
/*----------------------------------------------------------------------------*/
//--------début des implementation--------------------------------------------*/
/*----------------------------------------------------------------------------*/
template<class T>void MATRIX<T>::affiche() const
           {
           #ifdef _STDIO_H_
           for (int i=0;i<n;i++)
            {
            printf("| ");
            for (int j=0;j<m;j++)
               printf("%.1e | ",M[i+n*j]);
            printf("\n");
            }
           #else
           #if defined(_CPP_IOSTREAM) || defined(__IOSTREAM__)
           for (int i=0;i<n;i++)
               {
               std::cout<<"| ";
               for (int j=0;j<m;j++)
                   std::cout<<M[i+n*j]<<" | ";
               std::cout<<"\n";
               }
           #endif
           #endif
           }
/*------------------opération binaire sur les matrices------------------------*/
template <class T>MATRIX<T> operator+(const MATRIX<T>& A,const T B)
        {
        MATRIX<T> tmp;
        if (A.m==A.n)
           {
           for (int i=0;i<A.M.size();i++)
               tmp.push_back(A.M[i]+B);
           return tmp;
           }
        MATRIX<T> tmp1(A.n,A.m);
        for (int j=0;j<A.m;j++)
           for (int i=0;i<A.n;i++)
               tmp1(i,j)=A(i,j)+B;
        return tmp1;
        }
template <class T>MATRIX<T> operator+(const T B,const MATRIX<T>& A)
        {
        return A+B;
        }
template <class T>MATRIX<T> operator*(const MATRIX<T>& A,const T B)
        {
        MATRIX<T> tmp;
        if (A.m==A.n)
           {
           for (unsigned int i=0;i<A.M.size();i++)
               tmp.push_back(A.M[i]*B);
           return tmp;
           }
        MATRIX<T> tmp1(A.n,A.m);
        for (int j=0;j<A.m;j++)
           for (int i=0;i<A.n;i++)
               tmp1(i,j)=A(i,j)*B;
        return tmp1;
        }
template <class T>MATRIX<T> operator*(const T B,const MATRIX<T>& A)
        {
        return A*B;
        }
template <class T>MATRIX<T> operator+(MATRIX<T>& A,MATRIX<T>& B)
        {
        if (A.m==B.m & A.n==B.n)      {
        MATRIX<T> tmp;
        for (int i=0;i<A.M.size();i++)
            tmp.push_back(A.M[i]+B.M[i]);
        tmp.n=A.n;
        tmp.m=A.m;
        return tmp;                   }
        __ENVOI("problème de dimension (operator+)");
        return -1;
        }
template <class T>MATRIX<T> operator+(const MATRIX<T>& A,const MATRIX<T>& B)
         {
         MATRIX<T> a=A;
         MATRIX<T> b=B;
         return a+b;
         }
template <class T>MATRIX<T> operator-(MATRIX<T>& A,MATRIX<T>& B)
        {
        if (A.m==B.m & A.n==B.n)      {
        MATRIX<T> tmp;
        for (int i=0;i<A.M.size();i++)
            tmp.push_back(A.M[i]-B.M[i]);
        tmp.n=A.n;
        tmp.m=A.m;
        return tmp;                  }
        __ENVOI("problème de dimension (operator-)");
        return -1;
        }
template <class T>MATRIX<T> operator-(const MATRIX<T>& A,const MATRIX<T>& B)
         {
         MATRIX<T> a=A;
         MATRIX<T> b=B;
         return a-b;
         }
template <class T>MATRIX<T> operator&(MATRIX<T>& A,MATRIX<T>& B)
        {
        if (A.m==B.m & A.n==B.n)      {
        MATRIX<T> tmp;
        T a;
        for (int i=0;i<A.n*A.m;i++)
            {
            a=A.M[i]*B.M[i];
            tmp.push_back(a);
            }
        tmp.n=A.n;
        tmp.m=A.m;
        return tmp;                  }
        return -1;
        }
template <class T>MATRIX<T> operator&(const MATRIX<T>& A,const MATRIX<T>& B)
         {
         MATRIX<T> a=A;
         MATRIX<T> b=B;
         return a&b;
         }
template <class T>MATRIX<T> operator*(MATRIX<T>& A,MATRIX<T>& B)
        {
        if (A.m==B.n)
        {
        MATRIX<T> tmp;
        T a=0;
        for (int j=0;j<B.m;j++)
        for (int i=0;i<A.n;i++)
            {
            for (int k=0;k<A.m;k++)
                a=a+A(i,k)*B(k,j);
            tmp.push_back(a);
            a=0;
            }
        tmp.n=A.n;
        tmp.m=B.m;
        return (tmp);
        }
        if ((A.n==B.m)&(A.m==1)&(B.n==1)) //très con, reviens au test précedent Am=1 et Bn=1 => Am==Bn !
        {
        MATRIX<T> tmp;
        for (int i=0;i<A.n;i++)
        for (int j=0;j<B.m;j++)
            tmp.push_back(A(i,0)*B(0,j));
        return tmp;
        }
        if (A.n==1 & A.m==1)
                return A(0,0)*B;
        if (B.n==1 & B.m==1)
                return B*A;
        __ENVOI("multiplication de matrices de dimension incompatible\n");
        return -1;
        }
template <class T>MATRIX<T> operator*(const MATRIX<T>& A,const MATRIX<T>& B)
         {
         MATRIX<T> a=A;
         MATRIX<T> b=B;
         return a*b;
         }
template<class T> bool operator==(MATRIX<T> A,MATRIX<T> B)
    {
    if( (A.n!=B.n) | (A.m!=B.m) ) return false;
    for (int i=0;i<A.n*A.m;i++)
        if(A.M[i]!=B.M[i]) return false;
    return true;
    }
template<class T> bool operator!=(MATRIX<T> A,MATRIX<T> B){return !(A==B);}
/*------------------opération binaire  matrice vecteur------------------------*/
template <class T>std::vector<T> operator*(MATRIX<T> A,std::vector<T> v)
               {
               if (v.size()!=A.n)
               {
               __ENVOI("vecteur de mauvaise dimension\n");
               }
               std::vector<T> tmp;
               T a ;
               for (int i=0;i<A.n;i++)
                   {
                   a=0;
                   for (int j=0;j<A.m;j++)
                       a=a+A(i,j)*v[j];
                   tmp.push_back(a);
                   }
               return tmp;
               }
template <class T>std::vector<T> operator*(std::vector<T> v,MATRIX<T> A)
         {
         return A*v;
         }
/*---définitions des fonctions et constructeurs de la strut MATRIX------------*/
template <class T>MATRIX<T>::MATRIX<T>(int a,int b)
      {
             n=a;
             m=b;
             if (b==0) m=n;
             for (int i=0;i<n*m;i++)
                 M.push_back(T(0));
      }
template <class T>MATRIX<T>::MATRIX<T>(std::vector<T> v)
         {
         n=v.size();
         m=1;
         M=v;
         }
template <class T>MATRIX<T>::MATRIX<T>(T v)
         {
         n=1;
         m=1;
         M.push_back(v);
         }
template<class T> T& MATRIX<T>::operator()(int x,int y)
             {
             if ( (x<n)&(y<m) )
                return M[x+n*y];
             __ENVOI("indice trop grand [operator()]\n");
             exit(-1);
             }
template<class T> const T MATRIX<T>::operator()(int x,int y)const
             {
             if ( (x<n)&(y<m) )
                return M[x+n*y];
             __ENVOI("indice trop grand [operator()]\n");
             exit(-1);
             }
template<class T> MATRIX<T> MATRIX<T>::operator()(int a,int b,int c,int d)
    {
    if (a>b | b>n | c>d | d>m) {__ENVOI("erreur indice sous matrice");return 0;}
    MATRIX<T> tmp(b-a+1,d-c+1);
    for (int i=0;i<b-a+1;i++)
    for (int j=0;j<d-c+1;j++)
        tmp(i,j)=M[i+a+n*(j+c)];
    return tmp;
    }
template<class T> std::vector<T> MATRIX<T>::operator()(int x,char c)
               {
               std::vector<T> tmp;
               if (c=='c')
                  {
                  for (int i=0;i<n;i++)
                     tmp.push_back(M[i+n*x]);
                  return tmp;
                  }
               if (c=='r')
                  {
                  for (int i=0;i<m;i++)
                     tmp.push_back(M[x+n*i]);
                  return tmp;
                  }
               return tmp;
               }
/*template<class T> std::vector<T> MATRIX<T>::operator[](int r)
               {
               std::vector<T> tmp(m);
               for (int i=0;i<m;i++)
                   tmp[i]=M[i+n*r];
               return tmp;
               }*/
template<class T> T* MATRIX<T>::operator[](int r) //renvoie l'adresse T* du premier element de la colonne r, donc T[l] sera la ligne l de la colonne r
               {
               int nr=n*r;
               return &M[nr];
               }
template<class T>const T* MATRIX<T>::operator[](int r)const //renvoie l'adresse T* du premier element de la colonne r, donc T[l] sera la ligne l de la colonne r
               {
               int nr=n*r;
               return &M[nr];
               }
template <class T> MATRIX<T> MATRIX<T>::operator=(MATRIX S)
              {
                  M=S.M;
                  n=S.n;
                  m=S.m;
              return *this;
              }
template <class T> MATRIX<T> MATRIX<T>::operator+=(MATRIX S){*this=*this+S;return *this;}
template <class T> MATRIX<T> MATRIX<T>::operator-=(MATRIX S){*this=*this-S;return *this;}
template <class T> MATRIX<T> MATRIX<T>::operator*=(MATRIX S){*this=*this*S;return *this;}
template <class T> MATRIX<T> MATRIX<T>::operator-()const {return T(-1)**this;}             
template <class T> int MATRIX<T>::size() const
          {
          if (M.size()!=n*m) return -1;
          if (m==n) return n;
          return n*m;
          }
template <class T>MATRIX<T>& MATRIX<T>::Id()
        {
        if (m==n)
        for (int i=0;i<n;i++)
            M[i+n*i]=1.0;
        return *this;
        }
template <class T> void MATRIX<T>::clear()
           {
           M.clear();
           n=0;
           m=0;
           }
template <class T> T MATRIX<T>::Tr() const
             {
             T tr=0;
             for (int i=0;i<n;i++)
                 tr+=M[i*(1+n)];
             if (m==n)
             return tr;
             return 0;
             }
template <class T>MATRIX<T>& MATRIX<T>::inverse()
           {
           if (m!=n)
               {
               __ENVOI("inversion de matrices non carré non implémenter!\n");
               return *this;
               }
           MATRIX<T> MM(n),M1(n),Mi1(n),Minv(n);
           MM=*this;
           M1=MM;
           Mi1.Id();
           Minv=Mi1;
           for (int i=0;i<n;i++)
             {
             for (int j=0;j<n;j++)
                 {
                 if (MM(i,i)==0)
                    {
                    __ENVOI("inversion pivot GAUSS impossible : division par 0\n");
                    return *this;
                    }
                 M1(i,j)=MM(i,j)/MM(i,i);
                 Mi1(i,j)=Minv(i,j)/(MM(i,i));
                 }
             MM=M1;
             Minv=Mi1;
             for (int k=0;k<n;k++)// mise zéros de la colonne
             if (k!=i)
             for (int j=0;j<n;j++)
                 {
                 M1(k,j)=MM(k,j)-MM(i,j)*MM(k,i);
                 Mi1(k,j)=Minv(k,j)-Minv(i,j)*MM(k,i);
                 }
             MM=M1;
             Minv=Mi1;
             /*for (int k=0;k<n;k++)//mise à zéro de la ligne. ap essai : pas la peine
             if (k!=i)
             for (int j=0;j<n;j++)
                 {
                 M1(j,k)=MM(j,k)-MM(i,j)*MM(i,k);
                 Mi1(j,k)=Minv(j,k)-Minv(i,j)*MM(i,k);
                 }
             MM=M1;
             Minv=Mi1;*/
             }
           *this=Minv;
           return *this;
           }
template <class T> T MATRIX<T>::det() const
       {
           if (m!=n)
               {
               __ENVOI("déterminant de matrices non carré non implémenter!\n");
               return -1;
               }
       if (n==1) return M[0];
       //if (n==2) return M[0]*M[3]-M[1]*M[2];
       T d=0.0;
       MATRIX tmp;
       for (int i=0;i<n;i++)
           {
           for (int j=n;j<n*n;j++)
           if (int((j-i)*1.0/n)!=(j-i)*1.0/n)
               tmp.push_back(M[j]);
           d=d+pow(-1,double(i))*M[i]*tmp.det();
           tmp.clear();
           }
       return d;
       }
template <class T> MATRIX<T>& MATRIX<T>::transpose()
{
    MATRIX<T> P(m,n);
    for (int i=0; i<m; i++)
        for (int j=0; j<n; j++)
        {       P(i,j)=M[j+n*i];}
    *this=P;
    return *this;
}
template <class T> MATRIX<T>& MATRIX<T>::push_back(T x)
             {
             M.push_back(x);
             n=int(sqrt(double(M.size())));
             m=n;
             return *this;
             }
template <class T> int MATRIX<T>::genre() const
         {
         if ((m==1)&(n==1))
            return IS_SCALAR;
         if ((m==1)|(n==1))
            return IS_VECTOR;
         if ((m==0)&(n==0))
            return VIDE;
         if ((m>1)&(n>1))
            return IS_MATRICE;
         else
            return MATRICE_DIM_ERR;
         }
template<class T>void MATRIX<T>::resize(int l,int c)
    {
    if(l*c>m*n)
        {__ENVOI("changement de taille impossible : pas assez d'éléments\n");return;}
    if(l*c==m*n)
        {n=l;m=c;return;}
    for(int i=0;i<l*c-m*n;i++)
        M.pop_back();
    }
template<class T>M_taille MATRIX<T>::size_M() const
               {
               M_taille tmp;
               tmp.n=n;
               tmp.m=m;
               return tmp;
               }
template <class T>MATRIX<T> MATRIX<T>::fct(T t(T)) const
    {
    MATRIX<T> A=*this;
    for(int i=0;i<m*n;i++)
        A.M[i]=t(A.M[i]);
    return A;
    }
/*-----------fonction supplémentaire utile au calculs de matrices-------------*/
template<class T>std::vector<T> operator+(std::vector<T> x,std::vector<T>y)
               {
               std::vector<T> tmp;
               if ( x.size()==y.size() )
               for (int i=0;i<x.size();i++)
                   tmp.push_back(x[i]+y[i]);
               return tmp;
               }
template<class T>std::vector<T> operator-(std::vector<T> x,std::vector<T>y){return x+T(-1.)*y;}
template<class T>std::vector<T> operator*(T x,std::vector<T>y)
               {
               std::vector<T> tmp;
               for (int i=0;i<y.size();i++)
                   tmp.push_back(x*y[i]);
               return tmp;
               }
template<class T>std::vector<T> operator*(std::vector<T>y,T x)
               {
               return x*y;
               }
    // fonction transpose et inverse fonctionnant par copie volontairement
    // permettant de ne pas modifier A et de les utiliser avec une const MATRIX
template<class T>MATRIX<T> transpose(MATRIX<T> A)
        {return A.transpose();}
template<class T>MATRIX<T> inverse(MATRIX<T> A)
        {return A.inverse();}
template<class T>int TLU(const MATRIX<T>& A,MATRIX<T>& b,MATRIX<T>& c)
// fonction retournant une matrice triangulaire haute(c), et basse (b)
// tel que a=b×c
    {
    if (A.n!=A.m) return -1;
    //MATRIX<T> c(A.n),b(A.n);
    b=MATRIX<T>(A.n);c=b;
    T som=0;
    b(0,0)=A(0,0);
    c(0,0)=1;
    for (int j=1;j<A.n;j++)
        {
        c(0,j)=A(0,j)/c(0,0);
        b(j,0)=A(j,0)/b(0,0);
        }
    for (int k=1;k<A.n;k++)
    for (int kk=k;kk<A.n;kk++)
        {
        som=0;
        for (int l=0;l<kk;l++)
                som+=b(kk,l)*c(l,kk);
        b(kk,kk)=A(kk,kk)-som;
        c(kk,kk)=1;
        som=0;
        for (int l=0;l<k;l++)
                som+=b(k,l)*c(l,kk);
        c(k,kk)=(A(k,kk)-som)/b(kk,kk);
        b(kk,k)=(A(kk,k)-som)/c(k,k);
        }
    return 0;
    }
#if defined(_CPP_IOSTREAM) || defined(__IOSTREAM__)   
//template<class T>std::ostream& operator<<(std::ostream& s, MATRIX<T> A)
//        {
//        if(A.genre()==IS_SCALAR) {s<<A(0,0);return s;}
//        for (int i=0;i<A.n;i++)
//            {
//            s<<"\t|";
//            for (int j=0;j<A.m;j++)
//                s<<A(i,j)<<" | ";
//            s<<"\n";
//            }
//        return s;
//        }
//template<class T>MATRIX<T>& operator>>(std::istream& s, MATRIX<T>& A)
//        {
//        std::cout<<"Saisi de matrice:\nnombre de lignes\n";
//        s>>A.n;
//        std::cout<<"nombre de colonnes\n";
//        s>>A.m;
//        A=MATRIX<T>(A.n,A.m);
//        std::cout<<"Remplisage par ligne puis par colonne\n";
//        for (int i=0;i<A.n*A.m;i++)
//            s>>A.M[i];
//        //A.transpose();
//        return A;
//        }
std::ostream& operator<<(std::ostream& s, M_taille A)
        {
        s<<"| "<<A.n<<" | "<<A.m<<" | ";
        return s;
        }
#endif   // IOSTREAM
template <class T,typename T1> MATRIX<T>& operator>>(T1& s, MATRIX<T>& A)
        {
        __ENVOI("Saisi de matrice:\nnombre de lignes\n");
        s>>A.n;
        __ENVOI("nombre de colonnes\n");
        s>>A.m;
        A=MATRIX<T>(A.n,A.m);
        __ENVOI("Remplisage par ligne puis par colonne\n");
        for (int i=0;i<A.n*A.m;i++)
            s>>A.M[i];
        //A.transpose();
        return A;
        }
template <class T1,typename _T> T1& operator<<(T1& s, MATRIX<_T> A)
        {
        if(A.genre()==IS_SCALAR) {s<<A(0,0);return s;}
        for (int i=0;i<A.n;i++)
            {
            s<<"\t|";
            for (int j=0;j<A.m;j++)
                s<<A(i,j)<<" | ";
            s<<"\n";               
            }
        return s;
        }
template<class T>MATRIX<T> pow(const MATRIX<T>& r,const MATRIX<T>& e)// à généraliser pour les autres fonctions : sin cos exp ln etc...
    {
    MATRIX<T> tmp1=r;
    if ((r.genre()==IS_SCALAR) & (e.genre()==IS_SCALAR))
        return pow(r(0,0),e(0,0));
    else
        if(e.genre()==IS_SCALAR)
        if(e(0,0)>0)
                for(int i=1;i<e(0,0);i++)
                   tmp1=tmp1*r;
        else
                tmp1=inverse(pow(r,T(-1.0)*e));
    return tmp1;
    } 
template<class T>MATRIX<T> cos(const MATRIX<T>& r)// à généraliser pour les autres fonctions : sin cos exp ln etc...
    {
    MATRIX<T> tmp1=r;
    int n=tmp1.size_M().n,m=tmp1.size_M().m;
    for(int i=0;i<n;i++)
    for(int j=0;j<m;j++)
        tmp1(i,j)=cos(tmp1(i,j));
    return tmp1;
    } 
template<class T>MATRIX<T> log(const MATRIX<T>& r)// à généraliser pour les autres fonctions : sin cos exp ln etc...
    {
    MATRIX<T> tmp1=r;
    int n=tmp1.size_M().n,m=tmp1.size_M().m;
    for(int i=0;i<n;i++)
    for(int j=0;j<m;j++)
        tmp1(i,j)=log(tmp1(i,j));
    return tmp1;
    } 
template<class T>MATRIX<T> fabs(const MATRIX<T>& r)// à généraliser pour les autres fonctions : sin cos exp ln etc...
    {
    MATRIX<T> tmp1=r;
    int n=tmp1.size_M().n,m=tmp1.size_M().m;
    for(int i=0;i<n;i++)
    for(int j=0;j<m;j++)
        tmp1(i,j)=fabs(tmp1(i,j));
    return tmp1;
    } 
template<class T>MATRIX<T> acos(const MATRIX<T>& a){return a.fct(acos);}
template<class T>MATRIX<T> operator<(MATRIX<T> r,MATRIX<T> s)
// Compare les matrices. Tout d'abord le nombre d'éléments, si s a plus délément s que r
// Ensuite compare leur déterminant si elles sont carrées
// Pour finir leur ordre, ainsi 2×2>1×4
// le cas 1×4<4×1 renvoi zero car on ne peut comparer
    {
    MATRIX<T> tmp1(1);
    if(r.size()<s.size())
        {tmp1(0,0)=1;return tmp1;}
    else if(r.size()>s.size())
        {tmp1(0,0)=0;return tmp1;}
    else if((r.size_M().n==r.size_M().m)&(s.size_M().n==r.size_M().m))
        {tmp1(0,0)=(T)(r.det()<s.det());return tmp1;}
    int rn=r.size_M().n,rm=r.size_M().m,sn=s.size_M().n,sm=s.size_M().m;
    if(fmin(rn,rm)<fmin(sn,sm)){tmp1(0,0)=1;return tmp1;}
    tmp1(0,0)=0;
    return tmp1;
    } 
template<class T>MATRIX<T> operator>(MATRIX<T> r,MATRIX<T> s)
    {
    MATRIX<T> tmp1(1);
    if(r.size()>s.size())
        {tmp1(0,0)=1;return tmp1;}
    else if(r.size()<s.size())
        {tmp1(0,0)=0;return tmp1;}
    else if((r.size_M().n==r.size_M().m)&(s.size_M().n==r.size_M().m))
        {tmp1(0,0)=(T)(r.det()>s.det());return tmp1;}
    int rn=r.size_M().n,rm=r.size_M().m,sn=s.size_M().n,sm=s.size_M().m;
    if(fmin(rn,rm)>fmin(sn,sm)){tmp1(0,0)=1;return tmp1;}
    tmp1(0,0)=0;
    return tmp1;
    } 
//template<class T>MATRIX<T> operator==(MATRIX<T> r,MATRIX<T> s) //définie par défaut
//    {
//    MATRIX<T> tmp1(1);
//    int rn=r.size_M().n,rm=r.size_M().m,sn=s.size_M().n,sm=s.size_M().m;
//    if((rn!=sn)|(rm!=sm)){tmp1(0,0)=0;return tmp1;}
//    else
//    for(int i=0;i<rn;i++)
//    for(int j=0;j<rm;j++)
//        if(r(i,j)!=s(i,j)){tmp1(0,0)=0;return tmp1;}
//    tmp(0,0)=1;
//    return tmp1;
//    } 
#endif //ndef _matrice
/*int diagonalisation(MATRIX2,MATRIX2&,MATRIX2&,std::vector<double>&);
int diagonalisation(MATRIX2 A,MATRIX2& D,MATRIX2& P,std::vector<double>& Lambda)
    {
    double s,df;
    D=MATRIX2(A.n);P=MATRIX2(A.n);Lambda.clear();
    if (A.n!=A.m) {__ENVOI("imposible de diagonalise une matrice non carre!\n");exit(0);}
    double inter=0,tmp1,tmp2,tmp3;
    double pr=200;//precision
    for (int i=0;i<A.n*A.n;i++)inter+=fabs(A.M[i]);
    for (double l=-inter;l<inter;l+=inter/pr)
        / *{
        while(1)
                 {
                 df=((A-l*MATRIX2(A.n).Id()).det()-(A-(l-inter/(20*pr))*MATRIX2(A.n).Id()).det())/(inter/(20*pr));
                 s=l+(A-l*MATRIX2(A.n).Id()).det()/df;cout<<s-l<<endl;if(s<l) {l+=inter;continue;}
                 if (s-l!=0) l=s;
                 else break;if (!-inter<l) system("pause");
                 }
        Lambda.push_back(l);
        }//méthode de Newton(proto necessite d'avoir la dérivée du polynome caractéristique)* /
/ *        {
        tmp1=(A-l*MATRIX2(A.n).Id()).det();
        if (tmp1*tmp2<0 & l!=-inter)
                Lambda.push_back(l-inter*2.0/pr);
        if (tmp1==0)
                Lambda.push_back(l);
        tmp2=tmp1;
        }    //méthode perso, lente, possibilité de sauter des solutions, et impresice
    if (Lambda.size()!=A.n) {cout<<MATRIX2(Lambda);__ENVOI("probleme de nombre de VP\n");exit(0);}
    for (int i=0;i<A.n;i++)D(i,i)=Lambda[i];
    tmp2=0;for (int i=1;i<A.n;i++) {if (A(0,i)!=0)tmp2=1;}tmp3=tmp2;//traitement du cas où la 1ére coordonnée du VP est nulle.(proto)
    MATRIX2 AV(A.n-1),AVt;
    for(int i=0;i<A.n;i++)
        {
        for (int k=0;k<A.n-1;k++)for (int j=0;j<A.n-1;j++)AV(k,j)=A(k+1,j+1);
        AV=AV-Lambda[i]*MATRIX2(AV.n).Id();
        if(abs(A(0,0)-Lambda[i])<2*inter/pr) tmp3=1;else tmp3=tmp2;
        tmp1=AV.det();
        P(0,i)=tmp3;
        for(int j=1;j<A.n;j++)
                {
                AVt=AV;
                for (int k=0;k<AV.n;k++) AVt(k,j-1)=-A(k+1,0)/ **tmp3* /;
/ *                P(j,i)=AVt.det()/tmp1;
                }
        }
    return 0;
    }* /
/ * pour passer  n m et M en private, il faut definir les fonction qui les utilise en friend dans la declaration
 (si j'ai bien compris). et rajouter eventuellement quelques fonctions, comme resize pour aller avec la fonction 
 push_back, si l'on ne veut pas une matrice carré */