bessel.cpp

基于Visual C++环境的的Bessel函数程序

                sg = -sg;
            }
            v = an;
        }
        if( x<0.0 )
        {
            if( i%2!=0 )
            {
                sg = -sg;
            }
            x = -x;
        }
        if( v==0 )
        {
            result = besselj0(x);
            return result;
        }
        if( v==1 )
        {
            result = sg*besselj1(x);
            return result;
        }
    }
    if( x<0.0&&y!=an )
    {
        ap::ap_error::make_assertion(false);
        result = 0.0;
        return result;
    }
    y = fabs(x);
    if( y<ap::machineepsilon )
    {
        result = 0.0;
        return result;
    }
    k = 3.6*sqrt(y);
    t = 3.6*sqrt(an);
    if( y<t&&an>21.0 )
    {
        result = sg*besseljvs(v, x);
        return result;
    }
    if( an<k&&y>21.0 )
    {
        result = sg*besseljvhankel(v, x);
        return result;
    }
    if( an<500.0 )
    {
        if( nint!=0 )
        {
            k = 0.0;
            q = besseljvrecur(v, x, k, 1);
            if( k==0.0 )
            {
                result = sg*besselj0(x)/q;
                return result;
            }
            if( k==1.0 )
            {
                result = sg*besselj1(x)/q;
                return result;
            }
        }
        if( an>2.0*y||v>=0.0&&v<20.0&&y>6.0&&y<20.0 )
        {
            k = v;
            y = y+an+1.0;
            if( y<30.0 )
            {
                y = 30.0;
            }
            y = v+ap::ifloor(y-v);
            q = besseljvrecur(y, x, k, 0);
            result = sg*besseljvs(y, x)*q;
            return result;
        }
        if( k<=30.0 )
        {
            k = 2.0;
        }
        else
        {
            if( k<90.0 )
            {
                k = 3*k/4;
            }
        }
        if( an>k+3.0 )
        {
            if( v<0.0 )
            {
                k = -k;
            }
            q = v-ap::ifloor(v);
            k = ap::ifloor(k)+q;
            if( v>0.0 )
            {
                q = besseljvrecur(v, x, k, 1);
            }
            else
            {
                t = k;
                k = v;
                q = besseljvrecur(t, x, k, 1);
                k = t;
            }
            if( q==0.0 )
            {
                result = 0.0;
                return result;
            }
        }
        else
        {
            k = v;
            q = 1.0;
        }
        y = fabs(k);
        if( y<26.0 )
        {
            t = (0.0083*y+0.09)*y+12.9;
        }
        else
        {
            t = 0.9*y;
        }
        if( x>t )
        {
            y = besseljvhankel(k, x);
        }
        else
        {
            y = besseljvs(k, x);
        }
        if( v>0.0 )
        {
            y = y/q;
        }
        else
        {
            y = y*q;
        }
    }
    else
    {
        if( v<0.0 )
        {
            ap::ap_error::make_assertion(false);
            result = 0;
            return result;
        }
        t = x/v;
        t = t/v;
        if( t>0.3 )
        {
            y = besseljvhankel(v, x);
        }
        else
        {
            y = besseljnx(v, x);
        }
    }
    result = sg*y;
    return result;
}


/*************************************************************************
Modified Bessel function of noninteger order

Returns modified Bessel function of order v of the
argument.  If x is negative, v must be integer valued.

The function is defined as Iv(x) = Jv( ix ).  It is
here computed in terms of the confluent hypergeometric
function, according to the formula

             v  -x
Iv(x) = (x/2)  e   hyperg( v+0.5, 2v+1, 2x ) / gamma(v+1)

If v is a negative integer, then v is replaced by -v.


ACCURACY:

Tested at random points (v, x), with v between 0 and
30, x between 0 and 28.
                     Relative error:
arithmetic   domain     # trials      peak         rms
   DEC       0,30          2000      3.1e-15     5.4e-16
   IEEE      0,30         10000      1.7e-14     2.7e-15

Accuracy is diminished if v is near a negative integer.

Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1987, 1988, 2000 by Stephen L. Moshier
*************************************************************************/
BESSEL_API double besseliv(double v, double x)
{
    double result;
    int sg;
    double t;
    double ax;

    t = ap::ifloor(v);
    if( v<0.0 )
    {
        if( t==v )
        {
            v = -v;
            t = -t;
        }
    }
    sg = 1;
    if( x<0.0 )
    {
        ap::ap_error::make_assertion(t==v);
        if( v!=2.0*ap::ifloor(v/2.0) )
        {
            sg = -1;
        }
    }
    if( x==0.0 )
    {
        ap::ap_error::make_assertion(v>=0);
        if( v==0.0 )
        {
            result = 1;
            return result;
        }
        result = 0;
        return result;
    }
    ax = fabs(x);
    t = v*log(0.5*ax)-x;
    t = sg*exp(t)/gamma(v+1.0);
    ax = v+0.5;
    result = t*hypergeometric1f1(ax, 2.0*ax, 2.0*x);
    return result;
}


