ReactOS  0.4.13-dev-563-g0561610
j1_y1.c
Go to the documentation of this file.
1 /* @(#)e_j1.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_j1.c,v 1.8 1995/05/10 20:45:27 jtc Exp $";
18 #endif
19 
20 /* __ieee754_j1(x), __ieee754_y1(x)
21  * Bessel function of the first and second kinds of order zero.
22  * Method -- j1(x):
23  * 1. For tiny x, we use j1(x) = x/2 - x^3/16 + x^5/384 - ...
24  * 2. Reduce x to |x| since j1(x)=-j1(-x), and
25  * for x in (0,2)
26  * j1(x) = x/2 + x*z*R0/S0, where z = x*x;
27  * (precision: |j1/x - 1/2 - R0/S0 |<2**-61.51 )
28  * for x in (2,inf)
29  * j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1))
30  * y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
31  * where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
32  * as follow:
33  * cos(x1) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
34  * = 1/sqrt(2) * (sin(x) - cos(x))
35  * sin(x1) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
36  * = -1/sqrt(2) * (sin(x) + cos(x))
37  * (To avoid cancellation, use
38  * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
39  * to compute the worse one.)
40  *
41  * 3 Special cases
42  * j1(nan)= nan
43  * j1(0) = 0
44  * j1(inf) = 0
45  *
46  * Method -- y1(x):
47  * 1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN
48  * 2. For x<2.
49  * Since
50  * y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...)
51  * therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function.
52  * We use the following function to approximate y1,
53  * y1(x) = x*U(z)/V(z) + (2/pi)*(j1(x)*ln(x)-1/x), z= x^2
54  * where for x in [0,2] (abs err less than 2**-65.89)
55  * U(z) = U0[0] + U0[1]*z + ... + U0[4]*z^4
56  * V(z) = 1 + v0[0]*z + ... + v0[4]*z^5
57  * Note: For tiny x, 1/x dominate y1 and hence
58  * y1(tiny) = -2/pi/tiny, (choose tiny<2**-54)
59  * 3. For x>=2.
60  * y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
61  * where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
62  * by method mentioned above.
63  */
64 
65 #include "math.h"
66 #include "ieee754.h"
67 
68 #ifdef __STDC__
69 static double pone(double), qone(double);
70 #else
71 static double pone(), qone();
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] */
84 R[] = {-6.25000000000000000000e-02, /* 0xBFB00000, 0x00000000 */
85  1.40705666955189706048e-03, /* 0x3F570D9F, 0x98472C61 */
86  -1.59955631084035597520e-05, /* 0xBEF0C5C6, 0xBA169668 */
87  4.96727999609584448412e-08}, /* 0x3E6AAAFA, 0x46CA0BD9 */
88 S[] = {0.0, 1.91537599538363460805e-02, /* 0x3F939D0B, 0x12637E53 */
89  1.85946785588630915560e-04, /* 0x3F285F56, 0xB9CDF664 */
90  1.17718464042623683263e-06, /* 0x3EB3BFF8, 0x333F8498 */
91  5.04636257076217042715e-09, /* 0x3E35AC88, 0xC97DFF2C */
92  1.23542274426137913908e-11}; /* 0x3DAB2ACF, 0xCFB97ED8 */
93 
94 #ifdef __STDC__
95 static const double zero = 0.0;
96 #else
97 static double zero = 0.0;
98 #endif
99 
100 #ifdef __STDC__
101  double __ieee754_j1(double x)
102 #else
103  double __ieee754_j1(x)
104  double x;
105 #endif
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 }
155 
156 #ifdef __STDC__
157 static const double U0[5] = {
158 #else
159 static double U0[5] = {
160 #endif
161  -1.96057090646238940668e-01, /* 0xBFC91866, 0x143CBC8A */
162  5.04438716639811282616e-02, /* 0x3FA9D3C7, 0x76292CD1 */
163  -1.91256895875763547298e-03, /* 0xBF5F55E5, 0x4844F50F */
164  2.35252600561610495928e-05, /* 0x3EF8AB03, 0x8FA6B88E */
165  -9.19099158039878874504e-08, /* 0xBE78AC00, 0x569105B8 */
166 };
167 #ifdef __STDC__
168 static const double V0[5] = {
169 #else
170 static double V0[5] = {
171 #endif
172  1.99167318236649903973e-02, /* 0x3F94650D, 0x3F4DA9F0 */
173  2.02552581025135171496e-04, /* 0x3F2A8C89, 0x6C257764 */
174  1.35608801097516229404e-06, /* 0x3EB6C05A, 0x894E8CA6 */
175  6.22741452364621501295e-09, /* 0x3E3ABF1D, 0x5BA69A86 */
176  1.66559246207992079114e-11, /* 0x3DB25039, 0xDACA772A */
177 };
178 
179 #ifdef __STDC__
180  double __ieee754_y1(double x)
181 #else
182  double __ieee754_y1(x)
183  double x;
184 #endif
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 }
240 
241 /* For x >= 8, the asymptotic expansions of pone is
242  * 1 + 15/128 s^2 - 4725/2^15 s^4 - ..., where s = 1/x.
243  * We approximate pone by
244  * pone(x) = 1 + (R/S)
245  * where R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10
246  * S = 1 + ps0*s^2 + ... + ps4*s^10
247  * and
248  * | pone(x)-1-R/S | <= 2 ** ( -60.06)
249  */
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  1.17187499999988647970e-01, /* 0x3FBDFFFF, 0xFFFFFCCE */
258  1.32394806593073575129e+01, /* 0x402A7A9D, 0x357F7FCE */
259  4.12051854307378562225e+02, /* 0x4079C0D4, 0x652EA590 */
260  3.87474538913960532227e+03, /* 0x40AE457D, 0xA3A532CC */
261  7.91447954031891731574e+03, /* 0x40BEEA7A, 0xC32782DD */
262 };
263 #ifdef __STDC__
264 static const double ps8[5] = {
265 #else
266 static double ps8[5] = {
267 #endif
268  1.14207370375678408436e+02, /* 0x405C8D45, 0x8E656CAC */
269  3.65093083420853463394e+03, /* 0x40AC85DC, 0x964D274F */
270  3.69562060269033463555e+04, /* 0x40E20B86, 0x97C5BB7F */
271  9.76027935934950801311e+04, /* 0x40F7D42C, 0xB28F17BB */
272  3.08042720627888811578e+04, /* 0x40DE1511, 0x697A0B2D */
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.31990519556243522749e-11, /* 0x3DAD0667, 0xDAE1CA7D */
281  1.17187493190614097638e-01, /* 0x3FBDFFFF, 0xE2C10043 */
282  6.80275127868432871736e+00, /* 0x401B3604, 0x6E6315E3 */
283  1.08308182990189109773e+02, /* 0x405B13B9, 0x452602ED */
284  5.17636139533199752805e+02, /* 0x40802D16, 0xD052D649 */
285  5.28715201363337541807e+02, /* 0x408085B8, 0xBB7E0CB7 */
286 };
287 #ifdef __STDC__
288 static const double ps5[5] = {
289 #else
290 static double ps5[5] = {
291 #endif
292  5.92805987221131331921e+01, /* 0x404DA3EA, 0xA8AF633D */
293  9.91401418733614377743e+02, /* 0x408EFB36, 0x1B066701 */
294  5.35326695291487976647e+03, /* 0x40B4E944, 0x5706B6FB */
295  7.84469031749551231769e+03, /* 0x40BEA4B0, 0xB8A5BB15 */
296  1.50404688810361062679e+03, /* 0x40978030, 0x036F5E51 */
297 };
298 
299 #ifdef __STDC__
300 static const double pr3[6] = {
301 #else
302 static double pr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
303 #endif
304  3.02503916137373618024e-09, /* 0x3E29FC21, 0xA7AD9EDD */
305  1.17186865567253592491e-01, /* 0x3FBDFFF5, 0x5B21D17B */
306  3.93297750033315640650e+00, /* 0x400F76BC, 0xE85EAD8A */
307  3.51194035591636932736e+01, /* 0x40418F48, 0x9DA6D129 */
308  9.10550110750781271918e+01, /* 0x4056C385, 0x4D2C1837 */
309  4.85590685197364919645e+01, /* 0x4048478F, 0x8EA83EE5 */
310 };
311 #ifdef __STDC__
312 static const double ps3[5] = {
313 #else
314 static double ps3[5] = {
315 #endif
316  3.47913095001251519989e+01, /* 0x40416549, 0xA134069C */
317  3.36762458747825746741e+02, /* 0x40750C33, 0x07F1A75F */
318  1.04687139975775130551e+03, /* 0x40905B7C, 0x5037D523 */
319  8.90811346398256432622e+02, /* 0x408BD67D, 0xA32E31E9 */
320  1.03787932439639277504e+02, /* 0x4059F26D, 0x7C2EED53 */
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  1.07710830106873743082e-07, /* 0x3E7CE9D4, 0xF65544F4 */
329  1.17176219462683348094e-01, /* 0x3FBDFF42, 0xBE760D83 */
330  2.36851496667608785174e+00, /* 0x4002F2B7, 0xF98FAEC0 */
331  1.22426109148261232917e+01, /* 0x40287C37, 0x7F71A964 */
332  1.76939711271687727390e+01, /* 0x4031B1A8, 0x177F8EE2 */
333  5.07352312588818499250e+00, /* 0x40144B49, 0xA574C1FE */
334 };
335 #ifdef __STDC__
336 static const double ps2[5] = {
337 #else
338 static double ps2[5] = {
339 #endif
340  2.14364859363821409488e+01, /* 0x40356FBD, 0x8AD5ECDC */
341  1.25290227168402751090e+02, /* 0x405F5293, 0x14F92CD5 */
342  2.32276469057162813669e+02, /* 0x406D08D8, 0xD5A2DBD9 */
343  1.17679373287147100768e+02, /* 0x405D6B7A, 0xDA1884A9 */
344  8.36463893371618283368e+00, /* 0x4020BAB1, 0xF44E5192 */
345 };
346 
347 #ifdef __STDC__
348  static double pone(double x)
349 #else
350  static double pone(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,r1,r2,r3,s1,s2,s3,z2,z4;
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 qone is
386  * 3/8 s - 105/1024 s^3 - ..., where s = 1/x.
387  * We approximate pone by
388  * qone(x) = s*(0.375 + (R/S))
389  * where R = qr1*s^2 + qr2*s^4 + ... + qr5*s^10
390  * S = 1 + qs1*s^2 + ... + qs6*s^12
391  * and
392  * | qone(x)/s -0.375-R/S | <= 2 ** ( -61.13)
393  */
394 
395 #ifdef __STDC__
396 static const double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
397 #else
398 static double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
399 #endif
400  0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
401  -1.02539062499992714161e-01, /* 0xBFBA3FFF, 0xFFFFFDF3 */
402  -1.62717534544589987888e+01, /* 0xC0304591, 0xA26779F7 */
403  -7.59601722513950107896e+02, /* 0xC087BCD0, 0x53E4B576 */
404  -1.18498066702429587167e+04, /* 0xC0C724E7, 0x40F87415 */
405  -4.84385124285750353010e+04, /* 0xC0E7A6D0, 0x65D09C6A */
406 };
407 #ifdef __STDC__
408 static const double qs8[6] = {
409 #else
410 static double qs8[6] = {
411 #endif
412  1.61395369700722909556e+02, /* 0x40642CA6, 0xDE5BCDE5 */
413  7.82538599923348465381e+03, /* 0x40BE9162, 0xD0D88419 */
414  1.33875336287249578163e+05, /* 0x4100579A, 0xB0B75E98 */
415  7.19657723683240939863e+05, /* 0x4125F653, 0x72869C19 */
416  6.66601232617776375264e+05, /* 0x412457D2, 0x7719AD5C */
417  -2.94490264303834643215e+05, /* 0xC111F969, 0x0EA5AA18 */
418 };
419 
420 #ifdef __STDC__
421 static const double qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
422 #else
423 static double qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
424 #endif
425  -2.08979931141764104297e-11, /* 0xBDB6FA43, 0x1AA1A098 */
426  -1.02539050241375426231e-01, /* 0xBFBA3FFF, 0xCB597FEF */
427  -8.05644828123936029840e+00, /* 0xC0201CE6, 0xCA03AD4B */
428  -1.83669607474888380239e+02, /* 0xC066F56D, 0x6CA7B9B0 */
429  -1.37319376065508163265e+03, /* 0xC09574C6, 0x6931734F */
430  -2.61244440453215656817e+03, /* 0xC0A468E3, 0x88FDA79D */
431 };
432 #ifdef __STDC__
433 static const double qs5[6] = {
434 #else
435 static double qs5[6] = {
436 #endif
437  8.12765501384335777857e+01, /* 0x405451B2, 0xFF5A11B2 */
438  1.99179873460485964642e+03, /* 0x409F1F31, 0xE77BF839 */
439  1.74684851924908907677e+04, /* 0x40D10F1F, 0x0D64CE29 */
440  4.98514270910352279316e+04, /* 0x40E8576D, 0xAABAD197 */
441  2.79480751638918118260e+04, /* 0x40DB4B04, 0xCF7C364B */
442  -4.71918354795128470869e+03, /* 0xC0B26F2E, 0xFCFFA004 */
443 };
444 
445 #ifdef __STDC__
446 static const double qr3[6] = {
447 #else
448 static double qr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
449 #endif
450  -5.07831226461766561369e-09, /* 0xBE35CFA9, 0xD38FC84F */
451  -1.02537829820837089745e-01, /* 0xBFBA3FEB, 0x51AEED54 */
452  -4.61011581139473403113e+00, /* 0xC01270C2, 0x3302D9FF */
453  -5.78472216562783643212e+01, /* 0xC04CEC71, 0xC25D16DA */
454  -2.28244540737631695038e+02, /* 0xC06C87D3, 0x4718D55F */
455  -2.19210128478909325622e+02, /* 0xC06B66B9, 0x5F5C1BF6 */
456 };
457 #ifdef __STDC__
458 static const double qs3[6] = {
459 #else
460 static double qs3[6] = {
461 #endif
462  4.76651550323729509273e+01, /* 0x4047D523, 0xCCD367E4 */
463  6.73865112676699709482e+02, /* 0x40850EEB, 0xC031EE3E */
464  3.38015286679526343505e+03, /* 0x40AA684E, 0x448E7C9A */
465  5.54772909720722782367e+03, /* 0x40B5ABBA, 0xA61D54A6 */
466  1.90311919338810798763e+03, /* 0x409DBC7A, 0x0DD4DF4B */
467  -1.35201191444307340817e+02, /* 0xC060E670, 0x290A311F */
468 };
469 
470 #ifdef __STDC__
471 static const double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
472 #else
473 static double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
474 #endif
475  -1.78381727510958865572e-07, /* 0xBE87F126, 0x44C626D2 */
476  -1.02517042607985553460e-01, /* 0xBFBA3E8E, 0x9148B010 */
477  -2.75220568278187460720e+00, /* 0xC0060484, 0x69BB4EDA */
478  -1.96636162643703720221e+01, /* 0xC033A9E2, 0xC168907F */
479  -4.23253133372830490089e+01, /* 0xC04529A3, 0xDE104AAA */
480  -2.13719211703704061733e+01, /* 0xC0355F36, 0x39CF6E52 */
481 };
482 #ifdef __STDC__
483 static const double qs2[6] = {
484 #else
485 static double qs2[6] = {
486 #endif
487  2.95333629060523854548e+01, /* 0x403D888A, 0x78AE64FF */
488  2.52981549982190529136e+02, /* 0x406F9F68, 0xDB821CBA */
489  7.57502834868645436472e+02, /* 0x4087AC05, 0xCE49A0F7 */
490  7.39393205320467245656e+02, /* 0x40871B25, 0x48D4C029 */
491  1.55949003336666123687e+02, /* 0x40637E5E, 0x3C3ED8D4 */
492  -4.95949898822628210127e+00, /* 0xC013D686, 0xE71BE86B */
493 };
494 
495 #ifdef __STDC__
496  static double qone(double x)
497 #else
498  static double qone(x)
499  double x;
500 #endif
501 {
502 #ifdef __STDC__
503  const double *p = 0,*q = 0;
504 #else
505  double *p = 0,*q = 0;
506 #endif
507  double s,r,z,r1,r2,r3,s1,s2,s3,z2,z4,z6;
508  int32_t ix;
509  GET_HIGH_WORD(ix,x);
510  ix &= 0x7fffffff;
511  if(ix>=0x40200000) {p = qr8; q= qs8;}
512  else if(ix>=0x40122E8B){p = qr5; q= qs5;}
513  else if(ix>=0x4006DB6D){p = qr3; q= qs3;}
514  else if(ix>=0x40000000){p = qr2; q= qs2;}
515  z = one/(x*x);
516 #ifdef DO_NOT_USE_THIS
517  r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
518  s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
519 #else
520  r1 = p[0]+z*p[1]; z2=z*z;
521  r2 = p[2]+z*p[3]; z4=z2*z2;
522  r3 = p[4]+z*p[5]; z6=z4*z2;
523  r = r1 + z2*r2 + z4*r3;
524  s1 = one+z*q[0];
525  s2 = q[1]+z*q[2];
526  s3 = q[3]+z*q[4];
527  s = s1 + z2*s2 + z4*s3 + z6*q[5];
528 #endif
529  return (.375 + r/s)/x;
530 }
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
static double invsqrtpi
Definition: j1_y1.c:81
GLdouble GLdouble u2
Definition: glext.h:8308
struct S2 s2
static DNS_RECORDW r3
Definition: record.c:39
static double pr2[6]
Definition: j1_y1.c:326
static double zero
Definition: j1_y1.c:97
GLdouble GLdouble GLdouble r
Definition: gl.h:2055
static __inline double __ieee754_sqrt(double x)
Definition: ieee754.h:40
static DNS_RECORDW r1
Definition: record.c:37
GLint GLint GLint GLint GLint x
Definition: gl.h:1548
Definition: movable.cpp:7
GLdouble u1
Definition: glext.h:8308
static __inline double __cos(double x)
Definition: ieee754.h:42
static double tpi
Definition: j1_y1.c:82
static double ps3[5]
Definition: j1_y1.c:314
static double U0[5]
Definition: j1_y1.c:159
double __ieee754_y1(double x)
Definition: j1_y1.c:182
static double pr3[6]
Definition: j1_y1.c:302
#define GET_HIGH_WORD(i, d)
Definition: ieee754.h:26
static double qr5[6]
Definition: j1_y1.c:423
static double qr3[6]
Definition: j1_y1.c:448
GLdouble GLdouble z
Definition: glext.h:5874
double __ieee754_j1(double x)
Definition: j1_y1.c:103
static double ps5[5]
Definition: j1_y1.c:290
static double qs3[6]
Definition: j1_y1.c:460
static double qr2[6]
Definition: j1_y1.c:473
#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
GLdouble GLdouble GLdouble GLdouble q
Definition: gl.h:2063
static double qone()
static double V0[5]
Definition: j1_y1.c:170
GLdouble s
Definition: gl.h:2039
static double pone()
GLfloat GLfloat GLfloat GLfloat v3
Definition: glext.h:6064
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 double qs2[6]
Definition: j1_y1.c:485
static __inline void __sincos(double x, double *s, double *c)
Definition: ieee754.h:43
static double ps8[5]
Definition: j1_y1.c:266
const GLdouble * v
Definition: gl.h:2040
static double qr8[6]
Definition: j1_y1.c:398
static double ps2[5]
Definition: j1_y1.c:338
static double one
Definition: j1_y1.c:80
GLint GLint GLint GLint GLint GLint y
Definition: gl.h:1548
static double qs5[6]
Definition: j1_y1.c:435
#define HUGE_VAL
Definition: math.h:51
static double qs8[6]
Definition: j1_y1.c:410
#define c
Definition: ke_i.h:80
static __inline double __ieee754_log(double x)
Definition: ieee754.h:41
static double R[]
Definition: j1_y1.c:84
GLfloat GLfloat p
Definition: glext.h:8902
static double huge
Definition: j1_y1.c:79
#define ss
Definition: i386-dis.c:432
static double pr5[6]
Definition: j1_y1.c:278
static double pr8[6]
Definition: j1_y1.c:254
static TRIO_CONST char rcsid[]
Definition: trio.c:756
GLfloat GLfloat v1
Definition: glext.h:6062