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