📄 bessel.cpp
字号:
y = 8.0/x-2.0;
besselm1firstcheb(-5.75674448366501715755E-18, b0, b1, b2);
besselm1nextcheb(y, 1.79405087314755922667E-17, b0, b1, b2);
besselm1nextcheb(y, -5.68946255844285935196E-17, b0, b1, b2);
besselm1nextcheb(y, 1.83809354436663880070E-16, b0, b1, b2);
besselm1nextcheb(y, -6.05704724837331885336E-16, b0, b1, b2);
besselm1nextcheb(y, 2.03870316562433424052E-15, b0, b1, b2);
besselm1nextcheb(y, -7.01983709041831346144E-15, b0, b1, b2);
besselm1nextcheb(y, 2.47715442448130437068E-14, b0, b1, b2);
besselm1nextcheb(y, -8.97670518232499435011E-14, b0, b1, b2);
besselm1nextcheb(y, 3.34841966607842919884E-13, b0, b1, b2);
besselm1nextcheb(y, -1.28917396095102890680E-12, b0, b1, b2);
besselm1nextcheb(y, 5.13963967348173025100E-12, b0, b1, b2);
besselm1nextcheb(y, -2.12996783842756842877E-11, b0, b1, b2);
besselm1nextcheb(y, 9.21831518760500529508E-11, b0, b1, b2);
besselm1nextcheb(y, -4.19035475934189648750E-10, b0, b1, b2);
besselm1nextcheb(y, 2.01504975519703286596E-9, b0, b1, b2);
besselm1nextcheb(y, -1.03457624656780970260E-8, b0, b1, b2);
besselm1nextcheb(y, 5.74108412545004946722E-8, b0, b1, b2);
besselm1nextcheb(y, -3.50196060308781257119E-7, b0, b1, b2);
besselm1nextcheb(y, 2.40648494783721712015E-6, b0, b1, b2);
besselm1nextcheb(y, -1.93619797416608296024E-5, b0, b1, b2);
besselm1nextcheb(y, 1.95215518471351631108E-4, b0, b1, b2);
besselm1nextcheb(y, -2.85781685962277938680E-3, b0, b1, b2);
besselm1nextcheb(y, 1.03923736576817238437E-1, b0, b1, b2);
besselm1nextcheb(y, 2.72062619048444266945E0, b0, b1, b2);
v = 0.5*(b0-b2);
result = exp(-x)*v/sqrt(x);
}
return result;
}
/*************************************************************************
Modified Bessel function, second kind, integer order
Returns modified Bessel function of the second kind
of order n of the argument.
The range is partitioned into the two intervals [0,9.55] and
(9.55, infinity). An ascending power series is used in the
low range, and an asymptotic expansion in the high range.
ACCURACY:
Relative error:
arithmetic domain # trials peak rms
IEEE 0,30 90000 1.8e-8 3.0e-10
Error is high only near the crossover point x = 9.55
between the two expansions used.
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1987, 1988, 2000 by Stephen L. Moshier
*************************************************************************/
BESSEL_API double besselkn(int nn, double x)
{
double result;
double k;
double kf;
double nk1f;
double nkf;
double zn;
double t;
double s;
double z0;
double z;
double ans;
double fn;
double pn;
double pk;
double zmn;
double tlg;
double tox;
int i;
int n;
double eul;
eul = 5.772156649015328606065e-1;
if( nn<0 )
{
n = -nn;
}
else
{
n = nn;
}
ap::ap_error::make_assertion(n<=31);
ap::ap_error::make_assertion(x>0);
if( x<=9.55 )
{
ans = 0.0;
z0 = 0.25*x*x;
fn = 1.0;
pn = 0.0;
zmn = 1.0;
tox = 2.0/x;
if( n>0 )
{
pn = -eul;
k = 1.0;
for(i = 1; i <= n-1; i++)
{
pn = pn+1.0/k;
k = k+1.0;
fn = fn*k;
}
zmn = tox;
if( n==1 )
{
ans = 1.0/x;
}
else
{
nk1f = fn/n;
kf = 1.0;
s = nk1f;
z = -z0;
zn = 1.0;
for(i = 1; i <= n-1; i++)
{
nk1f = nk1f/(n-i);
kf = kf*i;
zn = zn*z;
t = nk1f*zn/kf;
s = s+t;
ap::ap_error::make_assertion(ap::maxrealnumber-fabs(t)>fabs(s));
ap::ap_error::make_assertion(!(tox>1.0&&ap::maxrealnumber/tox<zmn));
zmn = zmn*tox;
}
s = s*0.5;
t = fabs(s);
ap::ap_error::make_assertion(!(zmn>1.0&&ap::maxrealnumber/zmn<t));
ap::ap_error::make_assertion(!(t>1.0&&ap::maxrealnumber/t<zmn));
ans = s*zmn;
}
}
tlg = 2.0*log(0.5*x);
pk = -eul;
if( n==0 )
{
pn = pk;
t = 1.0;
}
else
{
pn = pn+1.0/n;
t = 1.0/fn;
}
s = (pk+pn-tlg)*t;
k = 1.0;
do
{
t = t*(z0/(k*(k+n)));
pk = pk+1.0/k;
pn = pn+1.0/(k+n);
s = s+(pk+pn-tlg)*t;
k = k+1.0;
}
while(fabs(t/s)>ap::machineepsilon);
s = 0.5*s/zmn;
if( n%2!=0 )
{
s = -s;
}
ans = ans+s;
result = ans;
return result;
}
if( x>log(ap::maxrealnumber) )
{
result = 0;
return result;
}
k = n;
pn = 4.0*k*k;
pk = 1.0;
z0 = 8.0*x;
fn = 1.0;
t = 1.0;
s = t;
nkf = ap::maxrealnumber;
i = 0;
do
{
z = pn-pk*pk;
t = t*z/(fn*z0);
nk1f = fabs(t);
if( i>=n&&nk1f>nkf )
{
break;
}
nkf = nk1f;
s = s+t;
fn = fn+1.0;
pk = pk+2.0;
i = i+1;
}
while(fabs(t/s)>ap::machineepsilon);
result = exp(-x)*sqrt(ap::pi()/(2.0*x))*s;
return result;
}
/*************************************************************************
Internal subroutine
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1987, 2000 by Stephen L. Moshier
*************************************************************************/
void besselmfirstcheb(double c, double& b0, double& b1, double& b2)
{
b0 = c;
b1 = 0.0;
b2 = 0.0;
}
/*************************************************************************
Internal subroutine
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1987, 2000 by Stephen L. Moshier
*************************************************************************/
void besselmnextcheb(double x, double c, double& b0, double& b1, double& b2)
{
b2 = b1;
b1 = b0;
b0 = x*b1-b2+c;
}
/*************************************************************************
Internal subroutine
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1987, 2000 by Stephen L. Moshier
*************************************************************************/
void besselm1firstcheb(double c, double& b0, double& b1, double& b2)
{
b0 = c;
b1 = 0.0;
b2 = 0.0;
}
/*************************************************************************
Internal subroutine
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1987, 2000 by Stephen L. Moshier
*************************************************************************/
void besselm1nextcheb(double x, double c, double& b0, double& b1, double& b2)
{
b2 = b1;
b1 = b0;
b0 = x*b1-b2+c;
}
void besselasympt0(double x, double& pzero, double& qzero)
{
double xsq;
double p2;
double q2;
double p3;
double q3;
xsq = 64.0/(x*x);
p2 = 0.0;
p2 = 2485.271928957404011288128951+xsq*p2;
p2 = 153982.6532623911470917825993+xsq*p2;
p2 = 2016135.283049983642487182349+xsq*p2;
p2 = 8413041.456550439208464315611+xsq*p2;
p2 = 12332384.76817638145232406055+xsq*p2;
p2 = 5393485.083869438325262122897+xsq*p2;
q2 = 1.0;
q2 = 2615.700736920839685159081813+xsq*q2;
q2 = 156001.7276940030940592769933+xsq*q2;
q2 = 2025066.801570134013891035236+xsq*q2;
q2 = 8426449.050629797331554404810+xsq*q2;
q2 = 12338310.22786324960844856182+xsq*q2;
q2 = 5393485.083869438325560444960+xsq*q2;
p3 = -0.0;
p3 = -4.887199395841261531199129300+xsq*p3;
p3 = -226.2630641933704113967255053+xsq*p3;
p3 = -2365.956170779108192723612816+xsq*p3;
p3 = -8239.066313485606568803548860+xsq*p3;
p3 = -10381.41698748464093880530341+xsq*p3;
p3 = -3984.617357595222463506790588+xsq*p3;
q3 = 1.0;
q3 = 408.7714673983499223402830260+xsq*q3;
q3 = 15704.89191515395519392882766+xsq*q3;
q3 = 156021.3206679291652539287109+xsq*q3;
q3 = 533291.3634216897168722255057+xsq*q3;
q3 = 666745.4239319826986004038103+xsq*q3;
q3 = 255015.5108860942382983170882+xsq*q3;
pzero = p2/q2;
qzero = 8*p3/q3/x;
}
void besselasympt1(double x, double& pzero, double& qzero)
{
double xsq;
double p2;
double q2;
double p3;
double q3;
xsq = 64.0/(x*x);
p2 = -1611.616644324610116477412898;
p2 = -109824.0554345934672737413139+xsq*p2;
p2 = -1523529.351181137383255105722+xsq*p2;
p2 = -6603373.248364939109255245434+xsq*p2;
p2 = -9942246.505077641195658377899+xsq*p2;
p2 = -4435757.816794127857114720794+xsq*p2;
q2 = 1.0;
q2 = -1455.009440190496182453565068+xsq*q2;
q2 = -107263.8599110382011903063867+xsq*q2;
q2 = -1511809.506634160881644546358+xsq*q2;
q2 = -6585339.479723087072826915069+xsq*q2;
q2 = -9934124.389934585658967556309+xsq*q2;
q2 = -4435757.816794127856828016962+xsq*q2;
p3 = 35.26513384663603218592175580;
p3 = 1706.375429020768002061283546+xsq*p3;
p3 = 18494.26287322386679652009819+xsq*p3;
p3 = 66178.83658127083517939992166+xsq*p3;
p3 = 85145.16067533570196555001171+xsq*p3;
p3 = 33220.91340985722351859704442+xsq*p3;
q3 = 1.0;
q3 = 863.8367769604990967475517183+xsq*q3;
q3 = 37890.22974577220264142952256+xsq*q3;
q3 = 400294.4358226697511708610813+xsq*q3;
q3 = 1419460.669603720892855755253+xsq*q3;
q3 = 1819458.042243997298924553839+xsq*q3;
q3 = 708712.8194102874357377502472+xsq*q3;
pzero = p2/q2;
qzero = 8*p3/q3/x;
}
/*************************************************************************
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.
*************************************************************************/
double besseljvrecur(double& n, double x, double& newn, int cancel);
double besseljvs(double n, double x);
double besseljvhankel(double n, double x);
double besseljnx(double n, double x);
double besseljnt(double n, double x);
double besseljvcbrt(double x);
double besselfrexp(double x, int& e);
/*************************************************************************
Bessel function of noninteger order
Returns Bessel function of order v of the argument,
where v is real. Negative x is allowed if v is an integer.
Several expansions are included: the ascending power
series, the Hankel expansion, and two transitional
expansions for large v. If v is not too large, it
is reduced by recurrence to a region of best accuracy.
The transitional expansions give 12D accuracy for v > 500.
ACCURACY:
Results for integer v are indicated by *, where x and v
both vary from -125 to +125. Otherwise,
x ranges from 0 to 125, v ranges as indicated by "domain."
Error criterion is absolute, except relative when |jv()| > 1.
arithmetic v domain x domain # trials peak rms
IEEE 0,125 0,125 100000 4.6e-15 2.2e-16
IEEE -125,0 0,125 40000 5.4e-11 3.7e-13
IEEE 0,500 0,500 20000 4.4e-15 4.0e-16
Integer v:
IEEE -125,125 -125,125 50000 3.5e-15* 1.9e-16*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
*************************************************************************/
BESSEL_API double besseljv(double v, double x)
{
double result;
double k;
double q;
double t;
double y;
double an;
int i;
int sg;
int nint;
nint = 0;
sg = 1;
an = fabs(v);
y = ap::ifloor(an);
if( y==an )
{
nint = 1;
i = ap::trunc(an-16384*ap::ifloor(an/16384.0));
if( v<0.0 )
{
if( i%2!=0 )
{
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -