#include <iostream>

#define Allocate(n,t) (t*)malloc((n)*sizeof(t))
#define AllocateC(n,t) (t*)calloc(n, sizeof(t))
#define Reallocate(p, n, t) ((p)=(t*)realloc((p),(n)*sizeof(t)))
#define AllocCopy(n,t,v) (t*)memcpy(Allocate(n,t),v,(n)*sizeof(t))

#define die(msg) {fprintf(stderr,"\n%s\n",msg);exit(-1);}

#ifndef max
#define max(a,b) ((a)>(b)?(a):(b))
#endif
#ifndef min
#define min(a,b) ((a)<(b)?(a):(b))
#endif


// Simple Vector data type.
typedef struct VectorRec {
  double *d;
  int len;
} *Vector;


#define VectorMove(v1, v2) memcpy((v1)->d, (v2)->d, (v2)->len*sizeof(double))

Vector VectorFromData(int length, double *data);

/* Returns the sum of all elements of v. */
double VectorSum(Vector v);

/* Returns a new, uninitialized Vector whose length is given. */
Vector VectorCreate(int length);

/* Compute vector dot product = VectorSum(VectorMultiply(v1,v2))
 * of v1 and v2. 
 */
double VectorDot(Vector v1, Vector v2);

/* Set all elements of v to x. */
void VectorSet(Vector v, double x);

/* Frees all storage associated with v. */
void VectorFree(Vector v);

/* Write a character representation of v to fp. */
void VectorPrint(FILE *fp, Vector v);
void VectorPrintNL(FILE *fp, Vector v);

/* If len < v->len, truncates v.
 * If len > v->len, adds garbage to the end of v. 
 */
void VectorResize(Vector v, int len);


// simple matrix
typedef struct MatrixRec {
  double *d;
  int len;
  int cols;
  double **d2;
  int rows;
} *Matrix;

/* m1 gets the contents of m2. m1 must be already allocated and have exactly
 * the same dimensions as m2.
 */
#define MatrixMove(m1, m2) memcpy((m1)->d, (m2)->d, (m2)->len*sizeof(double))

Matrix MatrixCopy(Matrix m);

/* Write a character representation of m to fp. This format was chosen so
 * that row Matrices print just like Vectors.
 */
void MatrixPrint(FILE *fp, Matrix m);

/* Returns a Matrix whose data is given. The array becomes property 
 * of the Matrix and will be freed when the Matrix is. The data must be
 * rows*cols elements long.
 */
Matrix MatrixFromData(int rows, int cols, double *data);

/* Returns a new, uninitialized Matrix whose dimensions are given. */
Matrix MatrixCreate(int rows, int cols);

/* Frees all storage associated with m. */
void MatrixFree(Matrix m);

void MatrixVectorMultiply(Vector dest, Matrix a, Vector v);

Matrix MatrixScan(FILE *fp);

Matrix MatrixReadBinary(FILE *fp);

void MatrixReshape(Matrix m, int rows, int cols);

void MatrixTranspose(Matrix m);

void MatrixMultiply(Matrix dest, Matrix a, Matrix b);

// Singular value decomposition
int MatrixSVD(Matrix A, Matrix U, Matrix V, Vector S);

int MatrixPseudoInvert(Matrix A);

int MatrixInvert(Matrix A);

int MatrixEigenvalues(Matrix A, double *vr, double *vi);

int MatrixEigenvectorsSym(Matrix A, double *v);

int MatrixEigenvectorsSymX(Matrix A, double *v, int count, Matrix V);


/* Volume implementation.
 * A Volume is a structure of {d, len, xdim, d2, zydim, ydim, d3, zdim}
 */

/*****************************************************************************/
typedef struct VolumeRec {
  double *d;
  int len;
  int xdim;
  double **d2;
  int zydim;
  int ydim;
  double ***d3;
  int zdim;
} *Volume;

/* v1 gets the contents of v2. v1 must be already allocated and have exactly
 * the same dimensions as v2.
 */
#define VolumeMove(v1, v2) memcpy((v1)->d, (v2)->d, (v2)->len*sizeof(double))

Volume VolumeCreate(int zdim, int ydim, int xdim);
Volume VolumeFromData(int zdim, int ydim, int xdim, double *data);
void VolumeFree(Volume v);
Volume VolumeCopy(Volume v);
Volume VolumeFromMatrix(Matrix m, int zdim);
Matrix MatrixFromVolume(Volume v, int z);

Matrix VolumeGray(Volume v);
void VolumePrint(FILE *fp, Volume v);

/* Returns a z slice of a volume. */
Matrix MatrixFromVolume(Volume v, int z);