/*************************************************************************
Reduce the order by backward recurrence.
AMS55 #9.1.27 and 9.1.73.
*************************************************************************/
double besseljvrecur(double& n, double x, double& newn, int cancel)
{
    double result;
    double pkm2;
    double pkm1;
    double pk;
    double qkm2;
    double qkm1;
    double k;
    double ans;
    double qk;
    double xk;
    double yk;
    double r;
    double t;
    double kf;
    double big;
    int nflag;
    int ctr;

    big = 1.44115188075855872E+17;
    if( n<0.0 )
    {
        nflag = 1;
    }
    else
    {
        nflag = 0;
    }
    while(true)
    {
        pkm2 = 0.0;
        qkm2 = 1.0;
        pkm1 = x;
        qkm1 = n+n;
        xk = -x*x;
        yk = qkm1;
        ans = 1.0;
        ctr = 0;
        do
        {
            yk = yk+2.0;
            pk = pkm1*yk+pkm2*xk;
            qk = qkm1*yk+qkm2*xk;
            pkm2 = pkm1;
            pkm1 = pk;
            qkm2 = qkm1;
            qkm1 = qk;
            if( qk!=0 )
            {
                r = pk/qk;
            }
            else
            {
                r = 0.0;
            }
            if( r!=0 )
            {
                t = fabs((ans-r)/r);
                ans = r;
            }
            else
            {
                t = 1.0;
            }
            ctr = ctr+1;
            if( ctr>1000 )
            {
                break;
            }
            if( t<ap::machineepsilon )
            {
                break;
            }
            if( fabs(pk)>big )
            {
                pkm2 = pkm2/big;
                pkm1 = pkm1/big;
                qkm2 = qkm2/big;
                qkm1 = qkm1/big;
            }
        }
        while(t>ap::machineepsilon);
        if( nflag>0 )
        {
            if( fabs(ans)<0.125 )
            {
                nflag = -1;
                n = n-1.0;
                continue;
            }
        }
        break;
    }
    kf = newn;
    
    //
    //backward recurrence
    //              2k
    //  J   (x)  =  --- J (x)  -  J   (x)
    //   k-1         x   k         k+1
    //
    pk = 1.0;
    pkm1 = 1.0/ans;
    k = n-1.0;
    r = 2*k;
    do
    {
        pkm2 = (pkm1*r-pk*x)/x;
        pk = pkm1;
        pkm1 = pkm2;
        r = r-2.0;
        k = k-1.0;
    }
    while(k>kf+0.5);
    if( cancel!=0 )
    {
        if( kf>=0.0&&fabs(pk)>fabs(pkm1) )
        {
            k = k+1.0;
            pkm2 = pk;
        }
    }
    newn = k;
    result = pkm2;
    return result;
}


/*************************************************************************
Ascending power series for Jv(x).
AMS55 #9.1.10.
*************************************************************************/
double besseljvs(double n, double x)
{
    double result;
    double sgngam;
    double t;
    double u;
    double y;
    double z;
    double k;
    int ex;

    z = -x*x/4.0;
    u = 1.0;
    y = u;
    k = 1.0;
    t = 1.0;
    while(t>ap::machineepsilon)
    {
        u = u*z/(k*(n+k));
        y = y+u;
        k = k+1.0;
        if( y!=0 )
        {
            t = fabs(u/y);
        }
    }
    t = besselfrexp(0.5*x, ex);
    ex = ap::trunc(ex*n);
    if( ex>-1023&&ex<1023&&n>0.0&&n<171.624376956302725-1.0 )
    {
        t = pow(0.5*x, n)/gamma(n+1.0);
        y = y*t;
    }
    else
    {
        t = n*log(0.5*x)-lngamma(n+1.0, sgngam);
        if( y<0 )
        {
            sgngam = -sgngam;
            y = -y;
        }
        t = t+log(y);
        if( t<-log(ap::maxrealnumber) )
        {
            result = 0;
            return result;
        }
        if( t>log(ap::maxrealnumber) )
        {
            ap::ap_error::make_assertion(false);
            result = ap::maxrealnumber;
            return result;
        }
        y = sgngam*exp(t);
    }
    result = y;
    return result;
}


/*************************************************************************
Hankel's asymptotic expansion
for large x.
AMS55 #9.2.5.