specialfunction.h

数学特殊函数算法库

/*! \file SpecialFunction.h 版权所有 (c) 2002-2008 , 北京迈达斯技术有限公司 
* 功能描述：特殊函数头文件                            \n 
* 编 制 者：赵志刚	 完 成  日 期：2007-6-21 16:21:47 \n 
* 修 改 者：	     最后修改日期： -  - \n 
* 历史记录：V 00.00.00(每次修改升级最后一个数字)  \n 
*/ 
#pragma once 
 
//#include 			//模板类输入输出流标准头文件 
//#include "Comm.h"			//公共头文件 
 
#include "MMathMacro.h" 
 
_MITC_MMATH_BEGIN 
//伽马函数 
 
/*! \fn 
* 函数功能：伽马函数  
* 算法原理：伽马函数  
* 输入参数：x：_Ty，自变量值。  
* 输出参数：  
* 返 回 值：_Ty类型  
*/ 
 
template  
_Ty GammaFunction(_Ty x) 
{ 
	_Ty t; 
	_Ty a[11] =  
	{ 
		0.0000677106, -0.0003442342,  0.0015397681, 
		-0.0024467480,  0.0109736958, -0.0002109075, 
		0.0742379071,  0.0815782188,  0.4118402518, 
		0.4227843370,  1.0 
	}; 
 
	if(x<0.0 || FloatEqual(x,0.0)) 
	{ 
		cout << "Error **x<=0!" << endl;  
		return(-1.0); 
	} 
	_Ty y = x; 
	if(y<1.0 || FloatEqual(y,1.0)) 
	{ 
		t=1.0/(y*(y+1.0));  
		y=y+2.0; 
	} 
	else 
		if(y<2.0 || FloatEqual(y,2.0)) 
		{ 
			t=1.0/y;  
			y=y+1.0; 
		} 
		else 
			if(y<3.0 || FloatEqual(y,3.0)) t=1.0; 
			else 
			{ 
				t=1.0; 
				while(y>3.0) 
				{ 
					y=y-1.0; 
					t=t*y; 
				} 
			} 
			_Ty s=a[0];  
			_Ty u = y - 2.0; 
			for(int i=1; i<11; i++)	s = s * u + a[i]; 
			s = s * t; 
			return(s); 
} 
 
//不完全伽马函数 
 
/*! \fn 
* 函数功能：不完全伽马函数  
* 算法原理：不完全伽马函数  
* 输入参数：a：_Ty类型类型变量 
            x：_Ty类型类型自变量 
* 输出参数：  
* 返 回 值：_Ty类型  
*/ 
 
