#include<stdlib.h>
#include<malloc.h>
#include<vector>
#include<iostream.h>

#define FREE_ARG char*


/* Vectors and matrices */
/*********************************************************************/


Number *vector(long nh)
/* allocate a Number vector with subscript range v[1..nh] */
{
	Number* v = new Number[nh+1];
	return v;

}


int *ivector(long nh)
/* allocate an int vector with subscript range v[1..nh] */
{
	int* v = new int[nh+1];
	return v;
}

unsigned long *lvector(long nh)
/* allocate an unsigned long vector with subscript range v[1..nh] */
{
	unsigned long* v = new unsigned long[nh+1];
	return v;
}

unsigned char *cvector(long nh)
/* allocate an unsigned char vector with subscript range v[1..nh] */
{
	unsigned char* v = new unsigned char[nh+1];
	return v;
}

Number **matrix(long nrl, long nrh, long ncl, long nch)
/* allocate a Number matrix with subscript range m[nrl..nrh][ncl..nch] */
{
	long i, nrow=nrh,ncol=nch;
	Number **m;

	/* allocate pointers to rows */
	m=(Number **) malloc((size_t)((nrow+1)*sizeof(Number*)));
	try{if (!m){
		throw "allocation failure 1 in matrix()";
	}
	}
	catch(char* message){
	cerr << message;
	}
	m += 1;
	m -= 1;

	/* allocate rows and set pointers to them */
	m[1]=(Number *) malloc((size_t)((nrow*ncol+1)*sizeof(Number)));
	if (!m[1]) cout << "allocation failure 2 in matrix()";
	m[1] += 1;
	m[1] -= 1;

	for(i=2;i<=nrh;i++) m[i]=m[i-1]+ncol;

	/* return pointer to array of pointers to rows */
	return m;
}


void free_matrix(Number **m, long nrl, long nrh, long ncl, long nch)
/* free a Number matrix allocated by matrix() */
{
	free((FREE_ARG) (m[nrl]+ncl-1));
	free((FREE_ARG) (m+nrl-1));
}


/* Other functions */
/*********************************************************************/
static Number maxarg1,maxarg2;
#define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
        (maxarg1) : (maxarg2))

static Number sqrarg;
#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)

static Number minarg1,minarg2;
#define FMIN(a,b) (minarg1=(a),minarg2=(b),(minarg1) < (minarg2) ?\
        (minarg1) : (minarg2))

static long lmaxarg1,lmaxarg2;
#define LMAX(a,b) (lmaxarg1=(a),lmaxarg2=(b),(lmaxarg1) > (lmaxarg2) ?\
        (lmaxarg1) : (lmaxarg2))

static long lminarg1,lminarg2;
#define LMIN(a,b) (lminarg1=(a),lminarg2=(b),(lminarg1) < (lminarg2) ?\
        (lminarg1) : (lminarg2))

static int imaxarg1,imaxarg2;
#define IMAX(a,b) (imaxarg1=(a),imaxarg2=(b),(imaxarg1) > (imaxarg2) ?\
        (imaxarg1) : (imaxarg2))

static int iminarg1,iminarg2;
#define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ?\
        (iminarg1) : (iminarg2))

#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))