#ifndef __MATLIB_H
#define __MATLIB_H


/*
 *    MATLIB header file
 *
 *    Copyright (c) Yaron Keren 1994-1996. All rights reserved.
 */

#include <stddef.h>
#include <math.h>

#define M_MATRIX 0
#define M_TEXT	 1

union i_o_t {
  int i;
  int o;
};

class Mr;

class M {
  char *name;
  int type;
  int m;
  int n;
  int p;
  double *pr;
  double *pi;

  void allocate_double();
  void delete_double();
public:
  M();
  M(int nrows,int ncols);
  M(int nrows,int ncols, char *mname);
  M(char *mname, double valr, double vali=0);
  M(char *mname, char *text);

  M(double val);
  M(char *text);
  M(const M& src);
  M(const Mr& src);
  ~M();

  void operator =(double val);
  void operator =(char *text);
  void operator =(const M& src);
  void operator =(const Mr& src);
  void deepcopy(const M& src);

  inline int rows() const { return(m); }
  inline int cols() const { return(n); }
  inline int size() const { return(m*n); }
  int gettype() const { return(type); }
  void settype(int t) { type=t; }
  char *getname() const { return(name); }
  int dirty() const;

  inline double *addr() const { return(pr); }
  inline double *addr(int row,int col) const { return(pr+row-1+(col-1)*m); }
  inline double *addi() const { return(pi); }
  inline double *addi(int row,int col) const { return(pi+row-1+(col-1)*m); }

  double& r() const;
  double& r(double idx) const;
  double& r(double row, double col) const;
  double& i() const;
  double& i(double idx) const;
  double& i(double row, double col) const;

  Mr operator ()(M src);
  Mr operator ()(M srcr, M srcc);

  void vwcopy(const M& src, const M& v);
  void vwcopy(const M& src, const M& v, const M& w);

  friend void rowcol(M& x);
};


class Mr {
  friend class M;

  int indexed;

  M& data;
  M idxr;
  M idxc;

public:
  Mr(M& src);
  Mr(M& src, const M& index);
  Mr(M& src, const M& indexc, const M& indexr);

  ~Mr();

  int rows();
  int cols();
  int size();
  int gettype();
  void settype(int t);
  char *getname();

  double *addr();
  double *addr(int row,int col);
  double *addi();
  double *addi(int row,int col);

  void operator =(double val);
  void operator =(const M& src);
  void operator =(const Mr& src);

  Mr operator ()(M& src);
  Mr operator ()(M& srcr, M& srcc);
};

M operator <(M x, M y);
M operator >(M x, M y);
M operator <=(M x, M y);
M operator >=(M x, M y);
M operator ==(M x, M y);
M operator !=(M x, M y);
M operator ||(M x, M y);
M operator &&(M x, M y);
M operator +(M x, M y);
M operator -(M x, M y);
M dot_add(M x, M y);
M dot_sub(M x, M y);
M dot_mul(M x,M y);
M dot_div(M x,M y);
M dot_invdiv(M x,M y);
M dot_pow(M x,M y);
M atan2(M x,M y);
M rem(M x,M y);
M min(M x,M y);
M min0(M x,M y);
M max(M x,M y);
M pow2(M f, M e);
M xorM(M x, M y);

void operator +=(M& x, M y);
void operator -=(M& x, M y);
void operator *=(M& x, M y);
void operator /=(M& x, M y);
void operator ++(M x);
void operator --(M x);

/*
M operator <(M x, double y);
M operator >(M x, double y);
M operator <=(M x, double y);
M operator >=(M x, double y);
M operator ==(M x, double y);
M operator !=(M x, double y);
M operator !=(M x, double y);
M operator ||(M x, double y);
M operator &&(M x, double y);
M operator +(M x, double y);
M operator -(M x, double y);
M dot_add(M x, double y);
M dot_sub(M x, double y);
M dot_mul(M x,double y);
M dot_div(M x,double y);
M dot_invdiv(M x,double y);
M dot_pow(M x,double y);

M operator <(double x, M y);
M operator >(double x, M y);
M operator <=(double x, M y);
M operator >=(double x, M y);
M operator ==(double x, M y);
M operator !=(double x, M y);
M operator !=(double x, M y);
M operator ||(double x, M y);
M operator &&(double x, M y);
M operator +(double x, M y);
M operator -(double x, M y);
M dot_add(double x, M y);
M dot_sub(double x, M y);
M dot_mul(double x,M y);
M dot_div(double x,M y);
M dot_invdiv(double x,M y);
M dot_pow(double x,M y);

M operator *(M x, double y);
M operator *(double x, M y);
*/

