#include <iostream> // cin, cout et cerr

using namespace std ;

inline void
fatal(const string & message)
{
	cerr << "Erreur fatale, message `" << message << "'." << endl ;
	exit(1) ;
}

class Matrice {

	const size_t lignes_ ;   // Nombre de lignes de la matrice.
	const size_t colonnes_ ; // Nombre de colonnes de la matrice.
	double * tableau ;       // Pointeur vers le début de l'emplacement alloué sur le tas.

	// Les deux méthodes suivantes sont `private', c'est-à-dire propres à la classe. On
	// ne peut pas les appeler depuis l'extérieur de la classe. Les éléments de la matrice
	// sont rangés en colonne comme en Fortran. Il serait ainsi possible d'utiliser des
	// sous-programmes de bibliothèques écrites en Fortran.
	// Référence (indices non contrôlés) à un élément de la matrice en lecture (référence
	// const à une matrice déclarée const).
	const double & imprudent(const size_t & i, const size_t & j) const {
		return tableau[j*lignes_ + i] ; }
	// Référence (indices non contrôlés) à un élément de la matrice en lecture (référence
	// non const à une matrice déclarée non const).
	double & imprudent(const size_t & i, const size_t & j) {
		return tableau[j*lignes_ + i] ; }

	public :
	// Constructeur normal.
	Matrice(const size_t & arg1, const size_t & arg2 = 1) :
		lignes_(arg1), colonnes_(arg2) {
		Matrice & A = *this ; // L'instance courante s'appelle `A'.
		if ( (A.lignes() == 0) || (A.colonnes() == 0) ) {
			fatal("Les arguments du constructeur normal sont invalides") ; }
		// Allocation sur le tas, à l'aide de l'opérateur new, d'un emplacement de taille
		// suffisante. L'opérateur new retourne 0 en cas d'échec.
		if ( (A.tableau = new double[A.lignes()*A.colonnes()]) == 0 ) {
			fatal("Échec de l'allocation") ; }
		// Initialisation à zéro de la matrice.
		for ( size_t i = 0 ; i < A.lignes() ; ++ i ) {
			for ( size_t j = 0 ; j < A.colonnes() ; ++ j ) {
				A.imprudent(i, j) = 0. ; } }
	}
	// Constructeur de copie.
	Matrice(const Matrice & B) : lignes_(B.lignes()), colonnes_(B.colonnes()) {
		Matrice & A = *this ; // L'instance courante s'appelle `A'.
		// Allocation sur le tas.
		if ( (tableau = new double[A.lignes()*A.colonnes()]) == 0 ) {
			fatal("Échec de l'allocation") ; }
		// Recopie de la matrice B.
		for ( size_t i = 0 ; i < A.lignes() ; ++ i ) {
			for ( size_t j = 0 ; j < A.colonnes() ; ++ j ) {
				A.imprudent(i, j) = B.imprudent(i, j) ; } }
	}
	// Destructeur.
	~Matrice() {
		if ( tableau == 0 ) { // Ceci ne devrait jamais arriver.
			fatal("Erreur interne à la bibliothèque, prévenir F. Legendre") ; }
		// La mémoire précédemment obtenue est retournée au système d'exploitation.
		delete [] tableau ;
	}
	// A.lignes() et A.colonnes() : nombre de lignes et de colonnes de A.
	size_t lignes() const { return lignes_ ; }
	size_t colonnes() const { return colonnes_ ; }

	// Opérateur = ; affectation : A = B.
	const Matrice & operator=(const Matrice & B) {
		Matrice & A = *this ; // L'instance courante s'appelle `A'.
		if ( (A.lignes() != B.lignes()) || (A.colonnes() != B.colonnes()) ) {
			cerr << "A.lignes() = " << A.lignes() << endl ;
			cerr << "A.colonnes() = " << A.colonnes() << endl ;
			cerr << "B.lignes() = " << B.lignes() << endl ;
			cerr << "B.colonnes() = " << B.colonnes() << endl ;
			fatal("Formats invalides pour l'opérateur `='") ; }
		for ( size_t i = 0 ; i < A.lignes() ; ++ i ) {
			for ( size_t j = 0 ; j < A.colonnes() ; ++ j ) {
				A.imprudent(i, j) = B.imprudent(i, j) ; } }
		return A ; }

