#ifndef _WIN32
#include "../include/lminmaxgd.h"
#else
#include "..\include\lminmaxgd.h"
using namespace std;
#endif
#include <time.h>




/////////////////////////////////////////////////////////////////////
/* Example :

*/
#include "..\src\gd.cpp"
#include "..\src\poly.cpp"
#include "..\src\serie.cpp"
#include "..\src\smatrix.cpp"
#include "..\src\tools.cpp"

smatrix removeline(smatrix & A, int  i);

int main()
{
try
{

   int n=6;
   int i=0,ok=0;
   smatrix   A(n,n); // state matrix 6x6


    A(0,2)=gd(2,3);
    A(1,3)=gd(2,4);
    A(2,0)=gd(1,5);
    A(2,1)=gd(2,6);
    A(3,1)=gd(2,3);
    A(4,2)=gd(1,7);
    A(4,3)=gd(1,4);
    A(5,3)=gd(0,14);

smatrix Ar,MN, M, N, Q;
int tab[n]; // a tab to memorize the removed lines
int num=0;
Q=A;
   for ( i=2;i<Q.getrow();i++) // start at 2 since we impose to keep the two first transtions x1 an x2
   {
 //cout<<"Q"<<Q<<endl;
 cout<<" i"<<i<<endl;
   N=removeline( Q,i);

   M=rfrac(Q,N);

  // cout<<"N"<<N<<endl; getchar();
//cout<<"N "<<N<<endl; // N composed of lines of A
  //  cout<<"M "<<M<<endl;
   // getchar();
   MN=otimes(M,N);

   if (MN==Q)
    {  if (num==0) tab[num]=i;
    else tab[num]=tab[num-1]+num;
      cout<<"N"<<N<<endl; getchar();
       cout<<" line removed (state) "<<tab[num]<<" ("<<tab[num]+1<<")"<<endl; // tab[num] is to memorize the removed state
       getchar();
        num++;
       Q=N;
       i--;
    }
   }
   N=Q;
   M=rfrac(A,N);
   MN=otimes(M,N);
   if(MN==A) cout<<" it is still ok, just to check that N and M are suitable, as expected"<<endl;

    Ar=otimes(N,M); // the dynamic of the equivalent system

    cout<<"A "<<A<<endl; // matrix of the system, it is equal to A=MN
    cout<<"Ar "<<Ar<<endl; // reduced matrix, it is equal to Ar=NM
    cout<<"N "<<N<<endl; // N composed of lines of A
    cout<<"M "<<M<<endl; // the transform from x to xr

     // checking properties;

     smatrix As, Ars;


     As=star(A);
     Ars=star(Ar);
     cout<<"Ars"<<Ars<<endl;

     // First property MAr*= A*M
     smatrix MArs, AsM;
     MArs=otimes(M,Ars);
     AsM=otimes(As,M);
     if(MArs==AsM)
     cout<<" Yes MAr*= A*M"<<endl;
     else
     cout<<" No MAr*= A*M"<<endl;


 // Second property NA*= Ar*N
     smatrix NAs, ArsN;
     NAs=otimes(N,As);
     ArsN=otimes(Ars,N);
     if(ArsN==NAs)
     cout<<" Yes Ar*N= NA*"<<endl;
     else
     cout<<" No Ar*N= NA*"<<endl;

// Third property A*=(Id + M Ar* N)

smatrix Id(n,n);
smatrix AMARN;
for(i=0;i<n;i++)
Id(i,i)=e;
AMARN=otimes(Ars,N);
AMARN=otimes(M,AMARN);
AMARN=oplus(Id,AMARN);
 if(AMARN==As)
     cout<<" Yes A*=(Id + M Ar* N)"<<endl;
     else
     cout<<" No A*=(Id + M Ar* N)"<<endl;

// Fourth property Ar*=(Id + N A* M)
smatrix Idr(Ar.getrow(),Ar.getcol());
smatrix ArNAM;
for(i=0;i<Ar.getrow();i++)
Idr(i,i)=e;
ArNAM=otimes(As,M);
ArNAM=otimes(N,ArNAM);
ArNAM=oplus(Idr,ArNAM);
 if(ArNAM==Ars)
     cout<<" Yes Ar*=(Id + N A* M)"<<endl;
     else
     cout<<" No Ar*=(Id + N A* M)"<<endl;
	return(0);
}

 catch(mem_limite l)
 {
	 cout<<"Exception : too many coefficent in polynom "<<l.memoire<<endl;
	 return(1);
 }

 catch(taille_incorrecte obj)
 { // 0 : r non causal
   // 1 : tentative d'accès à un element d'une matrice avec un indice incorrect
   // 2 : matrice de taille incompatible pour oplus, inf, otimes, rfrac, lfrac
   // 3 : etoile de matrice carré uniquement
	 cout<<"Exception  "<<obj.erreur<<endl;
	 return(1);
 }
}


smatrix removeline(smatrix & A, int nr)
{
    int nn=A.getrow();
    int i,j,k;
    smatrix N(A.getrow()-1,A.getcol());

    for (i=0,k=0;i<nn;i++) // we are removing the first line
    {
        if (i!=nr)
        {


        for(j=0;j<A.getcol();j++)
        {
            N(k,j)=A(i,j);

        }
            k++;
        }
    }



    return(N);

}
