ReactOS  0.4.14-dev-593-g1793dcc
jn_yn.c
Go to the documentation of this file.
1 /* @(#)e_jn.c 5.1 93/09/24 */
2 /*
3  * ====================================================
4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
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)
14 static 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__
47 static const double
48 #else
49 static double
50 #endif
51 invsqrtpi= 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
52 two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
53 one = 1.00000000000000000000e+00; /* 0x3FF00000, 0x00000000 */
54 
55 #ifdef __STDC__
56 static const double zero = 0.00000000000000000000e+00;
57 #else
58 static double zero = 0.00000000000000000000e+00;
59 #endif
60 
61 #ifdef __STDC__
62  double __ieee754_jn(int n, double x)
63 #else
64  double __ieee754_jn(n,x)
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
227  double __ieee754_yn(n,x)
228  int n; double x;
229 #endif
230 {
231  int32_t i,hx,ix,lx;
232  int32_t sign;
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 }
static size_t double int int int * sign
Definition: printf.c:69
double __ieee754_y0(double)
Definition: j0_y0.c:179
GLubyte GLubyte GLubyte GLubyte w
Definition: glext.h:6102
long sgn(REAL x)
Definition: varray.cc:48
static __inline double __ieee754_sqrt(double x)
Definition: ieee754.h:40
GLdouble n
Definition: glext.h:7729
GLdouble GLdouble t
Definition: gl.h:2047
GLint GLint GLint GLint GLint x
Definition: gl.h:1548
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
#define GET_HIGH_WORD(i, d)
Definition: ieee754.h:26
double __ieee754_yn(int n, double x)
Definition: jn_yn.c:227
#define a
Definition: ke_i.h:78
#define e
Definition: ke_i.h:82
static double two
Definition: jn_yn.c:52
double __ieee754_j0(double)
Definition: j0_y0.c:102
static const char mbstate_t *static wchar_t const char mbstate_t *static const wchar_t int *static double
Definition: string.c:80
GLdouble GLdouble z
Definition: glext.h:5874
#define b
Definition: ke_i.h:79
GLboolean GLboolean GLboolean b
Definition: glext.h:6204
#define EXTRACT_WORDS(ix0, ix1, d)
Definition: ieee754.h:16
const GLubyte * c
Definition: glext.h:8905
static double zero
Definition: jn_yn.c:58
double __ieee754_j1(double)
Definition: j1_y1.c:103
static double invsqrtpi
Definition: jn_yn.c:51
unsigned int u_int32_t
Definition: rosdhcp.h:35
GLdouble s
Definition: gl.h:2039
INT32 int32_t
Definition: types.h:71
_Check_return_ _CRT_JIT_INTRINSIC double __cdecl fabs(_In_ double x)
double __ieee754_y1(double)
Definition: j1_y1.c:182
static __inline void __sincos(double x, double *s, double *c)
Definition: ieee754.h:43
const GLdouble * v
Definition: gl.h:2040
static double one
Definition: jn_yn.c:53
static calc_node_t temp
Definition: rpn_ieee.c:38
#define HUGE_VAL
Definition: math.h:51
double __ieee754_jn(int n, double x)
Definition: jn_yn.c:64
#define c
Definition: ke_i.h:80
static __inline double __ieee754_log(double x)
Definition: ieee754.h:41
GLboolean GLboolean GLboolean GLboolean a
Definition: glext.h:6204
static TRIO_CONST char rcsid[]
Definition: trio.c:756
int k
Definition: mpi.c:3369