#ifndef FUNC_H
#define FUNC_H

#include <math.h>

#include "matrix.h"

#ifndef DOUBLE
#define DOUBLE double
#endif
#include "cmatrix.h"
class cmatrix;

DOUBLE gamma(DOUBLE x); // 計(jì)算伽馬函數(shù)
DOUBLE gamma2(DOUBLE a, DOUBLE x); // 計(jì)算不完全伽馬函數(shù)
DOUBLE erf(DOUBLE x); // 計(jì)算誤差函數(shù)
DOUBLE beta(DOUBLE a,DOUBLE b,DOUBLE x); // 計(jì)算不完全貝塔函數(shù)
DOUBLE gass(DOUBLE a,DOUBLE d,DOUBLE x); // 給定均值a，標(biāo)準(zhǔn)差d的正態(tài)分布函數(shù)
DOUBLE student(DOUBLE t, size_t n);	// t-分布函數(shù)
DOUBLE chii(DOUBLE x,size_t n); // X-方分布函數(shù)
DOUBLE fdisp(DOUBLE f,size_t n1,size_t n2); // F-分布函數(shù)
DOUBLE integral(DOUBLE (*f)(DOUBLE),DOUBLE a, DOUBLE b, DOUBLE eps=defaulterr);
 // 函數(shù)f,在(a,b)區(qū)間積分,采用勒讓德高斯求積法
DOUBLE sinn(DOUBLE x);	// 正弦積分函數(shù)
DOUBLE coss(DOUBLE x);	// 余弦積分函數(shù)
DOUBLE expp(DOUBLE x);	// 指數(shù)積分函數(shù)
DOUBLE getroot(DOUBLE (*f)(DOUBLE), DOUBLE x0=0.9, DOUBLE eps=defaulterr);
 // 求函數(shù)f(x)=0在x0附近的根，返回此根，如找不到根，則丟出例外
class algo;
DOUBLE getroot(algo & alg, DOUBLE x0=0.9, DOUBLE eps=defaulterr);
 // 同上，但函數(shù)為algo類變量


class algo	// 算法類
{
 public:
	DOUBLE yfactor;		// 乘因子，初始化為1
	DOUBLE xfactor;		// x軸放大因子，初始化為1
	DOUBLE addconst;	// 加和，初始化為0
	DOUBLE xshift;		// x平移量，初始化為0
	unsigned refnum;		// 引用數(shù),初始化為1
	algo():refnum(1),yfactor(1.0),xfactor(1.0),addconst(0.0),xshift(0.0){};
	 // 構(gòu)造函數(shù),產(chǎn)生y=x線性函數(shù)
	algo(DOUBLE xs, DOUBLE xf,DOUBLE adc=0, DOUBLE yf=1):refnum(1),yfactor(yf),
		addconst(adc),xshift(xs),xfactor(xf){};
	algo(DOUBLE a):refnum(1),yfactor(0.0),xfactor(1.0),addconst(a),xshift(0.0){};
	 // 常函數(shù)的構(gòu)造
	algo(algo & alg):yfactor(alg.yfactor),xfactor(alg.xfactor),
		addconst(alg.addconst),xshift(alg.xshift),refnum(1){}; // 拷貝構(gòu)造函數(shù)
	virtual ~algo(){}; // 虛析構(gòu)函數(shù)
	DOUBLE cal(DOUBLE x); // 計(jì)算算法值
	virtual DOUBLE calculate(DOUBLE x)
		{return x;}; // 本身算法,將被繼承子類改寫
	virtual algo * clone(); // 克隆自己，必須被繼承子類改寫
	algo * mul(DOUBLE a);	// 乘a
	algo * add(DOUBLE a);	// 加a
	algo * neg(); // 取負(fù)
	algo * setxfactor(DOUBLE x);		// 設(shè)置x軸因子
	algo * xroom(DOUBLE x);	// 將xfactor擴(kuò)大x倍
	algo * setxshift(DOUBLE x); // 設(shè)置xshift的值
	algo * xshiftas(DOUBLE x);	// 從當(dāng)前開始平移x
	DOUBLE integ(DOUBLE a, DOUBLE b, DOUBLE eps = defaulterr); // 在(a,b)區(qū)間積分
			// 采用勒讓德高斯求積法
};

enum	method {cadd,csub,cmul,cdiv,cpow,ccom}; // 枚舉加減乘除乘方復(fù)合這四種運(yùn)算

