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_j0(x), __ieee754_y0(x) |
34 | * Bessel function of the first and second kinds of order zero. |
35 | * Method -- j0(x): |
36 | * 1. For tiny x, we use j0(x) = 1 - x^2/4 + x^4/64 - ... |
37 | * 2. Reduce x to |x| since j0(x)=j0(-x), and |
38 | * for x in (0,2) |
39 | * j0(x) = 1 - z/4 + z^2*R0/S0, where z = x*x; |
40 | * for x in (2,inf) |
41 | * j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0)) |
42 | * where x0 = x-pi/4. It is better to compute sin(x0),cos(x0) |
43 | * as follow: |
44 | * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4) |
45 | * = 1/sqrt(2) * (cos(x) + sin(x)) |
46 | * sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4) |
47 | * = 1/sqrt(2) * (sin(x) - cos(x)) |
48 | * (To avoid cancellation, use |
49 | * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) |
50 | * to compute the worse one.) |
51 | * |
52 | * 3 Special cases |
53 | * j0(nan)= nan |
54 | * j0(0) = 1 |
55 | * j0(inf) = 0 |
56 | * |
57 | * Method -- y0(x): |
58 | * 1. For x<2. |
59 | * Since |
60 | * y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x^2/4 - ...) |
61 | * therefore y0(x)-2/pi*j0(x)*ln(x) is an even function. |
62 | * We use the following function to approximate y0, |
63 | * y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x^2 |
64 | * |
65 | * Note: For tiny x, U/V = u0 and j0(x)~1, hence |
66 | * y0(tiny) = u0 + (2/pi)*ln(tiny), (choose tiny<2**-27) |
67 | * 2. For x>=2. |
68 | * y0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)+q0(x)*sin(x0)) |
69 | * where x0 = x-pi/4. It is better to compute sin(x0),cos(x0) |
70 | * by the method mentioned above. |
71 | * 3. Special cases: y0(0)=-inf, y0(x<0)=NaN, y0(inf)=0. |
72 | */ |
73 | |
74 | #include <math.h> |
75 | #include <math_private.h> |
76 | |
77 | static long double pzero (long double), qzero (long double); |
78 | |
79 | static const long double |
80 | huge = 1e4930L, |
81 | one = 1.0L, |
82 | invsqrtpi = 5.6418958354775628694807945156077258584405e-1L, |
83 | tpi = 6.3661977236758134307553505349005744813784e-1L, |
84 | |
85 | /* J0(x) = 1 - x^2 / 4 + x^4 R0(x^2) / S0(x^2) |
86 | 0 <= x <= 2 |
87 | peak relative error 1.41e-22 */ |
88 | R[5] = { |
89 | 4.287176872744686992880841716723478740566E7L, |
90 | -6.652058897474241627570911531740907185772E5L, |
91 | 7.011848381719789863458364584613651091175E3L, |
92 | -3.168040850193372408702135490809516253693E1L, |
93 | 6.030778552661102450545394348845599300939E-2L, |
94 | }, |
95 | |
96 | S[4] = { |
97 | 2.743793198556599677955266341699130654342E9L, |
98 | 3.364330079384816249840086842058954076201E7L, |
99 | 1.924119649412510777584684927494642526573E5L, |
100 | 6.239282256012734914211715620088714856494E2L, |
101 | /* 1.000000000000000000000000000000000000000E0L,*/ |
102 | }; |
103 | |
104 | static const long double zero = 0.0; |
105 | |
106 | long double |
107 | __ieee754_j0l (long double x) |
108 | { |
109 | long double z, s, c, ss, cc, r, u, v; |
110 | int32_t ix; |
111 | uint32_t se; |
112 | |
113 | GET_LDOUBLE_EXP (se, x); |
114 | ix = se & 0x7fff; |
115 | if (__glibc_unlikely (ix >= 0x7fff)) |
116 | return one / (x * x); |
117 | x = fabsl (x); |
118 | if (ix >= 0x4000) /* |x| >= 2.0 */ |
119 | { |
120 | __sincosl (x, &s, &c); |
121 | ss = s - c; |
122 | cc = s + c; |
123 | if (ix < 0x7ffe) |
124 | { /* make sure x+x not overflow */ |
125 | z = -__cosl (x + x); |
126 | if ((s * c) < zero) |
127 | cc = z / ss; |
128 | else |
129 | ss = z / cc; |
130 | } |
131 | /* |
132 | * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) |
133 | * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) |
134 | */ |
135 | if (__glibc_unlikely (ix > 0x4080)) /* 2^129 */ |
136 | z = (invsqrtpi * cc) / __ieee754_sqrtl (x); |
137 | else |
138 | { |
139 | u = pzero (x); |
140 | v = qzero (x); |
141 | z = invsqrtpi * (u * cc - v * ss) / __ieee754_sqrtl (x); |
142 | } |
143 | return z; |
144 | } |
145 | if (__glibc_unlikely (ix < 0x3fef)) /* |x| < 2**-16 */ |
146 | { |
147 | /* raise inexact if x != 0 */ |
148 | math_force_eval (huge + x); |
149 | if (ix < 0x3fde) /* |x| < 2^-33 */ |
150 | return one; |
151 | else |
152 | return one - 0.25 * x * x; |
153 | } |
154 | z = x * x; |
155 | r = z * (R[0] + z * (R[1] + z * (R[2] + z * (R[3] + z * R[4])))); |
156 | s = S[0] + z * (S[1] + z * (S[2] + z * (S[3] + z))); |
157 | if (ix < 0x3fff) |
158 | { /* |x| < 1.00 */ |
159 | return (one - 0.25 * z + z * (r / s)); |
160 | } |
161 | else |
162 | { |
163 | u = 0.5 * x; |
164 | return ((one + u) * (one - u) + z * (r / s)); |
165 | } |
166 | } |
167 | strong_alias (__ieee754_j0l, __j0l_finite) |
168 | |
169 | |
170 | /* y0(x) = 2/pi ln(x) J0(x) + U(x^2)/V(x^2) |
171 | 0 < x <= 2 |
172 | peak relative error 1.7e-21 */ |
173 | static const long double |
174 | U[6] = { |
175 | -1.054912306975785573710813351985351350861E10L, |
176 | 2.520192609749295139432773849576523636127E10L, |
177 | -1.856426071075602001239955451329519093395E9L, |
178 | 4.079209129698891442683267466276785956784E7L, |
179 | -3.440684087134286610316661166492641011539E5L, |
180 | 1.005524356159130626192144663414848383774E3L, |
181 | }; |
182 | static const long double |
183 | V[5] = { |
184 | 1.429337283720789610137291929228082613676E11L, |
185 | 2.492593075325119157558811370165695013002E9L, |
186 | 2.186077620785925464237324417623665138376E7L, |
187 | 1.238407896366385175196515057064384929222E5L, |
188 | 4.693924035211032457494368947123233101664E2L, |
189 | /* 1.000000000000000000000000000000000000000E0L */ |
190 | }; |
191 | |
192 | long double |
193 | __ieee754_y0l (long double x) |
194 | { |
195 | long double z, s, c, ss, cc, u, v; |
196 | int32_t ix; |
197 | uint32_t se, i0, i1; |
198 | |
199 | GET_LDOUBLE_WORDS (se, i0, i1, x); |
200 | ix = se & 0x7fff; |
201 | /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */ |
202 | if (__glibc_unlikely (se & 0x8000)) |
203 | return zero / (zero * x); |
204 | if (__glibc_unlikely (ix >= 0x7fff)) |
205 | return one / (x + x * x); |
206 | if (__glibc_unlikely ((i0 | i1) == 0)) |
207 | return -HUGE_VALL + x; /* -inf and overflow exception. */ |
208 | if (ix >= 0x4000) |
209 | { /* |x| >= 2.0 */ |
210 | |
211 | /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0)) |
212 | * where x0 = x-pi/4 |
213 | * Better formula: |
214 | * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4) |
215 | * = 1/sqrt(2) * (sin(x) + cos(x)) |
216 | * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) |
217 | * = 1/sqrt(2) * (sin(x) - cos(x)) |
218 | * To avoid cancellation, use |
219 | * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) |
220 | * to compute the worse one. |
221 | */ |
222 | __sincosl (x, &s, &c); |
223 | ss = s - c; |
224 | cc = s + c; |
225 | /* |
226 | * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) |
227 | * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) |
228 | */ |
229 | if (ix < 0x7ffe) |
230 | { /* make sure x+x not overflow */ |
231 | z = -__cosl (x + x); |
232 | if ((s * c) < zero) |
233 | cc = z / ss; |
234 | else |
235 | ss = z / cc; |
236 | } |
237 | if (__glibc_unlikely (ix > 0x4080)) /* 1e39 */ |
238 | z = (invsqrtpi * ss) / __ieee754_sqrtl (x); |
239 | else |
240 | { |
241 | u = pzero (x); |
242 | v = qzero (x); |
243 | z = invsqrtpi * (u * ss + v * cc) / __ieee754_sqrtl (x); |
244 | } |
245 | return z; |
246 | } |
247 | if (__glibc_unlikely (ix <= 0x3fde)) /* x < 2^-33 */ |
248 | { |
249 | z = -7.380429510868722527629822444004602747322E-2L |
250 | + tpi * __ieee754_logl (x); |
251 | return z; |
252 | } |
253 | z = x * x; |
254 | u = U[0] + z * (U[1] + z * (U[2] + z * (U[3] + z * (U[4] + z * U[5])))); |
255 | v = V[0] + z * (V[1] + z * (V[2] + z * (V[3] + z * (V[4] + z)))); |
256 | return (u / v + tpi * (__ieee754_j0l (x) * __ieee754_logl (x))); |
257 | } |
258 | strong_alias (__ieee754_y0l, __y0l_finite) |
259 | |
260 | /* The asymptotic expansions of pzero is |
261 | * 1 - 9/128 s^2 + 11025/98304 s^4 - ..., where s = 1/x. |
262 | * For x >= 2, We approximate pzero by |
263 | * pzero(x) = 1 + s^2 R(s^2) / S(s^2) |
264 | */ |
265 | static const long double pR8[7] = { |
266 | /* 8 <= x <= inf |
267 | Peak relative error 4.62 */ |
268 | -4.094398895124198016684337960227780260127E-9L, |
269 | -8.929643669432412640061946338524096893089E-7L, |
270 | -6.281267456906136703868258380673108109256E-5L, |
271 | -1.736902783620362966354814353559382399665E-3L, |
272 | -1.831506216290984960532230842266070146847E-2L, |
273 | -5.827178869301452892963280214772398135283E-2L, |
274 | -2.087563267939546435460286895807046616992E-2L, |
275 | }; |
276 | static const long double pS8[6] = { |
277 | 5.823145095287749230197031108839653988393E-8L, |
278 | 1.279281986035060320477759999428992730280E-5L, |
279 | 9.132668954726626677174825517150228961304E-4L, |
280 | 2.606019379433060585351880541545146252534E-2L, |
281 | 2.956262215119520464228467583516287175244E-1L, |
282 | 1.149498145388256448535563278632697465675E0L, |
283 | /* 1.000000000000000000000000000000000000000E0L, */ |
284 | }; |
285 | |
286 | static const long double pR5[7] = { |
287 | /* 4.54541015625 <= x <= 8 |
288 | Peak relative error 6.51E-22 */ |
289 | -2.041226787870240954326915847282179737987E-7L, |
290 | -2.255373879859413325570636768224534428156E-5L, |
291 | -7.957485746440825353553537274569102059990E-4L, |
292 | -1.093205102486816696940149222095559439425E-2L, |
293 | -5.657957849316537477657603125260701114646E-2L, |
294 | -8.641175552716402616180994954177818461588E-2L, |
295 | -1.354654710097134007437166939230619726157E-2L, |
296 | }; |
297 | static const long double pS5[6] = { |
298 | 2.903078099681108697057258628212823545290E-6L, |
299 | 3.253948449946735405975737677123673867321E-4L, |
300 | 1.181269751723085006534147920481582279979E-2L, |
301 | 1.719212057790143888884745200257619469363E-1L, |
302 | 1.006306498779212467670654535430694221924E0L, |
303 | 2.069568808688074324555596301126375951502E0L, |
304 | /* 1.000000000000000000000000000000000000000E0L, */ |
305 | }; |
306 | |
307 | static const long double pR3[7] = { |
308 | /* 2.85711669921875 <= x <= 4.54541015625 |
309 | peak relative error 5.25e-21 */ |
310 | -5.755732156848468345557663552240816066802E-6L, |
311 | -3.703675625855715998827966962258113034767E-4L, |
312 | -7.390893350679637611641350096842846433236E-3L, |
313 | -5.571922144490038765024591058478043873253E-2L, |
314 | -1.531290690378157869291151002472627396088E-1L, |
315 | -1.193350853469302941921647487062620011042E-1L, |
316 | -8.567802507331578894302991505331963782905E-3L, |
317 | }; |
318 | static const long double pS3[6] = { |
319 | 8.185931139070086158103309281525036712419E-5L, |
320 | 5.398016943778891093520574483111255476787E-3L, |
321 | 1.130589193590489566669164765853409621081E-1L, |
322 | 9.358652328786413274673192987670237145071E-1L, |
323 | 3.091711512598349056276917907005098085273E0L, |
324 | 3.594602474737921977972586821673124231111E0L, |
325 | /* 1.000000000000000000000000000000000000000E0L, */ |
326 | }; |
327 | |
328 | static const long double pR2[7] = { |
329 | /* 2 <= x <= 2.85711669921875 |
330 | peak relative error 2.64e-21 */ |
331 | -1.219525235804532014243621104365384992623E-4L, |
332 | -4.838597135805578919601088680065298763049E-3L, |
333 | -5.732223181683569266223306197751407418301E-2L, |
334 | -2.472947430526425064982909699406646503758E-1L, |
335 | -3.753373645974077960207588073975976327695E-1L, |
336 | -1.556241316844728872406672349347137975495E-1L, |
337 | -5.355423239526452209595316733635519506958E-3L, |
338 | }; |
339 | static const long double pS2[6] = { |
340 | 1.734442793664291412489066256138894953823E-3L, |
341 | 7.158111826468626405416300895617986926008E-2L, |
342 | 9.153839713992138340197264669867993552641E-1L, |
343 | 4.539209519433011393525841956702487797582E0L, |
344 | 8.868932430625331650266067101752626253644E0L, |
345 | 6.067161890196324146320763844772857713502E0L, |
346 | /* 1.000000000000000000000000000000000000000E0L, */ |
347 | }; |
348 | |
349 | static long double |
350 | pzero (long double x) |
351 | { |
352 | const long double *p, *q; |
353 | long double z, r, s; |
354 | int32_t ix; |
355 | uint32_t se, i0, i1; |
356 | |
357 | GET_LDOUBLE_WORDS (se, i0, i1, x); |
358 | ix = se & 0x7fff; |
359 | /* ix >= 0x4000 for all calls to this function. */ |
360 | if (ix >= 0x4002) |
361 | { |
362 | p = pR8; |
363 | q = pS8; |
364 | } /* x >= 8 */ |
365 | else |
366 | { |
367 | i1 = (ix << 16) | (i0 >> 16); |
368 | if (i1 >= 0x40019174) /* x >= 4.54541015625 */ |
369 | { |
370 | p = pR5; |
371 | q = pS5; |
372 | } |
373 | else if (i1 >= 0x4000b6db) /* x >= 2.85711669921875 */ |
374 | { |
375 | p = pR3; |
376 | q = pS3; |
377 | } |
378 | else /* x >= 2 */ |
379 | { |
380 | p = pR2; |
381 | q = pS2; |
382 | } |
383 | } |
384 | z = one / (x * x); |
385 | r = |
386 | p[0] + z * (p[1] + |
387 | z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6]))))); |
388 | s = |
389 | q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * (q[5] + z))))); |
390 | return (one + z * r / s); |
391 | } |
392 | |
393 | |
394 | /* For x >= 8, the asymptotic expansions of qzero is |
395 | * -1/8 s + 75/1024 s^3 - ..., where s = 1/x. |
396 | * We approximate qzero by |
397 | * qzero(x) = s*(-.125 + R(s^2) / S(s^2)) |
398 | */ |
399 | static const long double qR8[7] = { |
400 | /* 8 <= x <= inf |
401 | peak relative error 2.23e-21 */ |
402 | 3.001267180483191397885272640777189348008E-10L, |
403 | 8.693186311430836495238494289942413810121E-8L, |
404 | 8.496875536711266039522937037850596580686E-6L, |
405 | 3.482702869915288984296602449543513958409E-4L, |
406 | 6.036378380706107692863811938221290851352E-3L, |
407 | 3.881970028476167836382607922840452192636E-2L, |
408 | 6.132191514516237371140841765561219149638E-2L, |
409 | }; |
410 | static const long double qS8[7] = { |
411 | 4.097730123753051126914971174076227600212E-9L, |
412 | 1.199615869122646109596153392152131139306E-6L, |
413 | 1.196337580514532207793107149088168946451E-4L, |
414 | 5.099074440112045094341500497767181211104E-3L, |
415 | 9.577420799632372483249761659674764460583E-2L, |
416 | 7.385243015344292267061953461563695918646E-1L, |
417 | 1.917266424391428937962682301561699055943E0L, |
418 | /* 1.000000000000000000000000000000000000000E0L, */ |
419 | }; |
420 | |
421 | static const long double qR5[7] = { |
422 | /* 4.54541015625 <= x <= 8 |
423 | peak relative error 1.03e-21 */ |
424 | 3.406256556438974327309660241748106352137E-8L, |
425 | 4.855492710552705436943630087976121021980E-6L, |
426 | 2.301011739663737780613356017352912281980E-4L, |
427 | 4.500470249273129953870234803596619899226E-3L, |
428 | 3.651376459725695502726921248173637054828E-2L, |
429 | 1.071578819056574524416060138514508609805E-1L, |
430 | 7.458950172851611673015774675225656063757E-2L, |
431 | }; |
432 | static const long double qS5[7] = { |
433 | 4.650675622764245276538207123618745150785E-7L, |
434 | 6.773573292521412265840260065635377164455E-5L, |
435 | 3.340711249876192721980146877577806687714E-3L, |
436 | 7.036218046856839214741678375536970613501E-2L, |
437 | 6.569599559163872573895171876511377891143E-1L, |
438 | 2.557525022583599204591036677199171155186E0L, |
439 | 3.457237396120935674982927714210361269133E0L, |
440 | /* 1.000000000000000000000000000000000000000E0L,*/ |
441 | }; |
442 | |
443 | static const long double qR3[7] = { |
444 | /* 2.85711669921875 <= x <= 4.54541015625 |
445 | peak relative error 5.24e-21 */ |
446 | 1.749459596550816915639829017724249805242E-6L, |
447 | 1.446252487543383683621692672078376929437E-4L, |
448 | 3.842084087362410664036704812125005761859E-3L, |
449 | 4.066369994699462547896426554180954233581E-2L, |
450 | 1.721093619117980251295234795188992722447E-1L, |
451 | 2.538595333972857367655146949093055405072E-1L, |
452 | 8.560591367256769038905328596020118877936E-2L, |
453 | }; |
454 | static const long double qS3[7] = { |
455 | 2.388596091707517488372313710647510488042E-5L, |
456 | 2.048679968058758616370095132104333998147E-3L, |
457 | 5.824663198201417760864458765259945181513E-2L, |
458 | 6.953906394693328750931617748038994763958E-1L, |
459 | 3.638186936390881159685868764832961092476E0L, |
460 | 7.900169524705757837298990558459547842607E0L, |
461 | 5.992718532451026507552820701127504582907E0L, |
462 | /* 1.000000000000000000000000000000000000000E0L, */ |
463 | }; |
464 | |
465 | static const long double qR2[7] = { |
466 | /* 2 <= x <= 2.85711669921875 |
467 | peak relative error 1.58e-21 */ |
468 | 6.306524405520048545426928892276696949540E-5L, |
469 | 3.209606155709930950935893996591576624054E-3L, |
470 | 5.027828775702022732912321378866797059604E-2L, |
471 | 3.012705561838718956481911477587757845163E-1L, |
472 | 6.960544893905752937420734884995688523815E-1L, |
473 | 5.431871999743531634887107835372232030655E-1L, |
474 | 9.447736151202905471899259026430157211949E-2L, |
475 | }; |
476 | static const long double qS2[7] = { |
477 | 8.610579901936193494609755345106129102676E-4L, |
478 | 4.649054352710496997203474853066665869047E-2L, |
479 | 8.104282924459837407218042945106320388339E-1L, |
480 | 5.807730930825886427048038146088828206852E0L, |
481 | 1.795310145936848873627710102199881642939E1L, |
482 | 2.281313316875375733663657188888110605044E1L, |
483 | 1.011242067883822301487154844458322200143E1L, |
484 | /* 1.000000000000000000000000000000000000000E0L, */ |
485 | }; |
486 | |
487 | static long double |
488 | qzero (long double x) |
489 | { |
490 | const long double *p, *q; |
491 | long double s, r, z; |
492 | int32_t ix; |
493 | uint32_t se, i0, i1; |
494 | |
495 | GET_LDOUBLE_WORDS (se, i0, i1, x); |
496 | ix = se & 0x7fff; |
497 | /* ix >= 0x4000 for all calls to this function. */ |
498 | if (ix >= 0x4002) /* x >= 8 */ |
499 | { |
500 | p = qR8; |
501 | q = qS8; |
502 | } |
503 | else |
504 | { |
505 | i1 = (ix << 16) | (i0 >> 16); |
506 | if (i1 >= 0x40019174) /* x >= 4.54541015625 */ |
507 | { |
508 | p = qR5; |
509 | q = qS5; |
510 | } |
511 | else if (i1 >= 0x4000b6db) /* x >= 2.85711669921875 */ |
512 | { |
513 | p = qR3; |
514 | q = qS3; |
515 | } |
516 | else /* x >= 2 */ |
517 | { |
518 | p = qR2; |
519 | q = qS2; |
520 | } |
521 | } |
522 | z = one / (x * x); |
523 | r = |
524 | p[0] + z * (p[1] + |
525 | z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6]))))); |
526 | s = |
527 | q[0] + z * (q[1] + |
528 | z * (q[2] + |
529 | z * (q[3] + z * (q[4] + z * (q[5] + z * (q[6] + z)))))); |
530 | return (-.125 + z * r / s) / x; |
531 | } |
532 | |