/*** *bessel.c - defines the bessel functions for C. * * Copyright (c) 1983-2001, 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 * 04-06-93 SKS Replace _CALLTYPE* with __cdecl * 08-28-96 JWM Disabled warning 4056. * *******************************************************************************/ /* * The functions sqrt, sin, cos, and log from the math library are used in * the computations of the bessel functions. */ #include #include #ifdef _M_IX86 #pragma warning(disable:4056) #endif #ifdef _X86SEG_ #include #define _CALLTYPE1 _PASCAL #else #include #endif #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 - double * p, q - 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. * * Author: written R.K. Wyss, Microsoft, Sept. 9, 1983 * * 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 - 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 * double precision accuracy. * * Side Effects: no global data other than the static coefficients above * is used or affected. * * Author: written R.K. Wyss, Microsoft, Sept. 9, 1983 * * History: */ #ifdef LD_VER D_TYPE _ cdecl _j0l( D_TYPE x ) #else D_TYPE __cdecl _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( sqrt(2.0/(PI*x))*(P0*cos(x-PI/4) - Q0*sin(x-PI/4)) ); } } /* * Function name: _j1 * * Arguments: x - 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. * * Author: written R.K. Wyss, Microsoft, Sept. 9, 1983 * * History: */ #ifdef LD_VER D_TYPE _cdecl _j1l( D_TYPE x ) #else D_TYPE __cdecl _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*sqrt(2.0/(PI*x))* ( P1*cos(x-3.0*PI/4.0) - Q1*sin(x-3.0*PI/4.0) ) ); } } /* * Function name: _y0 * * Arguments: x - 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 double * precision accuracy. * * Side Effects: no global data other than the static coefficients above * is used or affected. * * Author: written R.K. Wyss, Microsoft, Sept. 9, 1983 * * History: */ #ifdef LD_VER D_TYPE _cdecl _y0l( D_TYPE x ) #else D_TYPE __cdecl _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, D_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)*_j0(x)*log(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( sqrt(2.0/(PI*x))*(P0*sin(x-PI/4) + Q0*cos(x-PI/4)) ); } } /* * Function name: _y1 * * Arguments: x - 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 double * precision accuracy. * * Side Effects: no global data other than the static coefficients above * is used or affected. * * Author: written R.K. Wyss, Microsoft, Sept. 9, 1983 * * History: */ #ifdef LD_VER D_TYPE _cdecl _y1l( D_TYPE x ) #else D_TYPE __cdecl _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, D_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)*(_j1(x)*log(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( sqrt(2.0/(PI*x))* ( P1*sin(x-3.0*PI/4.0) + Q1*cos(x-3.0*PI/4.0) ) ); } } /* * Function name: _jn * * Arguments: n - integer * x - 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. * * Author: written R.K. Wyss, Microsoft, Sept. 9, 1983 * * History: * 07/29/85 Greg Whitten * rewrote _jn to use Hart suggested algorithm */ #ifdef LD_VER D_TYPE _cdecl _jnl( int n, D_TYPE x ) #else D_TYPE __cdecl _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 (_j0(x)); if (n == 1) return (_j1(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( (double)n < fabs(x) ) { /* stably compute _jn using forward recurrence above */ n <<= 1; /* n *= 2 (n is positive) */ jm1 = _j0(x); j = _j1(x); i = 2; for(;;) { hold = j; j = ((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/((double)(i) - hold); jnratio = j = x/((double)(n) - hold); jm1 = 1.0; /* have jn/jn-1 ratio - now use backward recurrence */ i = n-2; for (;;) { hold = jm1; jm1 = ((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(_j0(x)*jnratio/jm1); } } /* * Function name: _yn * * Arguments: n - integer * x - 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. * * Author: written R.K. Wyss, Microsoft, Sept. 9, 1983 * * 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 _cdecl _ynl( int n, D_TYPE x ) #else D_TYPE __cdecl _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, D_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*_y0(x) ); if (n == 1) return( sign*_y1(x) ); /* otherwise go ahead and compute the function by iteration */ yn0 = _y0(x); yn1 = _y1(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 double version not supported #endif uintptr_t savedcw; savedcw = _maskfp(); return _except1(FP_I, who, arg1, arg2, savedcw); }