POK
e_jn.c
1 /*
2  * POK header
3  *
4  * The following file is a part of the POK project. Any modification should
5  * made according to the POK licence. You CANNOT use this file or a part of
6  * this file is this part of a file for your own project
7  *
8  * For more information on the POK licence, please see our LICENCE FILE
9  *
10  * Please follow the coding guidelines described in doc/CODING_GUIDELINES
11  *
12  * Copyright (c) 2007-2009 POK team
13  *
14  * Created by julien on Fri Jan 30 14:41:34 2009
15  */
16 
17 /* @(#)e_jn.c 5.1 93/09/24 */
18 /*
19  * ====================================================
20  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
21  *
22  * Developed at SunPro, a Sun Microsystems, Inc. business.
23  * Permission to use, copy, modify, and distribute this
24  * software is freely granted, provided that this notice
25  * is preserved.
26  * ====================================================
27  */
28 
29 /*
30  * __ieee754_jn(n, x), __ieee754_yn(n, x)
31  * floating point Bessel's function of the 1st and 2nd kind
32  * of order n
33  *
34  * Special cases:
35  * y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
36  * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
37  * Note 2. About jn(n,x), yn(n,x)
38  * For n=0, j0(x) is called,
39  * for n=1, j1(x) is called,
40  * for n<x, forward recursion us used starting
41  * from values of j0(x) and j1(x).
42  * for n>x, a continued fraction approximation to
43  * j(n,x)/j(n-1,x) is evaluated and then backward
44  * recursion is used starting from a supposed value
45  * for j(n,x). The resulting value of j(0,x) is
46  * compared with the actual value to correct the
47  * supposed value of j(n,x).
48  *
49  * yn(n,x) is similar in all respects, except
50  * that forward recursion is used for all
51  * values of n>1.
52  *
53  */
54 
55 #ifdef POK_NEEDS_LIBMATH
56 
57 #include "namespace.h"
58 #include "math_private.h"
59 #include <libm.h>
60 
61 static const double
62 invsqrtpi= 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
63 two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
64 one = 1.00000000000000000000e+00; /* 0x3FF00000, 0x00000000 */
65 
66 static const double zero = 0.00000000000000000000e+00;
67 
68 double
69 __ieee754_jn(int n, double x)
70 {
71  int32_t i,hx,ix,lx, sgn;
72  double a, b, temp, di;
73  double z, w;
74 
75  temp = 0;
76  /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
77  * Thus, J(-n,x) = J(n,-x)
78  */
79  EXTRACT_WORDS(hx,lx,x);
80  ix = 0x7fffffff&hx;
81  /* if J(n,NaN) is NaN */
82  if((ix|((uint32_t)(lx|-lx))>>31)>0x7ff00000) return x+x;
83  if(n<0){
84  n = -n;
85  x = -x;
86  hx ^= 0x80000000;
87  }
88  if(n==0) return(__ieee754_j0(x));
89  if(n==1) return(__ieee754_j1(x));
90  sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */
91  x = fabs(x);
92  if((ix|lx)==0||ix>=0x7ff00000) /* if x is 0 or inf */
93  b = zero;
94  else if((double)n<=x) {
95  /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
96  if(ix>=0x52D00000) { /* x > 2**302 */
97  /* (x >> n**2)
98  * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
99  * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
100  * Let s=sin(x), c=cos(x),
101  * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
102  *
103  * n sin(xn)*sqt2 cos(xn)*sqt2
104  * ----------------------------------
105  * 0 s-c c+s
106  * 1 -s-c -c+s
107  * 2 -s+c -c-s
108  * 3 s+c c-s
109  */
110  switch(n&3) {
111  case 0: temp = cos(x)+sin(x); break;
112  case 1: temp = -cos(x)+sin(x); break;
113  case 2: temp = -cos(x)-sin(x); break;
114  case 3: temp = cos(x)-sin(x); break;
115  }
116  b = invsqrtpi*temp/sqrt(x);
117  } else {
118  a = __ieee754_j0(x);
119  b = __ieee754_j1(x);
120  for(i=1;i<n;i++){
121  temp = b;
122  b = b*((double)(i+i)/x) - a; /* avoid underflow */
123  a = temp;
124  }
125  }
126  } else {
127  if(ix<0x3e100000) { /* x < 2**-29 */
128  /* x is tiny, return the first Taylor expansion of J(n,x)
129  * J(n,x) = 1/n!*(x/2)^n - ...
130  */
131  if(n>33) /* underflow */
132  b = zero;
133  else {
134  temp = x*0.5; b = temp;
135  for (a=one,i=2;i<=n;i++) {
136  a *= (double)i; /* a = n! */
137  b *= temp; /* b = (x/2)^n */
138  }
139  b = b/a;
140  }
141  } else {
142  /* use backward recurrence */
143  /* x x^2 x^2
144  * J(n,x)/J(n-1,x) = ---- ------ ------ .....
145  * 2n - 2(n+1) - 2(n+2)
146  *
147  * 1 1 1
148  * (for large x) = ---- ------ ------ .....
149  * 2n 2(n+1) 2(n+2)
150  * -- - ------ - ------ -
151  * x x x
152  *
153  * Let w = 2n/x and h=2/x, then the above quotient
154  * is equal to the continued fraction:
155  * 1
156  * = -----------------------
157  * 1
158  * w - -----------------
159  * 1
160  * w+h - ---------
161  * w+2h - ...
162  *
163  * To determine how many terms needed, let
164  * Q(0) = w, Q(1) = w(w+h) - 1,
165  * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
166  * When Q(k) > 1e4 good for single
167  * When Q(k) > 1e9 good for double
168  * When Q(k) > 1e17 good for quadruple
169  */
170  /* determine k */
171  double t,v;
172  double q0,q1,h,tmp; int32_t k,m;
173  w = (n+n)/(double)x; h = 2.0/(double)x;
174  q0 = w; z = w+h; q1 = w*z - 1.0; k=1;
175  while(q1<1.0e9) {
176  k += 1; z += h;
177  tmp = z*q1 - q0;
178  q0 = q1;
179  q1 = tmp;
180  }
181  m = n+n;
182  for(t=zero, i = 2*(n+k); i>=m; i -= 2) t = one/(i/x-t);
183  a = t;
184  b = one;
185  /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
186  * Hence, if n*(log(2n/x)) > ...
187  * single 8.8722839355e+01
188  * double 7.09782712893383973096e+02
189  * long double 1.1356523406294143949491931077970765006170e+04
190  * then recurrent value may overflow and the result is
191  * likely underflow to zero
192  */
193  tmp = n;
194  v = two/x;
195  tmp = tmp*__ieee754_log(fabs(v*tmp));
196  if(tmp<7.09782712893383973096e+02) {
197  for(i=n-1,di=(double)(i+i);i>0;i--){
198  temp = b;
199  b *= di;
200  b = b/x - a;
201  a = temp;
202  di -= two;
203  }
204  } else {
205  for(i=n-1,di=(double)(i+i);i>0;i--){
206  temp = b;
207  b *= di;
208  b = b/x - a;
209  a = temp;
210  di -= two;
211  /* scale b to avoid spurious overflow */
212  if(b>1e100) {
213  a /= b;
214  t /= b;
215  b = one;
216  }
217  }
218  }
219  b = (t*__ieee754_j0(x)/b);
220  }
221  }
222  if(sgn==1) return -b; else return b;
223 }
224 
225 double
226 __ieee754_yn(int n, double x)
227 {
228  int32_t i,hx,ix,lx;
229  int32_t sign;
230  double a, b, temp;
231 
232  temp = 0;
233  EXTRACT_WORDS(hx,lx,x);
234  ix = 0x7fffffff&hx;
235  /* if Y(n,NaN) is NaN */
236  if((ix|((uint32_t)(lx|-lx))>>31)>0x7ff00000) return x+x;
237  if((ix|lx)==0) return -one/zero;
238  if(hx<0) return zero/zero;
239  sign = 1;
240  if(n<0){
241  n = -n;
242  sign = 1 - ((n&1)<<1);
243  }
244  if(n==0) return(__ieee754_y0(x));
245  if(n==1) return(sign*__ieee754_y1(x));
246  if(ix==0x7ff00000) return zero;
247  if(ix>=0x52D00000) { /* x > 2**302 */
248  /* (x >> n**2)
249  * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
250  * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
251  * Let s=sin(x), c=cos(x),
252  * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
253  *
254  * n sin(xn)*sqt2 cos(xn)*sqt2
255  * ----------------------------------
256  * 0 s-c c+s
257  * 1 -s-c -c+s
258  * 2 -s+c -c-s
259  * 3 s+c c-s
260  */
261  switch(n&3) {
262  case 0: temp = sin(x)-cos(x); break;
263  case 1: temp = -sin(x)-cos(x); break;
264  case 2: temp = -sin(x)+cos(x); break;
265  case 3: temp = sin(x)+cos(x); break;
266  }
267  b = invsqrtpi*temp/sqrt(x);
268  } else {
269  uint32_t high;
270  a = __ieee754_y0(x);
271  b = __ieee754_y1(x);
272  /* quit if b is -inf */
273  GET_HIGH_WORD(high,b);
274  for(i=1;i<n&&high!=0xfff00000;i++){
275  temp = b;
276  b = ((double)(i+i)/x)*b - a;
277  GET_HIGH_WORD(high,b);
278  a = temp;
279  }
280  }
281  if(sign>0) return b; else return -b;
282 }
283 #endif
284