1 | /* |
2 | * ==================================================== |
3 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
4 | * |
5 | * Developed at SunPro, a Sun Microsystems, Inc. business. |
6 | * Permission to use, copy, modify, and distribute this |
7 | * software is freely granted, provided that this notice |
8 | * is preserved. |
9 | * ==================================================== |
10 | */ |
11 | |
12 | /* Long double expansions are |
13 | Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov> |
14 | and are incorporated herein by permission of the author. The author |
15 | reserves the right to distribute this material elsewhere under different |
16 | copying permissions. These modifications are distributed here under |
17 | the following terms: |
18 | |
19 | This library is free software; you can redistribute it and/or |
20 | modify it under the terms of the GNU Lesser General Public |
21 | License as published by the Free Software Foundation; either |
22 | version 2.1 of the License, or (at your option) any later version. |
23 | |
24 | This library is distributed in the hope that it will be useful, |
25 | but WITHOUT ANY WARRANTY; without even the implied warranty of |
26 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
27 | Lesser General Public License for more details. |
28 | |
29 | You should have received a copy of the GNU Lesser General Public |
30 | License along with this library; if not, see |
31 | <http://www.gnu.org/licenses/>. */ |
32 | |
33 | /* __ieee754_j1(x), __ieee754_y1(x) |
34 | * Bessel function of the first and second kinds of order zero. |
35 | * Method -- j1(x): |
36 | * 1. For tiny x, we use j1(x) = x/2 - x^3/16 + x^5/384 - ... |
37 | * 2. Reduce x to |x| since j1(x)=-j1(-x), and |
38 | * for x in (0,2) |
39 | * j1(x) = x/2 + x*z*R0/S0, where z = x*x; |
40 | * for x in (2,inf) |
41 | * j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1)) |
42 | * y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1)) |
43 | * where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1) |
44 | * as follow: |
45 | * cos(x1) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4) |
46 | * = 1/sqrt(2) * (sin(x) - cos(x)) |
47 | * sin(x1) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) |
48 | * = -1/sqrt(2) * (sin(x) + cos(x)) |
49 | * (To avoid cancellation, use |
50 | * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) |
51 | * to compute the worse one.) |
52 | * |
53 | * 3 Special cases |
54 | * j1(nan)= nan |
55 | * j1(0) = 0 |
56 | * j1(inf) = 0 |
57 | * |
58 | * Method -- y1(x): |
59 | * 1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN |
60 | * 2. For x<2. |
61 | * Since |
62 | * y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...) |
63 | * therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function. |
64 | * We use the following function to approximate y1, |
65 | * y1(x) = x*U(z)/V(z) + (2/pi)*(j1(x)*ln(x)-1/x), z= x^2 |
66 | * Note: For tiny x, 1/x dominate y1 and hence |
67 | * y1(tiny) = -2/pi/tiny |
68 | * 3. For x>=2. |
69 | * y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1)) |
70 | * where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1) |
71 | * by method mentioned above. |
72 | */ |
73 | |
74 | #include <errno.h> |
75 | #include <float.h> |
76 | #include <math.h> |
77 | #include <math_private.h> |
78 | |
79 | static long double pone (long double), qone (long double); |
80 | |
81 | static const long double |
82 | huge = 1e4930L, |
83 | one = 1.0L, |
84 | invsqrtpi = 5.6418958354775628694807945156077258584405e-1L, |
85 | tpi = 6.3661977236758134307553505349005744813784e-1L, |
86 | |
87 | /* J1(x) = .5 x + x x^2 R(x^2) / S(x^2) |
88 | 0 <= x <= 2 |
89 | Peak relative error 4.5e-21 */ |
90 | R[5] = { |
91 | -9.647406112428107954753770469290757756814E7L, |
92 | 2.686288565865230690166454005558203955564E6L, |
93 | -3.689682683905671185891885948692283776081E4L, |
94 | 2.195031194229176602851429567792676658146E2L, |
95 | -5.124499848728030297902028238597308971319E-1L, |
96 | }, |
97 | |
98 | S[4] = |
99 | { |
100 | 1.543584977988497274437410333029029035089E9L, |
101 | 2.133542369567701244002565983150952549520E7L, |
102 | 1.394077011298227346483732156167414670520E5L, |
103 | 5.252401789085732428842871556112108446506E2L, |
104 | /* 1.000000000000000000000000000000000000000E0L, */ |
105 | }; |
106 | |
107 | static const long double zero = 0.0; |
108 | |
109 | |
110 | long double |
111 | __ieee754_j1l (long double x) |
112 | { |
113 | long double z, c, r, s, ss, cc, u, v, y; |
114 | int32_t ix; |
115 | uint32_t se; |
116 | |
117 | GET_LDOUBLE_EXP (se, x); |
118 | ix = se & 0x7fff; |
119 | if (__glibc_unlikely (ix >= 0x7fff)) |
120 | return one / x; |
121 | y = fabsl (x); |
122 | if (ix >= 0x4000) |
123 | { /* |x| >= 2.0 */ |
124 | __sincosl (y, &s, &c); |
125 | ss = -s - c; |
126 | cc = s - c; |
127 | if (ix < 0x7ffe) |
128 | { /* make sure y+y not overflow */ |
129 | z = __cosl (y + y); |
130 | if ((s * c) > zero) |
131 | cc = z / ss; |
132 | else |
133 | ss = z / cc; |
134 | } |
135 | /* |
136 | * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x) |
137 | * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x) |
138 | */ |
139 | if (__glibc_unlikely (ix > 0x4080)) |
140 | z = (invsqrtpi * cc) / __ieee754_sqrtl (y); |
141 | else |
142 | { |
143 | u = pone (y); |
144 | v = qone (y); |
145 | z = invsqrtpi * (u * cc - v * ss) / __ieee754_sqrtl (y); |
146 | } |
147 | if (se & 0x8000) |
148 | return -z; |
149 | else |
150 | return z; |
151 | } |
152 | if (__glibc_unlikely (ix < 0x3fde)) /* |x| < 2^-33 */ |
153 | { |
154 | if (huge + x > one) /* inexact if x!=0 necessary */ |
155 | { |
156 | long double ret = 0.5 * x; |
157 | math_check_force_underflow (ret); |
158 | if (ret == 0 && x != 0) |
159 | __set_errno (ERANGE); |
160 | return ret; |
161 | } |
162 | } |
163 | z = x * x; |
164 | r = z * (R[0] + z * (R[1]+ z * (R[2] + z * (R[3] + z * R[4])))); |
165 | s = S[0] + z * (S[1] + z * (S[2] + z * (S[3] + z))); |
166 | r *= x; |
167 | return (x * 0.5 + r / s); |
168 | } |
169 | strong_alias (__ieee754_j1l, __j1l_finite) |
170 | |
171 | |
172 | /* Y1(x) = 2/pi * (log(x) * j1(x) - 1/x) + x R(x^2) |
173 | 0 <= x <= 2 |
174 | Peak relative error 2.3e-23 */ |
175 | static const long double U0[6] = { |
176 | -5.908077186259914699178903164682444848615E10L, |
177 | 1.546219327181478013495975514375773435962E10L, |
178 | -6.438303331169223128870035584107053228235E8L, |
179 | 9.708540045657182600665968063824819371216E6L, |
180 | -6.138043997084355564619377183564196265471E4L, |
181 | 1.418503228220927321096904291501161800215E2L, |
182 | }; |
183 | static const long double V0[5] = { |
184 | 3.013447341682896694781964795373783679861E11L, |
185 | 4.669546565705981649470005402243136124523E9L, |
186 | 3.595056091631351184676890179233695857260E7L, |
187 | 1.761554028569108722903944659933744317994E5L, |
188 | 5.668480419646516568875555062047234534863E2L, |
189 | /* 1.000000000000000000000000000000000000000E0L, */ |
190 | }; |
191 | |
192 | |
193 | long double |
194 | __ieee754_y1l (long double x) |
195 | { |
196 | long double z, s, c, ss, cc, u, v; |
197 | int32_t ix; |
198 | uint32_t se, i0, i1; |
199 | |
200 | GET_LDOUBLE_WORDS (se, i0, i1, x); |
201 | ix = se & 0x7fff; |
202 | /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */ |
203 | if (__glibc_unlikely (se & 0x8000)) |
204 | return zero / (zero * x); |
205 | if (__glibc_unlikely (ix >= 0x7fff)) |
206 | return one / (x + x * x); |
207 | if (__glibc_unlikely ((i0 | i1) == 0)) |
208 | return -HUGE_VALL + x; /* -inf and overflow exception. */ |
209 | if (ix >= 0x4000) |
210 | { /* |x| >= 2.0 */ |
211 | __sincosl (x, &s, &c); |
212 | ss = -s - c; |
213 | cc = s - c; |
214 | if (ix < 0x7ffe) |
215 | { /* make sure x+x not overflow */ |
216 | z = __cosl (x + x); |
217 | if ((s * c) > zero) |
218 | cc = z / ss; |
219 | else |
220 | ss = z / cc; |
221 | } |
222 | /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0)) |
223 | * where x0 = x-3pi/4 |
224 | * Better formula: |
225 | * cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4) |
226 | * = 1/sqrt(2) * (sin(x) - cos(x)) |
227 | * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) |
228 | * = -1/sqrt(2) * (cos(x) + sin(x)) |
229 | * To avoid cancellation, use |
230 | * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) |
231 | * to compute the worse one. |
232 | */ |
233 | if (__glibc_unlikely (ix > 0x4080)) |
234 | z = (invsqrtpi * ss) / __ieee754_sqrtl (x); |
235 | else |
236 | { |
237 | u = pone (x); |
238 | v = qone (x); |
239 | z = invsqrtpi * (u * ss + v * cc) / __ieee754_sqrtl (x); |
240 | } |
241 | return z; |
242 | } |
243 | if (__glibc_unlikely (ix <= 0x3fbe)) |
244 | { /* x < 2**-65 */ |
245 | z = -tpi / x; |
246 | if (isinf (z)) |
247 | __set_errno (ERANGE); |
248 | return z; |
249 | } |
250 | z = x * x; |
251 | u = U0[0] + z * (U0[1] + z * (U0[2] + z * (U0[3] + z * (U0[4] + z * U0[5])))); |
252 | v = V0[0] + z * (V0[1] + z * (V0[2] + z * (V0[3] + z * (V0[4] + z)))); |
253 | return (x * (u / v) + |
254 | tpi * (__ieee754_j1l (x) * __ieee754_logl (x) - one / x)); |
255 | } |
256 | strong_alias (__ieee754_y1l, __y1l_finite) |
257 | |
258 | |
259 | /* For x >= 8, the asymptotic expansions of pone is |
260 | * 1 + 15/128 s^2 - 4725/2^15 s^4 - ..., where s = 1/x. |
261 | * We approximate pone by |
262 | * pone(x) = 1 + (R/S) |
263 | */ |
264 | |
265 | /* J1(x) cosX + Y1(x) sinX = sqrt( 2/(pi x)) P1(x) |
266 | P1(x) = 1 + z^2 R(z^2), z=1/x |
267 | 8 <= x <= inf (0 <= z <= 0.125) |
268 | Peak relative error 5.2e-22 */ |
269 | |
270 | static const long double pr8[7] = { |
271 | 8.402048819032978959298664869941375143163E-9L, |
272 | 1.813743245316438056192649247507255996036E-6L, |
273 | 1.260704554112906152344932388588243836276E-4L, |
274 | 3.439294839869103014614229832700986965110E-3L, |
275 | 3.576910849712074184504430254290179501209E-2L, |
276 | 1.131111483254318243139953003461511308672E-1L, |
277 | 4.480715825681029711521286449131671880953E-2L, |
278 | }; |
279 | static const long double ps8[6] = { |
280 | 7.169748325574809484893888315707824924354E-8L, |
281 | 1.556549720596672576431813934184403614817E-5L, |
282 | 1.094540125521337139209062035774174565882E-3L, |
283 | 3.060978962596642798560894375281428805840E-2L, |
284 | 3.374146536087205506032643098619414507024E-1L, |
285 | 1.253830208588979001991901126393231302559E0L, |
286 | /* 1.000000000000000000000000000000000000000E0L, */ |
287 | }; |
288 | |
289 | /* J1(x) cosX + Y1(x) sinX = sqrt( 2/(pi x)) P1(x) |
290 | P1(x) = 1 + z^2 R(z^2), z=1/x |
291 | 4.54541015625 <= x <= 8 |
292 | Peak relative error 7.7e-22 */ |
293 | static const long double pr5[7] = { |
294 | 4.318486887948814529950980396300969247900E-7L, |
295 | 4.715341880798817230333360497524173929315E-5L, |
296 | 1.642719430496086618401091544113220340094E-3L, |
297 | 2.228688005300803935928733750456396149104E-2L, |
298 | 1.142773760804150921573259605730018327162E-1L, |
299 | 1.755576530055079253910829652698703791957E-1L, |
300 | 3.218803858282095929559165965353784980613E-2L, |
301 | }; |
302 | static const long double ps5[6] = { |
303 | 3.685108812227721334719884358034713967557E-6L, |
304 | 4.069102509511177498808856515005792027639E-4L, |
305 | 1.449728676496155025507893322405597039816E-2L, |
306 | 2.058869213229520086582695850441194363103E-1L, |
307 | 1.164890985918737148968424972072751066553E0L, |
308 | 2.274776933457009446573027260373361586841E0L, |
309 | /* 1.000000000000000000000000000000000000000E0L,*/ |
310 | }; |
311 | |
312 | /* J1(x) cosX + Y1(x) sinX = sqrt( 2/(pi x)) P1(x) |
313 | P1(x) = 1 + z^2 R(z^2), z=1/x |
314 | 2.85711669921875 <= x <= 4.54541015625 |
315 | Peak relative error 6.5e-21 */ |
316 | static const long double pr3[7] = { |
317 | 1.265251153957366716825382654273326407972E-5L, |
318 | 8.031057269201324914127680782288352574567E-4L, |
319 | 1.581648121115028333661412169396282881035E-2L, |
320 | 1.179534658087796321928362981518645033967E-1L, |
321 | 3.227936912780465219246440724502790727866E-1L, |
322 | 2.559223765418386621748404398017602935764E-1L, |
323 | 2.277136933287817911091370397134882441046E-2L, |
324 | }; |
325 | static const long double ps3[6] = { |
326 | 1.079681071833391818661952793568345057548E-4L, |
327 | 6.986017817100477138417481463810841529026E-3L, |
328 | 1.429403701146942509913198539100230540503E-1L, |
329 | 1.148392024337075609460312658938700765074E0L, |
330 | 3.643663015091248720208251490291968840882E0L, |
331 | 3.990702269032018282145100741746633960737E0L, |
332 | /* 1.000000000000000000000000000000000000000E0L, */ |
333 | }; |
334 | |
335 | /* J1(x) cosX + Y1(x) sinX = sqrt( 2/(pi x)) P1(x) |
336 | P1(x) = 1 + z^2 R(z^2), z=1/x |
337 | 2 <= x <= 2.85711669921875 |
338 | Peak relative error 3.5e-21 */ |
339 | static const long double pr2[7] = { |
340 | 2.795623248568412225239401141338714516445E-4L, |
341 | 1.092578168441856711925254839815430061135E-2L, |
342 | 1.278024620468953761154963591853679640560E-1L, |
343 | 5.469680473691500673112904286228351988583E-1L, |
344 | 8.313769490922351300461498619045639016059E-1L, |
345 | 3.544176317308370086415403567097130611468E-1L, |
346 | 1.604142674802373041247957048801599740644E-2L, |
347 | }; |
348 | static const long double ps2[6] = { |
349 | 2.385605161555183386205027000675875235980E-3L, |
350 | 9.616778294482695283928617708206967248579E-2L, |
351 | 1.195215570959693572089824415393951258510E0L, |
352 | 5.718412857897054829999458736064922974662E0L, |
353 | 1.065626298505499086386584642761602177568E1L, |
354 | 6.809140730053382188468983548092322151791E0L, |
355 | /* 1.000000000000000000000000000000000000000E0L, */ |
356 | }; |
357 | |
358 | |
359 | static long double |
360 | pone (long double x) |
361 | { |
362 | const long double *p, *q; |
363 | long double z, r, s; |
364 | int32_t ix; |
365 | uint32_t se, i0, i1; |
366 | |
367 | GET_LDOUBLE_WORDS (se, i0, i1, x); |
368 | ix = se & 0x7fff; |
369 | /* ix >= 0x4000 for all calls to this function. */ |
370 | if (ix >= 0x4002) /* x >= 8 */ |
371 | { |
372 | p = pr8; |
373 | q = ps8; |
374 | } |
375 | else |
376 | { |
377 | i1 = (ix << 16) | (i0 >> 16); |
378 | if (i1 >= 0x40019174) /* x >= 4.54541015625 */ |
379 | { |
380 | p = pr5; |
381 | q = ps5; |
382 | } |
383 | else if (i1 >= 0x4000b6db) /* x >= 2.85711669921875 */ |
384 | { |
385 | p = pr3; |
386 | q = ps3; |
387 | } |
388 | else /* x >= 2 */ |
389 | { |
390 | p = pr2; |
391 | q = ps2; |
392 | } |
393 | } |
394 | z = one / (x * x); |
395 | r = p[0] + z * (p[1] + |
396 | z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6]))))); |
397 | s = q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * (q[5] + z))))); |
398 | return one + z * r / s; |
399 | } |
400 | |
401 | |
402 | /* For x >= 8, the asymptotic expansions of qone is |
403 | * 3/8 s - 105/1024 s^3 - ..., where s = 1/x. |
404 | * We approximate pone by |
405 | * qone(x) = s*(0.375 + (R/S)) |
406 | */ |
407 | |
408 | /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x), |
409 | Q1(x) = z(.375 + z^2 R(z^2)), z=1/x |
410 | 8 <= x <= inf |
411 | Peak relative error 8.3e-22 */ |
412 | |
413 | static const long double qr8[7] = { |
414 | -5.691925079044209246015366919809404457380E-10L, |
415 | -1.632587664706999307871963065396218379137E-7L, |
416 | -1.577424682764651970003637263552027114600E-5L, |
417 | -6.377627959241053914770158336842725291713E-4L, |
418 | -1.087408516779972735197277149494929568768E-2L, |
419 | -6.854943629378084419631926076882330494217E-2L, |
420 | -1.055448290469180032312893377152490183203E-1L, |
421 | }; |
422 | static const long double qs8[7] = { |
423 | 5.550982172325019811119223916998393907513E-9L, |
424 | 1.607188366646736068460131091130644192244E-6L, |
425 | 1.580792530091386496626494138334505893599E-4L, |
426 | 6.617859900815747303032860443855006056595E-3L, |
427 | 1.212840547336984859952597488863037659161E-1L, |
428 | 9.017885953937234900458186716154005541075E-1L, |
429 | 2.201114489712243262000939120146436167178E0L, |
430 | /* 1.000000000000000000000000000000000000000E0L, */ |
431 | }; |
432 | |
433 | /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x), |
434 | Q1(x) = z(.375 + z^2 R(z^2)), z=1/x |
435 | 4.54541015625 <= x <= 8 |
436 | Peak relative error 4.1e-22 */ |
437 | static const long double qr5[7] = { |
438 | -6.719134139179190546324213696633564965983E-8L, |
439 | -9.467871458774950479909851595678622044140E-6L, |
440 | -4.429341875348286176950914275723051452838E-4L, |
441 | -8.539898021757342531563866270278505014487E-3L, |
442 | -6.818691805848737010422337101409276287170E-2L, |
443 | -1.964432669771684034858848142418228214855E-1L, |
444 | -1.333896496989238600119596538299938520726E-1L, |
445 | }; |
446 | static const long double qs5[7] = { |
447 | 6.552755584474634766937589285426911075101E-7L, |
448 | 9.410814032118155978663509073200494000589E-5L, |
449 | 4.561677087286518359461609153655021253238E-3L, |
450 | 9.397742096177905170800336715661091535805E-2L, |
451 | 8.518538116671013902180962914473967738771E-1L, |
452 | 3.177729183645800174212539541058292579009E0L, |
453 | 4.006745668510308096259753538973038902990E0L, |
454 | /* 1.000000000000000000000000000000000000000E0L, */ |
455 | }; |
456 | |
457 | /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x), |
458 | Q1(x) = z(.375 + z^2 R(z^2)), z=1/x |
459 | 2.85711669921875 <= x <= 4.54541015625 |
460 | Peak relative error 2.2e-21 */ |
461 | static const long double qr3[7] = { |
462 | -3.618746299358445926506719188614570588404E-6L, |
463 | -2.951146018465419674063882650970344502798E-4L, |
464 | -7.728518171262562194043409753656506795258E-3L, |
465 | -8.058010968753999435006488158237984014883E-2L, |
466 | -3.356232856677966691703904770937143483472E-1L, |
467 | -4.858192581793118040782557808823460276452E-1L, |
468 | -1.592399251246473643510898335746432479373E-1L, |
469 | }; |
470 | static const long double qs3[7] = { |
471 | 3.529139957987837084554591421329876744262E-5L, |
472 | 2.973602667215766676998703687065066180115E-3L, |
473 | 8.273534546240864308494062287908662592100E-2L, |
474 | 9.613359842126507198241321110649974032726E-1L, |
475 | 4.853923697093974370118387947065402707519E0L, |
476 | 1.002671608961669247462020977417828796933E1L, |
477 | 7.028927383922483728931327850683151410267E0L, |
478 | /* 1.000000000000000000000000000000000000000E0L, */ |
479 | }; |
480 | |
481 | /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x), |
482 | Q1(x) = z(.375 + z^2 R(z^2)), z=1/x |
483 | 2 <= x <= 2.85711669921875 |
484 | Peak relative error 6.9e-22 */ |
485 | static const long double qr2[7] = { |
486 | -1.372751603025230017220666013816502528318E-4L, |
487 | -6.879190253347766576229143006767218972834E-3L, |
488 | -1.061253572090925414598304855316280077828E-1L, |
489 | -6.262164224345471241219408329354943337214E-1L, |
490 | -1.423149636514768476376254324731437473915E0L, |
491 | -1.087955310491078933531734062917489870754E0L, |
492 | -1.826821119773182847861406108689273719137E-1L, |
493 | }; |
494 | static const long double qs2[7] = { |
495 | 1.338768933634451601814048220627185324007E-3L, |
496 | 7.071099998918497559736318523932241901810E-2L, |
497 | 1.200511429784048632105295629933382142221E0L, |
498 | 8.327301713640367079030141077172031825276E0L, |
499 | 2.468479301872299311658145549931764426840E1L, |
500 | 2.961179686096262083509383820557051621644E1L, |
501 | 1.201402313144305153005639494661767354977E1L, |
502 | /* 1.000000000000000000000000000000000000000E0L, */ |
503 | }; |
504 | |
505 | |
506 | static long double |
507 | qone (long double x) |
508 | { |
509 | const long double *p, *q; |
510 | long double s, r, z; |
511 | int32_t ix; |
512 | uint32_t se, i0, i1; |
513 | |
514 | GET_LDOUBLE_WORDS (se, i0, i1, x); |
515 | ix = se & 0x7fff; |
516 | /* ix >= 0x4000 for all calls to this function. */ |
517 | if (ix >= 0x4002) /* x >= 8 */ |
518 | { |
519 | p = qr8; |
520 | q = qs8; |
521 | } |
522 | else |
523 | { |
524 | i1 = (ix << 16) | (i0 >> 16); |
525 | if (i1 >= 0x40019174) /* x >= 4.54541015625 */ |
526 | { |
527 | p = qr5; |
528 | q = qs5; |
529 | } |
530 | else if (i1 >= 0x4000b6db) /* x >= 2.85711669921875 */ |
531 | { |
532 | p = qr3; |
533 | q = qs3; |
534 | } |
535 | else /* x >= 2 */ |
536 | { |
537 | p = qr2; |
538 | q = qs2; |
539 | } |
540 | } |
541 | z = one / (x * x); |
542 | r = |
543 | p[0] + z * (p[1] + |
544 | z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6]))))); |
545 | s = |
546 | q[0] + z * (q[1] + |
547 | z * (q[2] + |
548 | z * (q[3] + z * (q[4] + z * (q[5] + z * (q[6] + z)))))); |
549 | return (.375 + z * r / s) / x; |
550 | } |
551 | |