bessel.cpp

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

// Bessel.cpp : Defines the entry point for the DLL application.
//

#include "stdafx.h"
#include "ap.h"
#include "Bessel.h"

#include "gammaf.h"
#include "hyper1f1.h"
#include "airyf.h"

BOOL APIENTRY DllMain( HANDLE hModule, 
                       DWORD  ul_reason_for_call, 
                       LPVOID lpReserved
					 )
{
    switch (ul_reason_for_call)
	{
		case DLL_PROCESS_ATTACH:
		case DLL_THREAD_ATTACH:
		case DLL_THREAD_DETACH:
		case DLL_PROCESS_DETACH:
			break;
    }
    return TRUE;
}


/*************************************************************************
Cephes Math Library Release 2.8:  June, 2000
Copyright by Stephen L. Moshier

Contributors:
    * Sergey Bochkanov (ALGLIB project). Translation from C to
      pseudocode.

See subroutines comments for additional copyrights.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:

- Redistributions of source code must retain the above copyright
  notice, this list of conditions and the following disclaimer.

- Redistributions in binary form must reproduce the above copyright
  notice, this list of conditions and the following disclaimer listed
  in this license in the documentation and/or other materials
  provided with the distribution.

- Neither the name of the copyright holders nor the names of its
  contributors may be used to endorse or promote products derived from
  this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*************************************************************************/

void besselmfirstcheb(double c, double& b0, double& b1, double& b2);
void besselmnextcheb(double x, double c, double& b0, double& b1, double& b2);
void besselm1firstcheb(double c, double& b0, double& b1, double& b2);
void besselm1nextcheb(double x, double c, double& b0, double& b1, double& b2);
void besselasympt0(double x, double& pzero, double& qzero);
void besselasympt1(double x, double& pzero, double& qzero);

/*************************************************************************
Bessel function of order zero

Returns Bessel function of order zero of the argument.

The domain is divided into the intervals [0, 5] and
(5, infinity). In the first interval the following rational
approximation is used:


       2         2
(w - r  ) (w - r  ) P (w) / Q (w)
      1         2    3       8

           2
where w = x  and the two r's are zeros of the function.

In the second interval, the Hankel asymptotic expansion
is employed with two rational functions of degree 6/6
and 7/7.

ACCURACY:

                     Absolute error:
arithmetic   domain     # trials      peak         rms
   IEEE      0, 30       60000       4.2e-16     1.1e-16

Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
*************************************************************************/
BESSEL_API double besselj0(double x)
{
    double result;
    double xsq;
    double nn;
    double pzero;
    double qzero;
    double p1;
    double q1;

    if( x<0 )
    {
        x = -x;
    }
    if( x>8.0 )
    {
        besselasympt0(x, pzero, qzero);
        nn = x-ap::pi()/4;
        result = sqrt(2/ap::pi()/x)*(pzero*cos(nn)-qzero*sin(nn));
        return result;
    }
    xsq = ap::sqr(x);
    p1 = 26857.86856980014981415848441;
    p1 = -40504123.71833132706360663322+xsq*p1;
    p1 = 25071582855.36881945555156435+xsq*p1;
    p1 = -8085222034853.793871199468171+xsq*p1;
    p1 = 1434354939140344.111664316553+xsq*p1;
    p1 = -136762035308817138.6865416609+xsq*p1;
    p1 = 6382059341072356562.289432465+xsq*p1;
    p1 = -117915762910761053603.8440800+xsq*p1;
    p1 = 493378725179413356181.6813446+xsq*p1;
    q1 = 1.0;
    q1 = 1363.063652328970604442810507+xsq*q1;
    q1 = 1114636.098462985378182402543+xsq*q1;
    q1 = 669998767.2982239671814028660+xsq*q1;
    q1 = 312304311494.1213172572469442+xsq*q1;
    q1 = 112775673967979.8507056031594+xsq*q1;
    q1 = 30246356167094626.98627330784+xsq*q1;
    q1 = 5428918384092285160.200195092+xsq*q1;
    q1 = 493378725179413356211.3278438+xsq*q1;
    result = p1/q1;
    return result;
}


