ReactOS 0.4.15-dev-7842-g558ab78
ieee754.h File Reference
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes

union  ieee_double_shape_type
 

Macros

#define EXTRACT_WORDS(ix0, ix1, d)
 
#define GET_HIGH_WORD(i, d)
 
#define GET_LOW_WORD(i, d)
 

Typedefs

typedef __int32 int32_t
 
typedef unsigned __int32 u_int32_t
 

Functions

static __inline double __ieee754_sqrt (double x)
 
static __inline double __ieee754_log (double x)
 
static __inline double __cos (double x)
 
static __inline void __sincos (double x, double *s, double *c)
 
double __ieee754_j0 (double)
 
double __ieee754_j1 (double)
 
double __ieee754_jn (int, double)
 
double __ieee754_y0 (double)
 
double __ieee754_y1 (double)
 
double __ieee754_yn (int, double)
 

Macro Definition Documentation

◆ EXTRACT_WORDS

#define EXTRACT_WORDS (   ix0,
  ix1,
  d 
)
Value:
do { \
ew_u.value = (d); \
(ix0) = ew_u.parts.msw; \
(ix1) = ew_u.parts.lsw; \
} while (0)
#define d
Definition: ke_i.h:81
struct ieee_double_shape_type::@4305 parts

Definition at line 16 of file ieee754.h.

◆ GET_HIGH_WORD

#define GET_HIGH_WORD (   i,
  d 
)
Value:
do { \
gh_u.value = (d); \
(i) = gh_u.parts.msw; \
} while (0)
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

Definition at line 26 of file ieee754.h.

◆ GET_LOW_WORD

#define GET_LOW_WORD (   i,
  d 
)
Value:
do { \
gl_u.value = (d); \
(i) = gl_u.parts.lsw; \
} while (0)

Definition at line 33 of file ieee754.h.

Typedef Documentation

◆ int32_t

typedef __int32 int32_t

Definition at line 3 of file ieee754.h.

◆ u_int32_t

Definition at line 4 of file ieee754.h.

Function Documentation

◆ __cos()

static __inline double __cos ( double  x)
static

Definition at line 42 of file ieee754.h.

42{return cos(x);}
_STLP_DECLSPEC complex< float > _STLP_CALL cos(const complex< float > &)
GLint GLint GLint GLint GLint x
Definition: gl.h:1548

Referenced by __ieee754_j0(), __ieee754_j1(), __ieee754_y0(), and __ieee754_y1().

◆ __ieee754_j0()

double __ieee754_j0 ( double  x)

Definition at line 102 of file j0_y0.c.

