ReactOS  0.4.13-dev-257-gfabbd7c
j0_y0.c
Go to the documentation of this file.
1 /* @(#)e_j0.c 5.1 93/09/24 */
2 /*
3  * ====================================================
4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
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)
17 static 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__
69 static double pzero(double), qzero(double);
70 #else
71 static double pzero(), qzero();
72 #endif
73 
74 #ifdef __STDC__
75 static const double
76 #else
77 static double
78 #endif
79 huge = 1e300,
80 one = 1.0,
81 invsqrtpi= 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
82 tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
83  /* R0/S0 on [0, 2.00] */
84 R[] = {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 */
88 S[] = {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__
94 static const double zero = 0.0;
95 #else
96 static double zero = 0.0;
97 #endif
98 
99 #ifdef __STDC__
100  double __ieee754_j0(double x)
101 #else
102  double __ieee754_j0(x)
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__
160 static const double
161 #else
162 static double
163 #endif
164 U[] = {-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 */
171 V[] = {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
179  double __ieee754_y0(x)
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__
252 static const double pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
253 #else
254 static 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__
264 static const double pS8[5] = {
265 #else
266 static 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__
276 static const double pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
277 #else
278 static 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__
288 static const double pS5[5] = {
289 #else
290 static 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__
300 static const double pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
301 #else
302 static 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__
312 static const double pS3[5] = {
313 #else
314 static 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__
324 static const double pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
325 #else
326 static 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__
336 static const double pS2[5] = {
337 #else
338 static 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__
395 static const double qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
396 #else
397 static 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__
407 static const double qS8[6] = {
408 #else
409 static 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__
420 static const double qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
421 #else
422 static 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__
432 static const double qS5[6] = {
433 #else
434 static 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__
445 static const double qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
446 #else
447 static 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__
457 static const double qS3[6] = {
458 #else
459 static 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__
470 static const double qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
471 #else
472 static 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__
482 static const double qS2[6] = {
483 #else
484 static 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 }
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
GLdouble GLdouble u2
Definition: glext.h:8308
static double huge
Definition: j0_y0.c:79
struct S2 s2
static DNS_RECORDW r3
Definition: record.c:39
static double pR3[6]
Definition: j0_y0.c:302
GLdouble GLdouble GLdouble r
Definition: gl.h:2055
static __inline double __ieee754_sqrt(double x)
Definition: ieee754.h:40
static double qS5[6]
Definition: j0_y0.c:434
static double qR5[6]
Definition: j0_y0.c:422
static DNS_RECORDW r1
Definition: record.c:37
static double pR2[6]
Definition: j0_y0.c:326
GLint GLint GLint GLint GLint x
Definition: gl.h:1548
double __ieee754_j0(double x)
Definition: j0_y0.c:102
Definition: movable.cpp:7
GLdouble u1
Definition: glext.h:8308
static __inline double __cos(double x)
Definition: ieee754.h:42
static double pR5[6]
Definition: j0_y0.c:278
static double one
Definition: j0_y0.c:80
#define GET_HIGH_WORD(i, d)
Definition: ieee754.h:26
static double pzero()
static double pS3[5]
Definition: j0_y0.c:314
double __ieee754_y0(double x)
Definition: j0_y0.c:179
GLdouble GLdouble z
Definition: glext.h:5874
static double qS2[6]
Definition: j0_y0.c:484
static double qR8[6]
Definition: j0_y0.c:397
static double qzero()
#define EXTRACT_WORDS(ix0, ix1, d)
Definition: ieee754.h:16
struct S1 s1
static DNS_RECORDW r2
Definition: record.c:38
const GLubyte * c
Definition: glext.h:8905
static double qR3[6]
Definition: j0_y0.c:447
GLdouble GLdouble GLdouble GLdouble q
Definition: gl.h:2063
static double zero
Definition: j0_y0.c:96
static double pS5[5]
Definition: j0_y0.c:290
GLdouble s
Definition: gl.h:2039
static double qS3[6]
Definition: j0_y0.c:459
static double invsqrtpi
Definition: j0_y0.c:81
INT32 int32_t
Definition: types.h:71
_Check_return_ _CRT_JIT_INTRINSIC double __cdecl fabs(_In_ double x)
GLfloat GLfloat GLfloat v2
Definition: glext.h:6063
uint32_t cc
Definition: isohybrid.c:75
static __inline void __sincos(double x, double *s, double *c)
Definition: ieee754.h:43
static double qS8[6]
Definition: j0_y0.c:409
static double tpi
Definition: j0_y0.c:82
const GLdouble * v
Definition: gl.h:2040
static double R[]
Definition: j0_y0.c:84
static double qR2[6]
Definition: j0_y0.c:472
#define HUGE_VAL
Definition: math.h:51
#define c
Definition: ke_i.h:80
static double pR8[6]
Definition: j0_y0.c:254
static double pS2[5]
Definition: j0_y0.c:338
static double pS8[5]
Definition: j0_y0.c:266
static __inline double __ieee754_log(double x)
Definition: ieee754.h:41
static BYTE u3[]
Definition: msg.c:580
GLfloat GLfloat p
Definition: glext.h:8902
#define ss
Definition: i386-dis.c:432
static TRIO_CONST char rcsid[]
Definition: trio.c:756
GLfloat GLfloat v1
Definition: glext.h:6062
static double U[]
Definition: j0_y0.c:164