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