ReactOS  0.4.13-dev-39-g8b6696f
jn_yn.c File Reference
#include "math.h"
#include "ieee754.h"
Include dependency graph for jn_yn.c:

Go to the source code of this file.

Functions

double __ieee754_jn (int n, double x)
 
double __ieee754_yn (int n, double x)
 

Variables

static double invsqrtpi = 5.64189583547756279280e-01
 
static double two = 2.00000000000000000000e+00
 
static double one = 1.00000000000000000000e+00
 
static double zero = 0.00000000000000000000e+00
 

Function Documentation

◆ __ieee754_jn()

double __ieee754_jn ( int  n,
double  x 
)

Definition at line 64 of file jn_yn.c.

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 }
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 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)
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 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
int k
Definition: mpi.c:3369

Referenced by _jn().

◆ __ieee754_yn()

double __ieee754_yn ( int  n,
double  x 
)

Definition at line 227 of file jn_yn.c.

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:64
double __ieee754_y0(double)
Definition: j0_y0.c:179
static __inline double __ieee754_sqrt(double x)
Definition: ieee754.h:40
GLdouble n
Definition: glext.h:7729
GLint GLint GLint GLint GLint x
Definition: gl.h:1548
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
#define a
Definition: ke_i.h:78
static const char mbstate_t *static wchar_t const char mbstate_t *static const wchar_t int *static double
Definition: string.c:80
#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
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
double __ieee754_y1(double)
Definition: j1_y1.c:182
static __inline void __sincos(double x, double *s, double *c)
Definition: ieee754.h:43
static calc_node_t temp
Definition: rpn_ieee.c:38
#define HUGE_VAL
Definition: math.h:51
#define c
Definition: ke_i.h:80
GLboolean GLboolean GLboolean GLboolean a
Definition: glext.h:6204

Referenced by _yn().

Variable Documentation

◆ invsqrtpi

double invsqrtpi = 5.64189583547756279280e-01
static

Definition at line 51 of file jn_yn.c.

Referenced by __ieee754_jn(), and __ieee754_yn().

◆ one

double one = 1.00000000000000000000e+00
static

Definition at line 53 of file jn_yn.c.

Referenced by __ieee754_jn().

◆ two

◆ zero

double zero = 0.00000000000000000000e+00
static

Definition at line 58 of file jn_yn.c.

Referenced by __ieee754_jn(), and __ieee754_yn().