?? specialfunction.inl

?? 主要是關于一些數值計算中所要用到的頭文件,希望對大家有所幫助
?? INL
?? 第 1 頁 / 共 2 頁
字號:
上一頁 12
			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);
}

//變形第二類整數階貝塞爾函數
template <class _Ty>
_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};
    _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<n; i++)
    { 
		p=b0+i*y*b1;
		b0=b1; 
		b1=p;
	}
    return(p);
}

//不完全貝塔函數
template <class _Ty>
_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 <class _Ty>
_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)<Abs(s1)*1.0e-07)	return(s1);
    }
    cout << "a or b too big !";
    return(s1);
}

//正態分布函數
template <class _Ty>
_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-分布函數
template <class _Ty>
_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-分布函數
template <class _Ty>
_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-分布函數
template <class _Ty>
_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);
}

//正弦積分
template <class _Ty>
_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);
}

//余弦積分
template <class _Ty>
_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);
}

//指數積分
template <class _Ty>
_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);
}

//第一類橢圓積分
template <class _Ty>
_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 <class _Ty>
_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);
}

//第二類橢圓積分
template <class _Ty>
_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 <class _Ty>
_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);
    
    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+q*c[j];
            }
            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);
}

#endif		// _SPECIALFUNCTION_INL
上一頁 12
?? 文件大小 76 K
?? 上傳用戶 qqwoshi
?? 所屬分類數學計算
??? 相關標簽

#數值計算 #頭文件 #家
?? 快捷鍵說明

復制代碼 Ctrl + C
搜索代碼 Ctrl + F
全屏模式 F11
切換主題 Ctrl + Shift + D
顯示快捷鍵 ?
增大字號 Ctrl + =
減小字號 Ctrl + -
亚洲欧美第一页_禁久久精品乱码_粉嫩av一区二区三区免费野_久草精品视频

?? specialfunction.inl

?? 快捷鍵說明
