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) */
107static const _Float128 ONEOSQPI = L(5.6418958354775628694807945156077258584405E-1);
108/* 2 / pi */
109static const _Float128 TWOOPI = L(6.3661977236758134307553505349005744813784E-1);
110static 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
116static 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
126static 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
141static 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
154static 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
172static 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
187static 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
206static 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
220static 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
239static 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
253static 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
271static 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
284static 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
302static 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
315static 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
332static 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
345static 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
362static 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
374static 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
392static 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
406static 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
425static 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
440static 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
461static 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
475static 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
495static 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
509static 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
528static 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
541static 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
560static 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
573static 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
592static 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
605static 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
623static 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
636static 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
652static _Float128
653neval (_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
670static _Float128
671deval (_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}
814libm_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
821static 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
832static 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 z = xx * xx;
873 p = xx * neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D);
874 p = -TWOOPI / xx + p;
875 p = TWOOPI * __ieee754_logl (x) * __ieee754_j1l (x) + p;
876 return p;
877 }
878
879 /* X = x - 3 pi/4
880 cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
881 = 1/sqrt(2) * (-cos(x) + sin(x))
882 sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
883 = -1/sqrt(2) * (sin(x) + cos(x))
884 cf. Fdlibm. */
885 __sincosl (xx, &s, &c);
886 ss = -s - c;
887 cc = s - c;
888 if (xx <= LDBL_MAX / 2)
889 {
890 z = __cosl (xx + xx);
891 if ((s * c) > 0)
892 cc = z / ss;
893 else
894 ss = z / cc;
895 }
896
897 if (xx > L(0x1p256))
898 return ONEOSQPI * ss / sqrtl (xx);
899
900 xinv = 1 / xx;
901 z = xinv * xinv;
902 if (xinv <= 0.25)
903 {
904 if (xinv <= 0.125)
905 {
906 if (xinv <= 0.0625)
907 {
908 p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID);
909 q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID);
910 }
911 else
912 {
913 p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D);
914 q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D);
915 }
916 }
917 else if (xinv <= 0.1875)
918 {
919 p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D);
920 q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D);
921 }
922 else
923 {
924 p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D);
925 q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D);
926 }
927 } /* .25 */
928 else /* if (xinv <= 0.5) */
929 {
930 if (xinv <= 0.375)
931 {
932 if (xinv <= 0.3125)
933 {
934 p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D);
935 q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D);
936 }
937 else
938 {
939 p = neval (z, P2r7_3r2N, NP2r7_3r2N)
940 / deval (z, P2r7_3r2D, NP2r7_3r2D);
941 q = neval (z, Q2r7_3r2N, NQ2r7_3r2N)
942 / deval (z, Q2r7_3r2D, NQ2r7_3r2D);
943 }
944 }
945 else if (xinv <= 0.4375)
946 {
947 p = neval (z, P2r3_2r7N, NP2r3_2r7N)
948 / deval (z, P2r3_2r7D, NP2r3_2r7D);
949 q = neval (z, Q2r3_2r7N, NQ2r3_2r7N)
950 / deval (z, Q2r3_2r7D, NQ2r3_2r7D);
951 }
952 else
953 {
954 p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
955 q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
956 }
957 }
958 p = 1 + z * p;
959 q = z * q;
960 q = q * xinv + L(0.375) * xinv;
961 z = ONEOSQPI * (p * ss + q * cc) / sqrtl (xx);
962 return z;
963}
964libm_alias_finite (__ieee754_y1l, __y1l)
965