105{
106 double z, s,c,ss,cc,r,u,v,r1,r2,s1,s2,z2,z4;
107 int32_t hx,ix;
108
109 GET_HIGH_WORD(hx,x);
110 ix = hx&0x7fffffff;
111 if(ix>=0x7ff00000) return one/(x*x);
112 x = fabs(x);
113 if(ix >= 0x40000000) { /* |x| >= 2.0 */
114 __sincos (x, &s, &c);
115 ss = s-c;
116 cc = s+c;
117 if(ix<0x7fe00000) { /* make sure x+x not overflow */
118 z = -__cos(x+x);
119 if ((s*c)<zero) cc = z/ss;
120 else ss = z/cc;
121 }
122 /*
123 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
124 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
125 */
126 if(ix>0x48000000) z = (invsqrtpi*cc)/__ieee754_sqrt(x);
127 else {
128 u = pzero(x); v = qzero(x);
130 }
131 return z;
132 }
133 if(ix<0x3f200000) { /* |x| < 2**-13 */
134 if(huge+x>one) { /* raise inexact if x != 0 */
135 if(ix<0x3e400000) return one; /* |x|<2**-27 */
136 else return one - 0.25*x*x;
137 }
138 }
139 z = x*x;
140#ifdef DO_NOT_USE_THIS
141 r = z*(R02+z*(R03+z*(R04+z*R05)));
142 s = one+z*(S01+z*(S02+z*(S03+z*S04)));
143#else
144 r1 = z*R[2]; z2=z*z;
145 r2 = R[3]+z*R[4]; z4=z2*z2;
146 r = r1 + z2*r2 + z4*R[5];
147 s1 = one+z*S[1];
148 s2 = S[2]+z*S[3];
149 s = s1 + z2*s2 + z4*S[4];
150#endif
151 if(ix < 0x3FF00000) { /* |x| < 1.00 */
152 return one + z*(-0.25+(r/s));
153 } else {
154 u = 0.5*x;
155 return((one+u)*(one-u)+z*(r/s));
156 }
157}
INT32 int32_t
Definition: types.h:71
const GLdouble * v
Definition: gl.h:2040
GLdouble s
Definition: gl.h:2039
GLdouble GLdouble GLdouble r
Definition: gl.h:2055
const GLubyte * c
Definition: glext.h:8905
GLdouble GLdouble z
Definition: glext.h:5874
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 * u
Definition: glfuncs.h:240
#define ss
Definition: i386-dis.c:441
static double huge
Definition: j0_y0.c:79
static double R[]
Definition: j0_y0.c:84
static double zero
Definition: j0_y0.c:96
static double invsqrtpi
Definition: j0_y0.c:81
static double qzero()
static double one
Definition: j0_y0.c:80
static double pzero()
static __inline void __sincos(double x, double *s, double *c)
Definition: ieee754.h:43
static __inline double __cos(double x)
Definition: ieee754.h:42
#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
uint32_t cc
Definition: isohybrid.c:75
#define c
Definition: ke_i.h:80
struct S1 s1
struct S2 s2
static DNS_RECORDW r1
Definition: record.c:37
static DNS_RECORDW r2
Definition: record.c:38
Definition: movable.cpp:9

Referenced by __ieee754_jn(), __ieee754_y0(), and _j0().

◆ __ieee754_j1()

double __ieee754_j1 ( double  x)

Definition at line 103 of file j1_y1.c.

106{
107 double z, s,c,ss,cc,r,u,v,y,r1,r2,s1,s2,s3,z2,z4;
108 int32_t hx,ix;
109
110 GET_HIGH_WORD(hx,x);
111 ix = hx&0x7fffffff;
112 if(ix>=0x7ff00000) return one/x;
113 y = fabs(x);
114 if(ix >= 0x40000000) { /* |x| >= 2.0 */
115 __sincos (y, &s, &c);
116 ss = -s-c;
117 cc = s-c;
118 if(ix<0x7fe00000) { /* make sure y+y not overflow */
119 z = __cos(y+y);
120 if ((s*c)>zero) cc = z/ss;
121 else ss = z/cc;
122 }
123 /*
124 * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
125 * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)
126 */
127 if(ix>0x48000000) z = (invsqrtpi*cc)/__ieee754_sqrt(y);
128 else {
129 u = pone(y); v = qone(y);
131 }
132 if(hx<0) return -z;
133 else return z;
134 }
135 if(ix<0x3e400000) { /* |x|<2**-27 */
136 if(huge+x>one) return 0.5*x;/* inexact if x!=0 necessary */
137 }
138 z = x*x;
139#ifdef DO_NOT_USE_THIS
140 r = z*(r00+z*(r01+z*(r02+z*r03)));
141 s = one+z*(s01+z*(s02+z*(s03+z*(s04+z*s05))));
142 r *= x;
143#else
144 r1 = z*R[0]; z2=z*z;
145 r2 = R[1]+z*R[2]; z4=z2*z2;
146 r = r1 + z2*r2 + z4*R[3];
147 r *= x;
148 s1 = one+z*S[1];
149 s2 = S[2]+z*S[3];
150 s3 = S[4]+z*S[5];
151 s = s1 + z2*s2 + z4*s3;
152#endif
153 return(x*0.5+r/s);
154}
GLint GLint GLint GLint GLint GLint y
Definition: gl.h:1548
static double huge
Definition: j1_y1.c:79
static double pone()
static double R[]
Definition: j1_y1.c:84
static double zero
Definition: j1_y1.c:97
static double invsqrtpi
Definition: j1_y1.c:81
static double qone()
static double one
Definition: j1_y1.c:80

