bessi.cpp

代码

double bessi(int n, double x)
{
	int j,m;
	double iacc,bigno,bigni,tox,bip,bi,temp;
	double bim;
    iacc = 40;
    bigno = 10000000000.0;
    bigni = 0.0000000001;
    if (n < 2)cout<<"bad argument n in bessi";
    tox = 2.0 / x;
    bip = 0.0;
    bi = 1.0;
    temp = 0.0;
    m = 2 * ((n + int(sqrt(iacc * n))));
    for( j = m;j>=1;j--)
	{
        bim = bip + j * tox * bi;
        bip = bi;
        bi = bim;
        if (fabs(bi) > bigno) 
		{
            temp = temp * bigni;
            bi = bi * bigni;
            bip = bip * bigni;
        }
        if( j == n) temp = bip;
	} 
    temp = temp * bessi0(x) / bi;
	return temp;
}