// Rascal, the Advanced Scientific CALculator
// Copyright (C) 2001, Sebastian Ritterbusch (Rascal@Ritterbusch.de)
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation; either version 2
// of the License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more detauls.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
//

#include <complex.hpp>
#include <rpeval.hpp> // Polynomial evaluation

valuematrix rpolyeval(const valuepoly &p,double t)
{
   RPolynomial rp(p.N);
   string err;
   
   int i;
   for(i=0;i<p.N;i++)
   {
      if(p.a[i].isDOUBLE())
         rp[i]=p.a[i].asDOUBLE();
      else if((p.a[i]*1.0).isDOUBLE())
         rp[i]=(p.a[i]*1.0).asDOUBLE();
      else
      {
         rp[i]=0;
         err=" Polynomial had non-real coefficients!";
      }
   }
   
   real     rt(t),ry;
   interval iy;
   int      No,Err;
   
   RPolyEval(rp,rt,ry,iy,No,Err);
   
   valuematrix erg(4,1);
   
   erg.a[0][0]=_double(ry);
   if(diam(iy)>0)
      erg.a[1][0]=iy;
   else
      erg.a[1][0]=_double(mid(iy));
   erg.a[2][0]=No;
   erg.a[3][0]=string(RPolyEvalErrMsg(Err)+err);
   
   return erg;
}

#include <cpzero.hpp>

valuematrix cpolyzero(valuepoly p,valuecomplex t)
{
   CPolynomial cp(p.N);
   string err;
   
   int i;
   for(i=0;i<p.N;i++)
   {
      if((p.a[i]*1.0).isDOUBLE())
         cp[i]=(p.a[i]*1.0).asDOUBLE();
      else if(p.a[i].isCOMPLEX())
      {
         valuecomplex q(p.a[i].asCOMPLEX());
         if((q.real()*1.0).isDOUBLE() && (q.imag()*1.0).isDOUBLE())
            cp[i]=cxsc::complex((q.real()*1.0).asDOUBLE(),(q.imag()*1.0).asDOUBLE());
         else
         {
            cp[i]=0;
            err=" Not all coefficients were usual complex coefficients!";
         }
      } else
      {
         cp[i]=0;
         err=" Not all coefficients were complex or real!";
      }
   }
   
   cxsc::complex ct;
   
   ct=0;
   
   if((t.imag()*1.0).isDOUBLE() && (t.real()*1.0).isDOUBLE())
      ct=cxsc::complex((t.real()*1.0).asDOUBLE(),(t.imag()*1.0).asDOUBLE());
   else
      err=" Approximation point is not complex or real!";

   CIPolynomial ci(p.N);
   int Err;
   cxsc::cinterval zz;

   CPolyZero(cp,ct,ci,zz,Err);
   
   valuematrix erg(3,1);
   
   erg.a[0][0]=valuecomplex(Re(zz),Im(zz));
   
   for(i=0;i<p.N;i++)
      if(diam(Re(ci[i]))>0 || diam(Im(ci[i]))>0)
         p.a[i]=valuecomplex(Re(ci[i]),Im(ci[i]));
      else
         p.a[i]=valuecomplex(_double(mid(Re(ci[i]))),_double(mid(Im(ci[i]))));

   erg.a[1][0]=p;
   
   erg.a[2][0]=CPolyZeroErrMsg(Err)+err;

   return erg;
}

#include <linsys.hpp>

