char *digit="0123456789";
void nrerror(char error_text[]);
char *cvector(long nl, long nh);
int *ivector(long nl, long nh);
long *lvector(long nl, long nh);
float *fvector(long nl, long nh);
double *dvector(long nl, long nh);
long double *ldvector(long nl, long nh);
/* double _huge *dvector(long nl, long nh); */
char **cmatrix(long nrl, long nrh, long ncl, long nch);
int **imatrix(long nrl, long nrh, long ncl, long nch);
double **dmatrix(long nrl, long nrh, long ncl, long nch);
float **fmatrix(long nrl, long nrh, long ncl, long nch);
long double **ldmatrix(long nrl, long nrh, long ncl, long nch);
int ***i3tensor(long nrl,long nrh,long ncl,long nch,long ndl,long ndh);
double ***d3tensor(long nrl,long nrh,long ncl,long nch,long ndl,long ndh);
long double ***ld3tensor(long nrl,long nrh,long ncl,long nch,long ndl,long ndh);
void forvect2mat(double *a, long nrl, long nrh, long ncl, long nch, double **m);
double *mat2forvect(double **m, long nrl, long nrh, long ncl, long nch);
double **foreigen2diagmat(double *a, long nrl, long nrh);

void free_ivector(int *v, long nl, long nh);
void free_lvector(long *v, long nl, long nh);
void free_dvector(double *v, long nl, long nh);
void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch);
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch);
void free_d3tensor(double ***t,long nrl,long nrh,long ncl,long nch,long ndl,long ndh);

#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#define NR_END 1
#define FREE_ARG char*


void nrerror(char error_text[]){
fprintf(stderr,"%s\n",error_text);
fprintf(stderr,"FATAL - execution terminated\n");
exit(1);
}


char *cvector(long nl, long nh){
char *v;
v=(char *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
if (!v) nrerror("allocation failure in ivector()");
return v-nl+NR_END;
}


int *ivector(long nl, long nh){
int *v;
v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
if (!v) nrerror("allocation failure in ivector()");
return v-nl+NR_END;
}

long *lvector(long nl, long nh){
long int *v;
v=(long int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(long int)));
if (!v) nrerror("allocation failure in lvector()");
return v-nl+NR_END;
}

float *fvector(long nl, long nh){
float *v;
v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float)));
if (!v) nrerror("allocation failure in fvector()");
return v-nl+NR_END;
}


double *dvector(long nl, long nh){
double *v;
v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
if (!v) nrerror("allocation failure in dvector()");
return v-nl+NR_END;
}



long double *ldvector(long nl, long nh){
long double *v;
v=(long double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(long double)));
if (!v) nrerror("allocation failure in ldvector()");
return v-nl+NR_END;
}

/*
double _huge *dvector(long nl, long nh){
double *v;
v=(double _huge *)halloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
if (!v) nrerror("allocation failure in dvector()");
return v-nl+NR_END;
}*/