class algojoin : public algo // 結(jié)合算法
{
 public:
	algo * leftalgo;	// 左算法，初始化為0
	algo * rightalgo; // 右算法，初始化為0
	method met;	// 指明算法
	algojoin(algo * l, algo * r, method m):leftalgo(l),
		rightalgo(r), met(m)
		{ if(leftalgo)
				leftalgo->refnum++;
		  if(rightalgo)
				rightalgo->refnum++;
		};
	algojoin(algojoin& alg):algo(alg),
		leftalgo(alg.leftalgo),rightalgo(alg.rightalgo),met(alg.met){
			if(leftalgo)
				leftalgo->refnum++;
			if(rightalgo)
				rightalgo->refnum++;};
			// 拷貝構(gòu)造函數(shù)
	virtual ~algojoin() {
		if(leftalgo) {	// 如左或者右算法已經(jīng)沒有被引用，則刪除
			leftalgo->refnum--;
			if(!leftalgo->refnum) delete leftalgo;
		}
		if(rightalgo) {
			rightalgo->refnum--;
			if(!rightalgo->refnum) delete rightalgo;
		}
	};
	virtual algo * clone(); // 克隆自己
	virtual DOUBLE calculate(DOUBLE x);	// 實(shí)施結(jié)合算法
};

class algofun : public algo	// 函數(shù)算法
{
 public:
	DOUBLE (*f)(DOUBLE);	// 函數(shù)指針
	algofun(DOUBLE (*fun)(DOUBLE)):f(fun){};	// 用函數(shù)指針進(jìn)行初始化
	algofun(algofun& alg):algo(alg),f(alg.f){}; // 拷貝構(gòu)造函數(shù)
	virtual DOUBLE calculate(DOUBLE x);	// 實(shí)施函數(shù)算法
	virtual algo * clone(); // 克隆自己
};

class algogass : public algo	// 正態(tài)分布函數(shù)
{
 public:
	algogass(DOUBLE a, DOUBLE d):algo(a,1.0/d){};	// 用均值a和標(biāo)準(zhǔn)差d進(jìn)行初始化
	algogass(algogass& alg):algo(alg){}; // 拷貝構(gòu)造函數(shù)
	virtual DOUBLE calculate(DOUBLE x);	// 實(shí)施函數(shù)算法
	virtual algo * clone(); // 克隆自己
};

class algoinverse : public algo	// 根據(jù)一個(gè)算法求反函數(shù)的算法
	// 通過不斷求f(x)-y=0在給定的各個(gè)y值下的根來進(jìn)行
{
 public:
	algo * al;
	algoinverse(algo * a):al(a){if(al) al->refnum++;};	// a 是要求反函數(shù)的算法
	algoinverse(algoinverse& alg):algo(alg),al(alg.al){
   	al->refnum++;	}; // 拷貝構(gòu)造函數(shù)
	~algoinverse() {
		if(al) {
			al->refnum--;
			if(!al->refnum)
				delete al;
		}
	};		// 析構(gòu)函數(shù)
	virtual DOUBLE calculate(DOUBLE x);	// 實(shí)施函數(shù)算法
	virtual algo * clone(); // 克隆自己
};

class algoregress : public algo // 線性回歸算法產(chǎn)生線性函數(shù)
{
 public:
	algoregress(matrix & xy, matrix* dt=0); // xy為n行2列的矩陣為n個(gè)樣本點(diǎn)
		// dt必須指向一個(gè)6列一行的矩陣向量，六個(gè)數(shù)依次為偏差平方和，平均標(biāo)準(zhǔn)
		// 偏差，回歸平方和，最大偏差，最小偏差，偏差平均值
};

class algopoly : public algo	// 多項(xiàng)式
{
 public:
	matrix data;	//	n乘1矩陣，存放n-1次多項(xiàng)式的系數(shù)a(0)到a(n-1)
	algopoly(){}; // 缺省構(gòu)造函數(shù)，給子類作調(diào)用
	algopoly(matrix& d):data(d){};	// 用矩陣構(gòu)造多項(xiàng)式
	algopoly(algopoly& alg):algo(alg),data(alg.data){}; // 拷貝構(gòu)造函數(shù)
	virtual DOUBLE calculate(DOUBLE x);	// 實(shí)施函數(shù)算法
	virtual algo * clone();	// 克隆自己
	cmatrix getroots();	// 求出此多項(xiàng)式的所有根
};