Referenced by __ieee754_jn(), __ieee754_y1(), and _j1().

◆ __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}
GLdouble GLdouble t
Definition: gl.h:2047
GLdouble n
Definition: glext.h:7729
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
const GLfloat * m
Definition: glext.h:10848
GLfloat GLfloat GLfloat GLfloat h
Definition: glext.h:7723
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
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 a
Definition: ke_i.h:78
#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_log()

static __inline double __ieee754_log ( double  x)
static

Definition at line 41 of file ieee754.h.

41{return log(x);}
#define log(outFile, fmt,...)
Definition: util.h:15

Referenced by __ieee754_jn(), __ieee754_y0(), and __ieee754_y1().

◆ __ieee754_sqrt()

static __inline double __ieee754_sqrt ( double  x)
static

Definition at line 40 of file ieee754.h.

40{return sqrt(x);}
_STLP_DECLSPEC complex< float > _STLP_CALL sqrt(const complex< float > &)
Definition: complex.cpp:188

Referenced by __ieee754_j0(), __ieee754_j1(), __ieee754_jn(), __ieee754_y0(), __ieee754_y1(), and __ieee754_yn().

◆ __ieee754_y0()

double __ieee754_y0 ( double  x)

Definition at line 179 of file j0_y0.c.

182{
183 double z, s,c,ss,cc,u,v,z2,z4,z6,u1,u2,u3,v1,v2;
184 int32_t hx,ix,lx;
185
186 EXTRACT_WORDS(hx,lx,x);
187 ix = 0x7fffffff&hx;
188 /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0, y0(0) is -inf. */
189 if(ix>=0x7ff00000) return one/(x+x*x);
190 if((ix|lx)==0) return -HUGE_VAL+x; /* -inf and overflow exception. */
191 if(hx<0) return zero/(zero*x);
192 if(ix >= 0x40000000) { /* |x| >= 2.0 */
193 /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
194 * where x0 = x-pi/4
195 * Better formula:
196 * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
197 * = 1/sqrt(2) * (sin(x) + cos(x))
198 * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
199 * = 1/sqrt(2) * (sin(x) - cos(x))
200 * To avoid cancellation, use
201 * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
202 * to compute the worse one.
203 */
204 __sincos (x, &s, &c);
205 ss = s-c;
206 cc = s+c;
207 /*
208 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
209 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
210 */
211 if(ix<0x7fe00000) { /* make sure x+x not overflow */
212 z = -__cos(x+x);
213 if ((s*c)<zero) cc = z/ss;
214 else ss = z/cc;
215 }
216 if(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrt(x);
217 else {
218 u = pzero(x); v = qzero(x);
220 }
221 return z;
222 }
223 if(ix<=0x3e400000) { /* x < 2**-27 */
224 return(U[0] + tpi*__ieee754_log(x));
225 }
226 z = x*x;
227#ifdef DO_NOT_USE_THIS
228 u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06)))));
229 v = one+z*(v01+z*(v02+z*(v03+z*v04)));
230#else
231 u1 = U[0]+z*U[1]; z2=z*z;
232 u2 = U[2]+z*U[3]; z4=z2*z2;
233 u3 = U[4]+z*U[5]; z6=z4*z2;
234 u = u1 + z2*u2 + z4*u3 + z6*U[6];
235 v1 = one+z*V[0];
236 v2 = V[1]+z*V[2];
237 v = v1 + z2*v2 + z4*V[3];
238#endif
239 return(u/v + tpi*(__ieee754_j0(x)*__ieee754_log(x)));
240}
GLdouble GLdouble u2
Definition: glext.h:8308
GLfloat GLfloat v1
Definition: glext.h:6062
GLfloat GLfloat GLfloat v2
Definition: glext.h:6063
GLdouble u1
Definition: glext.h:8308
static double U[]
Definition: j0_y0.c:164
static double tpi
Definition: j0_y0.c:82
double __ieee754_j0(double x)
Definition: j0_y0.c:102
#define HUGE_VAL
Definition: math.h:51
static BYTE u3[]
Definition: msg.c:580

