]>
Commit | Line | Data |
---|---|---|
b45e4f52 AJ |
1 | /* j1l.c |
2 | * | |
3 | * Bessel function of order one | |
4 | * | |
5 | * | |
6 | * | |
7 | * SYNOPSIS: | |
8 | * | |
9 | * long double x, y, j1l(); | |
10 | * | |
11 | * y = j1l( x ); | |
12 | * | |
13 | * | |
14 | * | |
15 | * DESCRIPTION: | |
16 | * | |
17 | * Returns Bessel function of first kind, order one of the argument. | |
18 | * | |
19 | * The domain is divided into two major intervals [0, 2] and | |
20 | * (2, infinity). In the first interval the rational approximation is | |
21 | * J1(x) = .5x + x x^2 R(x^2) | |
22 | * | |
23 | * The second interval is further partitioned into eight equal segments | |
24 | * of 1/x. | |
25 | * J1(x) = sqrt(2/(pi x)) (P1(x) cos(X) - Q1(x) sin(X)), | |
26 | * X = x - 3 pi / 4, | |
27 | * | |
28 | * and the auxiliary functions are given by | |
29 | * | |
30 | * J1(x)cos(X) + Y1(x)sin(X) = sqrt( 2/(pi x)) P1(x), | |
31 | * P1(x) = 1 + 1/x^2 R(1/x^2) | |
32 | * | |
33 | * Y1(x)cos(X) - J1(x)sin(X) = sqrt( 2/(pi x)) Q1(x), | |
34 | * Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)). | |
35 | * | |
36 | * | |
37 | * | |
38 | * ACCURACY: | |
39 | * | |
40 | * Absolute error: | |
41 | * arithmetic domain # trials peak rms | |
42 | * IEEE 0, 30 100000 2.8e-34 2.7e-35 | |
43 | * | |
44 | * | |
45 | */ | |
46 | ||
47 | /* y1l.c | |
48 | * | |
49 | * Bessel function of the second kind, order one | |
50 | * | |
51 | * | |
52 | * | |
53 | * SYNOPSIS: | |
54 | * | |
55 | * double x, y, y1l(); | |
56 | * | |
57 | * y = y1l( x ); | |
58 | * | |
59 | * | |
60 | * | |
61 | * DESCRIPTION: | |
62 | * | |
63 | * Returns Bessel function of the second kind, of order | |
64 | * one, of the argument. | |
65 | * | |
66 | * The domain is divided into two major intervals [0, 2] and | |
67 | * (2, infinity). In the first interval the rational approximation is | |
68 | * Y1(x) = 2/pi * (log(x) * J1(x) - 1/x) + x R(x^2) . | |
69 | * In the second interval the approximation is the same as for J1(x), and | |
70 | * Y1(x) = sqrt(2/(pi x)) (P1(x) sin(X) + Q1(x) cos(X)), | |
71 | * X = x - 3 pi / 4. | |
72 | * | |
73 | * ACCURACY: | |
74 | * | |
75 | * Absolute error, when y0(x) < 1; else relative error: | |
76 | * | |
77 | * arithmetic domain # trials peak rms | |
78 | * IEEE 0, 30 100000 2.7e-34 2.9e-35 | |
79 | * | |
80 | */ | |
81 | ||
cc7375ce RM |
82 | /* Copyright 2001 by Stephen L. Moshier (moshier@na-net.onrl.gov). |
83 | ||
84 | This library is free software; you can redistribute it and/or | |
85 | modify it under the terms of the GNU Lesser General Public | |
86 | License as published by the Free Software Foundation; either | |
87 | version 2.1 of the License, or (at your option) any later version. | |
88 | ||
89 | This library is distributed in the hope that it will be useful, | |
90 | but WITHOUT ANY WARRANTY; without even the implied warranty of | |
91 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
92 | Lesser General Public License for more details. | |
93 | ||
94 | You should have received a copy of the GNU Lesser General Public | |
59ba27a6 | 95 | License along with this library; if not, see |
5a82c748 | 96 | <https://www.gnu.org/licenses/>. */ |
b45e4f52 | 97 | |
46a3d3c7 | 98 | #include <errno.h> |
1ed0291c RH |
99 | #include <math.h> |
100 | #include <math_private.h> | |
70e2ba33 | 101 | #include <fenv_private.h> |
8f5b00d3 | 102 | #include <math-underflow.h> |
98c48fe5 | 103 | #include <float.h> |
b45e4f52 AJ |
104 | |
105 | /* 1 / sqrt(pi) */ | |
02bbfb41 | 106 | static const _Float128 ONEOSQPI = L(5.6418958354775628694807945156077258584405E-1); |
b45e4f52 | 107 | /* 2 / pi */ |
02bbfb41 PM |
108 | static const _Float128 TWOOPI = L(6.3661977236758134307553505349005744813784E-1); |
109 | static const _Float128 zero = 0; | |
b45e4f52 AJ |
110 | |
111 | /* J1(x) = .5x + x x^2 R(x^2) | |
112 | Peak relative error 1.9e-35 | |
113 | 0 <= x <= 2 */ | |
114 | #define NJ0_2N 6 | |
15089e04 | 115 | static const _Float128 J0_2N[NJ0_2N + 1] = { |
02bbfb41 PM |
116 | L(-5.943799577386942855938508697619735179660E16), |
117 | L(1.812087021305009192259946997014044074711E15), | |
118 | L(-2.761698314264509665075127515729146460895E13), | |
119 | L(2.091089497823600978949389109350658815972E11), | |
120 | L(-8.546413231387036372945453565654130054307E8), | |
121 | L(1.797229225249742247475464052741320612261E6), | |
122 | L(-1.559552840946694171346552770008812083969E3) | |
b45e4f52 AJ |
123 | }; |
124 | #define NJ0_2D 6 | |
15089e04 | 125 | static const _Float128 J0_2D[NJ0_2D + 1] = { |
02bbfb41 PM |
126 | L(9.510079323819108569501613916191477479397E17), |
127 | L(1.063193817503280529676423936545854693915E16), | |
128 | L(5.934143516050192600795972192791775226920E13), | |
129 | L(2.168000911950620999091479265214368352883E11), | |
130 | L(5.673775894803172808323058205986256928794E8), | |
131 | L(1.080329960080981204840966206372671147224E6), | |
132 | L(1.411951256636576283942477881535283304912E3), | |
b45e4f52 AJ |
133 | /* 1.000000000000000000000000000000000000000E0L */ |
134 | }; | |
135 | ||
136 | /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2), | |
137 | 0 <= 1/x <= .0625 | |
138 | Peak relative error 3.6e-36 */ | |
139 | #define NP16_IN 9 | |
15089e04 | 140 | static const _Float128 P16_IN[NP16_IN + 1] = { |
02bbfb41 PM |
141 | L(5.143674369359646114999545149085139822905E-16), |
142 | L(4.836645664124562546056389268546233577376E-13), | |
143 | L(1.730945562285804805325011561498453013673E-10), | |
144 | L(3.047976856147077889834905908605310585810E-8), | |
145 | L(2.855227609107969710407464739188141162386E-6), | |
146 | L(1.439362407936705484122143713643023998457E-4), | |
147 | L(3.774489768532936551500999699815873422073E-3), | |
148 | L(4.723962172984642566142399678920790598426E-2), | |
149 | L(2.359289678988743939925017240478818248735E-1), | |
150 | L(3.032580002220628812728954785118117124520E-1), | |
b45e4f52 AJ |
151 | }; |
152 | #define NP16_ID 9 | |
15089e04 | 153 | static const _Float128 P16_ID[NP16_ID + 1] = { |
02bbfb41 PM |
154 | L(4.389268795186898018132945193912677177553E-15), |
155 | L(4.132671824807454334388868363256830961655E-12), | |
156 | L(1.482133328179508835835963635130894413136E-9), | |
157 | L(2.618941412861122118906353737117067376236E-7), | |
158 | L(2.467854246740858470815714426201888034270E-5), | |
159 | L(1.257192927368839847825938545925340230490E-3), | |
160 | L(3.362739031941574274949719324644120720341E-2), | |
161 | L(4.384458231338934105875343439265370178858E-1), | |
162 | L(2.412830809841095249170909628197264854651E0), | |
163 | L(4.176078204111348059102962617368214856874E0), | |
b45e4f52 AJ |
164 | /* 1.000000000000000000000000000000000000000E0 */ |
165 | }; | |
166 | ||
167 | /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2), | |
168 | 0.0625 <= 1/x <= 0.125 | |
169 | Peak relative error 1.9e-36 */ | |
170 | #define NP8_16N 11 | |
15089e04 | 171 | static const _Float128 P8_16N[NP8_16N + 1] = { |
02bbfb41 PM |
172 | L(2.984612480763362345647303274082071598135E-16), |
173 | L(1.923651877544126103941232173085475682334E-13), | |
174 | L(4.881258879388869396043760693256024307743E-11), | |
175 | L(6.368866572475045408480898921866869811889E-9), | |
176 | L(4.684818344104910450523906967821090796737E-7), | |
177 | L(2.005177298271593587095982211091300382796E-5), | |
178 | L(4.979808067163957634120681477207147536182E-4), | |
179 | L(6.946005761642579085284689047091173581127E-3), | |
180 | L(5.074601112955765012750207555985299026204E-2), | |
181 | L(1.698599455896180893191766195194231825379E-1), | |
182 | L(1.957536905259237627737222775573623779638E-1), | |
183 | L(2.991314703282528370270179989044994319374E-2), | |
b45e4f52 AJ |
184 | }; |
185 | #define NP8_16D 10 | |
15089e04 | 186 | static const _Float128 P8_16D[NP8_16D + 1] = { |
02bbfb41 PM |
187 | L(2.546869316918069202079580939942463010937E-15), |
188 | L(1.644650111942455804019788382157745229955E-12), | |
189 | L(4.185430770291694079925607420808011147173E-10), | |
190 | L(5.485331966975218025368698195861074143153E-8), | |
191 | L(4.062884421686912042335466327098932678905E-6), | |
192 | L(1.758139661060905948870523641319556816772E-4), | |
193 | L(4.445143889306356207566032244985607493096E-3), | |
194 | L(6.391901016293512632765621532571159071158E-2), | |
195 | L(4.933040207519900471177016015718145795434E-1), | |
196 | L(1.839144086168947712971630337250761842976E0), | |
197 | L(2.715120873995490920415616716916149586579E0), | |
b45e4f52 AJ |
198 | /* 1.000000000000000000000000000000000000000E0 */ |
199 | }; | |
200 | ||
201 | /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2), | |
202 | 0.125 <= 1/x <= 0.1875 | |
203 | Peak relative error 1.3e-36 */ | |
204 | #define NP5_8N 10 | |
15089e04 | 205 | static const _Float128 P5_8N[NP5_8N + 1] = { |
02bbfb41 PM |
206 | L(2.837678373978003452653763806968237227234E-12), |
207 | L(9.726641165590364928442128579282742354806E-10), | |
208 | L(1.284408003604131382028112171490633956539E-7), | |
209 | L(8.524624695868291291250573339272194285008E-6), | |
210 | L(3.111516908953172249853673787748841282846E-4), | |
211 | L(6.423175156126364104172801983096596409176E-3), | |
212 | L(7.430220589989104581004416356260692450652E-2), | |
213 | L(4.608315409833682489016656279567605536619E-1), | |
214 | L(1.396870223510964882676225042258855977512E0), | |
215 | L(1.718500293904122365894630460672081526236E0), | |
216 | L(5.465927698800862172307352821870223855365E-1) | |
b45e4f52 AJ |
217 | }; |
218 | #define NP5_8D 10 | |
15089e04 | 219 | static const _Float128 P5_8D[NP5_8D + 1] = { |
02bbfb41 PM |
220 | L(2.421485545794616609951168511612060482715E-11), |
221 | L(8.329862750896452929030058039752327232310E-9), | |
222 | L(1.106137992233383429630592081375289010720E-6), | |
223 | L(7.405786153760681090127497796448503306939E-5), | |
224 | L(2.740364785433195322492093333127633465227E-3), | |
225 | L(5.781246470403095224872243564165254652198E-2), | |
226 | L(6.927711353039742469918754111511109983546E-1), | |
227 | L(4.558679283460430281188304515922826156690E0), | |
228 | L(1.534468499844879487013168065728837900009E1), | |
229 | L(2.313927430889218597919624843161569422745E1), | |
230 | L(1.194506341319498844336768473218382828637E1), | |
b45e4f52 AJ |
231 | /* 1.000000000000000000000000000000000000000E0 */ |
232 | }; | |
233 | ||
234 | /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2), | |
235 | Peak relative error 1.4e-36 | |
236 | 0.1875 <= 1/x <= 0.25 */ | |
237 | #define NP4_5N 10 | |
15089e04 | 238 | static const _Float128 P4_5N[NP4_5N + 1] = { |
02bbfb41 PM |
239 | L(1.846029078268368685834261260420933914621E-10), |
240 | L(3.916295939611376119377869680335444207768E-8), | |
241 | L(3.122158792018920627984597530935323997312E-6), | |
242 | L(1.218073444893078303994045653603392272450E-4), | |
243 | L(2.536420827983485448140477159977981844883E-3), | |
244 | L(2.883011322006690823959367922241169171315E-2), | |
245 | L(1.755255190734902907438042414495469810830E-1), | |
246 | L(5.379317079922628599870898285488723736599E-1), | |
247 | L(7.284904050194300773890303361501726561938E-1), | |
248 | L(3.270110346613085348094396323925000362813E-1), | |
249 | L(1.804473805689725610052078464951722064757E-2), | |
b45e4f52 AJ |
250 | }; |
251 | #define NP4_5D 9 | |
15089e04 | 252 | static const _Float128 P4_5D[NP4_5D + 1] = { |
02bbfb41 PM |
253 | L(1.575278146806816970152174364308980863569E-9), |
254 | L(3.361289173657099516191331123405675054321E-7), | |
255 | L(2.704692281550877810424745289838790693708E-5), | |
256 | L(1.070854930483999749316546199273521063543E-3), | |
257 | L(2.282373093495295842598097265627962125411E-2), | |
258 | L(2.692025460665354148328762368240343249830E-1), | |
259 | L(1.739892942593664447220951225734811133759E0), | |
260 | L(5.890727576752230385342377570386657229324E0), | |
261 | L(9.517442287057841500750256954117735128153E0), | |
262 | L(6.100616353935338240775363403030137736013E0), | |
b45e4f52 AJ |
263 | /* 1.000000000000000000000000000000000000000E0 */ |
264 | }; | |
265 | ||
266 | /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2), | |
267 | Peak relative error 3.0e-36 | |
268 | 0.25 <= 1/x <= 0.3125 */ | |
269 | #define NP3r2_4N 9 | |
15089e04 | 270 | static const _Float128 P3r2_4N[NP3r2_4N + 1] = { |
02bbfb41 PM |
271 | L(8.240803130988044478595580300846665863782E-8), |
272 | L(1.179418958381961224222969866406483744580E-5), | |
273 | L(6.179787320956386624336959112503824397755E-4), | |
274 | L(1.540270833608687596420595830747166658383E-2), | |
275 | L(1.983904219491512618376375619598837355076E-1), | |
276 | L(1.341465722692038870390470651608301155565E0), | |
277 | L(4.617865326696612898792238245990854646057E0), | |
278 | L(7.435574801812346424460233180412308000587E0), | |
279 | L(4.671327027414635292514599201278557680420E0), | |
280 | L(7.299530852495776936690976966995187714739E-1), | |
b45e4f52 AJ |
281 | }; |
282 | #define NP3r2_4D 9 | |
15089e04 | 283 | static const _Float128 P3r2_4D[NP3r2_4D + 1] = { |
02bbfb41 PM |
284 | L(7.032152009675729604487575753279187576521E-7), |
285 | L(1.015090352324577615777511269928856742848E-4), | |
286 | L(5.394262184808448484302067955186308730620E-3), | |
287 | L(1.375291438480256110455809354836988584325E-1), | |
288 | L(1.836247144461106304788160919310404376670E0), | |
289 | L(1.314378564254376655001094503090935880349E1), | |
290 | L(4.957184590465712006934452500894672343488E1), | |
291 | L(9.287394244300647738855415178790263465398E1), | |
292 | L(7.652563275535900609085229286020552768399E1), | |
293 | L(2.147042473003074533150718117770093209096E1), | |
b45e4f52 AJ |
294 | /* 1.000000000000000000000000000000000000000E0 */ |
295 | }; | |
296 | ||
297 | /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2), | |
298 | Peak relative error 1.0e-35 | |
299 | 0.3125 <= 1/x <= 0.375 */ | |
300 | #define NP2r7_3r2N 9 | |
15089e04 | 301 | static const _Float128 P2r7_3r2N[NP2r7_3r2N + 1] = { |
02bbfb41 PM |
302 | L(4.599033469240421554219816935160627085991E-7), |
303 | L(4.665724440345003914596647144630893997284E-5), | |
304 | L(1.684348845667764271596142716944374892756E-3), | |
305 | L(2.802446446884455707845985913454440176223E-2), | |
306 | L(2.321937586453963310008279956042545173930E-1), | |
307 | L(9.640277413988055668692438709376437553804E-1), | |
308 | L(1.911021064710270904508663334033003246028E0), | |
309 | L(1.600811610164341450262992138893970224971E0), | |
310 | L(4.266299218652587901171386591543457861138E-1), | |
311 | L(1.316470424456061252962568223251247207325E-2), | |
b45e4f52 AJ |
312 | }; |
313 | #define NP2r7_3r2D 8 | |
15089e04 | 314 | static const _Float128 P2r7_3r2D[NP2r7_3r2D + 1] = { |
02bbfb41 PM |
315 | L(3.924508608545520758883457108453520099610E-6), |
316 | L(4.029707889408829273226495756222078039823E-4), | |
317 | L(1.484629715787703260797886463307469600219E-2), | |
318 | L(2.553136379967180865331706538897231588685E-1), | |
319 | L(2.229457223891676394409880026887106228740E0), | |
320 | L(1.005708903856384091956550845198392117318E1), | |
321 | L(2.277082659664386953166629360352385889558E1), | |
322 | L(2.384726835193630788249826630376533988245E1), | |
323 | L(9.700989749041320895890113781610939632410E0), | |
b45e4f52 AJ |
324 | /* 1.000000000000000000000000000000000000000E0 */ |
325 | }; | |
326 | ||
327 | /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2), | |
328 | Peak relative error 1.7e-36 | |
329 | 0.3125 <= 1/x <= 0.4375 */ | |
330 | #define NP2r3_2r7N 9 | |
15089e04 | 331 | static const _Float128 P2r3_2r7N[NP2r3_2r7N + 1] = { |
02bbfb41 PM |
332 | L(3.916766777108274628543759603786857387402E-6), |
333 | L(3.212176636756546217390661984304645137013E-4), | |
334 | L(9.255768488524816445220126081207248947118E-3), | |
335 | L(1.214853146369078277453080641911700735354E-1), | |
336 | L(7.855163309847214136198449861311404633665E-1), | |
337 | L(2.520058073282978403655488662066019816540E0), | |
338 | L(3.825136484837545257209234285382183711466E0), | |
339 | L(2.432569427554248006229715163865569506873E0), | |
340 | L(4.877934835018231178495030117729800489743E-1), | |
341 | L(1.109902737860249670981355149101343427885E-2), | |
b45e4f52 AJ |
342 | }; |
343 | #define NP2r3_2r7D 8 | |
15089e04 | 344 | static const _Float128 P2r3_2r7D[NP2r3_2r7D + 1] = { |
02bbfb41 PM |
345 | L(3.342307880794065640312646341190547184461E-5), |
346 | L(2.782182891138893201544978009012096558265E-3), | |
347 | L(8.221304931614200702142049236141249929207E-2), | |
348 | L(1.123728246291165812392918571987858010949E0), | |
349 | L(7.740482453652715577233858317133423434590E0), | |
350 | L(2.737624677567945952953322566311201919139E1), | |
351 | L(4.837181477096062403118304137851260715475E1), | |
352 | L(3.941098643468580791437772701093795299274E1), | |
353 | L(1.245821247166544627558323920382547533630E1), | |
b45e4f52 AJ |
354 | /* 1.000000000000000000000000000000000000000E0 */ |
355 | }; | |
356 | ||
357 | /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2), | |
358 | Peak relative error 1.7e-35 | |
359 | 0.4375 <= 1/x <= 0.5 */ | |
360 | #define NP2_2r3N 8 | |
15089e04 | 361 | static const _Float128 P2_2r3N[NP2_2r3N + 1] = { |
02bbfb41 PM |
362 | L(3.397930802851248553545191160608731940751E-4), |
363 | L(2.104020902735482418784312825637833698217E-2), | |
364 | L(4.442291771608095963935342749477836181939E-1), | |
365 | L(4.131797328716583282869183304291833754967E0), | |
366 | L(1.819920169779026500146134832455189917589E1), | |
367 | L(3.781779616522937565300309684282401791291E1), | |
368 | L(3.459605449728864218972931220783543410347E1), | |
369 | L(1.173594248397603882049066603238568316561E1), | |
370 | L(9.455702270242780642835086549285560316461E-1), | |
b45e4f52 AJ |
371 | }; |
372 | #define NP2_2r3D 8 | |
15089e04 | 373 | static const _Float128 P2_2r3D[NP2_2r3D + 1] = { |
02bbfb41 PM |
374 | L(2.899568897241432883079888249845707400614E-3), |
375 | L(1.831107138190848460767699919531132426356E-1), | |
376 | L(3.999350044057883839080258832758908825165E0), | |
377 | L(3.929041535867957938340569419874195303712E1), | |
378 | L(1.884245613422523323068802689915538908291E2), | |
379 | L(4.461469948819229734353852978424629815929E2), | |
380 | L(5.004998753999796821224085972610636347903E2), | |
381 | L(2.386342520092608513170837883757163414100E2), | |
382 | L(3.791322528149347975999851588922424189957E1), | |
b45e4f52 AJ |
383 | /* 1.000000000000000000000000000000000000000E0 */ |
384 | }; | |
385 | ||
386 | /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x), | |
387 | Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)), | |
388 | Peak relative error 8.0e-36 | |
389 | 0 <= 1/x <= .0625 */ | |
390 | #define NQ16_IN 10 | |
15089e04 | 391 | static const _Float128 Q16_IN[NQ16_IN + 1] = { |
02bbfb41 PM |
392 | L(-3.917420835712508001321875734030357393421E-18), |
393 | L(-4.440311387483014485304387406538069930457E-15), | |
394 | L(-1.951635424076926487780929645954007139616E-12), | |
395 | L(-4.318256438421012555040546775651612810513E-10), | |
396 | L(-5.231244131926180765270446557146989238020E-8), | |
397 | L(-3.540072702902043752460711989234732357653E-6), | |
398 | L(-1.311017536555269966928228052917534882984E-4), | |
399 | L(-2.495184669674631806622008769674827575088E-3), | |
400 | L(-2.141868222987209028118086708697998506716E-2), | |
401 | L(-6.184031415202148901863605871197272650090E-2), | |
402 | L(-1.922298704033332356899546792898156493887E-2), | |
b45e4f52 AJ |
403 | }; |
404 | #define NQ16_ID 9 | |
15089e04 | 405 | static const _Float128 Q16_ID[NQ16_ID + 1] = { |
02bbfb41 PM |
406 | L(3.820418034066293517479619763498400162314E-17), |
407 | L(4.340702810799239909648911373329149354911E-14), | |
408 | L(1.914985356383416140706179933075303538524E-11), | |
409 | L(4.262333682610888819476498617261895474330E-9), | |
410 | L(5.213481314722233980346462747902942182792E-7), | |
411 | L(3.585741697694069399299005316809954590558E-5), | |
412 | L(1.366513429642842006385029778105539457546E-3), | |
413 | L(2.745282599850704662726337474371355160594E-2), | |
414 | L(2.637644521611867647651200098449903330074E-1), | |
415 | L(1.006953426110765984590782655598680488746E0), | |
b45e4f52 AJ |
416 | /* 1.000000000000000000000000000000000000000E0 */ |
417 | }; | |
418 | ||
419 | /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x), | |
420 | Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)), | |
421 | Peak relative error 1.9e-36 | |
422 | 0.0625 <= 1/x <= 0.125 */ | |
423 | #define NQ8_16N 11 | |
15089e04 | 424 | static const _Float128 Q8_16N[NQ8_16N + 1] = { |
02bbfb41 PM |
425 | L(-2.028630366670228670781362543615221542291E-17), |
426 | L(-1.519634620380959966438130374006858864624E-14), | |
427 | L(-4.540596528116104986388796594639405114524E-12), | |
428 | L(-7.085151756671466559280490913558388648274E-10), | |
429 | L(-6.351062671323970823761883833531546885452E-8), | |
430 | L(-3.390817171111032905297982523519503522491E-6), | |
431 | L(-1.082340897018886970282138836861233213972E-4), | |
432 | L(-2.020120801187226444822977006648252379508E-3), | |
433 | L(-2.093169910981725694937457070649605557555E-2), | |
434 | L(-1.092176538874275712359269481414448063393E-1), | |
435 | L(-2.374790947854765809203590474789108718733E-1), | |
436 | L(-1.365364204556573800719985118029601401323E-1), | |
b45e4f52 AJ |
437 | }; |
438 | #define NQ8_16D 11 | |
15089e04 | 439 | static const _Float128 Q8_16D[NQ8_16D + 1] = { |
02bbfb41 PM |
440 | L(1.978397614733632533581207058069628242280E-16), |
441 | L(1.487361156806202736877009608336766720560E-13), | |
442 | L(4.468041406888412086042576067133365913456E-11), | |
443 | L(7.027822074821007443672290507210594648877E-9), | |
444 | L(6.375740580686101224127290062867976007374E-7), | |
445 | L(3.466887658320002225888644977076410421940E-5), | |
446 | L(1.138625640905289601186353909213719596986E-3), | |
447 | L(2.224470799470414663443449818235008486439E-2), | |
448 | L(2.487052928527244907490589787691478482358E-1), | |
449 | L(1.483927406564349124649083853892380899217E0), | |
450 | L(4.182773513276056975777258788903489507705E0), | |
451 | L(4.419665392573449746043880892524360870944E0), | |
b45e4f52 AJ |
452 | /* 1.000000000000000000000000000000000000000E0 */ |
453 | }; | |
454 | ||
455 | /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x), | |
456 | Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)), | |
457 | Peak relative error 1.5e-35 | |
458 | 0.125 <= 1/x <= 0.1875 */ | |
459 | #define NQ5_8N 10 | |
15089e04 | 460 | static const _Float128 Q5_8N[NQ5_8N + 1] = { |
02bbfb41 PM |
461 | L(-3.656082407740970534915918390488336879763E-13), |
462 | L(-1.344660308497244804752334556734121771023E-10), | |
463 | L(-1.909765035234071738548629788698150760791E-8), | |
464 | L(-1.366668038160120210269389551283666716453E-6), | |
465 | L(-5.392327355984269366895210704976314135683E-5), | |
466 | L(-1.206268245713024564674432357634540343884E-3), | |
467 | L(-1.515456784370354374066417703736088291287E-2), | |
468 | L(-1.022454301137286306933217746545237098518E-1), | |
469 | L(-3.373438906472495080504907858424251082240E-1), | |
470 | L(-4.510782522110845697262323973549178453405E-1), | |
471 | L(-1.549000892545288676809660828213589804884E-1), | |
b45e4f52 AJ |
472 | }; |
473 | #define NQ5_8D 10 | |
15089e04 | 474 | static const _Float128 Q5_8D[NQ5_8D + 1] = { |
02bbfb41 PM |
475 | L(3.565550843359501079050699598913828460036E-12), |
476 | L(1.321016015556560621591847454285330528045E-9), | |
477 | L(1.897542728662346479999969679234270605975E-7), | |
478 | L(1.381720283068706710298734234287456219474E-5), | |
479 | L(5.599248147286524662305325795203422873725E-4), | |
480 | L(1.305442352653121436697064782499122164843E-2), | |
481 | L(1.750234079626943298160445750078631894985E-1), | |
482 | L(1.311420542073436520965439883806946678491E0), | |
483 | L(5.162757689856842406744504211089724926650E0), | |
484 | L(9.527760296384704425618556332087850581308E0), | |
485 | L(6.604648207463236667912921642545100248584E0), | |
b45e4f52 AJ |
486 | /* 1.000000000000000000000000000000000000000E0 */ |
487 | }; | |
488 | ||
489 | /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x), | |
490 | Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)), | |
491 | Peak relative error 1.3e-35 | |
492 | 0.1875 <= 1/x <= 0.25 */ | |
493 | #define NQ4_5N 10 | |
15089e04 | 494 | static const _Float128 Q4_5N[NQ4_5N + 1] = { |
02bbfb41 PM |
495 | L(-4.079513568708891749424783046520200903755E-11), |
496 | L(-9.326548104106791766891812583019664893311E-9), | |
497 | L(-8.016795121318423066292906123815687003356E-7), | |
498 | L(-3.372350544043594415609295225664186750995E-5), | |
499 | L(-7.566238665947967882207277686375417983917E-4), | |
500 | L(-9.248861580055565402130441618521591282617E-3), | |
501 | L(-6.033106131055851432267702948850231270338E-2), | |
502 | L(-1.966908754799996793730369265431584303447E-1), | |
503 | L(-2.791062741179964150755788226623462207560E-1), | |
504 | L(-1.255478605849190549914610121863534191666E-1), | |
505 | L(-4.320429862021265463213168186061696944062E-3), | |
b45e4f52 AJ |
506 | }; |
507 | #define NQ4_5D 9 | |
15089e04 | 508 | static const _Float128 Q4_5D[NQ4_5D + 1] = { |
02bbfb41 PM |
509 | L(3.978497042580921479003851216297330701056E-10), |
510 | L(9.203304163828145809278568906420772246666E-8), | |
511 | L(8.059685467088175644915010485174545743798E-6), | |
512 | L(3.490187375993956409171098277561669167446E-4), | |
513 | L(8.189109654456872150100501732073810028829E-3), | |
514 | L(1.072572867311023640958725265762483033769E-1), | |
515 | L(7.790606862409960053675717185714576937994E-1), | |
516 | L(3.016049768232011196434185423512777656328E0), | |
517 | L(5.722963851442769787733717162314477949360E0), | |
518 | L(4.510527838428473279647251350931380867663E0), | |
b45e4f52 AJ |
519 | /* 1.000000000000000000000000000000000000000E0 */ |
520 | }; | |
521 | ||
522 | /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x), | |
523 | Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)), | |
524 | Peak relative error 2.1e-35 | |
525 | 0.25 <= 1/x <= 0.3125 */ | |
526 | #define NQ3r2_4N 9 | |
15089e04 | 527 | static const _Float128 Q3r2_4N[NQ3r2_4N + 1] = { |
02bbfb41 PM |
528 | L(-1.087480809271383885936921889040388133627E-8), |
529 | L(-1.690067828697463740906962973479310170932E-6), | |
530 | L(-9.608064416995105532790745641974762550982E-5), | |
531 | L(-2.594198839156517191858208513873961837410E-3), | |
532 | L(-3.610954144421543968160459863048062977822E-2), | |
533 | L(-2.629866798251843212210482269563961685666E-1), | |
534 | L(-9.709186825881775885917984975685752956660E-1), | |
535 | L(-1.667521829918185121727268867619982417317E0), | |
536 | L(-1.109255082925540057138766105229900943501E0), | |
537 | L(-1.812932453006641348145049323713469043328E-1), | |
b45e4f52 AJ |
538 | }; |
539 | #define NQ3r2_4D 9 | |
15089e04 | 540 | static const _Float128 Q3r2_4D[NQ3r2_4D + 1] = { |
02bbfb41 PM |
541 | L(1.060552717496912381388763753841473407026E-7), |
542 | L(1.676928002024920520786883649102388708024E-5), | |
543 | L(9.803481712245420839301400601140812255737E-4), | |
544 | L(2.765559874262309494758505158089249012930E-2), | |
545 | L(4.117921827792571791298862613287549140706E-1), | |
546 | L(3.323769515244751267093378361930279161413E0), | |
547 | L(1.436602494405814164724810151689705353670E1), | |
548 | L(3.163087869617098638064881410646782408297E1), | |
549 | L(3.198181264977021649489103980298349589419E1), | |
550 | L(1.203649258862068431199471076202897823272E1), | |
b45e4f52 AJ |
551 | /* 1.000000000000000000000000000000000000000E0 */ |
552 | }; | |
553 | ||
554 | /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x), | |
555 | Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)), | |
556 | Peak relative error 1.6e-36 | |
557 | 0.3125 <= 1/x <= 0.375 */ | |
558 | #define NQ2r7_3r2N 9 | |
15089e04 | 559 | static const _Float128 Q2r7_3r2N[NQ2r7_3r2N + 1] = { |
02bbfb41 PM |
560 | L(-1.723405393982209853244278760171643219530E-7), |
561 | L(-2.090508758514655456365709712333460087442E-5), | |
562 | L(-9.140104013370974823232873472192719263019E-4), | |
563 | L(-1.871349499990714843332742160292474780128E-2), | |
564 | L(-1.948930738119938669637865956162512983416E-1), | |
565 | L(-1.048764684978978127908439526343174139788E0), | |
566 | L(-2.827714929925679500237476105843643064698E0), | |
567 | L(-3.508761569156476114276988181329773987314E0), | |
568 | L(-1.669332202790211090973255098624488308989E0), | |
569 | L(-1.930796319299022954013840684651016077770E-1), | |
b45e4f52 AJ |
570 | }; |
571 | #define NQ2r7_3r2D 9 | |
15089e04 | 572 | static const _Float128 Q2r7_3r2D[NQ2r7_3r2D + 1] = { |
02bbfb41 PM |
573 | L(1.680730662300831976234547482334347983474E-6), |
574 | L(2.084241442440551016475972218719621841120E-4), | |
575 | L(9.445316642108367479043541702688736295579E-3), | |
576 | L(2.044637889456631896650179477133252184672E-1), | |
577 | L(2.316091982244297350829522534435350078205E0), | |
578 | L(1.412031891783015085196708811890448488865E1), | |
579 | L(4.583830154673223384837091077279595496149E1), | |
580 | L(7.549520609270909439885998474045974122261E1), | |
581 | L(5.697605832808113367197494052388203310638E1), | |
582 | L(1.601496240876192444526383314589371686234E1), | |
b45e4f52 AJ |
583 | /* 1.000000000000000000000000000000000000000E0 */ |
584 | }; | |
585 | ||
586 | /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x), | |
587 | Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)), | |
588 | Peak relative error 9.5e-36 | |
589 | 0.375 <= 1/x <= 0.4375 */ | |
590 | #define NQ2r3_2r7N 9 | |
15089e04 | 591 | static const _Float128 Q2r3_2r7N[NQ2r3_2r7N + 1] = { |
02bbfb41 PM |
592 | L(-8.603042076329122085722385914954878953775E-7), |
593 | L(-7.701746260451647874214968882605186675720E-5), | |
594 | L(-2.407932004380727587382493696877569654271E-3), | |
595 | L(-3.403434217607634279028110636919987224188E-2), | |
596 | L(-2.348707332185238159192422084985713102877E-1), | |
597 | L(-7.957498841538254916147095255700637463207E-1), | |
598 | L(-1.258469078442635106431098063707934348577E0), | |
599 | L(-8.162415474676345812459353639449971369890E-1), | |
600 | L(-1.581783890269379690141513949609572806898E-1), | |
601 | L(-1.890595651683552228232308756569450822905E-3), | |
b45e4f52 AJ |
602 | }; |
603 | #define NQ2r3_2r7D 8 | |
15089e04 | 604 | static const _Float128 Q2r3_2r7D[NQ2r3_2r7D + 1] = { |
02bbfb41 PM |
605 | L(8.390017524798316921170710533381568175665E-6), |
606 | L(7.738148683730826286477254659973968763659E-4), | |
607 | L(2.541480810958665794368759558791634341779E-2), | |
608 | L(3.878879789711276799058486068562386244873E-1), | |
609 | L(3.003783779325811292142957336802456109333E0), | |
610 | L(1.206480374773322029883039064575464497400E1), | |
611 | L(2.458414064785315978408974662900438351782E1), | |
612 | L(2.367237826273668567199042088835448715228E1), | |
613 | L(9.231451197519171090875569102116321676763E0), | |
b45e4f52 AJ |
614 | /* 1.000000000000000000000000000000000000000E0 */ |
615 | }; | |
616 | ||
617 | /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x), | |
618 | Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)), | |
619 | Peak relative error 1.4e-36 | |
620 | 0.4375 <= 1/x <= 0.5 */ | |
621 | #define NQ2_2r3N 9 | |
15089e04 | 622 | static const _Float128 Q2_2r3N[NQ2_2r3N + 1] = { |
02bbfb41 PM |
623 | L(-5.552507516089087822166822364590806076174E-6), |
624 | L(-4.135067659799500521040944087433752970297E-4), | |
625 | L(-1.059928728869218962607068840646564457980E-2), | |
626 | L(-1.212070036005832342565792241385459023801E-1), | |
627 | L(-6.688350110633603958684302153362735625156E-1), | |
628 | L(-1.793587878197360221340277951304429821582E0), | |
629 | L(-2.225407682237197485644647380483725045326E0), | |
630 | L(-1.123402135458940189438898496348239744403E0), | |
631 | L(-1.679187241566347077204805190763597299805E-1), | |
632 | L(-1.458550613639093752909985189067233504148E-3), | |
b45e4f52 AJ |
633 | }; |
634 | #define NQ2_2r3D 8 | |
15089e04 | 635 | static const _Float128 Q2_2r3D[NQ2_2r3D + 1] = { |
02bbfb41 PM |
636 | L(5.415024336507980465169023996403597916115E-5), |
637 | L(4.179246497380453022046357404266022870788E-3), | |
638 | L(1.136306384261959483095442402929502368598E-1), | |
639 | L(1.422640343719842213484515445393284072830E0), | |
640 | L(8.968786703393158374728850922289204805764E0), | |
641 | L(2.914542473339246127533384118781216495934E1), | |
642 | L(4.781605421020380669870197378210457054685E1), | |
643 | L(3.693865837171883152382820584714795072937E1), | |
644 | L(1.153220502744204904763115556224395893076E1), | |
b45e4f52 AJ |
645 | /* 1.000000000000000000000000000000000000000E0 */ |
646 | }; | |
647 | ||
648 | ||
649 | /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */ | |
650 | ||
15089e04 PM |
651 | static _Float128 |
652 | neval (_Float128 x, const _Float128 *p, int n) | |
b45e4f52 | 653 | { |
15089e04 | 654 | _Float128 y; |
b45e4f52 AJ |
655 | |
656 | p += n; | |
657 | y = *p--; | |
658 | do | |
659 | { | |
660 | y = y * x + *p--; | |
661 | } | |
662 | while (--n > 0); | |
663 | return y; | |
664 | } | |
665 | ||
666 | ||
667 | /* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */ | |
668 | ||
15089e04 PM |
669 | static _Float128 |
670 | deval (_Float128 x, const _Float128 *p, int n) | |
b45e4f52 | 671 | { |
15089e04 | 672 | _Float128 y; |
b45e4f52 AJ |
673 | |
674 | p += n; | |
675 | y = x + *p--; | |
676 | do | |
677 | { | |
678 | y = y * x + *p--; | |
679 | } | |
680 | while (--n > 0); | |
681 | return y; | |
682 | } | |
683 | ||
684 | ||
685 | /* Bessel function of the first kind, order one. */ | |
686 | ||
15089e04 PM |
687 | _Float128 |
688 | __ieee754_j1l (_Float128 x) | |
b45e4f52 | 689 | { |
15089e04 | 690 | _Float128 xx, xinv, z, p, q, c, s, cc, ss; |
b45e4f52 | 691 | |
d81f90cc | 692 | if (! isfinite (x)) |
b45e4f52 AJ |
693 | { |
694 | if (x != x) | |
d73e7bdb | 695 | return x + x; |
b45e4f52 | 696 | else |
02bbfb41 | 697 | return 0; |
b45e4f52 | 698 | } |
02bbfb41 | 699 | if (x == 0) |
b45e4f52 AJ |
700 | return x; |
701 | xx = fabsl (x); | |
02bbfb41 | 702 | if (xx <= L(0x1p-58)) |
4aa10d01 | 703 | { |
02bbfb41 | 704 | _Float128 ret = x * L(0.5); |
d96164c3 | 705 | math_check_force_underflow (ret); |
c643db87 JM |
706 | if (ret == 0) |
707 | __set_errno (ERANGE); | |
4aa10d01 JM |
708 | return ret; |
709 | } | |
02bbfb41 | 710 | if (xx <= 2) |
b45e4f52 AJ |
711 | { |
712 | /* 0 <= x <= 2 */ | |
713 | z = xx * xx; | |
714 | p = xx * z * neval (z, J0_2N, NJ0_2N) / deval (z, J0_2D, NJ0_2D); | |
02bbfb41 | 715 | p += L(0.5) * xx; |
b45e4f52 AJ |
716 | if (x < 0) |
717 | p = -p; | |
718 | return p; | |
719 | } | |
720 | ||
2a185d32 JM |
721 | /* X = x - 3 pi/4 |
722 | cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4) | |
723 | = 1/sqrt(2) * (-cos(x) + sin(x)) | |
724 | sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4) | |
725 | = -1/sqrt(2) * (sin(x) + cos(x)) | |
726 | cf. Fdlibm. */ | |
727 | __sincosl (xx, &s, &c); | |
728 | ss = -s - c; | |
729 | cc = s - c; | |
02bbfb41 | 730 | if (xx <= LDBL_MAX / 2) |
98c48fe5 JM |
731 | { |
732 | z = __cosl (xx + xx); | |
733 | if ((s * c) > 0) | |
734 | cc = z / ss; | |
735 | else | |
736 | ss = z / cc; | |
737 | } | |
2a185d32 | 738 | |
02bbfb41 | 739 | if (xx > L(0x1p256)) |
2a185d32 | 740 | { |
f67a8147 | 741 | z = ONEOSQPI * cc / sqrtl (xx); |
2a185d32 JM |
742 | if (x < 0) |
743 | z = -z; | |
744 | return z; | |
745 | } | |
746 | ||
02bbfb41 | 747 | xinv = 1 / xx; |
b45e4f52 AJ |
748 | z = xinv * xinv; |
749 | if (xinv <= 0.25) | |
750 | { | |
751 | if (xinv <= 0.125) | |
752 | { | |
753 | if (xinv <= 0.0625) | |
754 | { | |
755 | p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID); | |
756 | q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID); | |
757 | } | |
758 | else | |
759 | { | |
760 | p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D); | |
761 | q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D); | |
762 | } | |
763 | } | |
764 | else if (xinv <= 0.1875) | |
765 | { | |
766 | p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D); | |
767 | q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D); | |
768 | } | |
769 | else | |
770 | { | |
771 | p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D); | |
772 | q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D); | |
773 | } | |
774 | } /* .25 */ | |
775 | else /* if (xinv <= 0.5) */ | |
776 | { | |
777 | if (xinv <= 0.375) | |
778 | { | |
779 | if (xinv <= 0.3125) | |
780 | { | |
781 | p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D); | |
782 | q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D); | |
783 | } | |
784 | else | |
785 | { | |
786 | p = neval (z, P2r7_3r2N, NP2r7_3r2N) | |
787 | / deval (z, P2r7_3r2D, NP2r7_3r2D); | |
788 | q = neval (z, Q2r7_3r2N, NQ2r7_3r2N) | |
789 | / deval (z, Q2r7_3r2D, NQ2r7_3r2D); | |
790 | } | |
791 | } | |
792 | else if (xinv <= 0.4375) | |
793 | { | |
794 | p = neval (z, P2r3_2r7N, NP2r3_2r7N) | |
795 | / deval (z, P2r3_2r7D, NP2r3_2r7D); | |
796 | q = neval (z, Q2r3_2r7N, NQ2r3_2r7N) | |
797 | / deval (z, Q2r3_2r7D, NQ2r3_2r7D); | |
798 | } | |
799 | else | |
800 | { | |
801 | p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D); | |
802 | q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D); | |
803 | } | |
804 | } | |
02bbfb41 | 805 | p = 1 + z * p; |
b45e4f52 | 806 | q = z * q; |
02bbfb41 | 807 | q = q * xinv + L(0.375) * xinv; |
f67a8147 | 808 | z = ONEOSQPI * (p * cc - q * ss) / sqrtl (xx); |
b45e4f52 AJ |
809 | if (x < 0) |
810 | z = -z; | |
811 | return z; | |
812 | } | |
0ac5ae23 | 813 | strong_alias (__ieee754_j1l, __j1l_finite) |
b45e4f52 AJ |
814 | |
815 | ||
816 | /* Y1(x) = 2/pi * (log(x) * J1(x) - 1/x) + x R(x^2) | |
817 | Peak relative error 6.2e-38 | |
818 | 0 <= x <= 2 */ | |
819 | #define NY0_2N 7 | |
53994f12 | 820 | static const _Float128 Y0_2N[NY0_2N + 1] = { |
02bbfb41 PM |
821 | L(-6.804415404830253804408698161694720833249E19), |
822 | L(1.805450517967019908027153056150465849237E19), | |
823 | L(-8.065747497063694098810419456383006737312E17), | |
824 | L(1.401336667383028259295830955439028236299E16), | |
825 | L(-1.171654432898137585000399489686629680230E14), | |
826 | L(5.061267920943853732895341125243428129150E11), | |
827 | L(-1.096677850566094204586208610960870217970E9), | |
828 | L(9.541172044989995856117187515882879304461E5), | |
b45e4f52 AJ |
829 | }; |
830 | #define NY0_2D 7 | |
53994f12 | 831 | static const _Float128 Y0_2D[NY0_2D + 1] = { |
02bbfb41 PM |
832 | L(3.470629591820267059538637461549677594549E20), |
833 | L(4.120796439009916326855848107545425217219E18), | |
834 | L(2.477653371652018249749350657387030814542E16), | |
835 | L(9.954678543353888958177169349272167762797E13), | |
836 | L(2.957927997613630118216218290262851197754E11), | |
837 | L(6.748421382188864486018861197614025972118E8), | |
838 | L(1.173453425218010888004562071020305709319E6), | |
839 | L(1.450335662961034949894009554536003377187E3), | |
b45e4f52 AJ |
840 | /* 1.000000000000000000000000000000000000000E0 */ |
841 | }; | |
842 | ||
843 | ||
844 | /* Bessel function of the second kind, order one. */ | |
845 | ||
15089e04 PM |
846 | _Float128 |
847 | __ieee754_y1l (_Float128 x) | |
b45e4f52 | 848 | { |
15089e04 | 849 | _Float128 xx, xinv, z, p, q, c, s, cc, ss; |
b45e4f52 | 850 | |
d81f90cc | 851 | if (! isfinite (x)) |
b987917e | 852 | return 1 / (x + x * x); |
02bbfb41 | 853 | if (x <= 0) |
b45e4f52 | 854 | { |
02bbfb41 | 855 | if (x < 0) |
1f510b3f | 856 | return (zero / (zero * x)); |
4918e5f4 | 857 | return -1 / zero; /* -inf and divide by zero exception. */ |
b45e4f52 AJ |
858 | } |
859 | xx = fabsl (x); | |
05b227bd | 860 | if (xx <= 0x1p-114) |
46a3d3c7 JM |
861 | { |
862 | z = -TWOOPI / x; | |
d81f90cc | 863 | if (isinf (z)) |
46a3d3c7 JM |
864 | __set_errno (ERANGE); |
865 | return z; | |
866 | } | |
02bbfb41 | 867 | if (xx <= 2) |
b45e4f52 AJ |
868 | { |
869 | /* 0 <= x <= 2 */ | |
8ea587db | 870 | SET_RESTORE_ROUNDL (FE_TONEAREST); |
b45e4f52 AJ |
871 | z = xx * xx; |
872 | p = xx * neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D); | |
873 | p = -TWOOPI / xx + p; | |
246ec411 | 874 | p = TWOOPI * __ieee754_logl (x) * __ieee754_j1l (x) + p; |
b45e4f52 AJ |
875 | return p; |
876 | } | |
877 | ||
2a185d32 JM |
878 | /* X = x - 3 pi/4 |
879 | cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4) | |
880 | = 1/sqrt(2) * (-cos(x) + sin(x)) | |
881 | sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4) | |
882 | = -1/sqrt(2) * (sin(x) + cos(x)) | |
883 | cf. Fdlibm. */ | |
884 | __sincosl (xx, &s, &c); | |
885 | ss = -s - c; | |
886 | cc = s - c; | |
02bbfb41 | 887 | if (xx <= LDBL_MAX / 2) |
98c48fe5 JM |
888 | { |
889 | z = __cosl (xx + xx); | |
890 | if ((s * c) > 0) | |
891 | cc = z / ss; | |
892 | else | |
893 | ss = z / cc; | |
894 | } | |
2a185d32 | 895 | |
02bbfb41 | 896 | if (xx > L(0x1p256)) |
f67a8147 | 897 | return ONEOSQPI * ss / sqrtl (xx); |
2a185d32 | 898 | |
02bbfb41 | 899 | xinv = 1 / xx; |
b45e4f52 AJ |
900 | z = xinv * xinv; |
901 | if (xinv <= 0.25) | |
902 | { | |
903 | if (xinv <= 0.125) | |
904 | { | |
905 | if (xinv <= 0.0625) | |
906 | { | |
907 | p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID); | |
908 | q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID); | |
909 | } | |
910 | else | |
911 | { | |
912 | p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D); | |
913 | q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D); | |
914 | } | |
915 | } | |
916 | else if (xinv <= 0.1875) | |
917 | { | |
918 | p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D); | |
919 | q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D); | |
920 | } | |
921 | else | |
922 | { | |
923 | p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D); | |
924 | q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D); | |
925 | } | |
926 | } /* .25 */ | |
927 | else /* if (xinv <= 0.5) */ | |
928 | { | |
929 | if (xinv <= 0.375) | |
930 | { | |
931 | if (xinv <= 0.3125) | |
932 | { | |
933 | p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D); | |
934 | q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D); | |
935 | } | |
936 | else | |
937 | { | |
938 | p = neval (z, P2r7_3r2N, NP2r7_3r2N) | |
939 | / deval (z, P2r7_3r2D, NP2r7_3r2D); | |
940 | q = neval (z, Q2r7_3r2N, NQ2r7_3r2N) | |
941 | / deval (z, Q2r7_3r2D, NQ2r7_3r2D); | |
942 | } | |
943 | } | |
944 | else if (xinv <= 0.4375) | |
945 | { | |
946 | p = neval (z, P2r3_2r7N, NP2r3_2r7N) | |
947 | / deval (z, P2r3_2r7D, NP2r3_2r7D); | |
948 | q = neval (z, Q2r3_2r7N, NQ2r3_2r7N) | |
949 | / deval (z, Q2r3_2r7D, NQ2r3_2r7D); | |
950 | } | |
951 | else | |
952 | { | |
953 | p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D); | |
954 | q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D); | |
955 | } | |
956 | } | |
02bbfb41 | 957 | p = 1 + z * p; |
b45e4f52 | 958 | q = z * q; |
02bbfb41 | 959 | q = q * xinv + L(0.375) * xinv; |
f67a8147 | 960 | z = ONEOSQPI * (p * ss + q * cc) / sqrtl (xx); |
b45e4f52 AJ |
961 | return z; |
962 | } | |
0ac5ae23 | 963 | strong_alias (__ieee754_y1l, __y1l_finite) |