ReactOS 0.4.15-dev-8393-g61b7fb9
j1_y1.c
Go to the documentation of this file.
1/* @(#)e_j1.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_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__
69static double pone(double), qone(double);
70#else
71static double pone(), qone();
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] */
84R[] = {-6.25000000000000000000e-02, /* 0xBFB00000, 0x00000000 */
85 1.40705666955189706048e-03, /* 0x3F570D9F, 0x98472C61 */
86 -1.59955631084035597520e-05, /* 0xBEF0C5C6, 0xBA169668 */
87 4.96727999609584448412e-08}, /* 0x3E6AAAFA, 0x46CA0BD9 */
88S[] = {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__
95static const double zero = 0.0;
96#else
97static double zero = 0.0;
98#endif
99
100#ifdef __STDC__
101 double __ieee754_j1(double x)
102#else
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__
157static const double U0[5] = {
158#else
159static 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__
168static const double V0[5] = {
169#else
170static 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
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__
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 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__
264static const double ps8[5] = {
265#else
266static 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__
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.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__
288static const double ps5[5] = {
289#else
290static 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__
300static const double pr3[6] = {
301#else
302static 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__
312static const double ps3[5] = {
313#else
314static 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__
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 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__
336static const double ps2[5] = {
337#else
338static 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__
396static const double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
397#else
398static 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__
408static const double qs8[6] = {
409#else
410static 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__
421static const double qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
422#else
423static 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__
433static const double qs5[6] = {
434#else
435static 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__
446static const double qr3[6] = {
447#else
448static 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__
458static const double qs3[6] = {
459#else
460static 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__
471static const double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
472#else
473static 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__
483static const double qs2[6] = {
484#else
485static 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}
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
GLint GLint GLint GLint GLint GLint y
Definition: gl.h:1548
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 GLfloat GLfloat v3
Definition: glext.h:6064
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 pr8[6]
Definition: j1_y1.c:254
static double U0[5]
Definition: j1_y1.c:159
static double pr2[6]
Definition: j1_y1.c:326
static double huge
Definition: j1_y1.c:79
static double qs5[6]
Definition: j1_y1.c:435
static double pone()
static double R[]
Definition: j1_y1.c:84
static double qs2[6]
Definition: j1_y1.c:485
static double qr5[6]
Definition: j1_y1.c:423
static double zero
Definition: j1_y1.c:97
static double V0[5]
Definition: j1_y1.c:170
static double invsqrtpi
Definition: j1_y1.c:81
static double qr2[6]
Definition: j1_y1.c:473
static double qr8[6]
Definition: j1_y1.c:398
static double tpi
Definition: j1_y1.c:82
static double ps2[5]
Definition: j1_y1.c:338
static double qone()
static double ps5[5]
Definition: j1_y1.c:290
static double qr3[6]
Definition: j1_y1.c:448
static double one
Definition: j1_y1.c:80
static double qs8[6]
Definition: j1_y1.c:410
double __ieee754_y1(double x)
Definition: j1_y1.c:182
static double pr5[6]
Definition: j1_y1.c:278
double __ieee754_j1(double x)
Definition: j1_y1.c:103
static double ps8[5]
Definition: j1_y1.c:266
static double qs3[6]
Definition: j1_y1.c:460
static double ps3[5]
Definition: j1_y1.c:314
static double pr3[6]
Definition: j1_y1.c:302
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 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