Referenced by __ieee754_yn(), and _y0().

◆ __ieee754_y1()

double __ieee754_y1 ( double  x)

Definition at line 182 of file j1_y1.c.

185{
186 double z, s,c,ss,cc,u,v,u1,u2,v1,v2,v3,z2,z4;
187 int32_t hx,ix,lx;
188
189 EXTRACT_WORDS(hx,lx,x);
190 ix = 0x7fffffff&hx;
191 /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
192 if(ix>=0x7ff00000) return one/(x+x*x);
193 if((ix|lx)==0) return -HUGE_VAL+x; /* -inf and overflow exception. */;
194 if(hx<0) return zero/(zero*x);
195 if(ix >= 0x40000000) { /* |x| >= 2.0 */
196 __sincos (x, &s, &c);
197 ss = -s-c;
198 cc = s-c;
199 if(ix<0x7fe00000) { /* make sure x+x not overflow */
200 z = __cos(x+x);
201 if ((s*c)>zero) cc = z/ss;
202 else ss = z/cc;
203 }
204 /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0))
205 * where x0 = x-3pi/4
206 * Better formula:
207 * cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
208 * = 1/sqrt(2) * (sin(x) - cos(x))
209 * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
210 * = -1/sqrt(2) * (cos(x) + sin(x))
211 * To avoid cancellation, use
212 * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
213 * to compute the worse one.
214 */
215 if(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrt(x);
216 else {
217 u = pone(x); v = qone(x);
219 }
220 return z;
221 }
222 if(ix<=0x3c900000) { /* x < 2**-54 */
223 return(-tpi/x);
224 }
225 z = x*x;
226#ifdef DO_NOT_USE_THIS
227 u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4])));
228 v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4]))));
229#else
230 u1 = U0[0]+z*U0[1];z2=z*z;
231 u2 = U0[2]+z*U0[3];z4=z2*z2;
232 u = u1 + z2*u2 + z4*U0[4];
233 v1 = one+z*V0[0];
234 v2 = V0[1]+z*V0[2];
235 v3 = V0[3]+z*V0[4];
236 v = v1 + z2*v2 + z4*v3;
237#endif
238 return(x*(u/v) + tpi*(__ieee754_j1(x)*__ieee754_log(x)-one/x));
239}
GLfloat GLfloat GLfloat GLfloat v3
Definition: glext.h:6064
static double U0[5]
Definition: j1_y1.c:159
static double V0[5]
Definition: j1_y1.c:170
static double tpi
Definition: j1_y1.c:82
double __ieee754_j1(double x)
Definition: j1_y1.c:103

Referenced by __ieee754_yn(), and _y1().

◆ __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 sign(x)
Definition: mapdesc.cc:613

Referenced by _yn().

◆ __sincos()

static __inline void __sincos ( double  x,
double s,
double c 
)
static

Definition at line 43 of file ieee754.h.

44{
45 *s = sin(x);
46 *c = cos(x);
47}
_STLP_DECLSPEC complex< float > _STLP_CALL sin(const complex< float > &)

Referenced by __ieee754_j0(), __ieee754_j1(), __ieee754_jn(), __ieee754_y0(), __ieee754_y1(), and __ieee754_yn().