ReactOS 0.4.15-dev-6680-g8c76870
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}
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
static double one
Definition: jn_yn.c:53
double __ieee754_j0(double)
Definition: j0_y0.c:102
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
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 a
Definition: ke_i.h:78
#define c
Definition: ke_i.h:80
#define b
Definition: ke_i.h:79
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
long sgn(REAL x)
Definition: varray.cc:48

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;
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}
double __ieee754_y0(double)
Definition: j0_y0.c:179
double __ieee754_y1(double)
Definition: j1_y1.c:182
#define GET_HIGH_WORD(i, d)
Definition: ieee754.h:26
#define HUGE_VAL
Definition: math.h:51
#define sign(x)
Definition: mapdesc.cc:613

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().