class algopair : public algopoly // 最小二乘擬合類
{
 public:
	algopair(matrix& xy, size_t m, DOUBLE & t0,DOUBLE &t1,DOUBLE &t2);
		// m為擬合多項(xiàng)式的項(xiàng)數(shù)，dt0為誤差平方和，dt1為誤差絕對值和，
		// dt2為最大誤差絕對值
	algopair(algopair& alg):algopoly(alg){}; // 拷貝構(gòu)造函數(shù)
};

class algoenter2 : public algo	// 一元全區(qū)間不等距插值
{
 public:
	matrix data;	//	n乘2矩陣，n個(gè)坐標(biāo)，先x后y,x必須從小到大
	algoenter2(matrix& d):data(d){};	// 用矩陣構(gòu)造多項(xiàng)式
	algoenter2(algoenter2& alg):algo(alg),data(alg.data){}; // 拷貝構(gòu)造函數(shù)
	virtual DOUBLE calculate(DOUBLE x);	// 實(shí)施函數(shù)算法
	virtual algo * clone();	// 克隆自己
};

class algoenter : public algo	// 一元全區(qū)間等距插值
{
 public:
	DOUBLE x0;		// 起始點(diǎn)
	DOUBLE h;		// 步長
	matrix data;	//	n乘1矩陣，n個(gè)函數(shù)值
	algoenter(matrix& d,DOUBLE xx0, DOUBLE hh):
		data(d),x0(xx0),h(hh){};
	algoenter(algoenter& alg):algo(alg),data(alg.data),
		x0(alg.x0),h(alg.h){}; // 拷貝構(gòu)造函數(shù)
	virtual DOUBLE calculate(DOUBLE x);	// 實(shí)施函數(shù)算法
	virtual algo * clone();	// 克隆自己
};

enum funckind {polyfunc, enter2func, gammafunc}; // 函數(shù)種類

class func {	// 函數(shù)類
 public:
	algo * alg;	// 決定函數(shù)的算法

	func();	// 缺省構(gòu)造函數(shù)
	func(DOUBLE a);	// 常函數(shù)的構(gòu)造函數(shù)
	func(DOUBLE (*fun)(DOUBLE));	// 函數(shù)指針的構(gòu)造函數(shù)
	func(func & fn);	// 拷貝構(gòu)造函數(shù)
	func(algo * a); // 算法構(gòu)造函數(shù)
   func(algo& a); // 算法構(gòu)造函數(shù)
	func(DOUBLE a, DOUBLE d);	// 產(chǎn)生正態(tài)分布函數(shù)a為均值，d為標(biāo)準(zhǔn)差
	func(matrix& m, funckind kind=polyfunc);	// 構(gòu)造數(shù)值相關(guān)函數(shù)，
				// 如kind = polyfunc, 則m為nX1矩陣，是n-1階多項(xiàng)式系數(shù)，
				// 其中m(0,0)為常數(shù)項(xiàng)，m(n-1,0)為n-1次項(xiàng)。
				// 如kind = enter2func, 則m為nX2矩陣，代表n個(gè)坐標(biāo)點(diǎn)
				// 其中第1列是由小到大排過序的各點(diǎn)的x坐標(biāo)
	func(matrix& m, DOUBLE x0, DOUBLE h); // 構(gòu)造等距插值函數(shù)
			// 其中m是nX1階矩陣，代表n個(gè)y值，x0是起始點(diǎn)，h是步長（采樣間隔）
	virtual ~func() {		// 析構(gòu)函數(shù)
		if(alg) {
			alg->refnum--;	// 引用數(shù)減一，如再無其它引用，則刪除算法
			if(!alg->refnum)
				delete alg;
		}
	};

	DOUBLE operator()(DOUBLE x){return alg->cal(x);}; // 計(jì)算x的函數(shù)值
	func& operator=(func& fn);	// 賦值運(yùn)算符
	func& operator=(DOUBLE (*fn)(DOUBLE)); // 用函數(shù)指針的賦值運(yùn)算符
	func& operator=(DOUBLE a); // 常函數(shù)的賦值運(yùn)算符

