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