specialfunction.h

数学特殊函数算法库

	_Ty c[7]={ 1.25331414,-0.07832358,0.02189568, 
		-0.01062446,0.00587872,-0.0025154,0.00053208}; 
	_Ty d[7]={ 1.25331414,0.23498619,-0.0365562, 
		0.01504268,-0.00780353,0.00325614,-0.00068245}; 
	if(n<0) n=-n; 
	if(x<0.0) x=-x; 
	if(FloatEqual(x,0.0)) return(1.0e+70); 
	if(n!=1) 
	{  
		if(x < 2.0 || FloatEqual(x, 2.0)) 
		{  
			y=x*x/4.0; p=a[6]; 
			for(i=5; i>=0; i--) p=p*y+a[i]; 
			p=p-TransformativeIntegerBessel1stFunction(0,x)*log(x/2.0); 
		} 
		else 
		{  
			y=2.0/x; p=c[6]; 
			for(i=5; i>=0; i--) p=p*y+c[i]; 
			p=p*exp(-x)/sqrt(x); 
		} 
	} 
	if(n==0) return(p); 
	b0=p; 
	if(x < 2.0 || FloatEqual(x, 2.0)) 
	{  
		y=x*x/4.0; 
		p=b[6]; 
		for(i=5; i>=0; i--) p=p*y+b[i]; 
		p=p/x+TransformativeIntegerBessel1stFunction(1,x)*log(x/2.0); 
	} 
	else 
	{  
		y=2.0/x; 
		p=d[6]; 
		for(i=5; i>=0; i--) p=p*y+d[i]; 
		p=p*exp(-x)/sqrt(x); 
	} 
	if(n==1) return(p); 
	b1=p; 
	y=2.0/x; 
	for(i=1; i 
_Ty IncompleteBetaFunction(_Ty a, _Ty b, _Ty x) 
{ 
	_Ty y; 
 
	if(a<0.0 || FloatEqual(a,0.0)) 
	{  
		cout << "err**a<=0!"; 
		return(-1.0); 
	} 
	if(b<0.0 || FloatEqual(b,0.0)) 
	{  
		cout << "err**b<=0!"; 
		return(-1.0); 
	} 
	if((x<0.0)||(x>1.0)) 
	{ 
		cout << "err**x<0 or x>1 !"; 
		return(1.0e+70); 
	} 
	if(FloatEqual(x,0.0)||FloatEqual(x,1.0)) y=0.0; 
	else 
	{  
		y=a*log(x)+b*log(1.0-x); 
		y=exp(y); 
		y=y*GammaFunction(a+b)/(GammaFunction(a)*GammaFunction(b)); 
	} 
	if(x<(a+1.0)/(a+b+2.0)) 
		y=y*beta(a,b,x)/a; 
	else 
		y=1.0-y*beta(b,a,1.0-x)/b; 
	return(y); 
} 
 
template  
_Ty beta(_Ty a, _Ty b, _Ty x) 
{ 
	_Ty p0(0.0), q0(1.0), p1(1.0), q1(1.0), s1; 
 
	for(int k=1; k<101; k++) 
	{ 
		_Ty d=(a+k)*(a+b+k)*x; 
		d=-d/((a+k+k)*(a+k+k+1.0)); 
		p0=p1+d*p0; 
		q0=q1+d*q0; 
		_Ty s0=p0/q0; 
		d=k*(b-k)*x; 
		d=d/((a+k+k-1.0)*(a+k+k)); 
		p1=p0+d*p1; 
		q1=q0+d*q1; 
		s1=p1/q1; 
		if(Abs(s1-s0) 
_Ty NormalDistributionFunction(_Ty a, _Ty d, _Ty x) 
{ 
	_Ty y; 
 
	if(d<0.0 || FloatEqual(d,0.0)) d=1.0e-10; 
	y=0.5+0.5*ErrorFunction((x-a)/(sqrt(2.0)*d)); 
	return(y); 
} 
//t-分布函数 
 
/*! \fn 
* 函数功能：t-分布函数  
* 算法原理：t-分布函数  
* 输入参数：t：_Ty类型变量 
            n：整型变量。自由度 
* 输出参数：  
* 返 回 值：_Ty类型  
*/ 
 
template  
_Ty tDistributionFunction(_Ty t, int n) 
{ 
	_Ty y; 
 
	if(t<0.0) t = -t; 
	y=1.0-IncompleteBetaFunction(n/2.0, 0.5, n/(n+t*t)); 
	return(y); 
} 
 
//X^2-分布函数 
 
/*! \fn 
* 函数功能：X^2-分布函数  
* 算法原理：X^2-分布函数  
* 输入参数：t：_Ty类型变量 
            n：整型变量。自由度  
* 输出参数：  
* 返 回 值：_Ty类型  
*/ 
 
template  
_Ty X2DistributionFunction(_Ty x, int n) 
{ 
	_Ty y; 
 
	if(x<0.0) x=-x; 
	y=IncompleteGammaFunction(n/2.0,x/2.0); 
	return(y); 
} 
//F-分布函数 
 
/*! \fn 
* 函数功能：F-分布函数  
* 算法原理：F-分布函数  
* 输入参数：t：_Ty类型变量。随机变量F的取值 
            n1：整型变量 自由度 
			n2：整型变量 自由度 
* 输出参数：  
* 返 回 值：_Ty类型 
*/ 
 
template  
_Ty FDistributionFunction(_Ty f, int n1, int n2) 
{ 
	_Ty y; 
	//extern double lbeta(); 
	if(f<0.0) f=-f; 
	y=IncompleteBetaFunction(n2/2.0,n1/2.0,n2/(n2+n1*f)); 
	return(y); 
} 
//正弦积分 
 
/*! \fn 
* 函数功能：正弦积分  
* 算法原理：正弦积分  
* 输入参数：x：_Ty类型变量自变量值x>0  
* 输出参数：  
* 返 回 值：_Ty类型  
*/ 
 
template  
_Ty SineIntegralFunction(_Ty x) 
{ 
	int m(1), i, j; 
	_Ty aa,bb,w,xx; 
	_Ty t[5]={-0.9061798459,-0.5384693101,0.0,0.5384693101,0.9061798459}; 
	_Ty c[5]={0.2369268851,0.4786286705,0.5688888889, 
		0.4786286705,0.2369268851}; 
 
	if(FloatEqual(x,0.0)) return(0.0); 
	_Ty h=Abs(x);  
	_Ty s=Abs(0.0001*h); 
	_Ty p(1.0e+35), ep(FLOATERROR), g(0.0); 
	while((ep>0.0000001 || FloatEqual(ep,0.0000001)) && (Abs(h)>s)) 
	{  
		g=0.0; 
		for(i=1;i<=m;i++) 
		{  
			aa=(i-1.0)*h;  
			bb=i*h; 
			w=0.0; 
			for(j=0;j<5;j++) 
			{ 
				xx=((bb-aa)*t[j]+(bb+aa))/2.0; 
				w=w+sin(xx)/xx*c[j]; 
			} 
			g=g+w; 
		} 
		g=g*h/2.0; 
		ep=Abs(g-p)/(1.0+Abs(g)); 
		p=g; 
		m=m+1; 
		h=Abs(x)/m; 
	} 
	return(g); 
} 
 
//余弦积分 
/*! \fn 
* 函数功能：余弦积分  
* 算法原理：余弦积分  
* 输入参数：x：_Ty类型变量自变量值x>0  
* 输出参数：  
* 返 回 值：_Ty类型 
template  
_Ty CosineIntegralFunction(_Ty x); 
*/ 
template  
_Ty CosineIntegralFunction(_Ty x) 
{ 
int m(1), i, j; 
_Ty ep(FLOATERROR),aa,bb,w,xx,g(0.0); 
_Ty t[5]={-0.9061798459,-0.5384693101,0.0,0.5384693101,0.9061798459}; 
_Ty c[5]={0.2369268851,0.4786286705,0.5688888889, 
0.4786286705,0.2369268851}; 
 
if(FloatEqual(x,0.0)) x = 1.0e-35; 
if(x<0.0) x = -x; 
_Ty r = 0.57721566490153286060651; 
_Ty q = r + log(x); 
_Ty h = x;  
_Ty s = Abs(0.0001*h); 
_Ty p=1.0e+35; 
while((ep>0.0000001||FloatEqual(ep,0.0000001))&&(Abs(h)>s)) 
{  
g=0.0; 
for(i=1;i<=m;i++) 
{  
aa=(i-1.0)*h;  
bb=i*h; 
w=0.0; 
for(j=0;j<5;j++) 
{ 
xx=((bb-aa)*t[j]+(bb+aa))/2.0; 
w=w+(1.0-cos(xx))/xx*c[j]; 
} 
g=g+w; 
} 
g=g*h/2.0; 
ep=Abs(g-p)/(1.0+Abs(g)); 
p=g; m=m+1; h=x/m; 
} 
g=q-g; 
return(g); 
} 
 
//指数积分 
/*! \fn 
* 函数功能：指数积分  
* 算法原理：指数积分  
* 输入参数：x：_Ty类型变量自变量值x>0  
* 输出参数：  
* 返 回 值：_Ty类型  
*/ 
 
// template  
// _Ty ExponentIntegralFunction(_Ty x); 
template  
_Ty ExponentIntegralFunction(_Ty x) 
{ 
	int m(1), i, j; 
	_Ty ep(FLOATERROR),aa,bb,w,xx,g(0.0); 
	_Ty t[5]={-0.9061798459,-0.5384693101,0.0,0.5384693101,0.9061798459}; 
	_Ty c[5]={0.2369268851,0.4786286705,0.5688888889, 
		0.4786286705,0.2369268851}; 
 
	if(FloatEqual(x,0.0)) x=1.0e-10; 
	if(x<0.0) x=-x; 
	_Ty r=0.57721566490153286060651; 
	_Ty q=r+log(x); 
	_Ty h=x; 
	_Ty s=Abs(0.0001*h); 
	_Ty p=1.0e+35; 
	while((ep>0.0000001||FloatEqual(ep,0.0000001))&&(Abs(h)>s)) 
	{  
		g=0.0; 
		for(i=1;i<=m;i++) 
		{ 
			aa=(i-1.0)*h; bb=i*h; 
			w=0.0; 
			for(j=0;j<5;j++) 
			{  
				xx=((bb-aa)*t[j]+(bb+aa))/2.0; 
				w=w+(exp(-xx)-1.0)/xx*c[j]; 
			} 
			g=g+w; 
		} 
		g=g*h/2.0; 
		ep=Abs(g-p)/(1.0+Abs(g)); 
		p=g; m=m+1; h=x/m; 
	} 
	g=q+g; 
	return(g); 
} 
 
//第一类椭圆积分 
 
/*! \fn 
* 函数功能：第一类椭圆积分  
* 算法原理：第一类椭圆积分  
* 输入参数：_Ty k, _Ty f  
* 输出参数：  
* 返 回 值：_Ty类型  
*/ 
 
template  
_Ty Ellipse1stIntegral(_Ty k, _Ty f) 
{ 
	int n(0); 
	_Ty pi(3.1415926), e(1.0), ff; 
	//_Ty fk(); 
	if(k<0.0) k=-k; 
	if(k>1.0) k=1.0/k; 
	_Ty y = Abs(f); 
 
	while(y>pi||FloatEqual(y,pi)) 
	{ 
		n=n+1; 
		y=y-pi; 
	} 
 
	if(y>pi/2.0||FloatEqual(y,pi/2.0)) 
	{ 
		n=n+1; 
		e=-e;  
		y=pi-y; 
	} 
	if(n==0) 
		ff=fk(k,y); 
	else 
	{  
		ff=fk(k,pi/2.0); 
		ff=2.0*n*ff+e*fk(k,y); 
	} 
	if(f<0.0) ff=-ff; 
	return(ff); 
} 
 
//第一类椭圆积分辅助函数 
template  
_Ty fk(_Ty k, _Ty f) 
{  
	int m(1), i, j; 
	_Ty p(1.0e+35),ep(FLOATERROR),aa,bb,w,xx,g(0.0),q; 
	_Ty t[5]={-0.9061798459,-0.5384693101,0.0,0.5384693101,0.9061798459}; 
	_Ty c[5]={0.2369268851,0.4786286705,0.5688888889, 
		0.4786286705,0.2369268851}; 
 
	_Ty h=Abs(f);  
	_Ty s=Abs(0.0001*h); 
	//p=1.0e+35; ep=0.000001; 
	while((ep>0.0000001||FloatEqual(ep,0.0000001))&&(Abs(h)>s)) 
	{  
		g=0.0; 
		for(i=1;i<=m;i++) 
		{ 
			aa=(i-1.0)*h;  
			bb=i*h; 
			w=0.0; 
			for(j=0;j<5;j++) 
			{  
				xx=((bb-aa)*t[j]+(bb+aa))/2.0; 
				q=sqrt(1.0-k*k*sin(xx)*sin(xx)); 
				w=w+c[j]/q; 
			} 
			g=g+w; 
		} 
		g=g*h/2.0; 
		ep=Abs(g-p)/(1.0+Abs(g)); 
		p=g;  
		m=m+m;  
		h=0.5*h; 
	} 
	return(g); 
} 
 
//第二类椭圆积分 
/*! \fn 
* 函数功能：第二类椭圆积分  
* 算法原理：第二类椭圆积分  
* 输入参数：_Ty k, _Ty f  
* 输出参数：  
* 返 回 值：_Ty类型  
*/ 
template  
_Ty Ellipse2ndIntegral(_Ty k, _Ty f) 
{ 
	int n(0); 
	_Ty pi(3.1415926), e(1.0), ff; 
	//_Ty ek(); 
	if(k<0.0) k=-k; 
	if(k>1.0) k=1.0/k; 
	_Ty y = Abs(f); 
 
	while(y>pi||FloatEqual(y,pi)) 
	{ 
		n=n+1;  
		y=y-pi; 
	} 
 
	if(y>pi/2.0||FloatEqual(y,pi/2.0)) 
	{ 
		n=n+1; 
		e=-e;  
		y=pi-y; 
	} 
	if(n==0) 
		ff=ek(k,y); 
	else 
	{  
		ff=ek(k,pi/2.0); 
		ff=2.0*n*ff+e*ek(k,y); 
	} 
	if(f<0.0) ff=-ff; 
	return(ff); 
} 
 
//第二类椭圆积分辅助函数 
template  
_Ty ek(_Ty k, _Ty f) 
{	int m(1), i, j; 
_Ty p(1.0e+35),ep(FLOATERROR),aa,bb,w,xx,g(0.0),q; 
_Ty t[5]={-0.9061798459,-0.5384693101,0.0,0.5384693101,0.9061798459}; 
_Ty c[5]={0.2369268851,0.4786286705,0.5688888889, 
0.4786286705,0.2369268851}; 
 
_Ty h=Abs(f); 
_Ty s=Abs(0.0001*h);