valuematrix linsolve(const valuematrix & a,const valuematrix & b)
{
   string err;
   int    n=a.N;
   
   if(a.N!=a.M)
   {
      err=" Argument matrix not square matrix!";
      if(a.M<n)
         n=a.M;
      if(b.N>b.M && b.N<n)
         n=b.N;
      if(b.M>b.N && b.M<n)
         n=b.M;
   }
   
   if(n>b.N && n>b.M)
   {
      err=" Argument vector too small!";
      if(b.N>b.M && b.N<n)
         n=b.N;
      if(b.M>b.N && b.M<n)
         n=b.M;
   }
   
   rmatrix  A(n,n);
   rvector  B(n);
   
   int i,j;
   value x;
   for(i=0;i<n;i++)
   {
      if(b.M>b.N)
         x=b.a[0][i];
      else
         x=b.a[i][0];

      if((x*1.0).isDOUBLE())
         B[i+1]=(x*1.0).asDOUBLE();
      else
      {
         B[i+1]=0; 
         err=" Not all vector entries are real!";
      }      
      
      for(j=0;j<n;j++)
      {
         x=a.a[i][j];
         if((x*1.0).isDOUBLE())
            A[i+1][j+1]=(x*1.0).asDOUBLE();
         else
         {
            A[i+1][j+1]=0; 
            err=" Not all matrix entries are real!";
         } 
      }
   }  
   ivector X(n);
   int     Err;
   real    Cond;
   
   
   LinSolve(A,B,X,Cond,Err);
   

   valuematrix erg(3,1);
   
   valuematrix ev(n,1);
   
   for(i=0;i<n;i++)
   {
      if(diam(X[i+1])>0)
         ev.a[i][0]=X[i+1];
      else
         ev.a[i][0]=_double(mid(X[i+1]));
   }
   erg.a[0][0]=ev;
   
   erg.a[1][0]=_double(Cond);
   
   erg.a[2][0]=LinSolveErrMsg(Err)+err;

   return erg;
}

valuematrix boothroyd(int n)
{
   valuematrix A(n,n);
   int i,j;
   for(i=0;i<n;i++)
      for(j=0;j<n;j++)
      {
         value I(i+1),J(j+1),N(n);
         A.a[i][j]=fac(N+I-1)/fac(I-1)/N/fac(N-J)/fac(J-1)*N/(I+J-1);
      }
   return A;
}

#include <rev_simp.hpp>

valuematrix revsimplex(const valuematrix &a)
{
   valuematrix iA(0,0),ib(0,0),ic(0,0);
   string err;
   
   if((a.N>=a.M && (a.N<3 || !a.a[0][0].isMATRIX() || !a.a[1][0].isMATRIX() || !a.a[2][0].isMATRIX()))
    ||(a.N<a.M && (a.M<3 || !a.a[0][0].isMATRIX() || !a.a[0][1].isMATRIX() || !a.a[0][2].isMATRIX())))
   {
      cerr << "Usage: RevSimplex(A,b,c)" << endl;
      return valuematrix(0,0);
   }   
   if(a.N>=a.M)
   {
      iA=a.a[0][0].asMATRIX();ib=a.a[1][0].asMATRIX();ic=a.a[2][0].asMATRIX();
   } else
   {
      iA=a.a[0][0].asMATRIX();ib=a.a[0][1].asMATRIX();ic=a.a[0][2].asMATRIX();
   }
   int n=iA.N;
   int m=iA.M;
   
   if(ib.M>ib.N)
      ib=transpose(ib);
   if(ic.M>ic.N)
      ic=transpose(ic);
      
   if(ib.N<n)
   {
      err=" Vector b too small!";
      n=ib.N;
   }
   if(ic.N<m)
   {
      err=" Vector c too small!";
      m=ic.M;
   }
   
   rmatrix A(n,m);
   rvector b(n),c(m);
   
   int i,j;
   for(i=0;i<n;i++)
      for(j=0;j<m;j++)
      {
         value x=iA.a[i][j];
         if((x*1.0).isDOUBLE())
            A[i+1][j+1]=(x*1.0).asDOUBLE();
         else
         {
            err=" Values in matrix are not real.";
            A[i+1][j+1]=0;
         }
      }
   for(i=0;i<m;i++)
   {
      value x=ic.a[i][0];
      if((x*1.0).isDOUBLE())
         c[i+1]=(x*1.0).asDOUBLE();
      else
      {
         err=" Values in vector c are not real.";
         c[i+1]=0;
      }
   }
   for(i=0;i<n;i++)
   {
      value x=ib.a[i][0];
      if((x*1.0).isDOUBLE())
         b[i+1]=(x*1.0).asDOUBLE();
      else
      {
         err=" Values in vector b are not real.";
         b[i+1]=0;
      }
   }
   
   rvector x;
   intvector v;
   real z;
   int Err;
   
   RevSimplex(A,b,c,x,v,z,Err);
   
   valuematrix erg(4,1);
   
   valuematrix oX(m,1);
   valuematrix oV(n,1);
   
   for(i=0;i<m;i++)
      oX.a[i][0]=_double(x[i+1]);
      
   erg.a[0][0]=oX;
     
   for(j=0;j<n;j++)
      oV.a[j][0]=(v[j+1]);
      
   erg.a[1][0]=oV;
   
   erg.a[2][0]=_double(z);
   
   erg.a[3][0]=RevSimplexErrMsg(Err)+err;

   return erg;
}