	func& operator+=(func& fn);	// 自身加一個(gè)函數(shù)
	func& operator+=(DOUBLE a){alg=alg->add(a);return (*this);};//自身加一個(gè)常數(shù)
	func& operator+=(DOUBLE (*f)(DOUBLE)); // 自身加一個(gè)函數(shù)指針
	func operator+(func& fn);	// 相加產(chǎn)生新函數(shù)
	func operator+(DOUBLE a);	// 與常數(shù)相加產(chǎn)生新函數(shù)
	friend func operator+(DOUBLE a, func& f); // 同上但常數(shù)在前
	func operator+(DOUBLE (*f)(DOUBLE)); // 加一個(gè)函數(shù)指針產(chǎn)生新函數(shù)
	friend func operator+(DOUBLE (*f)(DOUBLE),func& fn); // 同上但函數(shù)指針在前

	func& neg(); // 自身取負(fù)
	func operator-(); // 產(chǎn)生負(fù)函數(shù)

	func& operator-=(func& fn); // 自身減一個(gè)函數(shù)
	func& operator-=(DOUBLE a){alg=alg->add(-a);return (*this);};//自身減一個(gè)常數(shù)
	func& operator-=(DOUBLE (*f)(DOUBLE)); // 自身減一個(gè)函數(shù)指針
	func operator-(func& fn);	// 相減產(chǎn)生新函數(shù)
	func operator-(DOUBLE a);	// 與常數(shù)相減產(chǎn)生新函數(shù)
	friend func operator-(DOUBLE a, func& f); // 同上但常數(shù)在前
	func operator-(DOUBLE (*f)(DOUBLE)); // 減一個(gè)函數(shù)指針產(chǎn)生新函數(shù)
	friend func operator-(DOUBLE (*f)(DOUBLE),func& fn); // 函數(shù)指針減函數(shù)

	func& operator*=(func& fn);	// 自身乘一個(gè)函數(shù)
	func& operator*=(DOUBLE a){alg=alg->mul(a);return (*this);};//自身乘一個(gè)常數(shù)
	func& operator*=(DOUBLE (*f)(DOUBLE)); // 自身乘一個(gè)函數(shù)指針
	func operator*(func& fn);	// 相乘產(chǎn)生新函數(shù)
	func operator*(DOUBLE a);	// 與常數(shù)相乘產(chǎn)生新函數(shù)
	friend func operator*(DOUBLE a, func& f); // 同上但常數(shù)在前
	func operator*(DOUBLE (*f)(DOUBLE)); // 乘一個(gè)函數(shù)指針產(chǎn)生新函數(shù)
	friend func operator*(DOUBLE (*f)(DOUBLE),func& fn); // 函數(shù)指針乘函數(shù)

	func& operator/=(func& fn);	// 自身除以一個(gè)函數(shù)
	func& operator/=(DOUBLE a){alg=alg->mul(1.0/a);return (*this);
			};//自身除以常數(shù)
	func& operator/=(DOUBLE (*f)(DOUBLE)); // 自身除以一個(gè)函數(shù)指針
	func operator/(func& fn);	// 相除產(chǎn)生新函數(shù)
	func operator/(DOUBLE a);	// 與常數(shù)相除產(chǎn)生新函數(shù)
	friend func operator/(DOUBLE a, func& f); // 常數(shù)除以函數(shù)
	func operator/(DOUBLE (*f)(DOUBLE)); // 除以一個(gè)函數(shù)指針產(chǎn)生新函數(shù)
	friend func operator/(DOUBLE (*f)(DOUBLE),func& fn); // 函數(shù)指針除以函數(shù)

	void setxfactor(DOUBLE a);	// 設(shè)置x因子為a
	void xroom(DOUBLE a);	  // x方向擴(kuò)大a倍
	void setxshift(DOUBLE a);	// 設(shè)置函數(shù)沿x軸平移a
	void shiftxas(DOUBLE a); // 函數(shù)沿x軸右移a

	func& power(func& f);	// 函數(shù)的f次乘冪，函數(shù)自身改變
	func& power(DOUBLE a);	// 函數(shù)的a次冪，函數(shù)自身改變
	func operator^(func & fn);	// 函數(shù)的fn次乘冪，產(chǎn)生新函數(shù)，原函數(shù)不變
	func operator^(DOUBLE a);  // 函數(shù)的a次冪，產(chǎn)生新函數(shù)，原函數(shù)不變

	func operator()(func & fn);	// 復(fù)合函數(shù)，產(chǎn)生新的函數(shù)