/*************************************************************************
Bessel function of order one

Returns Bessel function of order one of the argument.

The domain is divided into the intervals [0, 8] and
(8, infinity). In the first interval a 24 term Chebyshev
expansion is used. In the second, the asymptotic
trigonometric representation is employed using two
rational functions of degree 5/5.

ACCURACY:

                     Absolute error:
arithmetic   domain      # trials      peak         rms
   IEEE      0, 30       30000       2.6e-16     1.1e-16

Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
*************************************************************************/
BESSEL_API double besselj1(double x)
{
    double result;
    double s;
    double xsq;
    double nn;
    double pzero;
    double qzero;
    double p1;
    double q1;

    s = ap::sign(x);
    if( x<0 )
    {
        x = -x;
    }
    if( x>8.0 )
    {
        besselasympt1(x, pzero, qzero);
        nn = x-3*ap::pi()/4;
        result = sqrt(2/ap::pi()/x)*(pzero*cos(nn)-qzero*sin(nn));
        if( s<0 )
        {
            result = -result;
        }
        return result;
    }
    xsq = ap::sqr(x);
    p1 = 2701.122710892323414856790990;
    p1 = -4695753.530642995859767162166+xsq*p1;
    p1 = 3413234182.301700539091292655+xsq*p1;
    p1 = -1322983480332.126453125473247+xsq*p1;
    p1 = 290879526383477.5409737601689+xsq*p1;
    p1 = -35888175699101060.50743641413+xsq*p1;
    p1 = 2316433580634002297.931815435+xsq*p1;
    p1 = -66721065689249162980.20941484+xsq*p1;
    p1 = 581199354001606143928.050809+xsq*p1;
    q1 = 1.0;
    q1 = 1606.931573481487801970916749+xsq*q1;
    q1 = 1501793.594998585505921097578+xsq*q1;
    q1 = 1013863514.358673989967045588+xsq*q1;
    q1 = 524371026216.7649715406728642+xsq*q1;
    q1 = 208166122130760.7351240184229+xsq*q1;
    q1 = 60920613989175217.46105196863+xsq*q1;
    q1 = 11857707121903209998.37113348+xsq*q1;
    q1 = 1162398708003212287858.529400+xsq*q1;
    result = s*x*p1/q1;
    return result;
}


/*************************************************************************
Bessel function of integer order

Returns Bessel function of order n, where n is a
(possibly negative) integer.

The ratio of jn(x) to j0(x) is computed by backward
recurrence.  First the ratio jn/jn-1 is found by a
continued fraction expansion.  Then the recurrence
relating successive orders is applied until j0 or j1 is
reached.

If n = 0 or 1 the routine for j0 or j1 is called
directly.

ACCURACY:

                     Absolute error:
arithmetic   range      # trials      peak         rms
   IEEE      0, 30        5000       4.4e-16     7.9e-17


Not suitable for large n or x. Use jv() (fractional order) instead.

Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1987, 2000 by Stephen L. Moshier
*************************************************************************/
BESSEL_API double besseljn(int n, double x)
{
    double result;
    double pkm2;
    double pkm1;
    double pk;
    double xk;
    double r;
    double ans;
    int k;
    int sg;

    if( n<0 )
    {
        n = -n;
        if( n%2==0 )
        {
            sg = 1;
        }
        else
        {
            sg = -1;
        }
    }
    else
    {
        sg = 1;
    }
    if( x<0 )
    {
        if( n%2!=0 )
        {
            sg = -sg;
        }
        x = -x;
    }
    if( n==0 )
    {
        result = sg*besselj0(x);
        return result;
    }
    if( n==1 )
    {
        result = sg*besselj1(x);
        return result;
    }
    if( n==2 )
    {
        if( x==0 )
        {
            result = 0;
        }
        else
        {
            result = sg*(2.0*besselj1(x)/x-besselj0(x));
        }
        return result;
    }
    if( x<ap::machineepsilon )
    {
        result = 0;
        return result;
    }
    k = 53;
    pk = 2*(n+k);
    ans = pk;
    xk = x*x;
    do
    {
        pk = pk-2.0;
        ans = pk-xk/ans;
        k = k-1;
    }
    while(k!=0);
    ans = x/ans;
    pk = 1.0;
    pkm1 = 1.0/ans;
    k = n-1;
    r = 2*k;
    do
    {
        pkm2 = (pkm1*r-pk*x)/x;
        pk = pkm1;
        pkm1 = pkm2;
        r = r-2.0;
        k = k-1;
    }
    while(k!=0);
    if( fabs(pk)>fabs(pkm1) )
    {
        ans = besselj1(x)/pk;
    }
    else
    {
        ans = besselj0(x)/pkm1;
    }
    result = sg*ans;
    return result;
}