#include <lop.hpp>

valuematrix linopt(const valuematrix &a)
{
   valuematrix iA(0,0),ib(0,0),ic(0,0),iBStart(0,0);
   string err;
   
   if((a.N>=a.M && (a.N<4 || !a.a[0][0].isMATRIX() || !a.a[1][0].isMATRIX() || !a.a[2][0].isMATRIX() || !a.a[3][0].isMATRIX()))
    ||(a.N<a.M && (a.M<4 || !a.a[0][0].isMATRIX() || !a.a[0][1].isMATRIX() || !a.a[0][2].isMATRIX() || !a.a[0][3].isMATRIX())))
   {
      cerr << "Usage: LinOpt(A,b,c,BStart)" << endl;
      return valuematrix(0,0);
   }   
   if(a.N>=a.M)
   {
      iA=a.a[0][0].asMATRIX();ib=a.a[1][0].asMATRIX();ic=a.a[2][0].asMATRIX();
      iBStart=a.a[3][0].asMATRIX();
   } else
   {
      iA=a.a[0][0].asMATRIX();ib=a.a[0][1].asMATRIX();ic=a.a[0][2].asMATRIX();
      iBStart=a.a[0][3].asMATRIX();
   }
   int n=iA.N;
   int m=iA.M;
   
   if(ib.M>ib.N)
      ib=transpose(ib);
   if(ic.M>ic.N)
      ic=transpose(ic);
   if(iBStart.M>iBStart.N)
      iBStart=transpose(iBStart);
      
   if(ib.N<n)
   {
      err=" Vector b too small!";
      n=ib.N;
   }
   if(iBStart.N<n)
   {
      err=" Vector BStart too small!";
   }
   if(ic.N<m)
   {
      err=" Vector c too small!";
      m=ic.M;
   }
   
   rmatrix A(n,m);
   rvector b(n),c(m);
   intvector BStart(n);
   
   int i,j;
   for(i=0;i<n;i++)
      for(j=0;j<m;j++)
      {
         value x=iA.a[i][j];
         if((x*1.0).isDOUBLE())
            A[i+1][j+1]=(x*1.0).asDOUBLE();
         else
         {
            err=" Values in matrix are not real.";
            A[i+1][j+1]=0;
         }
      }
   for(i=0;i<m;i++)
   {
      value x=ic.a[i][0];
      if((x*1.0).isDOUBLE())
         c[i+1]=(x*1.0).asDOUBLE();
      else
      {
         err=" Values in vector c are not real.";
         c[i+1]=0;
      }
   }
   for(i=0;i<n;i++)
   {
      value x=ib.a[i][0];
      if((x*1.0).isDOUBLE())
         b[i+1]=(x*1.0).asDOUBLE();
      else
      {
         err=" Values in vector b are not real.";
         b[i+1]=0;
      }
      x=toInteger(iBStart.a[i][0]);
      if(x.isINTEGER())
         BStart[i+1]=x.asINTEGER();
      else
      {
         BStart[i+1]=0;
         err=" Values in vector BStart are not integer.";
      }
   }
   
   interval  z;
   intmatrix V;
   imatrix   X;
   int       No;
   int       Err;
   
   LinOpt(A,b,c,BStart,z,V,X,No,Err);
   
   valuematrix erg(5,1);
   
   if(diam(z)==0)
      erg.a[0][0]=_double(mid(z));
   else
      erg.a[0][0]=z;
   
   valuematrix oV(No,n);
   
   for(i=0;i<No;i++)
      for(j=0;j<n;j++)
         oV.a[i][j]=V[i+1][j+1];
   
   erg.a[1][0]=oV;
   
   valuematrix oX(No,m);
   
   for(i=0;i<No;i++)
      for(j=0;j<m;j++)
         oX.a[i][j]=diam(X[i+1][j+1])==0?value(_double(mid(X[i+1][j+1]))):value(X[i+1][j+1]);
         
   erg.a[2][0]=oX;
   
   erg.a[3][0]=No;
   
   erg.a[4][0]=LinOptErrMsg(Err)+err;   
   
   return erg;
}