	// Opérateur (..., ...) ; les indices, qui commencent en 0, sont contrôlés. C'est une
	// référence sur un double qui est retournée : on peut ainsi écrire `A(0, 0) = 12 ;'.
	double & operator()(const size_t & i, const size_t & j = 0) {
		Matrice & A = *this ; // L'instance courante s'appelle `A'.
		if ( (i >= A.lignes()) || (j >= A.colonnes()) ) {
			cerr << "i = " << i << " A.lignes() = " << A.lignes() << endl ;
			cerr << "j = " << j << " A.colonnes() = " << A.colonnes() << endl ;
			fatal("Les indices sont invalides") ; }
		return A.imprudent(i, j) ; }
	// Même méthode mais qui retourne une référence const pour une const Matrice.
	const double & operator()(const size_t & i, const size_t & j = 0) const {
		const Matrice & A = *this ; // L'instance courante s'appelle `A'.
		if ( (i >= A.lignes()) || (j >= A.colonnes()) ) {
			cerr << "i = " << i << " A.lignes() = " << A.lignes() << endl ;
			cerr << "j = " << j << " A.colonnes() = " << A.colonnes() << endl ;
			fatal("Les indices sont invalides") ; }
		return A.imprudent(i, j) ; }

	// Opérateur * ; multiplication matricielle : C = A*B.
	Matrice operator*(const Matrice & B) const {
		const Matrice & A = *this ; // L'instance courante s'appelle `A'.
		if (A.colonnes() != B.lignes()) {
			cerr << "A.colonnes() = " << A.colonnes() << endl ;
			cerr << "B.lignes() = " << B.lignes() << endl ;
			fatal("Formats invalides pour l'opérateur `*'") ; }
		// Création d'une matrice temporaire, C, pour enregistrer le résultat.
		Matrice C(A.lignes(), B.colonnes()) ;
		for ( size_t i = 0 ; i < A.lignes() ; ++ i ) {
			for ( size_t j = 0 ; j < B.colonnes() ; ++ j ) {
				double somme = 0. ;
				for ( size_t k = 0 ; k < A.colonnes() ; ++ k ) {
					somme += A.imprudent(i, k) * B.imprudent(k, j) ; }
				C.imprudent(i, j) = somme ; } }
		return C ; }

	// Opérateur % ; multiplication terme à terme : C = A%B.
	Matrice operator%(const Matrice & B) const {
		const Matrice & A = *this ; // L'instance courante s'appelle `A'.
		if ( (A.lignes() != B.lignes()) || (A.colonnes() != B.colonnes()) ) {
			cerr << "A.lignes() = " << A.lignes() << endl ;
			cerr << "A.colonnes() = " << A.colonnes() << endl ;
			cerr << "B.lignes() = " << B.lignes() << endl ;
			cerr << "B.colonnes() = " << B.colonnes() << endl ;
			fatal("Formats invalides pour l'opérateur `%'") ; }
		Matrice C(A.lignes(), A.colonnes()) ; // Temporaire pour enregistrer le résultat.
		for ( size_t i = 0 ; i < A.lignes() ; ++ i ) {
			for ( size_t j = 0 ; j < A.colonnes() ; ++ j ) {
				C.imprudent(i, j) = A.imprudent(i, j) * B.imprudent(i, j) ; } }
		return C ; }

	// Opérateur ! ; transposition : B = !A.
	Matrice operator!() const {
		const Matrice & A = *this ;           // L'instance courante s'appelle `A'.
		Matrice B(A.colonnes(), A.lignes()) ; // Temporaire pour enregistrer le résultat.
		for ( size_t i = 0 ; i < A.lignes() ; ++ i ) {
			for ( size_t j = 0 ; j < A.colonnes() ; ++ j ) {
				B.imprudent(j, i) = A.imprudent(i, j) ; } }
		return B ; }

} ;

