ReactOS 0.4.15-dev-8510-g6f86c05
j0_y0.c
Go to the documentation of this file.
1/* @(#)e_j0.c 5.1 93/09/24 */
2/*
3 * ====================================================
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/* Modified by Naohiko Shimizu/Tokai University, Japan 1997/08/26,
13 for performance improvement on pipelined processors.
14*/
15
16#if defined(LIBM_SCCS) && !defined(lint)
17static char rcsid[] = "\$NetBSD: e_j0.c,v 1.8 1995/05/10 20:45:23 jtc Exp \$";
18#endif
19
20/* __ieee754_j0(x), __ieee754_y0(x)
21 * Bessel function of the first and second kinds of order zero.
22 * Method -- j0(x):
23 * 1. For tiny x, we use j0(x) = 1 - x^2/4 + x^4/64 - ...
24 * 2. Reduce x to |x| since j0(x)=j0(-x), and
25 * for x in (0,2)
26 * j0(x) = 1-z/4+ z^2*R0/S0, where z = x*x;
27 * (precision: |j0-1+z/4-z^2R0/S0 |<2**-63.67 )
28 * for x in (2,inf)
29 * j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
30 * where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
31 * as follow:
32 * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
33 * = 1/sqrt(2) * (cos(x) + sin(x))
34 * sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
35 * = 1/sqrt(2) * (sin(x) - cos(x))
36 * (To avoid cancellation, use
37 * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
38 * to compute the worse one.)
39 *
40 * 3 Special cases
41 * j0(nan)= nan
42 * j0(0) = 1
43 * j0(inf) = 0
44 *
45 * Method -- y0(x):
46 * 1. For x<2.
47 * Since
48 * y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x^2/4 - ...)
49 * therefore y0(x)-2/pi*j0(x)*ln(x) is an even function.
50 * We use the following function to approximate y0,
51 * y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x^2
52 * where
53 * U(z) = u00 + u01*z + ... + u06*z^6
54 * V(z) = 1 + v01*z + ... + v04*z^4
55 * with absolute approximation error bounded by 2**-72.
56 * Note: For tiny x, U/V = u0 and j0(x)~1, hence
57 * y0(tiny) = u0 + (2/pi)*ln(tiny), (choose tiny<2**-27)
58 * 2. For x>=2.
59 * y0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)+q0(x)*sin(x0))
60 * where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
61 * by the method mentioned above.
62 * 3. Special cases: y0(0)=-inf, y0(x<0)=NaN, y0(inf)=0.
63 */
64
65#include "math.h"
66#include "ieee754.h"
67
68#ifdef __STDC__
69static double pzero(double), qzero(double);
70#else
71static double pzero(), qzero();
72#endif
73
74#ifdef __STDC__
75static const double
76#else
77static double
78#endif
79huge = 1e300,
80one = 1.0,
81invsqrtpi= 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
82tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
83 /* R0/S0 on [0, 2.00] */
84R[] = {0.0, 0.0, 1.56249999999999947958e-02, /* 0x3F8FFFFF, 0xFFFFFFFD */
85 -1.89979294238854721751e-04, /* 0xBF28E6A5, 0xB61AC6E9 */
86 1.82954049532700665670e-06, /* 0x3EBEB1D1, 0x0C503919 */
87 -4.61832688532103189199e-09}, /* 0xBE33D5E7, 0x73D63FCE */
88S[] = {0.0, 1.56191029464890010492e-02, /* 0x3F8FFCE8, 0x82C8C2A4 */
89 1.16926784663337450260e-04, /* 0x3F1EA6D2, 0xDD57DBF4 */
90 5.13546550207318111446e-07, /* 0x3EA13B54, 0xCE84D5A9 */
91 1.16614003333790000205e-09}; /* 0x3E1408BC, 0xF4745D8F */
92
93#ifdef __STDC__
94static const double zero = 0.0;
95#else
96static double zero = 0.0;
97#endif
98
99#ifdef __STDC__
100 double __ieee754_j0(double x)
101#else
103 double x;
104#endif
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}
158
159#ifdef __STDC__
160static const double
161#else
162static double
163#endif
164U[] = {-7.38042951086872317523e-02, /* 0xBFB2E4D6, 0x99CBD01F */
165 1.76666452509181115538e-01, /* 0x3FC69D01, 0x9DE9E3FC */
166 -1.38185671945596898896e-02, /* 0xBF8C4CE8, 0xB16CFA97 */
167 3.47453432093683650238e-04, /* 0x3F36C54D, 0x20B29B6B */
168 -3.81407053724364161125e-06, /* 0xBECFFEA7, 0x73D25CAD */
169 1.95590137035022920206e-08, /* 0x3E550057, 0x3B4EABD4 */
170 -3.98205194132103398453e-11}, /* 0xBDC5E43D, 0x693FB3C8 */
171V[] = {1.27304834834123699328e-02, /* 0x3F8A1270, 0x91C9C71A */
172 7.60068627350353253702e-05, /* 0x3F13ECBB, 0xF578C6C1 */
173 2.59150851840457805467e-07, /* 0x3E91642D, 0x7FF202FD */
174 4.41110311332675467403e-10}; /* 0x3DFE5018, 0x3BD6D9EF */
175
176#ifdef __STDC__
177 double __ieee754_y0(double x)
178#else
180 double x;
181#endif
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}
241
242/* The asymptotic expansions of pzero is
243 * 1 - 9/128 s^2 + 11025/98304 s^4 - ..., where s = 1/x.
244 * For x >= 2, We approximate pzero by
245 * pzero(x) = 1 + (R/S)
246 * where R = pR0 + pR1*s^2 + pR2*s^4 + ... + pR5*s^10
247 * S = 1 + pS0*s^2 + ... + pS4*s^10
248 * and
249 * | pzero(x)-1-R/S | <= 2 ** ( -60.26)
250 */
251#ifdef __STDC__
252static const double pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
253#else
254static double pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
255#endif
256 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
257 -7.03124999999900357484e-02, /* 0xBFB1FFFF, 0xFFFFFD32 */
258 -8.08167041275349795626e+00, /* 0xC02029D0, 0xB44FA779 */
259 -2.57063105679704847262e+02, /* 0xC0701102, 0x7B19E863 */
260 -2.48521641009428822144e+03, /* 0xC0A36A6E, 0xCD4DCAFC */
261 -5.25304380490729545272e+03, /* 0xC0B4850B, 0x36CC643D */
262};
263#ifdef __STDC__
264static const double pS8[5] = {
265#else
266static double pS8[5] = {
267#endif
268 1.16534364619668181717e+02, /* 0x405D2233, 0x07A96751 */
269 3.83374475364121826715e+03, /* 0x40ADF37D, 0x50596938 */
270 4.05978572648472545552e+04, /* 0x40E3D2BB, 0x6EB6B05F */
271 1.16752972564375915681e+05, /* 0x40FC810F, 0x8F9FA9BD */
272 4.76277284146730962675e+04, /* 0x40E74177, 0x4F2C49DC */
273};
274
275#ifdef __STDC__
276static const double pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
277#else
278static double pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
279#endif
280 -1.14125464691894502584e-11, /* 0xBDA918B1, 0x47E495CC */
281 -7.03124940873599280078e-02, /* 0xBFB1FFFF, 0xE69AFBC6 */
282 -4.15961064470587782438e+00, /* 0xC010A370, 0xF90C6BBF */
283 -6.76747652265167261021e+01, /* 0xC050EB2F, 0x5A7D1783 */
284 -3.31231299649172967747e+02, /* 0xC074B3B3, 0x6742CC63 */
285 -3.46433388365604912451e+02, /* 0xC075A6EF, 0x28A38BD7 */
286};
287#ifdef __STDC__
288static const double pS5[5] = {
289#else
290static double pS5[5] = {
291#endif
292 6.07539382692300335975e+01, /* 0x404E6081, 0x0C98C5DE */
293 1.05125230595704579173e+03, /* 0x40906D02, 0x5C7E2864 */
294 5.97897094333855784498e+03, /* 0x40B75AF8, 0x8FBE1D60 */
295 9.62544514357774460223e+03, /* 0x40C2CCB8, 0xFA76FA38 */
296 2.40605815922939109441e+03, /* 0x40A2CC1D, 0xC70BE864 */
297};
298
299#ifdef __STDC__
300static const double pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
301#else
302static double pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
303#endif
304 -2.54704601771951915620e-09, /* 0xBE25E103, 0x6FE1AA86 */
305 -7.03119616381481654654e-02, /* 0xBFB1FFF6, 0xF7C0E24B */
306 -2.40903221549529611423e+00, /* 0xC00345B2, 0xAEA48074 */
307 -2.19659774734883086467e+01, /* 0xC035F74A, 0x4CB94E14 */
308 -5.80791704701737572236e+01, /* 0xC04D0A22, 0x420A1A45 */
309 -3.14479470594888503854e+01, /* 0xC03F72AC, 0xA892D80F */
310};
311#ifdef __STDC__
312static const double pS3[5] = {
313#else
314static double pS3[5] = {
315#endif
316 3.58560338055209726349e+01, /* 0x4041ED92, 0x84077DD3 */
317 3.61513983050303863820e+02, /* 0x40769839, 0x464A7C0E */
318 1.19360783792111533330e+03, /* 0x4092A66E, 0x6D1061D6 */
319 1.12799679856907414432e+03, /* 0x40919FFC, 0xB8C39B7E */
320 1.73580930813335754692e+02, /* 0x4065B296, 0xFC379081 */
321};
322
323#ifdef __STDC__
324static const double pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
325#else
326static double pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
327#endif
328 -8.87534333032526411254e-08, /* 0xBE77D316, 0xE927026D */
329 -7.03030995483624743247e-02, /* 0xBFB1FF62, 0x495E1E42 */
330 -1.45073846780952986357e+00, /* 0xBFF73639, 0x8A24A843 */
331 -7.63569613823527770791e+00, /* 0xC01E8AF3, 0xEDAFA7F3 */
332 -1.11931668860356747786e+01, /* 0xC02662E6, 0xC5246303 */
333 -3.23364579351335335033e+00, /* 0xC009DE81, 0xAF8FE70F */
334};
335#ifdef __STDC__
336static const double pS2[5] = {
337#else
338static double pS2[5] = {
339#endif
340 2.22202997532088808441e+01, /* 0x40363865, 0x908B5959 */
341 1.36206794218215208048e+02, /* 0x4061069E, 0x0EE8878F */
342 2.70470278658083486789e+02, /* 0x4070E786, 0x42EA079B */
343 1.53875394208320329881e+02, /* 0x40633C03, 0x3AB6FAFF */
344 1.46576176948256193810e+01, /* 0x402D50B3, 0x44391809 */
345};
346
347#ifdef __STDC__
348 static double pzero(double x)
349#else
350 static double pzero(x)
351 double x;
352#endif
353{
354#ifdef __STDC__
355 const double *p = 0,*q = 0;
356#else
357 double *p = 0,*q = 0;
358#endif
359 double z,r,s,z2,z4,r1,r2,r3,s1,s2,s3;
360 int32_t ix;
361 GET_HIGH_WORD(ix,x);
362 ix &= 0x7fffffff;
363 if(ix>=0x40200000) {p = pR8; q= pS8;}
364 else if(ix>=0x40122E8B){p = pR5; q= pS5;}
365 else if(ix>=0x4006DB6D){p = pR3; q= pS3;}
366 else if(ix>=0x40000000){p = pR2; q= pS2;}
367 z = one/(x*x);
368#ifdef DO_NOT_USE_THIS
369 r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
370 s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
371#else
372 r1 = p[0]+z*p[1]; z2=z*z;
373 r2 = p[2]+z*p[3]; z4=z2*z2;
374 r3 = p[4]+z*p[5];
375 r = r1 + z2*r2 + z4*r3;
376 s1 = one+z*q[0];
377 s2 = q[1]+z*q[2];
378 s3 = q[3]+z*q[4];
379 s = s1 + z2*s2 + z4*s3;
380#endif
381 return one+ r/s;
382}
383
384
385/* For x >= 8, the asymptotic expansions of qzero is
386 * -1/8 s + 75/1024 s^3 - ..., where s = 1/x.
387 * We approximate pzero by
388 * qzero(x) = s*(-1.25 + (R/S))
389 * where R = qR0 + qR1*s^2 + qR2*s^4 + ... + qR5*s^10
390 * S = 1 + qS0*s^2 + ... + qS5*s^12
391 * and
392 * | qzero(x)/s +1.25-R/S | <= 2 ** ( -61.22)
393 */
394#ifdef __STDC__
395static const double qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
396#else
397static double qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
398#endif
399 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
400 7.32421874999935051953e-02, /* 0x3FB2BFFF, 0xFFFFFE2C */
401 1.17682064682252693899e+01, /* 0x40278952, 0x5BB334D6 */
402 5.57673380256401856059e+02, /* 0x40816D63, 0x15301825 */
403 8.85919720756468632317e+03, /* 0x40C14D99, 0x3E18F46D */
404 3.70146267776887834771e+04, /* 0x40E212D4, 0x0E901566 */
405};
406#ifdef __STDC__
407static const double qS8[6] = {
408#else
409static double qS8[6] = {
410#endif
411 1.63776026895689824414e+02, /* 0x406478D5, 0x365B39BC */
412 8.09834494656449805916e+03, /* 0x40BFA258, 0x4E6B0563 */
413 1.42538291419120476348e+05, /* 0x41016652, 0x54D38C3F */
414 8.03309257119514397345e+05, /* 0x412883DA, 0x83A52B43 */
415 8.40501579819060512818e+05, /* 0x4129A66B, 0x28DE0B3D */
416 -3.43899293537866615225e+05, /* 0xC114FD6D, 0x2C9530C5 */
417};
418
419#ifdef __STDC__
420static const double qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
421#else
422static double qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
423#endif
424 1.84085963594515531381e-11, /* 0x3DB43D8F, 0x29CC8CD9 */
425 7.32421766612684765896e-02, /* 0x3FB2BFFF, 0xD172B04C */
426 5.83563508962056953777e+00, /* 0x401757B0, 0xB9953DD3 */
427 1.35111577286449829671e+02, /* 0x4060E392, 0x0A8788E9 */
428 1.02724376596164097464e+03, /* 0x40900CF9, 0x9DC8C481 */
429 1.98997785864605384631e+03, /* 0x409F17E9, 0x53C6E3A6 */
430};
431#ifdef __STDC__
432static const double qS5[6] = {
433#else
434static double qS5[6] = {
435#endif
436 8.27766102236537761883e+01, /* 0x4054B1B3, 0xFB5E1543 */
437 2.07781416421392987104e+03, /* 0x40A03BA0, 0xDA21C0CE */
438 1.88472887785718085070e+04, /* 0x40D267D2, 0x7B591E6D */
439 5.67511122894947329769e+04, /* 0x40EBB5E3, 0x97E02372 */
440 3.59767538425114471465e+04, /* 0x40E19118, 0x1F7A54A0 */
441 -5.35434275601944773371e+03, /* 0xC0B4EA57, 0xBEDBC609 */
442};
443
444#ifdef __STDC__
445static const double qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
446#else
447static double qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
448#endif
449 4.37741014089738620906e-09, /* 0x3E32CD03, 0x6ADECB82 */
450 7.32411180042911447163e-02, /* 0x3FB2BFEE, 0x0E8D0842 */
451 3.34423137516170720929e+00, /* 0x400AC0FC, 0x61149CF5 */
452 4.26218440745412650017e+01, /* 0x40454F98, 0x962DAEDD */
453 1.70808091340565596283e+02, /* 0x406559DB, 0xE25EFD1F */
454 1.66733948696651168575e+02, /* 0x4064D77C, 0x81FA21E0 */
455};
456#ifdef __STDC__
457static const double qS3[6] = {
458#else
459static double qS3[6] = {
460#endif
461 4.87588729724587182091e+01, /* 0x40486122, 0xBFE343A6 */
462 7.09689221056606015736e+02, /* 0x40862D83, 0x86544EB3 */
463 3.70414822620111362994e+03, /* 0x40ACF04B, 0xE44DFC63 */
464 6.46042516752568917582e+03, /* 0x40B93C6C, 0xD7C76A28 */
465 2.51633368920368957333e+03, /* 0x40A3A8AA, 0xD94FB1C0 */
466 -1.49247451836156386662e+02, /* 0xC062A7EB, 0x201CF40F */
467};
468
469#ifdef __STDC__
470static const double qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
471#else
472static double qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
473#endif
474 1.50444444886983272379e-07, /* 0x3E84313B, 0x54F76BDB */
475 7.32234265963079278272e-02, /* 0x3FB2BEC5, 0x3E883E34 */
476 1.99819174093815998816e+00, /* 0x3FFFF897, 0xE727779C */
477 1.44956029347885735348e+01, /* 0x402CFDBF, 0xAAF96FE5 */
478 3.16662317504781540833e+01, /* 0x403FAA8E, 0x29FBDC4A */
479 1.62527075710929267416e+01, /* 0x403040B1, 0x71814BB4 */
480};
481#ifdef __STDC__
482static const double qS2[6] = {
483#else
484static double qS2[6] = {
485#endif
486 3.03655848355219184498e+01, /* 0x403E5D96, 0xF7C07AED */
487 2.69348118608049844624e+02, /* 0x4070D591, 0xE4D14B40 */
488 8.44783757595320139444e+02, /* 0x408A6645, 0x22B3BF22 */
489 8.82935845112488550512e+02, /* 0x408B977C, 0x9C5CC214 */
490 2.12666388511798828631e+02, /* 0x406A9553, 0x0E001365 */
491 -5.31095493882666946917e+00, /* 0xC0153E6A, 0xF8B32931 */
492};
493
494#ifdef __STDC__
495 static double qzero(double x)
496#else
497 static double qzero(x)
498 double x;
499#endif
500{
501#ifdef __STDC__
502 const double *p = 0,*q = 0;
503#else
504 double *p = 0,*q = 0;
505#endif
506 double s,r,z,z2,z4,z6,r1,r2,r3,s1,s2,s3;
507 int32_t ix;
508 GET_HIGH_WORD(ix,x);
509 ix &= 0x7fffffff;
510 if(ix>=0x40200000) {p = qR8; q= qS8;}
511 else if(ix>=0x40122E8B){p = qR5; q= qS5;}
512 else if(ix>=0x4006DB6D){p = qR3; q= qS3;}
513 else if(ix>=0x40000000){p = qR2; q= qS2;}
514 z = one/(x*x);
515#ifdef DO_NOT_USE_THIS
516 r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
517 s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
518#else
519 r1 = p[0]+z*p[1]; z2=z*z;
520 r2 = p[2]+z*p[3]; z4=z2*z2;
521 r3 = p[4]+z*p[5]; z6=z4*z2;
522 r= r1 + z2*r2 + z4*r3;
523 s1 = one+z*q[0];
524 s2 = q[1]+z*q[2];
525 s3 = q[3]+z*q[4];
526 s = s1 + z2*s2 + z4*s3 +z6*q[5];
527#endif
528 return (-.125 + r/s)/x;
529}
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 GLdouble r
Definition: gl.h:2055
GLdouble GLdouble GLdouble GLdouble q
Definition: gl.h:2063
const GLubyte * c
Definition: glext.h:8905
GLdouble GLdouble u2
Definition: glext.h:8308
GLfloat GLfloat p
Definition: glext.h:8902
GLfloat GLfloat v1
Definition: glext.h:6062
GLfloat GLfloat GLfloat v2
Definition: glext.h:6063
GLdouble u1
Definition: glext.h:8308
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 pS3[5]
Definition: j0_y0.c:314
static double huge
Definition: j0_y0.c:79
static double R[]
Definition: j0_y0.c:84
static double qR5[6]
Definition: j0_y0.c:422
static double qR8[6]
Definition: j0_y0.c:397
static double pR2[6]
Definition: j0_y0.c:326
static double U[]
Definition: j0_y0.c:164
static double zero
Definition: j0_y0.c:96
static double pR5[6]
Definition: j0_y0.c:278
static double invsqrtpi
Definition: j0_y0.c:81
double __ieee754_y0(double x)
Definition: j0_y0.c:179
static double tpi
Definition: j0_y0.c:82
static double pS8[5]
Definition: j0_y0.c:266
static double pR8[6]
Definition: j0_y0.c:254
static double qS2[6]
Definition: j0_y0.c:484
static double qR2[6]
Definition: j0_y0.c:472
static double qzero()
static double one
Definition: j0_y0.c:80
static double qS3[6]
Definition: j0_y0.c:459
static double pS2[5]
Definition: j0_y0.c:338
static double qR3[6]
Definition: j0_y0.c:447
static double qS8[6]
Definition: j0_y0.c:409
static double pS5[5]
Definition: j0_y0.c:290
double __ieee754_j0(double x)
Definition: j0_y0.c:102
static double pR3[6]
Definition: j0_y0.c:302
static double qS5[6]
Definition: j0_y0.c:434
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 EXTRACT_WORDS(ix0, ix1, d)
Definition: ieee754.h:16
static __inline double __ieee754_log(double x)
Definition: ieee754.h:41
#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
#define HUGE_VAL
Definition: math.h:51
uint32_t cc
Definition: isohybrid.c:75
#define c
Definition: ke_i.h:80
struct S1 s1
struct S2 s2
static BYTE u3[]
Definition: msg.c:580
static DNS_RECORDW r3
Definition: record.c:39
static DNS_RECORDW r1
Definition: record.c:37
static DNS_RECORDW r2
Definition: record.c:38
Definition: movable.cpp:9
static TRIO_CONST char rcsid[]
Definition: trio.c:753