char **cmatrix(long nrl, long nrh, long ncl, long nch){
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
char **m;
m=(char **)malloc((size_t)((nrow+NR_END)*sizeof(char*)));
if (!m) nrerror("allocation failure 1 in cmatrix()");
m += NR_END;
m -= nrl;

m[nrl]=(char *)malloc((size_t)((nrow*ncol+NR_END)*sizeof(char)));
if (!m[nrl]) nrerror("allocation failure 2 in cmatrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;

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

return m;

}




int **imatrix(long nrl, long nrh, long ncl, long nch){
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
int **m;
m=(int **)malloc((size_t)((nrow+NR_END)*sizeof(int*)));
if (!m) nrerror("allocation failure 1 in imatrix()");
m += NR_END;
m -= nrl;

m[nrl]=(int *)malloc((size_t)((nrow*ncol+NR_END)*sizeof(int)));
if (!m[nrl]) nrerror("allocation failure 2 in imatrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;

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

return m;

}



double **dmatrix(long nrl, long nrh, long ncl, long nch){
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
double **m;
m=(double **)malloc((size_t)((nrow+NR_END)*sizeof(double*)));
if (!m) nrerror("allocation failure 1 in dmatrix()");
m += NR_END;
m -= nrl;

m[nrl]=(double *)malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
if (!m[nrl]) nrerror("allocation failure 2 in dmatrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;

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

return m;

}

float **fmatrix(long nrl, long nrh, long ncl, long nch){

long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
float **m;
m=(float **)malloc((size_t)((nrow+NR_END)*sizeof(float*)));
if (!m) nrerror("allocation failure 1 in fmatrix()");
m += NR_END;
m -= nrl;

m[nrl]=(float *)malloc((size_t)((nrow*ncol+NR_END)*sizeof(float)));
if (!m[nrl]) nrerror("allocation failure 2 in fmatrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;

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

return m;

}


long double **ldmatrix(long nrl, long nrh, long ncl, long nch){
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
long double **m;
m=(long double **)malloc((size_t)((nrow+NR_END)*sizeof(long double*)));
if (!m) nrerror("allocation failure 1 in dmatrix()");
m += NR_END;
m -= nrl;

m[nrl]=(long double *)malloc((size_t)((nrow*ncol+NR_END)*sizeof(long double)));
if (!m[nrl]) nrerror("allocation failure 2 in dmatrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;

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

return m;

}

double *mat2forvect(double **m, long nrl, long nrh, long ncl, long nch){
int i,j,k,n;
double *fa;
i=nrh-nrl+1;
j=nch-ncl+1;
n=i*j;
fa=dvector(0,n-1);
k=0;
for(j=ncl;j<=nch;++j)for(i=nrl;i<=nrh;++i){
	fa[k]=m[i][j];
	k++;
	}
return fa;
}


void forvect2mat(double *a, long nrl, long nrh, long ncl, long nch, double **m){
int i,j,i1,j1;
i1=nrh-nrl+1;
j1=nch-ncl+1;
for(i=1;i<=i1;++i)for(j=1;j<=j1;++j){
	m[i][j]=a[i1*(j-1)+i-1];
	}
}



double **foreigen2diagmat(double *a, long nrl, long nrh){
int i,j;
double **m;
m=dmatrix(nrl,nrh,nrl,nrh);
for(i=nrl;i<=nrh;++i)for(j=nrl;j<=nrh;++j){
        if(i==j)m[i][i]=a[i-nrl];
        else m[i][j]=0;
        }
return m;
}

int ***i3tensor(long nrl,long nrh,long ncl,long nch,long ndl,long ndh){
long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
int ***t;

t=(int ***) malloc((size_t)((nrow+NR_END)*sizeof(int**)));
if(!t)nrerror("allocation failure 1 in i3tensor()");
t += NR_END;
t -= nrl;

t[nrl]=(int **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(int*)));
if(!t[nrl])nrerror("allocation failure 2 in i3tensor()");
t[nrl] += NR_END;
t[nrl] -= ncl;

t[nrl][ncl]=(int *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(int)));
if(!t[nrl][ncl])nrerror("allocation failure 3 in i3tensor()");
t[nrl][ncl] += NR_END;
t[nrl][ncl] -= ndl;

for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
for(i=nrl+1;i<=nrh;i++){
	t[i]=t[i-1]+ncol;
	t[i][ncl]=t[i-1][ncl]+ncol*ndep;
	for(j=ncl+1;j<=nch;j++)t[i][j]=t[i][j-1]+ndep;
	}
return t;
}



double ***d3tensor(long nrl,long nrh,long ncl,long nch,long ndl,long ndh){
long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
double ***t;

t=(double ***) malloc((size_t)((nrow+NR_END)*sizeof(double**)));
if(!t)nrerror("allocation failure 1 in d3tensor()");
t += NR_END;
t -= nrl;

t[nrl]=(double **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double*)));
if(!t[nrl])nrerror("allocation failure 2 in d3tensor()");
t[nrl] += NR_END;
t[nrl] -= ncl;

t[nrl][ncl]=(double *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(double)));
if(!t[nrl][ncl])nrerror("allocation failure 3 in d3tensor()");
t[nrl][ncl] += NR_END;
t[nrl][ncl] -= ndl;

for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
for(i=nrl+1;i<=nrh;i++){
	t[i]=t[i-1]+ncol;
	t[i][ncl]=t[i-1][ncl]+ncol*ndep;
	for(j=ncl+1;j<=nch;j++)t[i][j]=t[i][j-1]+ndep;
	}
return t;
}


long double ***ld3tensor(long nrl,long nrh,long ncl,long nch,long ndl,long ndh){
long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
long double ***t;

t=(long double ***) malloc((size_t)((nrow+NR_END)*sizeof(long double**)));
if(!t)nrerror("allocation failure 1 in d3tensor()");
t += NR_END;
t -= nrl;

t[nrl]=(long double **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(long double*)));
if(!t[nrl])nrerror("allocation failure 2 in d3tensor()");
t[nrl] += NR_END;
t[nrl] -= ncl;

t[nrl][ncl]=(long double *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(long double)));
if(!t[nrl][ncl])nrerror("allocation failure 3 in d3tensor()");
t[nrl][ncl] += NR_END;
t[nrl][ncl] -= ndl;

for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
for(i=nrl+1;i<=nrh;i++){
	t[i]=t[i-1]+ncol;
	t[i][ncl]=t[i-1][ncl]+ncol*ndep;
	for(j=ncl+1;j<=nch;j++)t[i][j]=t[i][j-1]+ndep;
	}
return t;
}

void free_ivector(int *v, long nl, long nh){
free((FREE_ARG) (v+nl-NR_END));
}

void free_lvector(long *v, long nl, long nh){
free((FREE_ARG) (v+nl-NR_END));
}

void free_dvector(double *v, long nl, long nh){
free((FREE_ARG) (v+nl-NR_END));
}

void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch){
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}

void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch){
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}

void free_d3tensor(double ***t,long nrl,long nrh,long ncl,long nch,long ndl,long ndh){
free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
free((FREE_ARG) (t[nrl]+ncl-NR_END));
free((FREE_ARG) (t+nrl-NR_END));
}

