3 * Bessel function of order one
9 * long double x, y, j1l();
17 * Returns Bessel function of first kind, order one of the argument.
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)
23 * The second interval is further partitioned into eight equal segments
25 * J1(x) = sqrt(2/(pi x)) (P1(x) cos(X) - Q1(x) sin(X)),
28 * and the auxiliary functions are given by
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)
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)).
41 * arithmetic domain # trials peak rms
42 * IEEE 0, 30 100000 2.8e-34 2.7e-35
49 * Bessel function of the second kind, order one
63 * Returns Bessel function of the second kind, of order
64 * one, of the argument.
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)),
75 * Absolute error, when y0(x) < 1; else relative error:
77 * arithmetic domain # trials peak rms
78 * IEEE 0, 30 100000 2.7e-34 2.9e-35
82 /* Copyright 2001 by Stephen L. Moshier (moshier@na-net.onrl.gov).
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.
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.
94 You should have received a copy of the GNU Lesser General Public
95 License along with this library; if not, see
96 <https://www.gnu.org/licenses/>. */
100 #include <math_private.h>
101 #include <fenv_private.h>
102 #include <math-underflow.h>
104 #include <libm-alias-finite.h>
107 static const _Float128 ONEOSQPI
= L(5.6418958354775628694807945156077258584405E-1);
109 static const _Float128 TWOOPI
= L(6.3661977236758134307553505349005744813784E-1);
110 static const _Float128 zero
= 0;
112 /* J1(x) = .5x + x x^2 R(x^2)
113 Peak relative error 1.9e-35
116 static const _Float128 J0_2N
[NJ0_2N
+ 1] = {
117 L(-5.943799577386942855938508697619735179660E16
),
118 L(1.812087021305009192259946997014044074711E15
),
119 L(-2.761698314264509665075127515729146460895E13
),
120 L(2.091089497823600978949389109350658815972E11
),
121 L(-8.546413231387036372945453565654130054307E8
),
122 L(1.797229225249742247475464052741320612261E6
),
123 L(-1.559552840946694171346552770008812083969E3
)
126 static const _Float128 J0_2D
[NJ0_2D
+ 1] = {
127 L(9.510079323819108569501613916191477479397E17
),
128 L(1.063193817503280529676423936545854693915E16
),
129 L(5.934143516050192600795972192791775226920E13
),
130 L(2.168000911950620999091479265214368352883E11
),
131 L(5.673775894803172808323058205986256928794E8
),
132 L(1.080329960080981204840966206372671147224E6
),
133 L(1.411951256636576283942477881535283304912E3
),
134 /* 1.000000000000000000000000000000000000000E0L */
137 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
139 Peak relative error 3.6e-36 */
141 static const _Float128 P16_IN
[NP16_IN
+ 1] = {
142 L(5.143674369359646114999545149085139822905E-16),
143 L(4.836645664124562546056389268546233577376E-13),
144 L(1.730945562285804805325011561498453013673E-10),
145 L(3.047976856147077889834905908605310585810E-8),
146 L(2.855227609107969710407464739188141162386E-6),
147 L(1.439362407936705484122143713643023998457E-4),
148 L(3.774489768532936551500999699815873422073E-3),
149 L(4.723962172984642566142399678920790598426E-2),
150 L(2.359289678988743939925017240478818248735E-1),
151 L(3.032580002220628812728954785118117124520E-1),
154 static const _Float128 P16_ID
[NP16_ID
+ 1] = {
155 L(4.389268795186898018132945193912677177553E-15),
156 L(4.132671824807454334388868363256830961655E-12),
157 L(1.482133328179508835835963635130894413136E-9),
158 L(2.618941412861122118906353737117067376236E-7),
159 L(2.467854246740858470815714426201888034270E-5),
160 L(1.257192927368839847825938545925340230490E-3),
161 L(3.362739031941574274949719324644120720341E-2),
162 L(4.384458231338934105875343439265370178858E-1),
163 L(2.412830809841095249170909628197264854651E0
),
164 L(4.176078204111348059102962617368214856874E0
),
165 /* 1.000000000000000000000000000000000000000E0 */
168 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
169 0.0625 <= 1/x <= 0.125
170 Peak relative error 1.9e-36 */
172 static const _Float128 P8_16N
[NP8_16N
+ 1] = {
173 L(2.984612480763362345647303274082071598135E-16),
174 L(1.923651877544126103941232173085475682334E-13),
175 L(4.881258879388869396043760693256024307743E-11),
176 L(6.368866572475045408480898921866869811889E-9),
177 L(4.684818344104910450523906967821090796737E-7),
178 L(2.005177298271593587095982211091300382796E-5),
179 L(4.979808067163957634120681477207147536182E-4),
180 L(6.946005761642579085284689047091173581127E-3),
181 L(5.074601112955765012750207555985299026204E-2),
182 L(1.698599455896180893191766195194231825379E-1),
183 L(1.957536905259237627737222775573623779638E-1),
184 L(2.991314703282528370270179989044994319374E-2),
187 static const _Float128 P8_16D
[NP8_16D
+ 1] = {
188 L(2.546869316918069202079580939942463010937E-15),
189 L(1.644650111942455804019788382157745229955E-12),
190 L(4.185430770291694079925607420808011147173E-10),
191 L(5.485331966975218025368698195861074143153E-8),
192 L(4.062884421686912042335466327098932678905E-6),
193 L(1.758139661060905948870523641319556816772E-4),
194 L(4.445143889306356207566032244985607493096E-3),
195 L(6.391901016293512632765621532571159071158E-2),
196 L(4.933040207519900471177016015718145795434E-1),
197 L(1.839144086168947712971630337250761842976E0
),
198 L(2.715120873995490920415616716916149586579E0
),
199 /* 1.000000000000000000000000000000000000000E0 */
202 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
203 0.125 <= 1/x <= 0.1875
204 Peak relative error 1.3e-36 */
206 static const _Float128 P5_8N
[NP5_8N
+ 1] = {
207 L(2.837678373978003452653763806968237227234E-12),
208 L(9.726641165590364928442128579282742354806E-10),
209 L(1.284408003604131382028112171490633956539E-7),
210 L(8.524624695868291291250573339272194285008E-6),
211 L(3.111516908953172249853673787748841282846E-4),
212 L(6.423175156126364104172801983096596409176E-3),
213 L(7.430220589989104581004416356260692450652E-2),
214 L(4.608315409833682489016656279567605536619E-1),
215 L(1.396870223510964882676225042258855977512E0
),
216 L(1.718500293904122365894630460672081526236E0
),
217 L(5.465927698800862172307352821870223855365E-1)
220 static const _Float128 P5_8D
[NP5_8D
+ 1] = {
221 L(2.421485545794616609951168511612060482715E-11),
222 L(8.329862750896452929030058039752327232310E-9),
223 L(1.106137992233383429630592081375289010720E-6),
224 L(7.405786153760681090127497796448503306939E-5),
225 L(2.740364785433195322492093333127633465227E-3),
226 L(5.781246470403095224872243564165254652198E-2),
227 L(6.927711353039742469918754111511109983546E-1),
228 L(4.558679283460430281188304515922826156690E0
),
229 L(1.534468499844879487013168065728837900009E1
),
230 L(2.313927430889218597919624843161569422745E1
),
231 L(1.194506341319498844336768473218382828637E1
),
232 /* 1.000000000000000000000000000000000000000E0 */
235 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
236 Peak relative error 1.4e-36
237 0.1875 <= 1/x <= 0.25 */
239 static const _Float128 P4_5N
[NP4_5N
+ 1] = {
240 L(1.846029078268368685834261260420933914621E-10),
241 L(3.916295939611376119377869680335444207768E-8),
242 L(3.122158792018920627984597530935323997312E-6),
243 L(1.218073444893078303994045653603392272450E-4),
244 L(2.536420827983485448140477159977981844883E-3),
245 L(2.883011322006690823959367922241169171315E-2),
246 L(1.755255190734902907438042414495469810830E-1),
247 L(5.379317079922628599870898285488723736599E-1),
248 L(7.284904050194300773890303361501726561938E-1),
249 L(3.270110346613085348094396323925000362813E-1),
250 L(1.804473805689725610052078464951722064757E-2),
253 static const _Float128 P4_5D
[NP4_5D
+ 1] = {
254 L(1.575278146806816970152174364308980863569E-9),
255 L(3.361289173657099516191331123405675054321E-7),
256 L(2.704692281550877810424745289838790693708E-5),
257 L(1.070854930483999749316546199273521063543E-3),
258 L(2.282373093495295842598097265627962125411E-2),
259 L(2.692025460665354148328762368240343249830E-1),
260 L(1.739892942593664447220951225734811133759E0
),
261 L(5.890727576752230385342377570386657229324E0
),
262 L(9.517442287057841500750256954117735128153E0
),
263 L(6.100616353935338240775363403030137736013E0
),
264 /* 1.000000000000000000000000000000000000000E0 */
267 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
268 Peak relative error 3.0e-36
269 0.25 <= 1/x <= 0.3125 */
271 static const _Float128 P3r2_4N
[NP3r2_4N
+ 1] = {
272 L(8.240803130988044478595580300846665863782E-8),
273 L(1.179418958381961224222969866406483744580E-5),
274 L(6.179787320956386624336959112503824397755E-4),
275 L(1.540270833608687596420595830747166658383E-2),
276 L(1.983904219491512618376375619598837355076E-1),
277 L(1.341465722692038870390470651608301155565E0
),
278 L(4.617865326696612898792238245990854646057E0
),
279 L(7.435574801812346424460233180412308000587E0
),
280 L(4.671327027414635292514599201278557680420E0
),
281 L(7.299530852495776936690976966995187714739E-1),
284 static const _Float128 P3r2_4D
[NP3r2_4D
+ 1] = {
285 L(7.032152009675729604487575753279187576521E-7),
286 L(1.015090352324577615777511269928856742848E-4),
287 L(5.394262184808448484302067955186308730620E-3),
288 L(1.375291438480256110455809354836988584325E-1),
289 L(1.836247144461106304788160919310404376670E0
),
290 L(1.314378564254376655001094503090935880349E1
),
291 L(4.957184590465712006934452500894672343488E1
),
292 L(9.287394244300647738855415178790263465398E1
),
293 L(7.652563275535900609085229286020552768399E1
),
294 L(2.147042473003074533150718117770093209096E1
),
295 /* 1.000000000000000000000000000000000000000E0 */
298 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
299 Peak relative error 1.0e-35
300 0.3125 <= 1/x <= 0.375 */
302 static const _Float128 P2r7_3r2N
[NP2r7_3r2N
+ 1] = {
303 L(4.599033469240421554219816935160627085991E-7),
304 L(4.665724440345003914596647144630893997284E-5),
305 L(1.684348845667764271596142716944374892756E-3),
306 L(2.802446446884455707845985913454440176223E-2),
307 L(2.321937586453963310008279956042545173930E-1),
308 L(9.640277413988055668692438709376437553804E-1),
309 L(1.911021064710270904508663334033003246028E0
),
310 L(1.600811610164341450262992138893970224971E0
),
311 L(4.266299218652587901171386591543457861138E-1),
312 L(1.316470424456061252962568223251247207325E-2),
315 static const _Float128 P2r7_3r2D
[NP2r7_3r2D
+ 1] = {
316 L(3.924508608545520758883457108453520099610E-6),
317 L(4.029707889408829273226495756222078039823E-4),
318 L(1.484629715787703260797886463307469600219E-2),
319 L(2.553136379967180865331706538897231588685E-1),
320 L(2.229457223891676394409880026887106228740E0
),
321 L(1.005708903856384091956550845198392117318E1
),
322 L(2.277082659664386953166629360352385889558E1
),
323 L(2.384726835193630788249826630376533988245E1
),
324 L(9.700989749041320895890113781610939632410E0
),
325 /* 1.000000000000000000000000000000000000000E0 */
328 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
329 Peak relative error 1.7e-36
330 0.3125 <= 1/x <= 0.4375 */
332 static const _Float128 P2r3_2r7N
[NP2r3_2r7N
+ 1] = {
333 L(3.916766777108274628543759603786857387402E-6),
334 L(3.212176636756546217390661984304645137013E-4),
335 L(9.255768488524816445220126081207248947118E-3),
336 L(1.214853146369078277453080641911700735354E-1),
337 L(7.855163309847214136198449861311404633665E-1),
338 L(2.520058073282978403655488662066019816540E0
),
339 L(3.825136484837545257209234285382183711466E0
),
340 L(2.432569427554248006229715163865569506873E0
),
341 L(4.877934835018231178495030117729800489743E-1),
342 L(1.109902737860249670981355149101343427885E-2),
345 static const _Float128 P2r3_2r7D
[NP2r3_2r7D
+ 1] = {
346 L(3.342307880794065640312646341190547184461E-5),
347 L(2.782182891138893201544978009012096558265E-3),
348 L(8.221304931614200702142049236141249929207E-2),
349 L(1.123728246291165812392918571987858010949E0
),
350 L(7.740482453652715577233858317133423434590E0
),
351 L(2.737624677567945952953322566311201919139E1
),
352 L(4.837181477096062403118304137851260715475E1
),
353 L(3.941098643468580791437772701093795299274E1
),
354 L(1.245821247166544627558323920382547533630E1
),
355 /* 1.000000000000000000000000000000000000000E0 */
358 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
359 Peak relative error 1.7e-35
360 0.4375 <= 1/x <= 0.5 */
362 static const _Float128 P2_2r3N
[NP2_2r3N
+ 1] = {
363 L(3.397930802851248553545191160608731940751E-4),
364 L(2.104020902735482418784312825637833698217E-2),
365 L(4.442291771608095963935342749477836181939E-1),
366 L(4.131797328716583282869183304291833754967E0
),
367 L(1.819920169779026500146134832455189917589E1
),
368 L(3.781779616522937565300309684282401791291E1
),
369 L(3.459605449728864218972931220783543410347E1
),
370 L(1.173594248397603882049066603238568316561E1
),
371 L(9.455702270242780642835086549285560316461E-1),
374 static const _Float128 P2_2r3D
[NP2_2r3D
+ 1] = {
375 L(2.899568897241432883079888249845707400614E-3),
376 L(1.831107138190848460767699919531132426356E-1),
377 L(3.999350044057883839080258832758908825165E0
),
378 L(3.929041535867957938340569419874195303712E1
),
379 L(1.884245613422523323068802689915538908291E2
),
380 L(4.461469948819229734353852978424629815929E2
),
381 L(5.004998753999796821224085972610636347903E2
),
382 L(2.386342520092608513170837883757163414100E2
),
383 L(3.791322528149347975999851588922424189957E1
),
384 /* 1.000000000000000000000000000000000000000E0 */
387 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
388 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
389 Peak relative error 8.0e-36
392 static const _Float128 Q16_IN
[NQ16_IN
+ 1] = {
393 L(-3.917420835712508001321875734030357393421E-18),
394 L(-4.440311387483014485304387406538069930457E-15),
395 L(-1.951635424076926487780929645954007139616E-12),
396 L(-4.318256438421012555040546775651612810513E-10),
397 L(-5.231244131926180765270446557146989238020E-8),
398 L(-3.540072702902043752460711989234732357653E-6),
399 L(-1.311017536555269966928228052917534882984E-4),
400 L(-2.495184669674631806622008769674827575088E-3),
401 L(-2.141868222987209028118086708697998506716E-2),
402 L(-6.184031415202148901863605871197272650090E-2),
403 L(-1.922298704033332356899546792898156493887E-2),
406 static const _Float128 Q16_ID
[NQ16_ID
+ 1] = {
407 L(3.820418034066293517479619763498400162314E-17),
408 L(4.340702810799239909648911373329149354911E-14),
409 L(1.914985356383416140706179933075303538524E-11),
410 L(4.262333682610888819476498617261895474330E-9),
411 L(5.213481314722233980346462747902942182792E-7),
412 L(3.585741697694069399299005316809954590558E-5),
413 L(1.366513429642842006385029778105539457546E-3),
414 L(2.745282599850704662726337474371355160594E-2),
415 L(2.637644521611867647651200098449903330074E-1),
416 L(1.006953426110765984590782655598680488746E0
),
417 /* 1.000000000000000000000000000000000000000E0 */
420 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
421 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
422 Peak relative error 1.9e-36
423 0.0625 <= 1/x <= 0.125 */
425 static const _Float128 Q8_16N
[NQ8_16N
+ 1] = {
426 L(-2.028630366670228670781362543615221542291E-17),
427 L(-1.519634620380959966438130374006858864624E-14),
428 L(-4.540596528116104986388796594639405114524E-12),
429 L(-7.085151756671466559280490913558388648274E-10),
430 L(-6.351062671323970823761883833531546885452E-8),
431 L(-3.390817171111032905297982523519503522491E-6),
432 L(-1.082340897018886970282138836861233213972E-4),
433 L(-2.020120801187226444822977006648252379508E-3),
434 L(-2.093169910981725694937457070649605557555E-2),
435 L(-1.092176538874275712359269481414448063393E-1),
436 L(-2.374790947854765809203590474789108718733E-1),
437 L(-1.365364204556573800719985118029601401323E-1),
440 static const _Float128 Q8_16D
[NQ8_16D
+ 1] = {
441 L(1.978397614733632533581207058069628242280E-16),
442 L(1.487361156806202736877009608336766720560E-13),
443 L(4.468041406888412086042576067133365913456E-11),
444 L(7.027822074821007443672290507210594648877E-9),
445 L(6.375740580686101224127290062867976007374E-7),
446 L(3.466887658320002225888644977076410421940E-5),
447 L(1.138625640905289601186353909213719596986E-3),
448 L(2.224470799470414663443449818235008486439E-2),
449 L(2.487052928527244907490589787691478482358E-1),
450 L(1.483927406564349124649083853892380899217E0
),
451 L(4.182773513276056975777258788903489507705E0
),
452 L(4.419665392573449746043880892524360870944E0
),
453 /* 1.000000000000000000000000000000000000000E0 */
456 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
457 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
458 Peak relative error 1.5e-35
459 0.125 <= 1/x <= 0.1875 */
461 static const _Float128 Q5_8N
[NQ5_8N
+ 1] = {
462 L(-3.656082407740970534915918390488336879763E-13),
463 L(-1.344660308497244804752334556734121771023E-10),
464 L(-1.909765035234071738548629788698150760791E-8),
465 L(-1.366668038160120210269389551283666716453E-6),
466 L(-5.392327355984269366895210704976314135683E-5),
467 L(-1.206268245713024564674432357634540343884E-3),
468 L(-1.515456784370354374066417703736088291287E-2),
469 L(-1.022454301137286306933217746545237098518E-1),
470 L(-3.373438906472495080504907858424251082240E-1),
471 L(-4.510782522110845697262323973549178453405E-1),
472 L(-1.549000892545288676809660828213589804884E-1),
475 static const _Float128 Q5_8D
[NQ5_8D
+ 1] = {
476 L(3.565550843359501079050699598913828460036E-12),
477 L(1.321016015556560621591847454285330528045E-9),
478 L(1.897542728662346479999969679234270605975E-7),
479 L(1.381720283068706710298734234287456219474E-5),
480 L(5.599248147286524662305325795203422873725E-4),
481 L(1.305442352653121436697064782499122164843E-2),
482 L(1.750234079626943298160445750078631894985E-1),
483 L(1.311420542073436520965439883806946678491E0
),
484 L(5.162757689856842406744504211089724926650E0
),
485 L(9.527760296384704425618556332087850581308E0
),
486 L(6.604648207463236667912921642545100248584E0
),
487 /* 1.000000000000000000000000000000000000000E0 */
490 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
491 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
492 Peak relative error 1.3e-35
493 0.1875 <= 1/x <= 0.25 */
495 static const _Float128 Q4_5N
[NQ4_5N
+ 1] = {
496 L(-4.079513568708891749424783046520200903755E-11),
497 L(-9.326548104106791766891812583019664893311E-9),
498 L(-8.016795121318423066292906123815687003356E-7),
499 L(-3.372350544043594415609295225664186750995E-5),
500 L(-7.566238665947967882207277686375417983917E-4),
501 L(-9.248861580055565402130441618521591282617E-3),
502 L(-6.033106131055851432267702948850231270338E-2),
503 L(-1.966908754799996793730369265431584303447E-1),
504 L(-2.791062741179964150755788226623462207560E-1),
505 L(-1.255478605849190549914610121863534191666E-1),
506 L(-4.320429862021265463213168186061696944062E-3),
509 static const _Float128 Q4_5D
[NQ4_5D
+ 1] = {
510 L(3.978497042580921479003851216297330701056E-10),
511 L(9.203304163828145809278568906420772246666E-8),
512 L(8.059685467088175644915010485174545743798E-6),
513 L(3.490187375993956409171098277561669167446E-4),
514 L(8.189109654456872150100501732073810028829E-3),
515 L(1.072572867311023640958725265762483033769E-1),
516 L(7.790606862409960053675717185714576937994E-1),
517 L(3.016049768232011196434185423512777656328E0
),
518 L(5.722963851442769787733717162314477949360E0
),
519 L(4.510527838428473279647251350931380867663E0
),
520 /* 1.000000000000000000000000000000000000000E0 */
523 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
524 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
525 Peak relative error 2.1e-35
526 0.25 <= 1/x <= 0.3125 */
528 static const _Float128 Q3r2_4N
[NQ3r2_4N
+ 1] = {
529 L(-1.087480809271383885936921889040388133627E-8),
530 L(-1.690067828697463740906962973479310170932E-6),
531 L(-9.608064416995105532790745641974762550982E-5),
532 L(-2.594198839156517191858208513873961837410E-3),
533 L(-3.610954144421543968160459863048062977822E-2),
534 L(-2.629866798251843212210482269563961685666E-1),
535 L(-9.709186825881775885917984975685752956660E-1),
536 L(-1.667521829918185121727268867619982417317E0
),
537 L(-1.109255082925540057138766105229900943501E0
),
538 L(-1.812932453006641348145049323713469043328E-1),
541 static const _Float128 Q3r2_4D
[NQ3r2_4D
+ 1] = {
542 L(1.060552717496912381388763753841473407026E-7),
543 L(1.676928002024920520786883649102388708024E-5),
544 L(9.803481712245420839301400601140812255737E-4),
545 L(2.765559874262309494758505158089249012930E-2),
546 L(4.117921827792571791298862613287549140706E-1),
547 L(3.323769515244751267093378361930279161413E0
),
548 L(1.436602494405814164724810151689705353670E1
),
549 L(3.163087869617098638064881410646782408297E1
),
550 L(3.198181264977021649489103980298349589419E1
),
551 L(1.203649258862068431199471076202897823272E1
),
552 /* 1.000000000000000000000000000000000000000E0 */
555 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
556 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
557 Peak relative error 1.6e-36
558 0.3125 <= 1/x <= 0.375 */
560 static const _Float128 Q2r7_3r2N
[NQ2r7_3r2N
+ 1] = {
561 L(-1.723405393982209853244278760171643219530E-7),
562 L(-2.090508758514655456365709712333460087442E-5),
563 L(-9.140104013370974823232873472192719263019E-4),
564 L(-1.871349499990714843332742160292474780128E-2),
565 L(-1.948930738119938669637865956162512983416E-1),
566 L(-1.048764684978978127908439526343174139788E0
),
567 L(-2.827714929925679500237476105843643064698E0
),
568 L(-3.508761569156476114276988181329773987314E0
),
569 L(-1.669332202790211090973255098624488308989E0
),
570 L(-1.930796319299022954013840684651016077770E-1),
573 static const _Float128 Q2r7_3r2D
[NQ2r7_3r2D
+ 1] = {
574 L(1.680730662300831976234547482334347983474E-6),
575 L(2.084241442440551016475972218719621841120E-4),
576 L(9.445316642108367479043541702688736295579E-3),
577 L(2.044637889456631896650179477133252184672E-1),
578 L(2.316091982244297350829522534435350078205E0
),
579 L(1.412031891783015085196708811890448488865E1
),
580 L(4.583830154673223384837091077279595496149E1
),
581 L(7.549520609270909439885998474045974122261E1
),
582 L(5.697605832808113367197494052388203310638E1
),
583 L(1.601496240876192444526383314589371686234E1
),
584 /* 1.000000000000000000000000000000000000000E0 */
587 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
588 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
589 Peak relative error 9.5e-36
590 0.375 <= 1/x <= 0.4375 */
592 static const _Float128 Q2r3_2r7N
[NQ2r3_2r7N
+ 1] = {
593 L(-8.603042076329122085722385914954878953775E-7),
594 L(-7.701746260451647874214968882605186675720E-5),
595 L(-2.407932004380727587382493696877569654271E-3),
596 L(-3.403434217607634279028110636919987224188E-2),
597 L(-2.348707332185238159192422084985713102877E-1),
598 L(-7.957498841538254916147095255700637463207E-1),
599 L(-1.258469078442635106431098063707934348577E0
),
600 L(-8.162415474676345812459353639449971369890E-1),
601 L(-1.581783890269379690141513949609572806898E-1),
602 L(-1.890595651683552228232308756569450822905E-3),
605 static const _Float128 Q2r3_2r7D
[NQ2r3_2r7D
+ 1] = {
606 L(8.390017524798316921170710533381568175665E-6),
607 L(7.738148683730826286477254659973968763659E-4),
608 L(2.541480810958665794368759558791634341779E-2),
609 L(3.878879789711276799058486068562386244873E-1),
610 L(3.003783779325811292142957336802456109333E0
),
611 L(1.206480374773322029883039064575464497400E1
),
612 L(2.458414064785315978408974662900438351782E1
),
613 L(2.367237826273668567199042088835448715228E1
),
614 L(9.231451197519171090875569102116321676763E0
),
615 /* 1.000000000000000000000000000000000000000E0 */
618 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
619 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
620 Peak relative error 1.4e-36
621 0.4375 <= 1/x <= 0.5 */
623 static const _Float128 Q2_2r3N
[NQ2_2r3N
+ 1] = {
624 L(-5.552507516089087822166822364590806076174E-6),
625 L(-4.135067659799500521040944087433752970297E-4),
626 L(-1.059928728869218962607068840646564457980E-2),
627 L(-1.212070036005832342565792241385459023801E-1),
628 L(-6.688350110633603958684302153362735625156E-1),
629 L(-1.793587878197360221340277951304429821582E0
),
630 L(-2.225407682237197485644647380483725045326E0
),
631 L(-1.123402135458940189438898496348239744403E0
),
632 L(-1.679187241566347077204805190763597299805E-1),
633 L(-1.458550613639093752909985189067233504148E-3),
636 static const _Float128 Q2_2r3D
[NQ2_2r3D
+ 1] = {
637 L(5.415024336507980465169023996403597916115E-5),
638 L(4.179246497380453022046357404266022870788E-3),
639 L(1.136306384261959483095442402929502368598E-1),
640 L(1.422640343719842213484515445393284072830E0
),
641 L(8.968786703393158374728850922289204805764E0
),
642 L(2.914542473339246127533384118781216495934E1
),
643 L(4.781605421020380669870197378210457054685E1
),
644 L(3.693865837171883152382820584714795072937E1
),
645 L(1.153220502744204904763115556224395893076E1
),
646 /* 1.000000000000000000000000000000000000000E0 */
650 /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
653 neval (_Float128 x
, const _Float128
*p
, int n
)
668 /* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */
671 deval (_Float128 x
, const _Float128
*p
, int n
)
686 /* Bessel function of the first kind, order one. */
689 __ieee754_j1l (_Float128 x
)
691 _Float128 xx
, xinv
, z
, p
, q
, c
, s
, cc
, ss
;
703 if (xx
<= L(0x1p
-58))
705 _Float128 ret
= x
* L(0.5);
706 math_check_force_underflow (ret
);
708 __set_errno (ERANGE
);
715 p
= xx
* z
* neval (z
, J0_2N
, NJ0_2N
) / deval (z
, J0_2D
, NJ0_2D
);
723 cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
724 = 1/sqrt(2) * (-cos(x) + sin(x))
725 sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
726 = -1/sqrt(2) * (sin(x) + cos(x))
728 __sincosl (xx
, &s
, &c
);
731 if (xx
<= LDBL_MAX
/ 2)
733 z
= __cosl (xx
+ xx
);
742 z
= ONEOSQPI
* cc
/ sqrtl (xx
);
756 p
= neval (z
, P16_IN
, NP16_IN
) / deval (z
, P16_ID
, NP16_ID
);
757 q
= neval (z
, Q16_IN
, NQ16_IN
) / deval (z
, Q16_ID
, NQ16_ID
);
761 p
= neval (z
, P8_16N
, NP8_16N
) / deval (z
, P8_16D
, NP8_16D
);
762 q
= neval (z
, Q8_16N
, NQ8_16N
) / deval (z
, Q8_16D
, NQ8_16D
);
765 else if (xinv
<= 0.1875)
767 p
= neval (z
, P5_8N
, NP5_8N
) / deval (z
, P5_8D
, NP5_8D
);
768 q
= neval (z
, Q5_8N
, NQ5_8N
) / deval (z
, Q5_8D
, NQ5_8D
);
772 p
= neval (z
, P4_5N
, NP4_5N
) / deval (z
, P4_5D
, NP4_5D
);
773 q
= neval (z
, Q4_5N
, NQ4_5N
) / deval (z
, Q4_5D
, NQ4_5D
);
776 else /* if (xinv <= 0.5) */
782 p
= neval (z
, P3r2_4N
, NP3r2_4N
) / deval (z
, P3r2_4D
, NP3r2_4D
);
783 q
= neval (z
, Q3r2_4N
, NQ3r2_4N
) / deval (z
, Q3r2_4D
, NQ3r2_4D
);
787 p
= neval (z
, P2r7_3r2N
, NP2r7_3r2N
)
788 / deval (z
, P2r7_3r2D
, NP2r7_3r2D
);
789 q
= neval (z
, Q2r7_3r2N
, NQ2r7_3r2N
)
790 / deval (z
, Q2r7_3r2D
, NQ2r7_3r2D
);
793 else if (xinv
<= 0.4375)
795 p
= neval (z
, P2r3_2r7N
, NP2r3_2r7N
)
796 / deval (z
, P2r3_2r7D
, NP2r3_2r7D
);
797 q
= neval (z
, Q2r3_2r7N
, NQ2r3_2r7N
)
798 / deval (z
, Q2r3_2r7D
, NQ2r3_2r7D
);
802 p
= neval (z
, P2_2r3N
, NP2_2r3N
) / deval (z
, P2_2r3D
, NP2_2r3D
);
803 q
= neval (z
, Q2_2r3N
, NQ2_2r3N
) / deval (z
, Q2_2r3D
, NQ2_2r3D
);
808 q
= q
* xinv
+ L(0.375) * xinv
;
809 z
= ONEOSQPI
* (p
* cc
- q
* ss
) / sqrtl (xx
);
814 libm_alias_finite (__ieee754_j1l
, __j1l
)
817 /* Y1(x) = 2/pi * (log(x) * J1(x) - 1/x) + x R(x^2)
818 Peak relative error 6.2e-38
821 static const _Float128 Y0_2N
[NY0_2N
+ 1] = {
822 L(-6.804415404830253804408698161694720833249E19
),
823 L(1.805450517967019908027153056150465849237E19
),
824 L(-8.065747497063694098810419456383006737312E17
),
825 L(1.401336667383028259295830955439028236299E16
),
826 L(-1.171654432898137585000399489686629680230E14
),
827 L(5.061267920943853732895341125243428129150E11
),
828 L(-1.096677850566094204586208610960870217970E9
),
829 L(9.541172044989995856117187515882879304461E5
),
832 static const _Float128 Y0_2D
[NY0_2D
+ 1] = {
833 L(3.470629591820267059538637461549677594549E20
),
834 L(4.120796439009916326855848107545425217219E18
),
835 L(2.477653371652018249749350657387030814542E16
),
836 L(9.954678543353888958177169349272167762797E13
),
837 L(2.957927997613630118216218290262851197754E11
),
838 L(6.748421382188864486018861197614025972118E8
),
839 L(1.173453425218010888004562071020305709319E6
),
840 L(1.450335662961034949894009554536003377187E3
),
841 /* 1.000000000000000000000000000000000000000E0 */
845 /* Bessel function of the second kind, order one. */
848 __ieee754_y1l (_Float128 x
)
850 _Float128 xx
, xinv
, z
, p
, q
, c
, s
, cc
, ss
;
853 return 1 / (x
+ x
* x
);
857 return (zero
/ (zero
* x
));
858 return -1 / zero
; /* -inf and divide by zero exception. */
865 __set_errno (ERANGE
);
871 SET_RESTORE_ROUNDL (FE_TONEAREST
);
872 xx
= math_opt_barrier (xx
);
873 x
= math_opt_barrier (x
);
875 p
= xx
* neval (z
, Y0_2N
, NY0_2N
) / deval (z
, Y0_2D
, NY0_2D
);
876 p
= -TWOOPI
/ xx
+ p
;
877 p
= TWOOPI
* __ieee754_logl (x
) * __ieee754_j1l (x
) + p
;
883 cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
884 = 1/sqrt(2) * (-cos(x) + sin(x))
885 sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
886 = -1/sqrt(2) * (sin(x) + cos(x))
888 __sincosl (xx
, &s
, &c
);
891 if (xx
<= LDBL_MAX
/ 2)
893 z
= __cosl (xx
+ xx
);
901 return ONEOSQPI
* ss
/ sqrtl (xx
);
911 p
= neval (z
, P16_IN
, NP16_IN
) / deval (z
, P16_ID
, NP16_ID
);
912 q
= neval (z
, Q16_IN
, NQ16_IN
) / deval (z
, Q16_ID
, NQ16_ID
);
916 p
= neval (z
, P8_16N
, NP8_16N
) / deval (z
, P8_16D
, NP8_16D
);
917 q
= neval (z
, Q8_16N
, NQ8_16N
) / deval (z
, Q8_16D
, NQ8_16D
);
920 else if (xinv
<= 0.1875)
922 p
= neval (z
, P5_8N
, NP5_8N
) / deval (z
, P5_8D
, NP5_8D
);
923 q
= neval (z
, Q5_8N
, NQ5_8N
) / deval (z
, Q5_8D
, NQ5_8D
);
927 p
= neval (z
, P4_5N
, NP4_5N
) / deval (z
, P4_5D
, NP4_5D
);
928 q
= neval (z
, Q4_5N
, NQ4_5N
) / deval (z
, Q4_5D
, NQ4_5D
);
931 else /* if (xinv <= 0.5) */
937 p
= neval (z
, P3r2_4N
, NP3r2_4N
) / deval (z
, P3r2_4D
, NP3r2_4D
);
938 q
= neval (z
, Q3r2_4N
, NQ3r2_4N
) / deval (z
, Q3r2_4D
, NQ3r2_4D
);
942 p
= neval (z
, P2r7_3r2N
, NP2r7_3r2N
)
943 / deval (z
, P2r7_3r2D
, NP2r7_3r2D
);
944 q
= neval (z
, Q2r7_3r2N
, NQ2r7_3r2N
)
945 / deval (z
, Q2r7_3r2D
, NQ2r7_3r2D
);
948 else if (xinv
<= 0.4375)
950 p
= neval (z
, P2r3_2r7N
, NP2r3_2r7N
)
951 / deval (z
, P2r3_2r7D
, NP2r3_2r7D
);
952 q
= neval (z
, Q2r3_2r7N
, NQ2r3_2r7N
)
953 / deval (z
, Q2r3_2r7D
, NQ2r3_2r7D
);
957 p
= neval (z
, P2_2r3N
, NP2_2r3N
) / deval (z
, P2_2r3D
, NP2_2r3D
);
958 q
= neval (z
, Q2_2r3N
, NQ2_2r3N
) / deval (z
, Q2_2r3D
, NQ2_2r3D
);
963 q
= q
* xinv
+ L(0.375) * xinv
;
964 z
= ONEOSQPI
* (p
* ss
+ q
* cc
) / sqrtl (xx
);
967 libm_alias_finite (__ieee754_y1l
, __y1l
)