M all(M x);
M any(M x);
M max(M x);
M min(M x);
M sum(M x);
M sumabs(M x);
M prod(M x);

M colon(double j, double i, double k);
M colon(M j_Mat, M i_Mat, M k_Mat);

extern M nargin, nargout;
extern int nargin_set, nargout_set;
extern M TICTOC, eps, version, Inf, NaN, nop_M;
extern M i,j,pi,zero_M,one_M,two_M,junk_M, sixpack_M, lasterr, mfilename;
extern M stdinM, stdoutM, stderrM;
extern M s_M, l_M, c_p;
extern char *s_s, *l_s;
extern double s_d, l_d;
extern double *s_Mp, *l_Mp, *x_Mp, *c_pp;
extern i_o_t i_o;

typedef M func1M(M);
typedef M func2M(M, M);
typedef M func3M(M, M, M);
typedef M func4M(M, M, M, M);

M junk_func1M(M);
M junk_func2M(M,M);
M junk_func3M(M,M,M);
M junk_func4M(M,M,M,M);

class D {
  double d;
  char *name;
public:
  D(double v) { name=NULL; d=v; }
  D(int, int, char *dname) { name=dname; d=0; }
  D(const M& src);
  operator double() const { return d; }
  char* getname() const { return name; }

  double operator +=(double v) { d+=v; return d; }
  double operator -=(double v) { d-=v; return d; }
  double operator *=(double v) { d*=v; return d; }
  double operator /=(double v) { d/=v; return d; }
  void operator ++() { d++; }
  void operator --() { d--; }
};

/*
D operator <(D x, D y);
D operator >(D x, D y);
D operator <=(D x, D y);
D operator >=(D x, D y);
D operator ==(D x, D y);
D operator !=(D x, D y);
D operator ||(D x, D y);
D operator &&(D x, D y);
D operator +(D x, D y);
D operator -(D x, D y);
D operator *(D x, D y);
D operator /(D x, D y);
D invdiv(D x, D y);
D pow(D x, D y);
*/
M diag(M x);
M diag(M x, M kk);

void memreport();
double *newdouble(int size);
char *newchar(int size);
char **newcharptr(int size);
int *newint(int size);
long *newlong(int size);
void deletedouble(double *p);
void deletechar(char *p);
void deletecharptr(char **p);
void deleteint(int *p);
void deletelong(long *p);

typedef double	*double_p;
typedef int	*int_p;

void copy_double(double *dest, double *src, int n);
void fill_double(double *dest, int n, double val);
void zero_double(double *dest, int n);
void zero_char(char *dest, int n);

M error(char *s);
M error(char *s, char *name1);
M error(char *s, char *name1, char *name2);
M error(char *s, char *name1, char *name2, char *name3);
M error(char *s, char *name1, char *name2, char *name3, char *name4);
M error(char *s, char *name1, char *name2, char *name3, char *name4, char *name5);
M warning(char *s);
M warning(char *s, char *name1);
M warning(char *s, char *name1, char *name2);

M eye(M x);
M eye(M x, M y);
M ones(double m);
M ones(double m, double n);
M ones(M x);
M ones(M x, M y);
M randM();
M randM(M x);
M randM(M x, M y);
M randn();
M randn(M x);
M randn(M x, M y);
M zeros(double m);
M zeros(double m, double n);
M zeros(M x);
M zeros(M x, M y);

M fft(M x);
M fft(M x, M N);
M ifft(M x);
M ifft(M x, M N);
M dft(M x);


