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