/*************************************************************************
Bessel function of the second kind, order zero

Returns Bessel function of the second kind, of order
zero, of the argument.

The domain is divided into the intervals [0, 5] and
(5, infinity). In the first interval a rational approximation
R(x) is employed to compute
  y0(x)  = R(x)  +   2 * log(x) * j0(x) / PI.
Thus a call to j0() is required.

In the second interval, the Hankel asymptotic expansion
is employed with two rational functions of degree 6/6
and 7/7.



ACCURACY:

 Absolute error, when y0(x) < 1; else relative error:

arithmetic   domain     # trials      peak         rms
   IEEE      0, 30       30000       1.3e-15     1.6e-16

Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
*************************************************************************/
BESSEL_API double bessely0(double x)
{
    double result;
    double nn;
    double xsq;
    double pzero;
    double qzero;
    double p4;
    double q4;

    if( x>8.0 )
    {
        besselasympt0(x, pzero, qzero);
        nn = x-ap::pi()/4;
        result = sqrt(2/ap::pi()/x)*(pzero*sin(nn)+qzero*cos(nn));
        return result;
    }
    xsq = ap::sqr(x);
    p4 = -41370.35497933148554125235152;
    p4 = 59152134.65686889654273830069+xsq*p4;
    p4 = -34363712229.79040378171030138+xsq*p4;
    p4 = 10255208596863.94284509167421+xsq*p4;
    p4 = -1648605817185729.473122082537+xsq*p4;
    p4 = 137562431639934407.8571335453+xsq*p4;
    p4 = -5247065581112764941.297350814+xsq*p4;
    p4 = 65874732757195549259.99402049+xsq*p4;
    p4 = -27502866786291095837.01933175+xsq*p4;
    q4 = 1.0;
    q4 = 1282.452772478993804176329391+xsq*q4;
    q4 = 1001702.641288906265666651753+xsq*q4;
    q4 = 579512264.0700729537480087915+xsq*q4;
    q4 = 261306575504.1081249568482092+xsq*q4;
    q4 = 91620380340751.85262489147968+xsq*q4;
    q4 = 23928830434997818.57439356652+xsq*q4;
    q4 = 4192417043410839973.904769661+xsq*q4;
    q4 = 372645883898616588198.9980+xsq*q4;
    result = p4/q4+2/ap::pi()*besselj0(x)*log(x);
    return result;
}


/*************************************************************************
Bessel function of second kind of order one

Returns Bessel function of the second kind of order one
of the argument.

The domain is divided into the intervals [0, 8] and
(8, infinity). In the first interval a 25 term Chebyshev
expansion is used, and a call to j1() is required.
In the second, the asymptotic trigonometric representation
is employed using two rational functions of degree 5/5.

ACCURACY:

                     Absolute error:
arithmetic   domain      # trials      peak         rms
   IEEE      0, 30       30000       1.0e-15     1.3e-16

Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
*************************************************************************/
BESSEL_API double bessely1(double x)
{
    double result;
    double nn;
    double xsq;
    double pzero;