	DOUBLE integ(DOUBLE a, DOUBLE b, DOUBLE eps=defaulterr);
		// 從a到b計(jì)算函數(shù)的定積分
	DOUBLE singleroot(DOUBLE x=0.9, DOUBLE eps = defaulterr);
		// 計(jì)算函數(shù)f(x)=0的一個(gè)單根
	func inverse();	// 產(chǎn)生反函數(shù)

	void getab(DOUBLE& a,DOUBLE& b) {	// 主要被線性回歸子類用，返回線性因子
   												// 和加常數(shù)
		a = alg->yfactor;
		b = alg->addconst;
	};
};

inline func operator+(DOUBLE a, func& f) // 常數(shù)加函數(shù)
{	return f+a; }

inline func operator+(DOUBLE (*f)(DOUBLE),func& fn) // 函數(shù)指針加函數(shù)
{	return fn+f;}

func operator-(DOUBLE a, func& f); // 常數(shù)減函數(shù)
func operator-(DOUBLE (*f)(DOUBLE),func& fn); // 函數(shù)指針減函數(shù)

inline func operator*(DOUBLE a, func& f) // 常數(shù)乘函數(shù)
{	return f*a; }
inline func operator*(DOUBLE (*f)(DOUBLE),func& fn) // 函數(shù)指針乘函數(shù)
{	return fn*f;}

func operator/(DOUBLE a, func& f); // 常數(shù)除以函數(shù)
func operator/(DOUBLE (*f)(DOUBLE),func& fn); // 函數(shù)指針除以函數(shù)

class funcgass : public func // 正態(tài)分布函數(shù)
{
 public:
	funcgass(DOUBLE a=0.0, DOUBLE d=1.0):func(new algogass(a,d)){};
	void setmandd(DOUBLE a, DOUBLE d){
		alg->xshift = a;
		alg->xfactor = 1.0/d; }; // 設(shè)置均值和標(biāo)準(zhǔn)差
	void getmandd(DOUBLE &a, DOUBLE &d) {
		a = alg->xshift; d=1.0/alg->xfactor;}; // 獲得均值和標(biāo)準(zhǔn)差
};

class funcpoly : public func	// 多項(xiàng)式函數(shù)
{
 public:
	funcpoly(matrix& d):func(new algopoly(d)){};
	matrix& getdata(){return ((algopoly*)alg)->data;};
	cmatrix getroots(){return ((algopoly*)alg)->getroots();};
};

class funcenter2 : public func // 一元全區(qū)間不等距插值
{
 public:
	funcenter2(matrix& d):func(new algoenter2(d)){};
	matrix& getdata(){return ((algoenter2*)alg)->data;};
};

class funcenter : public func // 一元全區(qū)間等距插值
{
 public:
	funcenter(matrix& d, DOUBLE x0, DOUBLE h):
		func(new algoenter(d,x0,h)){};
	matrix& getdata(DOUBLE &xx0, DOUBLE &hh){
		xx0 = ((algoenter*)alg)->x0; hh = ((algoenter*)alg)->h;
		return ((algoenter*)alg)->data; };
};

class funcpair : public func // 最小二乘擬合函數(shù)
{
 public:
	funcpair(matrix& xy, size_t m, DOUBLE & t0,DOUBLE &t1,DOUBLE &t2);
		// xy為n行2列數(shù)組，存放n個(gè)樣本點(diǎn)
			// m為擬合多項(xiàng)式的項(xiàng)數(shù)，dt0為誤差平方和，dt1為誤差絕對值和，
		// dt2為最大誤差絕對值
	matrix & getdata(){return ((algopair*)alg)->data;};
	 // 返回結(jié)果擬合系數(shù)一維矩陣
	cmatrix getroots(){return ((algopoly*)alg)->getroots();};
		// 求出所有根
};

class funcregress : public func // 線性回歸函數(shù)，其實(shí)就是線性函數(shù)ax+b
	// 但由數(shù)據(jù)樣本點(diǎn)構(gòu)成a與b
{
 public:
	funcregress(matrix & xy, matrix* dt=0);  // xy為n行2列的矩陣為n個(gè)樣本點(diǎn)
		// dt必須指向一個(gè)6列一行的矩陣向量，六個(gè)數(shù)依次為偏差平方和，平均標(biāo)準(zhǔn)
		// 偏差，回歸平方和，最大偏差，最小偏差，偏差平均值
};

#endif // FUNC_H