ReactOS 0.4.15-dev-8614-gbc76250
jn_yn.c
Go to the documentation of this file.
1/* @(#)e_jn.c 5.1 93/09/24 */
2/*
3 * ====================================================
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13#if defined(LIBM_SCCS) && !defined(lint)
14static char rcsid[] = "\$NetBSD: e_jn.c,v 1.9 1995/05/10 20:45:34 jtc Exp \$";
15#endif
16
17/*
18 * __ieee754_jn(n, x), __ieee754_yn(n, x)
19 * floating point Bessel's function of the 1st and 2nd kind
20 * of order n
21 *
22 * Special cases:
23 * y0(0)=y1(0)=yn(n,0) = -inf with overflow signal;
24 * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
25 * Note 2. About jn(n,x), yn(n,x)
26 * For n=0, j0(x) is called,
27 * for n=1, j1(x) is called,
28 * for n<x, forward recursion us used starting
29 * from values of j0(x) and j1(x).
30 * for n>x, a continued fraction approximation to
31 * j(n,x)/j(n-1,x) is evaluated and then backward
32 * recursion is used starting from a supposed value
33 * for j(n,x). The resulting value of j(0,x) is
34 * compared with the actual value to correct the
35 * supposed value of j(n,x).
36 *
37 * yn(n,x) is similar in all respects, except
38 * that forward recursion is used for all
39 * values of n>1.
40 *
41 */
42
43#include "math.h"
44#include "ieee754.h"
45
46#ifdef __STDC__
47static const double
48#else
49static double
50#endif
51invsqrtpi= 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
52two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
53one = 1.00000000000000000000e+00; /* 0x3FF00000, 0x00000000 */
54
55#ifdef __STDC__
56static const double zero = 0.00000000000000000000e+00;
57#else
58static double zero = 0.00000000000000000000e+00;
59#endif
60
61#ifdef __STDC__
62 double __ieee754_jn(int n, double x)
63#else
65 int n; double x;
66#endif
67{
68 int32_t i,hx,ix,lx, sgn;
69 double a, b, temp, di;
70 double z, w;
71
72 /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
73 * Thus, J(-n,x) = J(n,-x)
74 */
75 EXTRACT_WORDS(hx,lx,x);
76 ix = 0x7fffffff&hx;
77 /* if J(n,NaN) is NaN */
78 if((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000) return x+x;
79 if(n<0){
80 n = -n;
81 x = -x;
82 hx ^= 0x80000000;
83 }
84 if(n==0) return(__ieee754_j0(x));
85 if(n==1) return(__ieee754_j1(x));
86 sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */
87 x = fabs(x);
88 if((ix|lx)==0||ix>=0x7ff00000) /* if x is 0 or inf */
89 b = zero;
90 else if((double)n<=x) {
91 /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
92 if(ix>=0x52D00000) { /* x > 2**302 */
93 /* (x >> n**2)
94 * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
95 * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
96 * Let s=sin(x), c=cos(x),
97 * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
98 *
99 * n sin(xn)*sqt2 cos(xn)*sqt2
100 * ----------------------------------
101 * 0 s-c c+s
102 * 1 -s-c -c+s
103 * 2 -s+c -c-s
104 * 3 s+c c-s
105 */
106 double s;
107 double c;
108 __sincos (x, &s, &c);
109 switch(n&3) {
110 case 0: temp = c + s; break;
111 case 1: temp = -c + s; break;
112 case 2: temp = -c - s; break;
113 case 3: temp = c - s; break;
114 }
116 } else {
117 a = __ieee754_j0(x);
118 b = __ieee754_j1(x);
119 for(i=1;i<n;i++){
120 temp = b;
121 b = b*((double)(i+i)/x) - a; /* avoid underflow */
122 a = temp;
123 }
124 }
125 } else {
126 if(ix<0x3e100000) { /* x < 2**-29 */
127 /* x is tiny, return the first Taylor expansion of J(n,x)
128 * J(n,x) = 1/n!*(x/2)^n - ...
129 */
130 if(n>33) /* underflow */
131 b = zero;
132 else {
133 temp = x*0.5; b = temp;
134 for (a=one,i=2;i<=n;i++) {
135 a *= (double)i; /* a = n! */
136 b *= temp; /* b = (x/2)^n */
137 }
138 b = b/a;
139 }
140 } else {
141 /* use backward recurrence */
142 /* x x^2 x^2
143 * J(n,x)/J(n-1,x) = ---- ------ ------ .....
144 * 2n - 2(n+1) - 2(n+2)
145 *
146 * 1 1 1
147 * (for large x) = ---- ------ ------ .....
148 * 2n 2(n+1) 2(n+2)
149 * -- - ------ - ------ -
150 * x x x
151 *
152 * Let w = 2n/x and h=2/x, then the above quotient
153 * is equal to the continued fraction:
154 * 1
155 * = -----------------------
156 * 1
157 * w - -----------------
158 * 1
159 * w+h - ---------
160 * w+2h - ...
161 *
162 * To determine how many terms needed, let
163 * Q(0) = w, Q(1) = w(w+h) - 1,
164 * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
165 * When Q(k) > 1e4 good for single
166 * When Q(k) > 1e9 good for double
167 * When Q(k) > 1e17 good for quadruple
168 */
169 /* determine k */
170 double t,v;
171 double q0,q1,h,tmp; int32_t k,m;
172 w = (n+n)/(double)x; h = 2.0/(double)x;
173 q0 = w; z = w+h; q1 = w*z - 1.0; k=1;
174 while(q1<1.0e9) {
175 k += 1; z += h;
176 tmp = z*q1 - q0;
177 q0 = q1;
178 q1 = tmp;
179 }
180 m = n+n;
181 for(t=zero, i = 2*(n+k); i>=m; i -= 2) t = one/(i/x-t);
182 a = t;
183 b = one;
184 /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
185 * Hence, if n*(log(2n/x)) > ...
186 * single 8.8722839355e+01
187 * double 7.09782712893383973096e+02
188 * long double 1.1356523406294143949491931077970765006170e+04
189 * then recurrent value may overflow and the result is
190 * likely underflow to zero
191 */
192 tmp = n;
193 v = two/x;
194 tmp = tmp*__ieee754_log(fabs(v*tmp));
195 if(tmp<7.09782712893383973096e+02) {
196 for(i=n-1,di=(double)(i+i);i>0;i--){
197 temp = b;
198 b *= di;
199 b = b/x - a;
200 a = temp;
201 di -= two;
202 }
203 } else {
204 for(i=n-1,di=(double)(i+i);i>0;i--){
205 temp = b;
206 b *= di;
207 b = b/x - a;
208 a = temp;
209 di -= two;
210 /* scale b to avoid spurious overflow */
211 if(b>1e100) {
212 a /= b;
213 t /= b;
214 b = one;
215 }
216 }
217 }
218 b = (t*__ieee754_j0(x)/b);
219 }
220 }
221 if(sgn==1) return -b; else return b;
222}
223
224#ifdef __STDC__
225 double __ieee754_yn(int n, double x)
226#else
228 int n; double x;
229#endif
230{
231 int32_t i,hx,ix,lx;
233 double a, b, temp;
234
235 EXTRACT_WORDS(hx,lx,x);
236 ix = 0x7fffffff&hx;
237 /* if Y(n,NaN) is NaN */
238 if((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000) return x+x;
239 if((ix|lx)==0) return -HUGE_VAL+x; /* -inf and overflow exception. */;
240 if(hx<0) return zero/(zero*x);
241 sign = 1;
242 if(n<0){
243 n = -n;
244 sign = 1 - ((n&1)<<1);
245 }
246 if(n==0) return(__ieee754_y0(x));
247 if(n==1) return(sign*__ieee754_y1(x));
248 if(ix==0x7ff00000) return zero;
249 if(ix>=0x52D00000) { /* x > 2**302 */
250 /* (x >> n**2)
251 * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
252 * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
253 * Let s=sin(x), c=cos(x),
254 * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
255 *
256 * n sin(xn)*sqt2 cos(xn)*sqt2
257 * ----------------------------------
258 * 0 s-c c+s
259 * 1 -s-c -c+s
260 * 2 -s+c -c-s
261 * 3 s+c c-s
262 */
263 double c;
264 double s;
265 __sincos (x, &s, &c);
266 switch(n&3) {
267 case 0: temp = s - c; break;
268 case 1: temp = -s - c; break;
269 case 2: temp = -s + c; break;
270 case 3: temp = s + c; break;
271 }
273 } else {
274 u_int32_t high;
275 a = __ieee754_y0(x);
276 b = __ieee754_y1(x);
277 /* quit if b is -inf */
278 GET_HIGH_WORD(high,b);
279 for(i=1;i<n&&high!=0xfff00000;i++){
280 temp = b;
281 b = ((double)(i+i)/x)*b - a;
282 GET_HIGH_WORD(high,b);
283 a = temp;
284 }
285 }
286 if(sign>0) return b; else return -b;
287}
INT32 int32_t
Definition: types.h:71
GLint GLint GLint GLint GLint x
Definition: gl.h:1548
const GLdouble * v
Definition: gl.h:2040
GLdouble s
Definition: gl.h:2039
GLdouble GLdouble t
Definition: gl.h:2047
GLdouble n
Definition: glext.h:7729
const GLubyte * c
Definition: glext.h:8905
GLboolean GLboolean GLboolean b
Definition: glext.h:6204
GLboolean GLboolean GLboolean GLboolean a
Definition: glext.h:6204
GLubyte GLubyte GLubyte GLubyte w
Definition: glext.h:6102
GLdouble GLdouble z
Definition: glext.h:5874
const GLfloat * m
Definition: glext.h:10848
GLfloat GLfloat GLfloat GLfloat h
Definition: glext.h:7723
GLsizei GLenum const GLvoid GLsizei GLenum GLbyte GLbyte GLbyte GLdouble GLdouble GLdouble GLfloat GLfloat GLfloat GLint GLint GLint GLshort GLshort GLshort GLubyte GLubyte GLubyte GLuint GLuint GLuint GLushort GLushort GLushort GLbyte GLbyte GLbyte GLbyte GLdouble GLdouble GLdouble GLdouble GLfloat GLfloat GLfloat GLfloat GLint GLint GLint GLint GLshort GLshort GLshort GLshort GLubyte GLubyte GLubyte GLubyte GLuint GLuint GLuint GLuint GLushort GLushort GLushort GLushort GLboolean const GLdouble const GLfloat const GLint const GLshort const GLbyte const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLdouble const GLfloat const GLfloat const GLint const GLint const GLshort const GLshort const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLfloat const GLint const GLshort GLenum GLenum GLenum GLfloat GLenum GLint GLenum GLenum GLenum GLfloat GLenum GLenum GLint GLenum GLfloat GLenum GLint GLint GLushort GLenum GLenum GLfloat GLenum GLenum GLint GLfloat const GLubyte GLenum GLenum GLenum const GLfloat GLenum GLenum const GLint GLenum GLint GLint GLsizei GLsizei GLint GLenum GLenum const GLvoid GLenum GLenum const GLfloat GLenum GLenum const GLint GLenum GLenum const GLdouble GLenum GLenum const GLfloat GLenum GLenum const GLint GLsizei GLuint GLfloat GLuint GLbitfield GLfloat GLint GLuint GLboolean GLenum GLfloat GLenum GLbitfield GLenum GLfloat GLfloat GLint GLint const GLfloat GLenum GLfloat GLfloat GLint GLint GLfloat GLfloat GLint GLint const GLfloat GLint GLfloat GLfloat GLint GLfloat GLfloat GLint GLfloat GLfloat const GLdouble const GLfloat const GLdouble const GLfloat GLint i
Definition: glfuncs.h:248
static double two
Definition: jn_yn.c:52
static double zero
Definition: jn_yn.c:58
static double invsqrtpi
Definition: jn_yn.c:51
double __ieee754_yn(int n, double x)
Definition: jn_yn.c:227
static double one
Definition: jn_yn.c:53
double __ieee754_jn(int n, double x)
Definition: jn_yn.c:64
double __ieee754_j0(double)
Definition: j0_y0.c:102
double __ieee754_y0(double)
Definition: j0_y0.c:179
double __ieee754_y1(double)
Definition: j1_y1.c:182
static __inline void __sincos(double x, double *s, double *c)
Definition: ieee754.h:43
double __ieee754_j1(double)
Definition: j1_y1.c:103
#define EXTRACT_WORDS(ix0, ix1, d)
Definition: ieee754.h:16
static __inline double __ieee754_log(double x)
Definition: ieee754.h:41
#define GET_HIGH_WORD(i, d)
Definition: ieee754.h:26
static __inline double __ieee754_sqrt(double x)
Definition: ieee754.h:40
_Check_return_ _CRT_JIT_INTRINSIC double __cdecl fabs(_In_ double x)
Definition: fabs.c:17
#define HUGE_VAL
Definition: math.h:51
#define a
Definition: ke_i.h:78
#define c
Definition: ke_i.h:80
#define b
Definition: ke_i.h:79
#define sign(x)
Definition: mapdesc.cc:613
static const char mbstate_t *static wchar_t const char mbstate_t *static const wchar_t int *static double
Definition: string.c:80
int k
Definition: mpi.c:3369
unsigned int u_int32_t
Definition: rosdhcp.h:35
static calc_node_t temp
Definition: rpn_ieee.c:38
static TRIO_CONST char rcsid[]
Definition: trio.c:753
long sgn(REAL x)
Definition: varray.cc:48