#ifndef _NR_H_
#define _NR_H_
#include <fstream>
#include <complex>
#include "nrutil.h"
#include "nrtypes.h"
using namespace std;

namespace NR {

void addint(Mat_O_DP &uf, Mat_I_DP &uc, Mat_O_DP &res);
void airy(const DP x, DP &ai, DP &bi, DP &aip, DP &bip);
void amebsa(Mat_IO_DP &p, Vec_IO_DP &y, Vec_O_DP &pb, DP &yb, const DP ftol,
	DP funk(Vec_I_DP &), int &iter, const DP temptr);
void amoeba(Mat_IO_DP &p, Vec_IO_DP &y, const DP ftol, DP funk(Vec_I_DP &),
	int &nfunk);
DP amotry(Mat_IO_DP &p, Vec_O_DP &y, Vec_IO_DP &psum, DP funk(Vec_I_DP &),
	const int ihi, const DP fac);
DP amotsa(Mat_IO_DP &p, Vec_O_DP &y, Vec_IO_DP &psum, Vec_O_DP &pb, DP &yb,
	DP funk(Vec_I_DP &), const int ihi, DP &yhi, const DP fac);
void anneal(Vec_I_DP &x, Vec_I_DP &y, Vec_IO_INT &iorder);
DP anorm2(Mat_I_DP &a);
void arcmak(Vec_I_ULNG &nfreq, unsigned long nchh, unsigned long nradd,
	arithcode &acode);
void arcode(unsigned long &ich, string &code, unsigned long &lcd,
	const int isign, arithcode &acode);
void arcsum(Vec_I_ULNG &iin, Vec_O_ULNG &iout, unsigned long ja,
	const int nwk, const unsigned long nrad, const unsigned long nc);
void asolve(Vec_I_DP &b, Vec_O_DP &x, const int itrnsp);
void atimes(Vec_I_DP &x, Vec_O_DP &r, const int itrnsp);
void avevar(Vec_I_DP &data, DP &ave, DP &var);
void balanc(Mat_IO_DP &a);
void banbks(Mat_I_DP &a, const int m1, const int m2, Mat_I_DP &al,
	Vec_I_INT &indx, Vec_IO_DP &b);
void bandec(Mat_IO_DP &a, const int m1, const int m2, Mat_O_DP &al,
	Vec_O_INT &indx, DP &d);
void banmul(Mat_I_DP &a, const int m1, const int m2, Vec_I_DP &x,
	Vec_O_DP &b);
void bcucof(Vec_I_DP &y, Vec_I_DP &y1, Vec_I_DP &y2, Vec_I_DP &y12,
	const DP d1, const DP d2, Mat_O_DP &c);
void bcuint(Vec_I_DP &y, Vec_I_DP &y1, Vec_I_DP &y2, Vec_I_DP &y12,
	const DP x1l, const DP x1u, const DP x2l, const DP x2u,
	const DP x1, const DP x2, DP &ansy, DP &ansy1, DP &ansy2);
void beschb(const DP x, DP &gam1, DP &gam2, DP &gampl, DP &gammi);
DP bessi(const int n, const DP x);
DP bessi0(const DP x);
DP bessi1(const DP x);
void bessik(const DP x, const DP xnu, DP &ri, DP &rk, DP &rip, DP &rkp);
DP bessj(const int n, const DP x);
DP bessj0(const DP x);
DP bessj1(const DP x);
void bessjy(const DP x, const DP xnu, DP &rj, DP &ry, DP &rjp, DP &ryp);
DP bessk(const int n, const DP x);
DP bessk0(const DP x);
DP bessk1(const DP x);
DP bessy(const int n, const DP x);
DP bessy0(const DP x);
DP bessy1(const DP x);
DP beta(const DP z, const DP w);
DP betacf(const DP a, const DP b, const DP x);
DP betai(const DP a, const DP b, const DP x);
DP bico(const int n, const int k);
void bksub(const int ne, const int nb, const int jf, const int k1,
	const int k2, Mat3D_IO_DP &c);
DP bnldev(const DP pp, const int n, int &idum);
DP brent(const DP ax, const DP bx, const DP cx, DP f(const DP),
	const DP tol, DP &xmin);
void broydn(Vec_IO_DP &x, bool &check, void vecfunc(Vec_I_DP &, Vec_O_DP &));
void bsstep(Vec_IO_DP &y, Vec_IO_DP &dydx, DP &xx, const DP htry,
	const DP eps, Vec_I_DP &yscal, DP &hdid, DP &hnext,
	void derivs(const DP, Vec_I_DP &, Vec_O_DP &));
void caldat(const int julian, int &mm, int &id, int &iyyy);
void chder(const DP a, const DP b, Vec_I_DP &c, Vec_O_DP &cder, const int n);
DP chebev(const DP a, const DP b, Vec_I_DP &c, const int m, const DP x);
void chebft(const DP a, const DP b, Vec_O_DP &c, DP func(const DP));
void chebpc(Vec_I_DP &c, Vec_O_DP &d);
void chint(const DP a, const DP b, Vec_I_DP &c, Vec_O_DP &cint, const int n);
DP chixy(const DP bang);
void choldc(Mat_IO_DP &a, Vec_O_DP &p);
void cholsl(Mat_I_DP &a, Vec_I_DP &p, Vec_I_DP &b, Vec_O_DP &x);
void chsone(Vec_I_DP &bins, Vec_I_DP &ebins, const int knstrn, DP &df,
	DP &chsq, DP &prob);
void chstwo(Vec_I_DP &bins1, Vec_I_DP &bins2, const int knstrn, DP &df,
	DP &chsq, DP &prob);
void cisi(const DP x, complex<DP> &cs);
void cntab1(Mat_I_INT &nn, DP &chisq, DP &df, DP &prob, DP &cramrv, DP &ccc);
void cntab2(Mat_I_INT &nn, DP &h, DP &hx, DP &hy, DP &hygx, DP &hxgy,
	DP &uygx, DP &uxgy, DP &uxy);
void convlv(Vec_I_DP &data, Vec_I_DP &respns, const int isign,
	Vec_O_DP &ans);
void copy(Mat_O_DP &aout, Mat_I_DP &ain);
void correl(Vec_I_DP &data1, Vec_I_DP &data2, Vec_O_DP &ans);
void cosft1(Vec_IO_DP &y);
void cosft2(Vec_IO_DP &y, const int isign);
void covsrt(Mat_IO_DP &covar, Vec_I_BOOL &ia, const int mfit);
void crank(Vec_IO_DP &w, DP &s);
void cyclic(Vec_I_DP &a, Vec_I_DP &b, Vec_I_DP &c, const DP alpha,
	const DP beta, Vec_I_DP &r, Vec_O_DP &x);
void daub4(Vec_IO_DP &a, const int n, const int isign);
DP dawson(const DP x);
DP dbrent(const DP ax, const DP bx, const DP cx, DP f(const DP),
	DP df(const DP), const DP tol, DP &xmin);
void ddpoly(Vec_I_DP &c, const DP x, Vec_O_DP &pd);
bool decchk(string str, char &ch);
void derivs_s(const DP x, Vec_I_DP &y, Vec_O_DP &dydx);
DP df1dim(const DP x);
void dfpmin(Vec_IO_DP &p, const DP gtol, int &iter, DP &fret,
	DP func(Vec_I_DP &), void dfunc(Vec_I_DP &, Vec_O_DP &));
DP dfridr(DP func(const DP), const DP x, const DP h, DP &err);
void dftcor(const DP w, const DP delta, const DP a, const DP b,
	Vec_I_DP &endpts, DP &corre, DP &corim, DP &corfac);
void dftint(DP func(const DP), const DP a, const DP b, const DP w,
	DP &cosint, DP &sinint);
void difeq(const int k, const int k1, const int k2, const int jsf,
	const int is1, const int isf, Vec_I_INT &indexv, Mat_O_DP &s,
	Mat_I_DP &y);
void dlinmin(Vec_IO_DP &p, Vec_IO_DP &xi, DP &fret, DP func(Vec_I_DP &),
	void dfunc(Vec_I_DP &, Vec_O_DP &));
void eclass(Vec_O_INT &nf, Vec_I_INT &lista, Vec_I_INT &listb);
void eclazz(Vec_O_INT &nf, bool equiv(const int, const int));
DP ei(const DP x);
void eigsrt(Vec_IO_DP &d, Mat_IO_DP &v);
DP elle(const DP phi, const DP ak);
DP ellf(const DP phi, const DP ak);
DP ellpi(const DP phi, const DP en, const DP ak);
void elmhes(Mat_IO_DP &a);
DP erfcc(const DP x);
DP erff(const DP x);
DP erffc(const DP x);
void eulsum(DP &sum, const DP term, const int jterm, Vec_IO_DP &wksp);
DP evlmem(const DP fdt, Vec_I_DP &d, const DP xms);
DP expdev(int &idum);
DP expint(const int n, const DP x);
DP f1dim(const DP x);
DP factln(const int n);
DP factrl(const int n);
void fasper(Vec_I_DP &x, Vec_I_DP &y, const DP ofac, const DP hifac,
	Vec_O_DP &wk1, Vec_O_DP &wk2, int &nout, int &jmax, DP &prob);
void fdjac(Vec_IO_DP &x, Vec_I_DP &fvec, Mat_O_DP &df,
	void vecfunc(Vec_I_DP &, Vec_O_DP &));
void fgauss(const DP x, Vec_I_DP &a, DP &y, Vec_O_DP &dyda);
void fit(Vec_I_DP &x, Vec_I_DP &y, Vec_I_DP &sig, const bool mwt, DP &a,
	DP &b, DP &siga, DP &sigb, DP &chi2, DP &q);
void fitexy(Vec_I_DP &x, Vec_I_DP &y, Vec_I_DP &sigx, Vec_I_DP &sigy,
	DP &a, DP &b, DP &siga, DP &sigb, DP &chi2, DP &q);
void fixrts(Vec_IO_DP &d);
void fleg(const DP x, Vec_O_DP &pl);
void flmoon(const int n, const int nph, int &jd, DP &frac);
DP fmin(Vec_I_DP &x);
void four1(Vec_IO_DP &data, const int isign);
void fourew(Vec_FSTREAM_p &file, int &na, int &nb, int &nc, int &nd);
void fourfs(Vec_FSTREAM_p &file, Vec_I_INT &nn, const int isign);
void fourn(Vec_IO_DP &data, Vec_I_INT &nn, const int isign);
void fpoly(const DP x, Vec_O_DP &p);
void fred2(const DP a, const DP b, Vec_O_DP &t, Vec_O_DP &f, Vec_O_DP &w,
	DP g(const DP), DP ak(const DP, const DP));
DP fredin(const DP x, const DP a, const DP b, Vec_I_DP &t, Vec_I_DP &f,
	Vec_I_DP &w, DP g(const DP), DP ak(const DP, const DP));
void frenel(const DP x, complex<DP> &cs);
void frprmn(Vec_IO_DP &p, const DP ftol, int &iter, DP &fret,
	DP func(Vec_I_DP &), void dfunc(Vec_I_DP &, Vec_O_DP &));
void ftest(Vec_I_DP &data1, Vec_I_DP &data2, DP &f, DP &prob);
DP gamdev(const int ia, int &idum);
DP gammln(const DP xx);
DP gammp(const DP a, const DP x);
DP gammq(const DP a, const DP x);
DP gasdev(int &idum);
void gaucof(Vec_IO_DP &a, Vec_IO_DP &b, const DP amu0, Vec_O_DP &x,
	Vec_O_DP &w);
void gauher(Vec_O_DP &x, Vec_O_DP &w);
void gaujac(Vec_O_DP &x, Vec_O_DP &w, const DP alf, const DP bet);
void gaulag(Vec_O_DP &x, Vec_O_DP &w, const DP alf);
void gauleg(const DP x1, const DP x2, Vec_O_DP &x, Vec_O_DP &w);
void gaussj(Mat_IO_DP &a, Mat_IO_DP &b);
void gcf(DP &gammcf, const DP a, const DP x, DP &gln);
DP golden(const DP ax, const DP bx, const DP cx, DP f(const DP),
	const DP tol, DP &xmin);
void gser(DP &gamser, const DP a, const DP x, DP &gln);
void hpsel(Vec_I_DP &arr, Vec_O_DP &heap);
void hpsort(Vec_IO_DP &ra);
void hqr(Mat_IO_DP &a, Vec_O_CPLX_DP &wri);
void hufapp(Vec_IO_ULNG &index, Vec_I_ULNG &nprob, const unsigned long n,
	const unsigned long m);
void hufdec(unsigned long &ich, string &code, const unsigned long lcode,
	unsigned long &nb, huffcode &hcode);
void hufenc(const unsigned long ich, string &code, unsigned long &nb,
	huffcode &hcode);
void hufmak(Vec_I_ULNG &nfreq, const unsigned long nchin,
	unsigned long &ilong, unsigned long &nlong, huffcode &hcode);
void hunt(Vec_I_DP &xx, const DP x, int &jlo);
void hypdrv(const DP s, Vec_I_DP &yy, Vec_O_DP &dyyds);
complex<DP> hypgeo(const complex<DP> &a, const complex<DP> &b,
	const complex<DP> &c, const complex<DP> &z);
void hypser(const complex<DP> &a, const complex<DP> &b,
	const complex<DP> &c, const complex<DP> &z,
	complex<DP> &series, complex<DP> &deriv);
unsigned short icrc(const unsigned short crc, const string &bufptr,
	const short jinit, const int jrev);
unsigned short icrc1(const unsigned short crc, const unsigned char onech);
unsigned long igray(const unsigned long n, const int is);
void indexx(Vec_I_DP &arr, Vec_O_INT &indx);
void indexx(Vec_I_INT &arr, Vec_O_INT &indx);
void interp(Mat_O_DP &uf, Mat_I_DP &uc);
int irbit1(unsigned long &iseed);
int irbit2(unsigned long &iseed);
void jacobi(Mat_IO_DP &a, Vec_O_DP &d, Mat_O_DP &v, int &nrot);
void jacobn_s(const DP x, Vec_I_DP &y, Vec_O_DP &dfdx, Mat_O_DP &dfdy);
int julday(const int mm, const int id, const int iyyy);
void kendl1(Vec_I_DP &data1, Vec_I_DP &data2, DP &tau, DP &z, DP &prob);
void kendl2(Mat_I_DP &tab, DP &tau, DP &z, DP &prob);
void kermom(Vec_O_DP &w, const DP y);
void ks2d1s(Vec_I_DP &x1, Vec_I_DP &y1, void quadvl(const DP, const DP,
	DP &, DP &, DP &, DP &), DP &d1, DP &prob);
void ks2d2s(Vec_I_DP &x1, Vec_I_DP &y1, Vec_I_DP &x2, Vec_I_DP &y2, DP &d,
	DP &prob);
void ksone(Vec_IO_DP &data, DP func(const DP), DP &d, DP &prob);
void kstwo(Vec_IO_DP &data1, Vec_IO_DP &data2, DP &d, DP &prob);
void laguer(Vec_I_CPLX_DP &a, complex<DP> &x, int &its);
void lfit(Vec_I_DP &x, Vec_I_DP &y, Vec_I_DP &sig, Vec_IO_DP &a,
	Vec_I_BOOL &ia, Mat_O_DP &covar, DP &chisq,
	void funcs(const DP, Vec_O_DP &));
void linbcg(Vec_I_DP &b, Vec_IO_DP &x, const int itol, const DP tol,
	const int itmax, int &iter, DP &err);
void linmin(Vec_IO_DP &p, Vec_IO_DP &xi, DP &fret, DP func(Vec_I_DP &));
void lnsrch(Vec_I_DP &xold, const DP fold, Vec_I_DP &g, Vec_IO_DP &p,
	Vec_O_DP &x, DP &f, const DP stpmax, bool &check, DP func(Vec_I_DP &));
void locate(Vec_I_DP &xx, const DP x, int &j);
void lop(Mat_O_DP &out, Mat_I_DP &u);
void lubksb(Mat_I_DP &a, Vec_I_INT &indx, Vec_IO_DP &b);
void ludcmp(Mat_IO_DP &a, Vec_O_INT &indx, DP &d);
void machar(int &ibeta, int &it, int &irnd, int &ngrd, int &machep,
	int &negep, int &iexp, int &minexp, int &maxexp, DP &eps, DP &epsneg,
	DP &xmin, DP &xmax);
void matadd(Mat_I_DP &a, Mat_I_DP &b, Mat_O_DP &c);
void matsub(Mat_I_DP &a, Mat_I_DP &b, Mat_O_DP &c);
void medfit(Vec_I_DP &x, Vec_I_DP &y, DP &a, DP &b, DP &abdev);
void memcof(Vec_I_DP &data, DP &xms, Vec_O_DP &d);
bool metrop(const DP de, const DP t);
void mgfas(Mat_IO_DP &u, const int maxcyc);
void mglin(Mat_IO_DP &u, const int ncycle);