M fcloseM(M fid);
M feofM(M fid);
M ferrorM(M fid, M clear=zero_M);
M ferrorM(M fid, i_o_t, Mr msg, Mr errnum);
M ferrorM(M fid, M clear, i_o_t, Mr msg, Mr errnum);
M fopenM(M filename, M permission="r", M machine="n");
M fopenM(M filename, i_o_t, Mr fid, Mr msg);
M fopenM(M filename, M permission, i_o_t, Mr fid, Mr msg);
M fopenM(M filename, M permission, M machine, i_o_t, Mr fid, Mr msg);
M fprintfM(M string);
M fprintfM(M fid, M format, M a=l_M, M b=l_M, M c=l_M, M d=l_M, M e=l_M, M f=l_M);
M printfM(M format, M a=l_M, M b=l_M, M c=l_M, M d=l_M, M e=l_M, M f=l_M);
M fread(M fid, M size=Inf, M precision="uchar", M skip=zero_M);
M fread(M fid, i_o_t, Mr A, Mr count);
M fread(M fid, M size, i_o_t, Mr A, Mr count);
M fread(M fid, M size, M precision, i_o_t, Mr A, Mr count);
M fread(M fid, M size, M precision, M skip, i_o_t, Mr A, Mr count);
M frewind(M fid);
M fscanfM(M fid, M format);
M fscanfM(M fid, M format, M size);
M fscanfM(M fid, M format, i_o_t, Mr A, Mr count);
M fscanfM(M fid, M format, M size, i_o_t, Mr A, Mr count);
M fscanfM(M fid, M format, M size, i_o_t, Mr A, Mr count, Mr errmsg);
M fseek(M fid, M offset, M origin);
M ftell(M fid);
M fwrite(M fid, M A, M precision="uchar", M skip=zero_M);

M ltitr(M a, M b, M u, M x0);
M ltifr(M a, M b, M s);
M filter(M b, M a, M x);
M filter(M b, M a, M x, M zi);
M filter(M b, M a, M x, i_o_t, Mr y, Mr zf);
M filter(M b, M a, M x, M zi, i_o_t, Mr y, Mr zf);

M axis();
M axis(M a);
M contour(M Z);
M contour(M X, M Y, M Z);
M contour(M Z, M N);
M contour(M X, M Y, M Z, M N);
M drawnow();
M errorbar(M y, M l);
M errorbar(M x, M y, M l);
M errorbar(M x, M y, M l, M u);
M grid();
M grid(M a);
M hold(M a);
M ishold();

M loglog(M y);
M loglog(M a1	 ,  M a2    ,  M a3=l_M,  M a4=l_M,  M a5=l_M,
	 M a6=l_M,  M a7=l_M,  M a8=l_M,  M a9=l_M,  M a10=l_M,
	 M a11=l_M, M a12=l_M, M a13=l_M, M a14=l_M, M a15=l_M);
M mesh(M Z);
M mesh(M x, M y, M z);
M meshc(M Z);
M meshc(M x, M y, M z);
M meshdom(M x, M y, i_o_t, Mr X, Mr Y);
M plot(M y);
M plot(M a1    ,  M a2	  ,  M a3=l_M,	M a4=l_M,  M a5=l_M,
       M a6=l_M,  M a7=l_M,  M a8=l_M,	M a9=l_M,  M a10=l_M,
       M a11=l_M, M a12=l_M, M a13=l_M, M a14=l_M, M a15=l_M);
M polar(M y);
M polar(M a1	,  M a2    ,  M a3=l_M,  M a4=l_M,  M a5=l_M,
	M a6=l_M,  M a7=l_M,  M a8=l_M,  M a9=l_M,  M a10=l_M,
	M a11=l_M, M a12=l_M, M a13=l_M, M a14=l_M, M a15=l_M);

M semilogx(M y);
M semilogx(M a1    ,  M a2    ,  M a3=l_M,  M a4=l_M,  M a5=l_M,
	  M a6=l_M,   M a7=l_M,  M a8=l_M,  M a9=l_M,  M a10=l_M,
	  M a11=l_M,  M a12=l_M, M a13=l_M, M a14=l_M, M a15=l_M);
M semilogy(M y);
M semilogy(M a1    ,  M a2    ,  M a3=l_M,  M a4=l_M,  M a5=l_M,
	  M a6=l_M,   M a7=l_M,  M a8=l_M,  M a9=l_M,  M a10=l_M,
	  M a11=l_M,  M a12=l_M, M a13=l_M, M a14=l_M, M a15=l_M);

M quiver(M x);
M quiver(M x, M y);
M quiver(M x, M y, M s);
M text(M x, M y, M a);
M title(M a,...);
M shg();
M view(M a);
M view(M az, M el);
M xlabel(M a,...);
M ylabel(M a,...);
M zlabel(M a,...);

