// 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;
}
syntax highlighted by Code2HTML, v. 0.9.1