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
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
95 License along with this library; if not, see
96 <https://www.gnu.org/licenses/>. */
97
98 #include <errno.h>
99 #include <math.h>
100 #include <math_private.h>
101 #include <fenv_private.h>
102 #include <math-underflow.h>
103 #include <float.h>
104 #include <libm-alias-finite.h>
105
106 /* 1 / sqrt(pi) */
107 static const _Float128 ONEOSQPI = L(5.6418958354775628694807945156077258584405E-1);
108 /* 2 / pi */
109 static const _Float128 TWOOPI = L(6.3661977236758134307553505349005744813784E-1);
110 static const _Float128 zero = 0;
111
112 /* J1(x) = .5x + x x^2 R(x^2)
113 Peak relative error 1.9e-35
114 0 <= x <= 2 */
115 #define NJ0_2N 6
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)
124 };
125 #define NJ0_2D 6
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 */
135 };
136
137 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
138 0 <= 1/x <= .0625
139 Peak relative error 3.6e-36 */
140 #define NP16_IN 9
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),
152 };
153 #define NP16_ID 9
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 */
166 };
167
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 */
171 #define NP8_16N 11
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),
185 };
186 #define NP8_16D 10
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 */
200 };
201
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 */
205 #define NP5_8N 10
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)
218 };
219 #define NP5_8D 10
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 */
233 };
234
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 */
238 #define NP4_5N 10
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),
251 };
252 #define NP4_5D 9
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 */
265 };
266
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 */
270 #define NP3r2_4N 9
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),
282 };
283 #define NP3r2_4D 9
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 */
296 };
297
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 */
301 #define NP2r7_3r2N 9
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),
313 };
314 #define NP2r7_3r2D 8
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 */
326 };
327
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 */
331 #define NP2r3_2r7N 9
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),
343 };
344 #define NP2r3_2r7D 8
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 */
356 };
357
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 */
361 #define NP2_2r3N 8
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),
372 };
373 #define NP2_2r3D 8
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 */
385 };
386
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
390 0 <= 1/x <= .0625 */
391 #define NQ16_IN 10
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),
404 };
405 #define NQ16_ID 9
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 */
418 };
419
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 */
424 #define NQ8_16N 11
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),
438 };
439 #define NQ8_16D 11
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 */
454 };
455
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 */
460 #define NQ5_8N 10
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),
473 };
474 #define NQ5_8D 10
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 */
488 };
489
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 */
494 #define NQ4_5N 10
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),
507 };
508 #define NQ4_5D 9
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 */
521 };
522
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 */
527 #define NQ3r2_4N 9
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),
539 };
540 #define NQ3r2_4D 9
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 */
553 };
554
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 */
559 #define NQ2r7_3r2N 9
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),
571 };
572 #define NQ2r7_3r2D 9
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 */
585 };
586
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 */
591 #define NQ2r3_2r7N 9
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),
603 };
604 #define NQ2r3_2r7D 8
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 */
616 };
617
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 */
622 #define NQ2_2r3N 9
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),
634 };
635 #define NQ2_2r3D 8
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 */
647 };
648
649
650 /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
651
652 static _Float128
653 neval (_Float128 x, const _Float128 *p, int n)
654 {
655 _Float128 y;
656
657 p += n;
658 y = *p--;
659 do
660 {
661 y = y * x + *p--;
662 }
663 while (--n > 0);
664 return y;
665 }
666
667
668 /* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */
669
670 static _Float128
671 deval (_Float128 x, const _Float128 *p, int n)
672 {
673 _Float128 y;
674
675 p += n;
676 y = x + *p--;
677 do
678 {
679 y = y * x + *p--;
680 }
681 while (--n > 0);
682 return y;
683 }
684
685
686 /* Bessel function of the first kind, order one. */
687
688 _Float128
689 __ieee754_j1l (_Float128 x)
690 {
691 _Float128 xx, xinv, z, p, q, c, s, cc, ss;
692
693 if (! isfinite (x))
694 {
695 if (x != x)
696 return x + x;
697 else
698 return 0;
699 }
700 if (x == 0)
701 return x;
702 xx = fabsl (x);
703 if (xx <= L(0x1p-58))
704 {
705 _Float128 ret = x * L(0.5);
706 math_check_force_underflow (ret);
707 if (ret == 0)
708 __set_errno (ERANGE);
709 return ret;
710 }
711 if (xx <= 2)
712 {
713 /* 0 <= x <= 2 */
714 z = xx * xx;
715 p = xx * z * neval (z, J0_2N, NJ0_2N) / deval (z, J0_2D, NJ0_2D);
716 p += L(0.5) * xx;
717 if (x < 0)
718 p = -p;
719 return p;
720 }
721
722 /* X = x - 3 pi/4
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))
727 cf. Fdlibm. */
728 __sincosl (xx, &s, &c);
729 ss = -s - c;
730 cc = s - c;
731 if (xx <= LDBL_MAX / 2)
732 {
733 z = __cosl (xx + xx);
734 if ((s * c) > 0)
735 cc = z / ss;
736 else
737 ss = z / cc;
738 }
739
740 if (xx > L(0x1p256))
741 {
742 z = ONEOSQPI * cc / sqrtl (xx);
743 if (x < 0)
744 z = -z;
745 return z;
746 }
747
748 xinv = 1 / xx;
749 z = xinv * xinv;
750 if (xinv <= 0.25)
751 {
752 if (xinv <= 0.125)
753 {
754 if (xinv <= 0.0625)
755 {
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);
758 }
759 else
760 {
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);
763 }
764 }
765 else if (xinv <= 0.1875)
766 {
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);
769 }
770 else
771 {
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);
774 }
775 } /* .25 */
776 else /* if (xinv <= 0.5) */
777 {
778 if (xinv <= 0.375)
779 {
780 if (xinv <= 0.3125)
781 {
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);
784 }
785 else
786 {
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);
791 }
792 }
793 else if (xinv <= 0.4375)
794 {
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);
799 }
800 else
801 {
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);
804 }
805 }
806 p = 1 + z * p;
807 q = z * q;
808 q = q * xinv + L(0.375) * xinv;
809 z = ONEOSQPI * (p * cc - q * ss) / sqrtl (xx);
810 if (x < 0)
811 z = -z;
812 return z;
813 }
814 libm_alias_finite (__ieee754_j1l, __j1l)
815
816
817 /* Y1(x) = 2/pi * (log(x) * J1(x) - 1/x) + x R(x^2)
818 Peak relative error 6.2e-38
819 0 <= x <= 2 */
820 #define NY0_2N 7
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),
830 };
831 #define NY0_2D 7
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 */
842 };
843
844
845 /* Bessel function of the second kind, order one. */
846
847 _Float128
848 __ieee754_y1l (_Float128 x)
849 {
850 _Float128 xx, xinv, z, p, q, c, s, cc, ss;
851
852 if (! isfinite (x))
853 return 1 / (x + x * x);
854 if (x <= 0)
855 {
856 if (x < 0)
857 return (zero / (zero * x));
858 return -1 / zero; /* -inf and divide by zero exception. */
859 }
860 xx = fabsl (x);
861 if (xx <= 0x1p-114)
862 {
863 z = -TWOOPI / x;
864 if (isinf (z))
865 __set_errno (ERANGE);
866 return z;
867 }
868 if (xx <= 2)
869 {
870 /* 0 <= x <= 2 */
871 SET_RESTORE_ROUNDL (FE_TONEAREST);
872 xx = math_opt_barrier (xx);
873 x = math_opt_barrier (x);
874 z = xx * xx;
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;
878 math_force_eval (p);
879 return p;
880 }
881
882 /* X = x - 3 pi/4
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))
887 cf. Fdlibm. */
888 __sincosl (xx, &s, &c);
889 ss = -s - c;
890 cc = s - c;
891 if (xx <= LDBL_MAX / 2)
892 {
893 z = __cosl (xx + xx);
894 if ((s * c) > 0)
895 cc = z / ss;
896 else
897 ss = z / cc;
898 }
899
900 if (xx > L(0x1p256))
901 return ONEOSQPI * ss / sqrtl (xx);
902
903 xinv = 1 / xx;
904 z = xinv * xinv;
905 if (xinv <= 0.25)
906 {
907 if (xinv <= 0.125)
908 {
909 if (xinv <= 0.0625)
910 {
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);
913 }
914 else
915 {
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);
918 }
919 }
920 else if (xinv <= 0.1875)
921 {
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);
924 }
925 else
926 {
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);
929 }
930 } /* .25 */
931 else /* if (xinv <= 0.5) */
932 {
933 if (xinv <= 0.375)
934 {
935 if (xinv <= 0.3125)
936 {
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);
939 }
940 else
941 {
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);
946 }
947 }
948 else if (xinv <= 0.4375)
949 {
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);
954 }
955 else
956 {
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);
959 }
960 }
961 p = 1 + z * p;
962 q = z * q;
963 q = q * xinv + L(0.375) * xinv;
964 z = ONEOSQPI * (p * ss + q * cc) / sqrtl (xx);
965 return z;
966 }
967 libm_alias_finite (__ieee754_y1l, __y1l)