template  
_Ty IncompleteGammaFunction(_Ty a, _Ty x) 
{ 
	int n; 
	_Ty p,d,s,s1,p0,q0,p1,q1; 
	if((a<0.0 || FloatEqual(a,0.0)) || (x<0.0)) 
	{  
		if(a<0.0 || FloatEqual(a,0.0)) cout <<"Error **a<=0!" <1.0e+35) return(1.0); 
	_Ty q = log(x);  
	q = a * q;  
	_Ty qq = exp(q); 
	if(x<1.0+a) 
	{ 
		p=a;  
		d=1.0/a;  
		s=d; 
		for(n=1; n<101; n++) 
		{ 
			p=1.0+p; 
			d=d*x/p; 
			s=s+d; 
			if(Abs(d) 

_Ty ErrorFunction(_Ty x) 
{ 
	_Ty y; 
 
	if(x>0.0 || FloatEqual(x,0.0)) 
		y=IncompleteGammaFunction(0.5,x*x); 
	else 
		y=-IncompleteGammaFunction(0.5,x*x); 
	return(y); 
} 
 
//第一类整数阶贝塞尔函数 
 
/*! \fn 
* 函数功能：第一类整数阶贝塞尔函数  
* 算法原理：第一类整数阶贝塞尔函数  
* 输入参数：n：整型变量。 
            x：_Ty类型变量。自变量值 
* 输出参数：  
* 返 回 值：_Ty类型  
*/ 
 
template  
_Ty IntegerBessel1stFunction(int n, _Ty x) 
{ 
	int i,m; 
	_Ty y,z,p,q,s,b0,b1; 
	_Ty a[6]={ 57568490574.0,-13362590354.0, 
		651619640.7,-11214424.18,77392.33017,-184.9052456}; 
	_Ty b[6]={ 57568490411.0,1029532985.0, 
		9494680.718,59272.64853,267.8532712,1.0}; 
	_Ty c[6]={ 72362614232.0,-7895059235.0, 
		242396853.1,-2972611.439,15704.4826,-30.16036606}; 
	_Ty d[6]={ 144725228443.0,2300535178.0, 
		18583304.74,99447.43394,376.9991397,1.0}; 
	_Ty e[5]={ 1.0,-0.1098628627e-02, 
		0.2734510407e-04,-0.2073370639e-05,0.2093887211e-06}; 
	_Ty f[5]={ -0.1562499995e-01, 
		0.1430488765e-03,-0.6911147651e-05, 
		0.7621095161e-06,-0.934935152e-07}; 
	_Ty g[5]={ 1.0,0.183105e-02, 
		-0.3516396496e-04,0.2457520174e-05, -0.240337019e-06}; 
	_Ty h[5]={ 0.4687499995e-01, 
		-0.2002690873e-03,0.8449199096e-05, 
		-0.88228987e-06,0.105787412e-06}; 
	_Ty t=Abs(x); 
	if(n<0) n=-n; 
	if(n!=1) 
	{  
		if(t<8.0) 
		{ 
			y=t*t; 
			p=a[5]; 
			q=b[5]; 
			for(i=4; i>=0; i--) 
			{  
				p=p*y+a[i]; 
				q=q*y+b[i]; 
			} 
			p=p/q; 
		} 
		else 
		{  
			z=8.0/t; 
			y=z*z; 
			p=e[4]; 
			q=f[4]; 
			for(i=3; i>=0; i--) 
			{ 
				p=p*y+e[i]; 
				q=q*y+f[i]; 
			} 
			s=t-0.785398164; 
			p=p*cos(s)-z*q*sin(s); 
			p=p*sqrt(0.636619772/t); 
		} 
	} 
	if(n==0) return(p); 
	b0=p; 
	if(t<8.0) 
	{  
		y=t*t;  
		p=c[5]; 
		q=d[5]; 
		for(i=4; i>=0; i--) 
		{ 
			p=p*y+c[i]; 
			q=q*y+d[i]; 
		} 
		p=x*p/q; 
	} 
	else 
	{ 
		z=8.0/t;  
		y=z*z; 
		p=g[4]; 
		q=h[4]; 
		for(i=3; i>=0; i--) 
		{ 
			p=p*y+g[i]; 
			q=q*y+h[i]; 
		} 
		s=t-2.356194491; 
		p=p*cos(s)-z*q*sin(s); 
		p=p*x*sqrt(0.636619772/t)/t; 
	} 
	if(n==1) return(p); 
	b1=p; 
	if(x==0.0) return(0.0); 
	s=2.0/t; 
	if(t>1.0*n) 
	{ 
		if(x<0.0) b1=-b1; 
		for(i=1; i=0; i--) 
		{  
			t=s*(i+1)*b0-b1; 
			b1=b0; 
			b0=t; 
			if(Abs(b0)>1.0e+10) 
			{ 
				b0=b0*1.0e-10;  
				b1=b1*1.0e-10; 
				p=p*1.0e-10;  
				q=q*1.0e-10; 
			} 
			if((i+2)%2==0) q=q+b0; 
			if((i+1)==n) p=b1; 
		} 
		q=2.0*q-b0; 
		p=p/q; 
	} 
	if((x<0.0)&&(n%2==1)) p=-p; 
	return(p); 
} 
 
//第二类整数阶贝塞尔函数 
 
/*! \fn 
* 函数功能：第二类整数阶贝塞尔函数  
* 算法原理：第二类整数阶贝塞尔函数  
* 输入参数：n：整型变量。 
            x：_Ty类型变量。自变量值  
* 输出参数：  
* 返 回 值：_Ty类型  
*/ 
 
template  
_Ty IntegerBessel2ndFunction(int n, _Ty x) 
{ 
	int i; 
	_Ty y,z,p,q,s,b0,b1; 
 
	_Ty a[6]={ -2.957821389e+9,7.062834065e+9, 
		-5.123598036e+8,1.087988129e+7,-8.632792757e+4, 
		2.284622733e+2}; 
	_Ty b[6]={ 4.0076544269e+10,7.452499648e+8, 
		7.189466438e+6,4.74472647e+4,2.261030244e+2,1.0}; 
	_Ty c[6]={ -4.900604943e+12,1.27527439e+12, 
		-5.153438139e+10,7.349264551e+8,-4.237922726e+6, 
		8.511937935e+3}; 
	_Ty d[7]={ 2.49958057e+13,4.244419664e+11, 
		3.733650367e+9,2.245904002e+7,1.02042605e+5, 
		3.549632885e+2,1.0}; 
	_Ty e[5]={ 1.0,-0.1098628627e-02, 
		0.2734510407e-04,-0.2073370639e-05, 
		0.2093887211e-06}; 
	_Ty f[5]={ -0.1562499995e-01, 
		0.1430488765e-03,-0.6911147651e-05, 
		0.7621095161e-06,-0.934935152e-07}; 
	_Ty g[5]={ 1.0,0.183105e-02, 
		-0.3516396496e-04,0.2457520174e-05, 
		-0.240337019e-06}; 
	_Ty h[5]={ 0.4687499995e-01, 
		-0.2002690873e-03,0.8449199096e-05, 
		-0.88228987e-06,0.105787412e-06}; 
	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<8.0) 
		{ 
			y=x*x;  
			p=a[5];  
			q=b[5]; 
			for(i=4; i>=0; i--) 
			{ 
				p=p*y+a[i]; 
				q=q*y+b[i];} 
			p=p/q+0.636619772*IntegerBessel1stFunction(0,x)*log(x); 
		} 
		else 
		{  
			z=8.0/x; 
			y=z*z; 
			p=e[4]; 
			q=f[4]; 
			for(i=3; i>=0; i--) 
			{ 
				p=p*y+e[i]; 
				q=q*y+f[i]; 
			} 
			s=x-0.785398164; 
			p=p*sin(s)+z*q*cos(s); 
			p=p*sqrt(0.636619772/x); 
		} 
	} 
	if(n==0) return(p); 
	b0=p; 
	if(x<8.0) 
	{ 
		y=x*x;  
		p=c[5]; 
		q=d[6]; 
		for(i=4; i>=0; i--) 
		{ 
			p=p*y+c[i]; 
			q=q*y+d[i+1]; 
		} 
		q=q*y+d[0]; 
		p=x*p/q+0.636619772*(IntegerBessel1stFunction(1,x)*log(x)-1.0/x);; 
	} 
	else 
	{  
		z=8.0/x; 
		y=z*z; 
		p=g[4];  
		q=h[4]; 
		for(i=3; i>=0; i--) 
		{ 
			p=p*y+g[i]; 
			q=q*y+h[i]; 
		} 
		s=x-2.356194491; 
		p=p*sin(s)+z*q*cos(s); 
		p=p*sqrt(0.636619772/x); 
	} 
	if(n==1) return(p); 
	b1=p; 
	s=2.0/x; 
	for(i=1; i 
_Ty TransformativeIntegerBessel1stFunction(int n,_Ty x) 
{ 
	int i,m; 
	_Ty t,y,p,b0,b1,q; 
	_Ty a[7]={ 1.0,3.5156229,3.0899424,1.2067492, 
		0.2659732,0.0360768,0.0045813}; 
	_Ty b[7]={ 0.5,0.87890594,0.51498869, 
		0.15084934,0.02658773,0.00301532,0.00032411}; 
	_Ty c[9]={ 0.39894228,0.01328592,0.00225319, 
		-0.00157565,0.00916281,-0.02057706, 
		0.02635537,-0.01647633,0.00392377}; 
	_Ty d[9]={ 0.39894228,-0.03988024,-0.00362018, 
		0.00163801,-0.01031555,0.02282967, 
		-0.02895312,0.01787654,-0.00420059}; 
	if(n<0) n=-n; 
	t=Abs(x); 
	if(n!=1) 
	{  
		if(t<3.75) 
		{ 
			y=(x/3.75)*(x/3.75); 
			p=a[6]; 
			for(i=5; i>=0; i--)	p=p*y+a[i]; 
		} 
		else 
		{  
			y=3.75/t; 
			p=c[8]; 
			for(i=7; i>=0; i--)	p=p*y+c[i]; 
			p=p*exp(t)/sqrt(t); 
		} 
	} 
	if(n==0) return(p); 
	q=p; 
	if(t<3.75) 
	{ 
		y=(x/3.75)*(x/3.75);  
		p=b[6]; 
		for(i=5; i>=0; i--) p=p*y+b[i]; 
		p=p*t; 
	} 
	else 
	{  
		y=3.75/t; 
		p=d[8]; 
		for(i=7; i>=0; i--) p=p*y+d[i]; 
		p=p*exp(t)/sqrt(t); 
	} 
	if(x<0.0) p=-p; 
	if(n==1) return(p); 
	if(FloatEqual(x,0.0)) return(0.0); 
	y=2.0/t; 
	t = b0 = 0.0; 
	b1=1.0; 
	m=n+(int)sqrt(40.0*n); 
	m=2*m; 
	for(i=m; i>0; i--) 
	{  
		p=b0+i*y*b1;  
		b0=b1; 
		b1=p; 
		if(Abs(b1)>1.0e+10) 
		{ 
			t=t*1.0e-10;  
			b0=b0*1.0e-10; 
			b1=b1*1.0e-10; 
		} 
		if(i==n) t=b0; 
	} 
	p=t*q/b1; 
	if((x<0.0)&&(n%2==1)) p=-p; 
	return(p); 
} 
 
//变形第二类整数阶贝塞尔函数 
 
/*! \fn 
* 函数功能：变形第二类整数阶贝塞尔函数  
* 算法原理：变形第二类整数阶贝塞尔函数  
* 输入参数：n：整型变量。 
            x：_Ty类型变量。自变量值  
* 输出参数：  
* 返 回 值：_Ty类型  
*/ 
 
 
template  
_Ty TransformativeIntegerBessel2ndFunction(int n, _Ty x) 
{ 
	int i; 
	_Ty y,p,b0,b1; 
 
	_Ty a[7]={ -0.57721566,0.4227842,0.23069756, 
		0.0348859,0.00262698,0.0001075,0.0000074}; 
	_Ty b[7]={ 1.0,0.15443144,-0.67278579, 
		-0.18156897,-0.01919402,-0.00110404,-0.00004686};