// Représentation externe d'une matrice sur un flux de sortie. La matrice est affichée
// par ligne alors qu'elle est enregistrée par colonne. Le caractère `;' sépare les dif-
// férentes lignes. La suite de nombres est encadrée par les caractères `{' et `}'.
ostream & operator<<(ostream & arg1, const Matrice & A)
{
	arg1 << "{ " ;
	for ( size_t i = 0 ; i < A.lignes() ; ++ i ) {
		for ( size_t j = 0 ; j < A.colonnes() ; ++ j ) {
			arg1 << A(i, j) << ' ' ;  }
		if ( i != (A.lignes()-1) ) {
			arg1 << "; " ; } }
	arg1 << '}' ;
	return arg1 ;
}

Matrice
inv(const Matrice & A)
{
	if ( A.lignes() != A.colonnes() ) {
		cerr << " A.lignes() = " << A.lignes() << endl ;
		cerr << " A.colonnes() = " << A.colonnes() << endl ;
		fatal("La matrice n'est pas carrée dans `inv(...)'") ; }
	Matrice B(A.lignes(), A.lignes()) ;
	double det ;
	switch ( A.lignes() ) {

		case 1 : B(0, 0) = 1 / A(0, 0) ;
		break ;

		case 2 : det = A(0, 0)*A(1, 1) - A(1, 0)*A(0, 1) ;
		B(0, 0) = A(1, 1) / det ;
		B(1, 1) = A(0, 0) / det ;
		B(1, 0) = -A(1, 0) / det ;
		B(0, 1) = -A(0, 1) / det ;
		break ;

		default : cout << "L'inversion d'une matrice de taille > 2 n'est pas encore "
			"implémentée, se plaindre auprès de F. Legendre" << endl ;
		break ;
	}
	return B ;
}


int
main()
{
	Matrice A(2, 3) ; // Constructeur normal.
	// Pour montrer que les éléments de la matrice sont initialisés à 0.
	cout << "A = " << A << endl ;
	// Affectation individuelle des éléments de la matrice.
	A(0, 0) = 1 ; A(1, 0) = 2 ; A(0, 1) = 3 ; A(1, 1) = 4 ;  A(0, 2) = 5 ; A(1, 2) = 6 ; 
	// Impression de contrôle.
	cout << "A = " << A << endl ;

	Matrice B(A) ; // Constructeur de copie.
	cout << "B = " << B << endl ;
	cout << "A%B = " << A%B << endl ; // Test de l'opérateur `%'.
	cout << "!A = " << !A << endl ;   // Test de l'opérateur `!'.
	Matrice C(!A) ;                   // Constructeur de copie à partir d'un temporaire.
	cout << "A*C = " << A*C << endl ; // Test de l'opérateur `*'.

	// Test de l'opérateur d'affectation.
	Matrice D(2, 3) ;
	Matrice E(2, 3) ;
	E = D = A%B ;
	cout << "D = " << D << endl ;
	cout << "E = " << E << endl ;

	// Test de l'inversion des matrices.
	Matrice F(A*!A) ;
	cout << "F*inv(F) = " << F*inv(F) << endl ;
	cout << "inv(F)*F = " << inv(F)*F << endl ;

	// Programmation des moindres carrés ordinaires.
	Matrice X(5, 2) ;
	X(0, 0) = -1 ; X(0, 1) = 1 ;
	X(1, 0) = 3 ; X(1, 1) = 1 ;
	X(2, 0) = -4 ; X(2, 1) = 1 ;
	X(3, 0) = 2 ; X(3, 1) = 1 ;
	X(4, 0) = -1 ; X(4, 1) = 1 ;
	Matrice y(5) ;
	for ( size_t i = 0 ; i < X.lignes() ; ++ i ) {
		y(i) = 3*X(i, 0) + 4 + i%2-.5 ; } // `i%2-.5' est égal à -.5, .5, -.5, ...
	Matrice a_chapeau( inv(!X*X) * (!X*y) ) ;
	cout << "â = " << a_chapeau << endl ;

	return 0 ; // Retour au système d'exploitation avec le code `0'.
}