M bmpwrite(M img, M map, M filename);
M image(M img);
M colormap(M map);
M black(M l=64);
M gray(M l=64);
M grey(M l=64);
M white(M l=64);

M axes(...);
M cla(...);
M close(...);
M clf(...);
M figure(...);
M gca(...);
M gcf(...);
M get(...);
M line(...);
M newplot(...);
M reset(...);
M set(...);
M subplot(...);
M surf(...);
M surface(...);
M waitforbuttonpress(...);

M wavwrite(M wavedata, M rate, M filename);
M soundM(M y, M rate=8192.0);

int isempty(M x);
M   isglobal(M);
int isreal(M x);
int issame(M x, M y);
M   issparse(M);
int isstr(M x);
int istrue(M x);
int istrue(double x);

M balance(M a);
M balance(M a, i_o_t, Mr t, Mr b);
M chol(M b);
M chol(M x, i_o_t, Mr R, Mr p);
M cond(M x);
M det(M x);
M dongarra(M x);
M eig(M x);
M eig(M x, i_o_t, Mr v, Mr d);
M eig(M a, M b);
M eig(M a, M b, i_o_t, Mr v, Mr d);
M expm(M A);
M hess(M a);
M hess(M a, i_o_t, Mr p, Mr h);
M inv(M x);
M invdiv(M a, M b);
M lu(M x);
M lu(M x, i_o_t, Mr l, Mr u);
M lu(M x, i_o_t, Mr l, Mr u, Mr p);
M norm(M x, M P=two_M);
M operator /(M x, M y);
M powM(M x,M y);
M qr(M x);
M qr(M x, M);
M qr(M x, i_o_t, Mr q, Mr r);
M qr(M x, M, i_o_t, Mr q, Mr r);
M qr(M x, i_o_t, Mr q, Mr r, Mr e);
M qr(M x, M, i_o_t, Mr q, Mr r, Mr e);
M qz(M a, M b, i_o_t, Mr aa, Mr bb);
M qz(M a, M b, i_o_t, Mr aa, Mr bb, Mr q);
M qz(M a, M b, i_o_t, Mr aa, Mr bb, Mr q, Mr z);
M qz(M a, M b, i_o_t, Mr aa, Mr bb, Mr q, Mr z, Mr v);
M rank(M x);
M rank(M x, M tol);
M rcond(M x);
M schur(M x);
M schur(M x, i_o_t, Mr u, Mr t);
M svd(M x);
M svd(M x, M z);
M svd(M x, i_o_t, Mr U, Mr S, Mr V);
M svd(M x, M z, i_o_t, Mr U, Mr S, Mr V);


M hilb(M n);
M invhilb(M n);
M magic(M n);
M pascalM(M n, M r=zero_M);


M clear(Mr x);
M computer(i_o_t, Mr c, Mr maxsize);
M computer();
M cputime();
M diff(M x, M order=one_M);
M echo(char* st);
M etime(M t1, M t0);
M exit();
M eval(...);
M exist(M x);
M feval(...);
M find(M x);
M find(M x, i_o_t, Mr I, Mr J);
M find(M x, i_o_t, Mr I, Mr J, Mr v);
int forsize(M& x);
void forelem(M& var, M& x, int index);
M full(M x);
M flops(M x);
M flops();
M getenv(M x);
M isieee();
M joinrows(M a, M b);
M joincols(M a, M b);
M linspace(M x1, M x2, M n=100);
M logspace(M x1, M x2, M n=50);
M matlabpath();
M matlabpath(M x);
M mean(M x);
M nnz(M x);
M nonzeros(M x);
M quit();
M setstr(M x);
M sparse(...);
M spdiags(...);
M spparms(...);
M tic();
M toc();
M then();

#ifdef MATLAB_MEX_FILE
extern "C" {
}

M Matrix2M(const Matrix& mat);
void M2Matrix(const M& X, Matrix & mat);

#endif // MATLAB_MEX_FILE
M min(M x, i_o_t, Mr y, Mr i);
M max(M x, i_o_t, Mr y, Mr i);

#include <iostream.h>

extern istream *matlibin_p;
extern ostream *matlibout_p, *matliberr_p;

#define matlibin (*matlibin_p)
#define matlibout (*matlibout_p)
#define matliberr (*matliberr_p)

M operator *(M x, M y);


M cos(M x);
M sin(M x);
M tan(M x);
M sec(M x);
M csc(M x);
M cot(M x);

M acos(M x);
M asin(M x);
M atan(M x);
M asec(M x);
M acsc(M x);
M acot(M x);

M cosh(M x);
M sinh(M x);
M tanh(M x);
M sech(M x);
M csch(M x);
M coth(M x);

M acosh(M x);
M asinh(M x);
M atanh(M x);
M asech(M x);
M acsch(M x);
M acoth(M x);

M operator !(M x);
M operator -(M x);
M operator +(M x);
M abs(M x);
M angle(M x);
M arg(M x);
M ceil(M x);
M cumprod(M x);
M cumsum(M x);
M conj(M x);
M exp(M x);
M finite(M x);
M fix(M x);
M floor(M x);
M imag(M x);
M isinf(M x);
M isnan(M x);
M log(M x);
M log10(M x);
M real(M x);
M round(M x);
M sign(M x);
M sqrtM(M x);
M tril(M x);
M triu(M x);
M tril(M x, M y);
M triu(M x, M y);
M log2(M x);
M log2(M x, i_o_t, Mr f, Mr e);
M pow2(M x);

double abs(double x);
double abs(double xr, double xi);
double max(double x, double y);
double max(double x, double y, double z);
double min(double x, double y);
double rem(double x, double y);
double round(double x);

M length(M x);
M size(M x);
M size(M x, M type);
M size(M x, i_o_t, Mr m_Mat, Mr n_Mat);

M sort(M x);
M sort(M x, i_o_t, Mr S, Mr I);


M int2str(M x);
M num2str(M x, M prec=4);
M sscanfM(M s, M format);
M sscanfM(M s, M format, M size);
M sscanfM(M s, M format, M size, i_o_t, Mr A, Mr count);
M sscanfM(M s, M format, M size, i_o_t, Mr A, Mr count, Mr errmsg);

M sprintfM(M string);
M sprintfM(M format, M a, M b=l_M, M c=l_M, M d=l_M, M e=l_M, M f=l_M);

M str2mat(M a1=l_M,  M a2=l_M,	M a3=l_M,  M a4=l_M,  M a5=l_M,
	  M a6=l_M,  M a7=l_M,	M a8=l_M,  M a9=l_M,  M a10=l_M,
	  M a11=l_M, M a12=l_M, M a13=l_M, M a14=l_M, M a15=l_M);

M strcmp(M a, M b);
void initM(istream& stin, ostream& stout, ostream& sterr);
void exitM(void);
M brackets(int type,...);
M cd();
M cd(M dir);
M clockM();
M deleteM(M filename);
M dir();
M dos(M command);
M load(char *fname, M& a=l_M, M& b=l_M, M& c=l_M, M& d=l_M, M& e=l_M, M& f=l_M, M& g=l_M, M& h=l_M, M& i=l_M, M& j=l_M);
M loadascii(char *fname, M& a);
M save(char *fname, M a=l_M, M b=l_M, M c=l_M, M d=l_M, M e=l_M, M f=l_M, M g=l_M, M h=l_M, M i=l_M, M j=l_M);
M saveascii(char *fname, M a=l_M, M b=l_M, M c=l_M, M d=l_M, M e=l_M, M f=l_M);
M type(M fname);
M systemM(M cmd);
M system_dependent(...);
M whos();


M clc();
M disp(M x,ostream& outstream=matlibout);
M error(M x);
M format(M x="short");
M home();
M input(M x);
M input(M x, M s);
M menu(M title,
	  M a1=l_M,  M a2=l_M,	M a3=l_M,  M a4=l_M,  M a5=l_M,
	  M a6=l_M,  M a7=l_M,	M a8=l_M,  M a9=l_M,  M a10=l_M,
	  M a11=l_M, M a12=l_M, M a13=l_M, M a14=l_M, M a15=l_M);
ostream& operator<<(ostream& outstream, const M& x);
ostream& operator<<(ostream& outstream, const D& x);
M pauseM();
M pauseM(M x);
M smiley(...);
M warning(M x);
void whos_line(M& x);

M dot_transpose(M x);
M fliplr(M x);
M flipud(M x);
M reshape(M x, M m, M n);
M rot90(M x, M k=one_M);
M transpose(M x, int d=0);


#endif

