#ifndef __MATLIB_H
#define __MATLIB_H


/*
 *    MATLIB header file
 *
 *    Copyright (c) Yaron Keren 1995. All rights reserved.
 */

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

#define M_MATRIX 0
#define M_TEXT   1

union c_p_t {
  int c;
  int p;
};

union i_o_t {
  int i;
  int o;
};

class M {
  friend class Mr;
private:
  char *name;
  int type;
  int m;
  int n;
  double *pr;
  double *pi;

  void allocate_double();
  void delete_double();
public:

  M();
  M(int nrows,int ncols);
  M(char *mname, double valr, double vali=0);
  M(char *mname, char *text);

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

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

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

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

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

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

  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, M& index);
  Mr(M& src, c_p_t);
  Mr(M& src, M& indexc, M& indexr);
  Mr(M& src, M& indexc, c_p_t);
  Mr(M& src, c_p_t, M& indexr);
  Mr(M& src, c_p_t, c_p_t);

  ~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 =(M& src);
  void operator =(Mr& src);

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

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 xor(M 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);

extern double 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 s_M, l_M;
extern char *s_s, *l_s;
extern double s_d, l_d;
extern double *s_Mp, *l_Mp, *x_Mp;
extern c_p_t c_p;
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);
typedef M func5M(M, M, M, M, M);
typedef M func6M(M, M, M, M, M, M);
typedef M func7M(M, M, M, M, M, M, M);
typedef M func8M(M, M, M, M, M, M, M, M);

M warning(M x);
M error(M x);
M error(char *s);
M warning(char *s);

M eye(M x);
M eye(M x, M y);
M ones(int m);
M ones(int m, int 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(int m);
M zeros(int m, int 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 fclose(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="l");
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 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 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 i_o, Mr y, Mr zf);
M filter(M b, M a, M x, M zi, i_o_t i_o, 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 get(...);
M line(...);
M set(...);
M subplot(...);
M surf(...);

M wavwrite(M wavedata, M rate, M filename);

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

M qz(M a, M b);
M qz(M a, M b, i_o_t, Mr aa);
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 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 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 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, i_o_t, Mr U, Mr S, Mr V);


M angle(M x);
M clear(Mr x);
M colon(double j, double i, double k);
M colon(M j_Mat, M i_Mat, M k_Mat);
M computer(i_o_t, Mr c, Mr maxsize);
M computer();
M cputime();
M diag(M x);
M diag(M x, M kk);
M diff(M x, M order=one_M);
M dot_transpose(M x);
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);
M fliplr(M x);
M flipud(M x);
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 length(M x);
M linspace(M x1, M x2);
M linspace(M x1, M x2, M n);
M logspace(M x1, M x2);
M logspace(M x1, M x2, M n);
M log2(M x, i_o_t, Mr f, Mr e);
M magic(M n);
M matlabpath(M x);
M mean(M x);
M nnz(M x);
M nonzeros(M x);
M quit();
M reshape(M x, M m, M n);
M rot90(M x);
M rot90(M x, M k);
M setstr(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 sparse(...);
M spdiags(...);
M spparms(...);
M tic();
M toc();
M transpose(M x, int d=0);

M min(M x, i_o_t, Mr y, Mr i);
M max(M x, i_o_t, Mr y, Mr i);

M operator !(M x);
M operator -(M x);
M abs(M x);
M acos(M x);
M acot(M x);
M acsc(M x);
M asin(M x);
M asec(M x);
M atan(M x);
M ceil(M x);
M cos(M x);
M cot(M x);
M csc(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 sec(M x);
M sign(M x);
M sin(M x);
M sqrtM(M x);
M tan(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 pow2(M x);

double scalar(M& x);
double scalar(M& x, int idx);
double scalar(M& x,int row,int col);
double scalai(M& x);
double scalai(M& x, int idx);
double scalai(M& x,int row,int col);

void setr(M& x, double r);
void setr(M& x, double r, int idx);
void setr(M& x, double r, int row, int col);
void seti(M& x, double r);
void seti(M& x, double r, int idx);
void seti(M& x, double r, int row, int col);

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


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 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 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(void);
void exitM(void);
M brackets(char type,...);
M cd();
M cd(M dir);
M clockM();
M deleteM(M filename);
M dir(M );
M dos(M command);
M load(char *fname);
M loadascii(char *fname);
M save(char *fname, M& x);
M saveascii(char *fname, M& x);
M type(M fname);
M system_dependent(...);
M whos();

#include <iostream.h>

M clc();
M disp(M x,ostream& outstream=cout);
M home();
M input(M x);
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, M& x);
M pauseM();
M pauseM(M x);
M smiley(...);


#endif

