/***
*bessel.c - defines the bessel functions for C.
*
* Copyright (c) 1983-1989, Microsoft Corporation. All rights reserved.
*
*Purpose:
*
* This is a collection of routines for computing the bessel functions j0, j1,
* y0, y1, jn and yn. The approximations used for j0, j1, y0, and y1 are
* from the approximations listed in Hart, Computer Approximations, 1978.
* For these functions, a rational approximation with 18 places of accuracy
* after the decimal point has been selected. jn and yn are computed using
* the recursive formula that the bessel functions satisfy. Using these
* formulas their values can be computed from the values of the bessel
* functions of order 0 and 1. In the case of jn, the recursive formula
*
* jn(n-1,x) = (2.0*n/x)*jn(n,x) - jn(n+1,x)
*
* is used to stabily compute in the downward direction, normalizing in the
* the end by j0(x) in the usual manner. In the case of yn, the recursive
* formula
*
* yn(n+1,x) = (2.0*n/x)*yn(n,x) - yn(n-1,x)
*
* is used to stably compute the functions in the forward direction.
*
*
* Note: upon testing and experimentation the low range approximations were
* found to have an error on the order of 1.0e-14 in the neighborhood of
* 8.0. Moving the boundary point between the low range and high
* range approximations down to 7.5 reduced this error to less than
* 1.0e-14. This is not suprising. The high range asymptotoic is
* likely to have greater precision in the neighborhood of 8.0.
*
*Revision History:
*
* 06/05/89 WAJ Added this header. Made changes for C6 and -W3
* 06/06/89 WAJ Moved some of the routines into _RTEXT if MTHREAD.
* 08/17/90 WAJ Now uses _stdcall.
* 01/13/92 GDP changed domain_err. No full IEEE support yet
* 06-04-92 XY change to long double version
*
*******************************************************************************/
/*
* The functions sqrt, sin, cos, and log from the math library are used in
* the computations of the bessel functions.
*/
#include <math.h>
#include <transl.h>
#ifdef _X86SEG_
#include <os2supp.h>
#define _CALLTYPE1 _PASCAL
#else
#include <cruntime.h>
#endif
#define LD_VER
#ifdef LD_VER
#define D_TYPE long double
#else
#define D_TYPE double
#endif
static D_TYPE domain_err( int who, D_TYPE arg1, D_TYPE arg2 ); /* error routine for y0, y1, yn */
static D_TYPE evaluate( D_TYPE x, D_TYPE p[], int n1, D_TYPE q[], int n2 );
#ifdef FAR_CODE
#ifdef LD_VER
#pragma alloc_text( _RTEXT, _y0l, _y1l, _ynl, _j0l, _j1l, _jnl )
#else
#pragma alloc_text( _RTEXT, _y0, _y1, _yn, _j0, _j1, _jn )
#endif
#endif
/*
* Following are the constants needed for the computations of the bessel
* functions as in Hart.
*/
#define PI 3.14159265358979323846264338327950288
/* coefficients for Hart JZERO 5848, the low range approximation for _j0 */
static D_TYPE J0p[12] = {
0.1208181340866561224763662419e+12 ,
-0.2956513002312076810191727211e+11 ,
0.1729413174598080383355729444e+10 ,
-0.4281611621547871420502838045e+08 ,
0.5645169313685735094277826749e+06 ,
-0.4471963251278787165486324342e+04 ,
0.2281027164345610253338043760e+02 ,
-0.7777570245675629906097285039e-01 ,
0.1792464784997734953753734861e-03 ,
-0.2735011670747987792661294323e-06 ,
0.2553996162031530552738418047e-09 ,
-0.1135416951138795305302383379e-12
};
static D_TYPE J0q[5] = {
0.1208181340866561225104607422e+12 ,
0.6394034985432622416780183619e+09 ,
0.1480704129894421521840387092e+07 ,
0.1806405145147135549477896097e+04 ,
0.1e+01
};
/* coefficients for Hart 6548, P0 of the high range approximation for j0
and _y0 */
static D_TYPE P0p[6] = {
0.2277909019730468430227002627e+05 ,
0.4134538663958076579678016384e+05 ,
0.2117052338086494432193395727e+05 ,
0.3480648644324927034744531110e+04 ,
0.1537620190900835429577172500e+03 ,
0.8896154842421045523607480000e+00
};
static D_TYPE P0q[6] = {
0.2277909019730468431768423768e+05 ,
0.4137041249551041663989198384e+05 ,
0.2121535056188011573042256764e+05 ,
0.3502873513823560820735614230e+04 ,
0.1571115985808089364906848200e+03 ,
0.1e+01
};
/* coefficients for Hart 6948, Q0 of the high range approximation for _j0
and _y0 */
static D_TYPE Q0p[6] = {
-0.8922660020080009409846916000e+02 ,
-0.1859195364434299380025216900e+03 ,
-0.1118342992048273761126212300e+03 ,
-0.2230026166621419847169915000e+02 ,
-0.1244102674583563845913790000e+01 ,
-0.8803330304868075181663000000e-02
};
static D_TYPE Q0q[6] = {
0.5710502412851206190524764590e+04 ,
0.1195113154343461364695265329e+05 ,
0.7264278016921101883691345060e+04 ,
0.1488723123228375658161346980e+04 ,
0.9059376959499312585881878000e+02 ,
0.1e+01
};
/* coefficients for Hart JONE 6047, the low range approximation for _j1 */
static D_TYPE J1p[11] = {
0.4276440148317146125749678272e+11 ,
-0.5101551390663600782363700742e+10 ,
0.1928444249391651825203957853e+09 ,
-0.3445216851469225845312168656e+07 ,
0.3461845033978656620861683039e+05 ,
-0.2147334276854853222870548439e+03 ,
0.8645934990693258061130801001e+00 ,
-0.2302415336775925186376173217e-02 ,
0.3991878933072250766608485041e-05 ,
-0.4179409142757237977587032616e-08 ,
0.2060434024597835939153003596e-11
};
static D_TYPE J1q[5] = {
0.8552880296634292263013618479e+11 ,
0.4879975894656629161544052051e+09 ,
0.1226033111836540909388789681e+07 ,
0.1635396109098603257687643236e+04 ,
0.1e+01
};
/* coefficients for Hart PONE 6749, P1 of the high range approximation for
_j1 and y1 */
static D_TYPE P1p[6] = {
0.3522466491336797983417243730e+05 ,
0.6275884524716128126900567500e+05 ,
0.3135396311091595742386698880e+05 ,
0.4985483206059433843450045500e+04 ,
0.2111529182853962382105718000e+03 ,
0.1257171692914534155849500000e+01
};
static D_TYPE P1q[6] = {
0.3522466491336797980683904310e+05 ,
0.6269434695935605118888337310e+05 ,
0.3124040638190410399230157030e+05 ,
0.4930396490181088978386097000e+04 ,
0.2030775189134759322293574000e+03 ,
0.1e+01
};
/* coefficients for Hart QONE 7149, Q1 of the high range approximation for _j1
and y1 */
static D_TYPE Q1p[6] = {
0.3511751914303552822533318000e+03 ,
0.7210391804904475039280863000e+03 ,
0.4259873011654442389886993000e+03 ,
0.8318989576738508273252260000e+02 ,
0.4568171629551226706440500000e+01 ,
0.3532840052740123642735000000e-01
};
static D_TYPE Q1q[6] = {
0.7491737417180912771451950500e+04 ,
0.1541417733926509704998480510e+05 ,
0.9152231701516992270590472700e+04 ,
0.1811186700552351350672415800e+04 ,
0.1038187587462133728776636000e+03 ,
0.1e+01
};
/* coeffiecients for Hart YZERO 6245, the low range approximation for y0 */
static D_TYPE Y0p[9] = {
-0.2750286678629109583701933175e+20 ,
0.6587473275719554925999402049e+20 ,
-0.5247065581112764941297350814e+19 ,
0.1375624316399344078571335453e+18 ,
-0.1648605817185729473122082537e+16 ,
0.1025520859686394284509167421e+14 ,
-0.3436371222979040378171030138e+11 ,
0.5915213465686889654273830069e+08 ,
-0.4137035497933148554125235152e+05
};
static D_TYPE Y0q[9] = {
0.3726458838986165881989980739e+21 ,
0.4192417043410839973904769661e+19 ,
0.2392883043499781857439356652e+17 ,
0.9162038034075185262489147968e+14 ,
0.2613065755041081249568482092e+12 ,
0.5795122640700729537480087915e+09 ,
0.1001702641288906265666651753e+07 ,
0.1282452772478993804176329391e+04 ,
0.1e+01
};
/* coefficients for Hart YONE 6444, the low range approximation for y1 */
static D_TYPE Y1p[8] = {
-0.2923821961532962543101048748e+20 ,
0.7748520682186839645088094202e+19 ,
-0.3441048063084114446185461344e+18 ,
0.5915160760490070618496315281e+16 ,
-0.4863316942567175074828129117e+14 ,
0.2049696673745662182619800495e+12 ,
-0.4289471968855248801821819588e+09 ,
0.3556924009830526056691325215e+06
};
static D_TYPE Y1q[9] = {
0.1491311511302920350174081355e+21 ,
0.1818662841706134986885065935e+19 ,
0.1131639382698884526905082830e+17 ,
0.4755173588888137713092774006e+14 ,
0.1500221699156708987166369115e+12 ,
0.3716660798621930285596927703e+09 ,
0.7269147307198884569801913150e+06 ,
0.1072696143778925523322126700e+04 ,
0.1e+01
};
/*
* Function name: evaluate
*
* Arguments: x - long double
* p, q - long double arrays of coefficients
* n1, n2 - the order of the numerator and denominator
* polynomials
*
* Description: evaluate is meant strictly as a helper routine for the
* bessel function routines to evaluate the rational polynomial
* aproximations appearing in _j0, _j1, y0, and y1. Given the
* coefficient arrays in p and q, it evaluates the numerator
* and denominator polynomials through orders n1 and n2
* respectively, returning p(x)/q(x). This routine is not
* available to the user of the bessel function routines.
*
* Side Effects: evaluate uses the global data stored in the coefficients
* above. No other global data is used or affected.
*
* Copyright: written R.K. Wyss, Microsoft, Sept. 9, 1983
* copyright (c) Microsoft Corp. 1984
*
* History:
*/
static D_TYPE evaluate( D_TYPE x, D_TYPE p[], int n1, D_TYPE q[], int n2 )
{
D_TYPE numerator, denominator;
int i;
numerator = x*p[n1];
for ( i = n1-1 ; i > 0 ; i-- )
numerator = x*(p[i] + numerator);
numerator += p[0];
denominator = x*q[n2];
for ( i = n2-1 ; i > 0 ; i-- )
denominator = x*(q[i] + denominator);
denominator += q[0];
return( numerator/denominator );
}
/*
* Function name: _j0
*
* Arguments: x - long double
*
* Description: _j0 computes the bessel function of the first kind of zero
* order for real values of its argument x, where x can range
* from - infinity to + infinity. The algorithm is taken
* from Hart, Computer Approximations, 1978, and yields full
* long double precision accuracy.
*
* Side Effects: no global data other than the static coefficients above
* is used or affected.
*
* Copyright: written R.K. Wyss, Microsoft, Sept. 9, 1983
* copyright (c) Microsoft Corp. 1984
*
* History:
*/
#ifdef LD_VER
D_TYPE _j0l( D_TYPE x )
#else
D_TYPE _CALLTYPE1 _j0( D_TYPE x )
#endif
{
D_TYPE z, P0, Q0;
/* if the argument is negative, take the absolute value */
if ( x < 0.0 )
x = - x;
/* if x <= 7.5 use Hart JZERO 5847 */
if ( x <= 7.5 )
return( evaluate( x*x, J0p, 11, J0q, 4) );
/* else if x >= 7.5 use Hart PZERO 6548 and QZERO 6948, the high range
approximation */
else {
z = 8.0/x;
P0 = evaluate( z*z, P0p, 5, P0q, 5);
Q0 = z*evaluate( z*z, Q0p, 5, Q0q, 5);
return( sqrtl(2.0/(PI*x))*(P0*cosl(x-PI/4) - Q0*sinl(x-PI/4)) );
}
}
/*
* Function name: _j1
*
* Arguments: x - long double
*
* Description: _j1 computes the bessel function of the first kind of the
* first order for real values of its argument x, where x can
* range from - infinity to + infinity. The algorithm is taken
* from Hart, Computer Approximations, 1978, and yields full
* D_TYPE precision accuracy.
*
* Side Effects: no global data other than the static coefficients above
* is used or affected.
*
* Copyright: written R.K. Wyss, Microsoft, Sept. 9, 1983
* copyright (c) Microsoft Corp. 1984
*
* History:
*/
#ifdef LD_VER
D_TYPE _j1l( D_TYPE x )
#else
D_TYPE _CALLTYPE1 _j1( D_TYPE x )
#endif
{
D_TYPE z, P1, Q1;
int sign;
/* if the argument is negative, take the absolute value and set sign */
sign = 1;
if( x < 0.0 ){
x = -x;
sign = -1;
}
/* if x <= 7.5 use Hart JONE 6047 */
if ( x <= 7.5 )
return( sign*x*evaluate( x*x, J1p, 10, J1q, 4) );
/* else if x > 7.5 use Hart PONE 6749 and QONE 7149, the high range
approximation */
else {
z = 8.0/x;
P1 = evaluate( z*z, P1p, 5, P1q, 5);
Q1 = z*evaluate( z*z, Q1p, 5, Q1q, 5);
return( sign*sqrtl(2.0/(PI*x))*
( P1*cosl(x-3.0*PI/4.0) - Q1*sinl(x-3.0*PI/4.0) ) );
}
}
/*
* Function name: _y0
*
* Arguments: x - long double
*
* Description: y0 computes the bessel function of the second kind of zero
* order for real values of its argument x, where x can range
* from 0 to + infinity. The algorithm is taken from Hart,
* Computer Approximations, 1978, and yields full long double
* precision accuracy.
*
* Side Effects: no global data other than the static coefficients above
* is used or affected.
*
* Copyright: written R.K. Wyss, Microsoft, Sept. 9, 1983
* copyright (c) Microsoft Corp. 1984
*
* History:
*/
#ifdef LD_VER
D_TYPE _y0l( D_TYPE x )
#else
D_TYPE _CALLTYPE1 _y0( D_TYPE x )
#endif
{
D_TYPE z, P0, Q0;
/* if the argument is negative, set EDOM error, print an error message,
* and return -HUGE
*/
if (x < 0.0)
return( domain_err(OP_Y0 , x, LD_IND) );
/* if x <= 7.5 use Hart YZERO 6245, the low range approximation */
if ( x <= 7.5 )
return( evaluate( x*x, Y0p, 8, Y0q, 8) + (2.0/PI)*_j0l(x)*logl(x) );
/* else if x > 7.5 use Hart PZERO 6548 and QZERO 6948, the high range
approximation */
else {
z = 8.0/x;
P0 = evaluate( z*z, P0p, 5, P0q, 5);
Q0 = z*evaluate( z*z, Q0p, 5, Q0q, 5);
return( sqrtl(2.0/(PI*x))*(P0*sinl(x-PI/4) + Q0*cosl(x-PI/4)) );
}
}
/*
* Function name: _y1
*
* Arguments: x - long double
*
* Description: y1 computes the bessel function of the second kind of first
* order for real values of its argument x, where x can range
* from 0 to + infinity. The algorithm is taken from Hart,
* Computer Approximations, 1978, and yields full long double
* precision accuracy.
*
* Side Effects: no global data other than the static coefficients above
* is used or affected.
*
* Copyright: written R.K. Wyss, Microsoft, Sept. 9, 1983
* copyright (c) Microsoft Corp. 1984
*
* History:
*/
#ifdef LD_VER
D_TYPE _y1l( D_TYPE x )
#else
D_TYPE _CALLTYPE1 _y1( D_TYPE x )
#endif
{
D_TYPE z, P1, Q1;
/* if the argument is negative, set EDOM error, print an error message,
* and return -HUGE
*/
if (x < 0.0)
return( domain_err(OP_Y1, x, LD_IND) );
/* if x <= 7.5 use Hart YONE 6444, the low range approximation */
if ( x <= 7.5 )
return( x*evaluate( x*x, Y1p, 7, Y1q, 8)
+ (2.0/PI)*(_j1l(x)*logl(x) - 1.0/x) );
/* else if x > 7.5 use Hart PONE 6749 and QONE 7149, the high range
approximation */
else {
z = 8.0/x;
P1 = evaluate( z*z, P1p, 5, P1q, 5);
Q1 = z*evaluate( z*z, Q1p, 5, Q1q, 5);
return( sqrtl(2.0/(PI*x))*
( P1*sinl(x-3.0*PI/4.0) + Q1*cosl(x-3.0*PI/4.0) ) );
}
}
/*
* Function name: _jn
*
* Arguments: n - integer
* x - long double
*
* Description: _jn computes the bessel function of the first kind of order
* n for real values of its argument, where x can range from
* - infinity to + infinity, and n can range over the integers
* from - infinity to + infinity. The function is computed
* by recursion, using the formula
*
* _jn(n-1,x) = (2.0*n/x)*_jn(n,x) - _jn(n+1,x)
*
* stabilly in the downward direction, normalizing by _j0(x)
* in the end in the usual manner.
*
* Side Effects: the routines _j0, y0, and yn are called during the
* execution of this routine.
*
* Copyright: written R.K. Wyss, Microsoft, Sept. 9, 1983
* copyright (c) Microsoft Corp. 1984
*
* History:
* 07/29/85 Greg Whitten
* rewrote _jn to use Hart suggested algorithm
*/
#ifdef LD_VER
D_TYPE _jnl( int n, D_TYPE x )
#else
D_TYPE _CALLTYPE1 _jn( int n, D_TYPE x )
#endif
{
int i;
D_TYPE x2, jm1, j, jnratio, hold;
/* use symmetry relationships: _j(-n,x) = _j(n,-x) */
if( n < 0 ){
n = -n;
x = -x;
}
/* if n = 0 use _j0(x) and if n = 1 use _j1(x) functions */
if (n == 0)
return (_j0l(x));
if (n == 1)
return (_j1l(x));
/* if x = 0.0 then _j(n,0.0) = 0.0 for n > 0 (_j(0,x) = 1.0) */
if (x == 0.0)
return (0.0);
/* otherwise - must use the recurrence relation
*
* _jn(n+1,x) = (2.0*n/x)*_jn(n,x) - _jn(n-1,x) forward
* _jn(n-1,x) = (2.0*n/x)*_jn(n,x) - _jn(n+1,x) backward
*/
if( (long double)n < fabsl(x) ) {
/* stably compute _jn using forward recurrence above */
n <<= 1; /* n *= 2 (n is positive) */
jm1 = _j0l(x);
j = _j1l(x);
i = 2;
for(;;) {
hold = j;
j = ((long double)(i))*j/x - jm1;
i += 2;
if (i == n)
return (j);
jm1 = hold;
}
}
else {
/* stably compute _jn using backward recurrence above */
/* use Hart continued fraction formula for j(n,x)/j(n-1,x)
* so that we can compute a normalization factor
*/
n <<= 1; /* n *= 2 (n is positive) */
x2 = x*x;
hold = 0.0; /* initial continued fraction tail value */
for (i=n+36; i>n; i-=2)
hold = x2/((long double)(i) - hold);
jnratio = j = x/((long double)(n) - hold);
jm1 = 1.0;
/* have jn/jn-1 ratio - now use backward recurrence */
i = n-2;
for (;;) {
hold = jm1;
jm1 = ((long double)(i))*jm1/x - j;
i -= 2;
if (i == 0)
break;
j = hold;
}
/* jm1 is relative j0(x) so normalize it for final result
*
* jnratio = K*j(n,x) and jm1 = K*_j0(x)
*/
return(_j0l(x)*jnratio/jm1);
}
}
/*
* Function name: _yn
*
* Arguments: n - integer
* x - long double
*
* Description: yn computes the bessel function of the second kind of order
* n for real values of its argument x, where x can range from
* 0 to + infinity, and n can range over the integers from
* - infinity to + infinity. The function is computed by
* recursion from y0 and y1, using the recursive formula
*
* yn(n+1,x) = (2.0*n/x)*yn(n,x) - yn(n-1,x)
*
* in the forward direction.
*
* Side Effects: the routines y0 and y1 are called during the execution
* of this routine.
*
* Copyright: written R.K. Wyss, Microsoft, Sept. 9, 1983
* copyright (c) Microsoft Corp. 1984
*
* History:
* 08/09/85 Greg Whitten
* added check for n==0 and n==1
* 04/20/87 Barry McCord
* eliminated use of "const" as an identifier for ANSI conformance
*/
#ifdef LD_VER
D_TYPE _ynl( int n, D_TYPE x )
#else
D_TYPE _CALLTYPE1 _yn( int n, D_TYPE x )
#endif
{
int i;
int sign;
D_TYPE constant, yn2, yn1, yn0;
/* if the argument is negative, set EDOM error, print an error message,
* and return -HUGE
*/
if (x < 0.0)
return(domain_err(OP_YN, x, LD_IND));
/* take the absolute value of n, and set sign accordingly */
sign = 1;
if( n < 0 ){
n = -n;
if( n&1 )
sign = -1;
}
if( n == 0 )
return( sign*_y0l(x) );
if (n == 1)
return( sign*_y1l(x) );
/* otherwise go ahead and compute the function by iteration */
yn0 = _y0l(x);
yn1 = _y1l(x);
constant = 2.0/x;
for( i = 1 ; i < n ; i++ ){
yn2 = constant*i*yn1 - yn0;
yn0 = yn1;
yn1 = yn2;
}
return( sign*yn2 );
}
static D_TYPE domain_err( int who, D_TYPE arg1, D_TYPE arg2 )
{
/*
#ifdef LD_VER
#error long long double version not supported
#endif
*/
unsigned int savedcw;
savedcw = _maskfp();
return _except1(FP_I, who, arg1, arg2, savedcw);
}