1 /* j0l.c
2 *
3 * Bessel function of order zero
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * long double x, y, j0l();
10 *
11 * y = j0l( x );
12 *
13 *
14 *
15 * DESCRIPTION:
16 *
17 * Returns Bessel function of first kind, order zero 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
21 * is J0(x) = 1 - x^2 / 4 + x^4 R(x^2)
22 * The second interval is further partitioned into eight equal segments
23 * of 1/x.
24 *
25 * J0(x) = sqrt(2/(pi x)) (P0(x) cos(X) - Q0(x) sin(X)),
26 * X = x - pi/4,
27 *
28 * and the auxiliary functions are given by
29 *
30 * J0(x)cos(X) + Y0(x)sin(X) = sqrt( 2/(pi x)) P0(x),
31 * P0(x) = 1 + 1/x^2 R(1/x^2)
32 *
33 * Y0(x)cos(X) - J0(x)sin(X) = sqrt( 2/(pi x)) Q0(x),
34 * Q0(x) = 1/x (-.125 + 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 1.7e-34 2.4e-35
43 *
44 *
45 */
46
47 /* y0l.c
48 *
49 * Bessel function of the second kind, order zero
50 *
51 *
52 *
53 * SYNOPSIS:
54 *
55 * double x, y, y0l();
56 *
57 * y = y0l( x );
58 *
59 *
60 *
61 * DESCRIPTION:
62 *
63 * Returns Bessel function of the second kind, of order
64 * zero, of the argument.
65 *
66 * The approximation is the same as for J0(x), and
67 * Y0(x) = sqrt(2/(pi x)) (P0(x) sin(X) + Q0(x) cos(X)).
68 *
69 * ACCURACY:
70 *
71 * Absolute error, when y0(x) < 1; else relative error:
72 *
73 * arithmetic domain # trials peak rms
74 * IEEE 0, 30 100000 3.0e-34 2.7e-35
75 *
76 */
77
78 /* Copyright 2001 by Stephen L. Moshier (moshier@na-net.ornl.gov).
79
80 This library is free software; you can redistribute it and/or
81 modify it under the terms of the GNU Lesser General Public
82 License as published by the Free Software Foundation; either
83 version 2.1 of the License, or (at your option) any later version.
84
85 This library is distributed in the hope that it will be useful,
86 but WITHOUT ANY WARRANTY; without even the implied warranty of
87 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
88 Lesser General Public License for more details.
89
90 You should have received a copy of the GNU Lesser General Public
91 License along with this library; if not, see
92 <https://www.gnu.org/licenses/>. */
93
94 #include <math.h>
95 #include <math_private.h>
96 #include <float.h>
97 #include <libm-alias-finite.h>
98
99 /* 1 / sqrt(pi) */
100 static const _Float128 ONEOSQPI = L(5.6418958354775628694807945156077258584405E-1);
101 /* 2 / pi */
102 static const _Float128 TWOOPI = L(6.3661977236758134307553505349005744813784E-1);
103 static const _Float128 zero = 0;
104
105 /* J0(x) = 1 - x^2/4 + x^2 x^2 R(x^2)
106 Peak relative error 3.4e-37
107 0 <= x <= 2 */
108 #define NJ0_2N 6
109 static const _Float128 J0_2N[NJ0_2N + 1] = {
110 L(3.133239376997663645548490085151484674892E16),
111 L(-5.479944965767990821079467311839107722107E14),
112 L(6.290828903904724265980249871997551894090E12),
113 L(-3.633750176832769659849028554429106299915E10),
114 L(1.207743757532429576399485415069244807022E8),
115 L(-2.107485999925074577174305650549367415465E5),
116 L(1.562826808020631846245296572935547005859E2),
117 };
118 #define NJ0_2D 6
119 static const _Float128 J0_2D[NJ0_2D + 1] = {
120 L(2.005273201278504733151033654496928968261E18),
121 L(2.063038558793221244373123294054149790864E16),
122 L(1.053350447931127971406896594022010524994E14),
123 L(3.496556557558702583143527876385508882310E11),
124 L(8.249114511878616075860654484367133976306E8),
125 L(1.402965782449571800199759247964242790589E6),
126 L(1.619910762853439600957801751815074787351E3),
127 /* 1.000000000000000000000000000000000000000E0 */
128 };
129
130 /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2),
131 0 <= 1/x <= .0625
132 Peak relative error 3.3e-36 */
133 #define NP16_IN 9
134 static const _Float128 P16_IN[NP16_IN + 1] = {
135 L(-1.901689868258117463979611259731176301065E-16),
136 L(-1.798743043824071514483008340803573980931E-13),
137 L(-6.481746687115262291873324132944647438959E-11),
138 L(-1.150651553745409037257197798528294248012E-8),
139 L(-1.088408467297401082271185599507222695995E-6),
140 L(-5.551996725183495852661022587879817546508E-5),
141 L(-1.477286941214245433866838787454880214736E-3),
142 L(-1.882877976157714592017345347609200402472E-2),
143 L(-9.620983176855405325086530374317855880515E-2),
144 L(-1.271468546258855781530458854476627766233E-1),
145 };
146 #define NP16_ID 9
147 static const _Float128 P16_ID[NP16_ID + 1] = {
148 L(2.704625590411544837659891569420764475007E-15),
149 L(2.562526347676857624104306349421985403573E-12),
150 L(9.259137589952741054108665570122085036246E-10),
151 L(1.651044705794378365237454962653430805272E-7),
152 L(1.573561544138733044977714063100859136660E-5),
153 L(8.134482112334882274688298469629884804056E-4),
154 L(2.219259239404080863919375103673593571689E-2),
155 L(2.976990606226596289580242451096393862792E-1),
156 L(1.713895630454693931742734911930937246254E0),
157 L(3.231552290717904041465898249160757368855E0),
158 /* 1.000000000000000000000000000000000000000E0 */
159 };
160
161 /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2)
162 0.0625 <= 1/x <= 0.125
163 Peak relative error 2.4e-35 */
164 #define NP8_16N 10
165 static const _Float128 P8_16N[NP8_16N + 1] = {
166 L(-2.335166846111159458466553806683579003632E-15),
167 L(-1.382763674252402720401020004169367089975E-12),
168 L(-3.192160804534716696058987967592784857907E-10),
169 L(-3.744199606283752333686144670572632116899E-8),
170 L(-2.439161236879511162078619292571922772224E-6),
171 L(-9.068436986859420951664151060267045346549E-5),
172 L(-1.905407090637058116299757292660002697359E-3),
173 L(-2.164456143936718388053842376884252978872E-2),
174 L(-1.212178415116411222341491717748696499966E-1),
175 L(-2.782433626588541494473277445959593334494E-1),
176 L(-1.670703190068873186016102289227646035035E-1),
177 };
178 #define NP8_16D 10
179 static const _Float128 P8_16D[NP8_16D + 1] = {
180 L(3.321126181135871232648331450082662856743E-14),
181 L(1.971894594837650840586859228510007703641E-11),
182 L(4.571144364787008285981633719513897281690E-9),
183 L(5.396419143536287457142904742849052402103E-7),
184 L(3.551548222385845912370226756036899901549E-5),
185 L(1.342353874566932014705609788054598013516E-3),
186 L(2.899133293006771317589357444614157734385E-2),
187 L(3.455374978185770197704507681491574261545E-1),
188 L(2.116616964297512311314454834712634820514E0),
189 L(5.850768316827915470087758636881584174432E0),
190 L(5.655273858938766830855753983631132928968E0),
191 /* 1.000000000000000000000000000000000000000E0 */
192 };
193
194 /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2)
195 0.125 <= 1/x <= 0.1875
196 Peak relative error 2.7e-35 */
197 #define NP5_8N 10
198 static const _Float128 P5_8N[NP5_8N + 1] = {
199 L(-1.270478335089770355749591358934012019596E-12),
200 L(-4.007588712145412921057254992155810347245E-10),
201 L(-4.815187822989597568124520080486652009281E-8),
202 L(-2.867070063972764880024598300408284868021E-6),
203 L(-9.218742195161302204046454768106063638006E-5),
204 L(-1.635746821447052827526320629828043529997E-3),
205 L(-1.570376886640308408247709616497261011707E-2),
206 L(-7.656484795303305596941813361786219477807E-2),
207 L(-1.659371030767513274944805479908858628053E-1),
208 L(-1.185340550030955660015841796219919804915E-1),
209 L(-8.920026499909994671248893388013790366712E-3),
210 };
211 #define NP5_8D 9
212 static const _Float128 P5_8D[NP5_8D + 1] = {
213 L(1.806902521016705225778045904631543990314E-11),
214 L(5.728502760243502431663549179135868966031E-9),
215 L(6.938168504826004255287618819550667978450E-7),
216 L(4.183769964807453250763325026573037785902E-5),
217 L(1.372660678476925468014882230851637878587E-3),
218 L(2.516452105242920335873286419212708961771E-2),
219 L(2.550502712902647803796267951846557316182E-1),
220 L(1.365861559418983216913629123778747617072E0),
221 L(3.523825618308783966723472468855042541407E0),
222 L(3.656365803506136165615111349150536282434E0),
223 /* 1.000000000000000000000000000000000000000E0 */
224 };
225
226 /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2)
227 Peak relative error 3.5e-35
228 0.1875 <= 1/x <= 0.25 */
229 #define NP4_5N 9
230 static const _Float128 P4_5N[NP4_5N + 1] = {
231 L(-9.791405771694098960254468859195175708252E-10),
232 L(-1.917193059944531970421626610188102836352E-7),
233 L(-1.393597539508855262243816152893982002084E-5),
234 L(-4.881863490846771259880606911667479860077E-4),
235 L(-8.946571245022470127331892085881699269853E-3),
236 L(-8.707474232568097513415336886103899434251E-2),
237 L(-4.362042697474650737898551272505525973766E-1),
238 L(-1.032712171267523975431451359962375617386E0),
239 L(-9.630502683169895107062182070514713702346E-1),
240 L(-2.251804386252969656586810309252357233320E-1),
241 };
242 #define NP4_5D 9
243 static const _Float128 P4_5D[NP4_5D + 1] = {
244 L(1.392555487577717669739688337895791213139E-8),
245 L(2.748886559120659027172816051276451376854E-6),
246 L(2.024717710644378047477189849678576659290E-4),
247 L(7.244868609350416002930624752604670292469E-3),
248 L(1.373631762292244371102989739300382152416E-1),
249 L(1.412298581400224267910294815260613240668E0),
250 L(7.742495637843445079276397723849017617210E0),
251 L(2.138429269198406512028307045259503811861E1),
252 L(2.651547684548423476506826951831712762610E1),
253 L(1.167499382465291931571685222882909166935E1),
254 /* 1.000000000000000000000000000000000000000E0 */
255 };
256
257 /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2)
258 Peak relative error 2.3e-36
259 0.25 <= 1/x <= 0.3125 */
260 #define NP3r2_4N 9
261 static const _Float128 P3r2_4N[NP3r2_4N + 1] = {
262 L(-2.589155123706348361249809342508270121788E-8),
263 L(-3.746254369796115441118148490849195516593E-6),
264 L(-1.985595497390808544622893738135529701062E-4),
265 L(-5.008253705202932091290132760394976551426E-3),
266 L(-6.529469780539591572179155511840853077232E-2),
267 L(-4.468736064761814602927408833818990271514E-1),
268 L(-1.556391252586395038089729428444444823380E0),
269 L(-2.533135309840530224072920725976994981638E0),
270 L(-1.605509621731068453869408718565392869560E0),
271 L(-2.518966692256192789269859830255724429375E-1),
272 };
273 #define NP3r2_4D 9
274 static const _Float128 P3r2_4D[NP3r2_4D + 1] = {
275 L(3.682353957237979993646169732962573930237E-7),
276 L(5.386741661883067824698973455566332102029E-5),
277 L(2.906881154171822780345134853794241037053E-3),
278 L(7.545832595801289519475806339863492074126E-2),
279 L(1.029405357245594877344360389469584526654E0),
280 L(7.565706120589873131187989560509757626725E0),
281 L(2.951172890699569545357692207898667665796E1),
282 L(5.785723537170311456298467310529815457536E1),
283 L(5.095621464598267889126015412522773474467E1),
284 L(1.602958484169953109437547474953308401442E1),
285 /* 1.000000000000000000000000000000000000000E0 */
286 };
287
288 /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2)
289 Peak relative error 1.0e-35
290 0.3125 <= 1/x <= 0.375 */
291 #define NP2r7_3r2N 9
292 static const _Float128 P2r7_3r2N[NP2r7_3r2N + 1] = {
293 L(-1.917322340814391131073820537027234322550E-7),
294 L(-1.966595744473227183846019639723259011906E-5),
295 L(-7.177081163619679403212623526632690465290E-4),
296 L(-1.206467373860974695661544653741899755695E-2),
297 L(-1.008656452188539812154551482286328107316E-1),
298 L(-4.216016116408810856620947307438823892707E-1),
299 L(-8.378631013025721741744285026537009814161E-1),
300 L(-6.973895635309960850033762745957946272579E-1),
301 L(-1.797864718878320770670740413285763554812E-1),
302 L(-4.098025357743657347681137871388402849581E-3),
303 };
304 #define NP2r7_3r2D 8
305 static const _Float128 P2r7_3r2D[NP2r7_3r2D + 1] = {
306 L(2.726858489303036441686496086962545034018E-6),
307 L(2.840430827557109238386808968234848081424E-4),
308 L(1.063826772041781947891481054529454088832E-2),
309 L(1.864775537138364773178044431045514405468E-1),
310 L(1.665660052857205170440952607701728254211E0),
311 L(7.723745889544331153080842168958348568395E0),
312 L(1.810726427571829798856428548102077799835E1),
313 L(1.986460672157794440666187503833545388527E1),
314 L(8.645503204552282306364296517220055815488E0),
315 /* 1.000000000000000000000000000000000000000E0 */
316 };
317
318 /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2)
319 Peak relative error 1.3e-36
320 0.3125 <= 1/x <= 0.4375 */
321 #define NP2r3_2r7N 9
322 static const _Float128 P2r3_2r7N[NP2r3_2r7N + 1] = {
323 L(-1.594642785584856746358609622003310312622E-6),
324 L(-1.323238196302221554194031733595194539794E-4),
325 L(-3.856087818696874802689922536987100372345E-3),
326 L(-5.113241710697777193011470733601522047399E-2),
327 L(-3.334229537209911914449990372942022350558E-1),
328 L(-1.075703518198127096179198549659283422832E0),
329 L(-1.634174803414062725476343124267110981807E0),
330 L(-1.030133247434119595616826842367268304880E0),
331 L(-1.989811539080358501229347481000707289391E-1),
332 L(-3.246859189246653459359775001466924610236E-3),
333 };
334 #define NP2r3_2r7D 8
335 static const _Float128 P2r3_2r7D[NP2r3_2r7D + 1] = {
336 L(2.267936634217251403663034189684284173018E-5),
337 L(1.918112982168673386858072491437971732237E-3),
338 L(5.771704085468423159125856786653868219522E-2),
339 L(8.056124451167969333717642810661498890507E-1),
340 L(5.687897967531010276788680634413789328776E0),
341 L(2.072596760717695491085444438270778394421E1),
342 L(3.801722099819929988585197088613160496684E1),
343 L(3.254620235902912339534998592085115836829E1),
344 L(1.104847772130720331801884344645060675036E1),
345 /* 1.000000000000000000000000000000000000000E0 */
346 };
347
348 /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2)
349 Peak relative error 1.2e-35
350 0.4375 <= 1/x <= 0.5 */
351 #define NP2_2r3N 8
352 static const _Float128 P2_2r3N[NP2_2r3N + 1] = {
353 L(-1.001042324337684297465071506097365389123E-4),
354 L(-6.289034524673365824853547252689991418981E-3),
355 L(-1.346527918018624234373664526930736205806E-1),
356 L(-1.268808313614288355444506172560463315102E0),
357 L(-5.654126123607146048354132115649177406163E0),
358 L(-1.186649511267312652171775803270911971693E1),
359 L(-1.094032424931998612551588246779200724257E1),
360 L(-3.728792136814520055025256353193674625267E0),
361 L(-3.000348318524471807839934764596331810608E-1),
362 };
363 #define NP2_2r3D 8
364 static const _Float128 P2_2r3D[NP2_2r3D + 1] = {
365 L(1.423705538269770974803901422532055612980E-3),
366 L(9.171476630091439978533535167485230575894E-2),
367 L(2.049776318166637248868444600215942828537E0),
368 L(2.068970329743769804547326701946144899583E1),
369 L(1.025103500560831035592731539565060347709E2),
370 L(2.528088049697570728252145557167066708284E2),
371 L(2.992160327587558573740271294804830114205E2),
372 L(1.540193761146551025832707739468679973036E2),
373 L(2.779516701986912132637672140709452502650E1),
374 /* 1.000000000000000000000000000000000000000E0 */
375 };
376
377 /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x),
378 Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2))
379 Peak relative error 2.2e-35
380 0 <= 1/x <= .0625 */
381 #define NQ16_IN 10
382 static const _Float128 Q16_IN[NQ16_IN + 1] = {
383 L(2.343640834407975740545326632205999437469E-18),
384 L(2.667978112927811452221176781536278257448E-15),
385 L(1.178415018484555397390098879501969116536E-12),
386 L(2.622049767502719728905924701288614016597E-10),
387 L(3.196908059607618864801313380896308968673E-8),
388 L(2.179466154171673958770030655199434798494E-6),
389 L(8.139959091628545225221976413795645177291E-5),
390 L(1.563900725721039825236927137885747138654E-3),
391 L(1.355172364265825167113562519307194840307E-2),
392 L(3.928058355906967977269780046844768588532E-2),
393 L(1.107891967702173292405380993183694932208E-2),
394 };
395 #define NQ16_ID 9
396 static const _Float128 Q16_ID[NQ16_ID + 1] = {
397 L(3.199850952578356211091219295199301766718E-17),
398 L(3.652601488020654842194486058637953363918E-14),
399 L(1.620179741394865258354608590461839031281E-11),
400 L(3.629359209474609630056463248923684371426E-9),
401 L(4.473680923894354600193264347733477363305E-7),
402 L(3.106368086644715743265603656011050476736E-5),
403 L(1.198239259946770604954664925153424252622E-3),
404 L(2.446041004004283102372887804475767568272E-2),
405 L(2.403235525011860603014707768815113698768E-1),
406 L(9.491006790682158612266270665136910927149E-1),
407 /* 1.000000000000000000000000000000000000000E0 */
408 };
409
410 /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x),
411 Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2))
412 Peak relative error 5.1e-36
413 0.0625 <= 1/x <= 0.125 */
414 #define NQ8_16N 11
415 static const _Float128 Q8_16N[NQ8_16N + 1] = {
416 L(1.001954266485599464105669390693597125904E-17),
417 L(7.545499865295034556206475956620160007849E-15),
418 L(2.267838684785673931024792538193202559922E-12),
419 L(3.561909705814420373609574999542459912419E-10),
420 L(3.216201422768092505214730633842924944671E-8),
421 L(1.731194793857907454569364622452058554314E-6),
422 L(5.576944613034537050396518509871004586039E-5),
423 L(1.051787760316848982655967052985391418146E-3),
424 L(1.102852974036687441600678598019883746959E-2),
425 L(5.834647019292460494254225988766702933571E-2),
426 L(1.290281921604364618912425380717127576529E-1),
427 L(7.598886310387075708640370806458926458301E-2),
428 };
429 #define NQ8_16D 11
430 static const _Float128 Q8_16D[NQ8_16D + 1] = {
431 L(1.368001558508338469503329967729951830843E-16),
432 L(1.034454121857542147020549303317348297289E-13),
433 L(3.128109209247090744354764050629381674436E-11),
434 L(4.957795214328501986562102573522064468671E-9),
435 L(4.537872468606711261992676606899273588899E-7),
436 L(2.493639207101727713192687060517509774182E-5),
437 L(8.294957278145328349785532236663051405805E-4),
438 L(1.646471258966713577374948205279380115839E-2),
439 L(1.878910092770966718491814497982191447073E-1),
440 L(1.152641605706170353727903052525652504075E0),
441 L(3.383550240669773485412333679367792932235E0),
442 L(3.823875252882035706910024716609908473970E0),
443 /* 1.000000000000000000000000000000000000000E0 */
444 };
445
446 /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x),
447 Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2))
448 Peak relative error 3.9e-35
449 0.125 <= 1/x <= 0.1875 */
450 #define NQ5_8N 10
451 static const _Float128 Q5_8N[NQ5_8N + 1] = {
452 L(1.750399094021293722243426623211733898747E-13),
453 L(6.483426211748008735242909236490115050294E-11),
454 L(9.279430665656575457141747875716899958373E-9),
455 L(6.696634968526907231258534757736576340266E-7),
456 L(2.666560823798895649685231292142838188061E-5),
457 L(6.025087697259436271271562769707550594540E-4),
458 L(7.652807734168613251901945778921336353485E-3),
459 L(5.226269002589406461622551452343519078905E-2),
460 L(1.748390159751117658969324896330142895079E-1),
461 L(2.378188719097006494782174902213083589660E-1),
462 L(8.383984859679804095463699702165659216831E-2),
463 };
464 #define NQ5_8D 10
465 static const _Float128 Q5_8D[NQ5_8D + 1] = {
466 L(2.389878229704327939008104855942987615715E-12),
467 L(8.926142817142546018703814194987786425099E-10),
468 L(1.294065862406745901206588525833274399038E-7),
469 L(9.524139899457666250828752185212769682191E-6),
470 L(3.908332488377770886091936221573123353489E-4),
471 L(9.250427033957236609624199884089916836748E-3),
472 L(1.263420066165922645975830877751588421451E-1),
473 L(9.692527053860420229711317379861733180654E-1),
474 L(3.937813834630430172221329298841520707954E0),
475 L(7.603126427436356534498908111445191312181E0),
476 L(5.670677653334105479259958485084550934305E0),
477 /* 1.000000000000000000000000000000000000000E0 */
478 };
479
480 /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x),
481 Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2))
482 Peak relative error 3.2e-35
483 0.1875 <= 1/x <= 0.25 */
484 #define NQ4_5N 10
485 static const _Float128 Q4_5N[NQ4_5N + 1] = {
486 L(2.233870042925895644234072357400122854086E-11),
487 L(5.146223225761993222808463878999151699792E-9),
488 L(4.459114531468296461688753521109797474523E-7),
489 L(1.891397692931537975547242165291668056276E-5),
490 L(4.279519145911541776938964806470674565504E-4),
491 L(5.275239415656560634702073291768904783989E-3),
492 L(3.468698403240744801278238473898432608887E-2),
493 L(1.138773146337708415188856882915457888274E-1),
494 L(1.622717518946443013587108598334636458955E-1),
495 L(7.249040006390586123760992346453034628227E-2),
496 L(1.941595365256460232175236758506411486667E-3),
497 };
498 #define NQ4_5D 9
499 static const _Float128 Q4_5D[NQ4_5D + 1] = {
500 L(3.049977232266999249626430127217988047453E-10),
501 L(7.120883230531035857746096928889676144099E-8),
502 L(6.301786064753734446784637919554359588859E-6),
503 L(2.762010530095069598480766869426308077192E-4),
504 L(6.572163250572867859316828886203406361251E-3),
505 L(8.752566114841221958200215255461843397776E-2),
506 L(6.487654992874805093499285311075289932664E-1),
507 L(2.576550017826654579451615283022812801435E0),
508 L(5.056392229924022835364779562707348096036E0),
509 L(4.179770081068251464907531367859072157773E0),
510 /* 1.000000000000000000000000000000000000000E0 */
511 };
512
513 /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x),
514 Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2))
515 Peak relative error 1.4e-36
516 0.25 <= 1/x <= 0.3125 */
517 #define NQ3r2_4N 10
518 static const _Float128 Q3r2_4N[NQ3r2_4N + 1] = {
519 L(6.126167301024815034423262653066023684411E-10),
520 L(1.043969327113173261820028225053598975128E-7),
521 L(6.592927270288697027757438170153763220190E-6),
522 L(2.009103660938497963095652951912071336730E-4),
523 L(3.220543385492643525985862356352195896964E-3),
524 L(2.774405975730545157543417650436941650990E-2),
525 L(1.258114008023826384487378016636555041129E-1),
526 L(2.811724258266902502344701449984698323860E-1),
527 L(2.691837665193548059322831687432415014067E-1),
528 L(7.949087384900985370683770525312735605034E-2),
529 L(1.229509543620976530030153018986910810747E-3),
530 };
531 #define NQ3r2_4D 9
532 static const _Float128 Q3r2_4D[NQ3r2_4D + 1] = {
533 L(8.364260446128475461539941389210166156568E-9),
534 L(1.451301850638956578622154585560759862764E-6),
535 L(9.431830010924603664244578867057141839463E-5),
536 L(3.004105101667433434196388593004526182741E-3),
537 L(5.148157397848271739710011717102773780221E-2),
538 L(4.901089301726939576055285374953887874895E-1),
539 L(2.581760991981709901216967665934142240346E0),
540 L(7.257105880775059281391729708630912791847E0),
541 L(1.006014717326362868007913423810737369312E1),
542 L(5.879416600465399514404064187445293212470E0),
543 /* 1.000000000000000000000000000000000000000E0*/
544 };
545
546 /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x),
547 Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2))
548 Peak relative error 3.8e-36
549 0.3125 <= 1/x <= 0.375 */
550 #define NQ2r7_3r2N 9
551 static const _Float128 Q2r7_3r2N[NQ2r7_3r2N + 1] = {
552 L(7.584861620402450302063691901886141875454E-8),
553 L(9.300939338814216296064659459966041794591E-6),
554 L(4.112108906197521696032158235392604947895E-4),
555 L(8.515168851578898791897038357239630654431E-3),
556 L(8.971286321017307400142720556749573229058E-2),
557 L(4.885856732902956303343015636331874194498E-1),
558 L(1.334506268733103291656253500506406045846E0),
559 L(1.681207956863028164179042145803851824654E0),
560 L(8.165042692571721959157677701625853772271E-1),
561 L(9.805848115375053300608712721986235900715E-2),
562 };
563 #define NQ2r7_3r2D 9
564 static const _Float128 Q2r7_3r2D[NQ2r7_3r2D + 1] = {
565 L(1.035586492113036586458163971239438078160E-6),
566 L(1.301999337731768381683593636500979713689E-4),
567 L(5.993695702564527062553071126719088859654E-3),
568 L(1.321184892887881883489141186815457808785E-1),
569 L(1.528766555485015021144963194165165083312E0),
570 L(9.561463309176490874525827051566494939295E0),
571 L(3.203719484883967351729513662089163356911E1),
572 L(5.497294687660930446641539152123568668447E1),
573 L(4.391158169390578768508675452986948391118E1),
574 L(1.347836630730048077907818943625789418378E1),
575 /* 1.000000000000000000000000000000000000000E0 */
576 };
577
578 /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x),
579 Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2))
580 Peak relative error 2.2e-35
581 0.375 <= 1/x <= 0.4375 */
582 #define NQ2r3_2r7N 9
583 static const _Float128 Q2r3_2r7N[NQ2r3_2r7N + 1] = {
584 L(4.455027774980750211349941766420190722088E-7),
585 L(4.031998274578520170631601850866780366466E-5),
586 L(1.273987274325947007856695677491340636339E-3),
587 L(1.818754543377448509897226554179659122873E-2),
588 L(1.266748858326568264126353051352269875352E-1),
589 L(4.327578594728723821137731555139472880414E-1),
590 L(6.892532471436503074928194969154192615359E-1),
591 L(4.490775818438716873422163588640262036506E-1),
592 L(8.649615949297322440032000346117031581572E-2),
593 L(7.261345286655345047417257611469066147561E-4),
594 };
595 #define NQ2r3_2r7D 8
596 static const _Float128 Q2r3_2r7D[NQ2r3_2r7D + 1] = {
597 L(6.082600739680555266312417978064954793142E-6),
598 L(5.693622538165494742945717226571441747567E-4),
599 L(1.901625907009092204458328768129666975975E-2),
600 L(2.958689532697857335456896889409923371570E-1),
601 L(2.343124711045660081603809437993368799568E0),
602 L(9.665894032187458293568704885528192804376E0),
603 L(2.035273104990617136065743426322454881353E1),
604 L(2.044102010478792896815088858740075165531E1),
605 L(8.445937177863155827844146643468706599304E0),
606 /* 1.000000000000000000000000000000000000000E0 */
607 };
608
609 /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x),
610 Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2))
611 Peak relative error 3.1e-36
612 0.4375 <= 1/x <= 0.5 */
613 #define NQ2_2r3N 9
614 static const _Float128 Q2_2r3N[NQ2_2r3N + 1] = {
615 L(2.817566786579768804844367382809101929314E-6),
616 L(2.122772176396691634147024348373539744935E-4),
617 L(5.501378031780457828919593905395747517585E-3),
618 L(6.355374424341762686099147452020466524659E-2),
619 L(3.539652320122661637429658698954748337223E-1),
620 L(9.571721066119617436343740541777014319695E-1),
621 L(1.196258777828426399432550698612171955305E0),
622 L(6.069388659458926158392384709893753793967E-1),
623 L(9.026746127269713176512359976978248763621E-2),
624 L(5.317668723070450235320878117210807236375E-4),
625 };
626 #define NQ2_2r3D 8
627 static const _Float128 Q2_2r3D[NQ2_2r3D + 1] = {
628 L(3.846924354014260866793741072933159380158E-5),
629 L(3.017562820057704325510067178327449946763E-3),
630 L(8.356305620686867949798885808540444210935E-2),
631 L(1.068314930499906838814019619594424586273E0),
632 L(6.900279623894821067017966573640732685233E0),
633 L(2.307667390886377924509090271780839563141E1),
634 L(3.921043465412723970791036825401273528513E1),
635 L(3.167569478939719383241775717095729233436E1),
636 L(1.051023841699200920276198346301543665909E1),
637 /* 1.000000000000000000000000000000000000000E0*/
638 };
639
640
641 /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
642
643 static _Float128
644 neval (_Float128 x, const _Float128 *p, int n)
645 {
646 _Float128 y;
647
648 p += n;
649 y = *p--;
650 do
651 {
652 y = y * x + *p--;
653 }
654 while (--n > 0);
655 return y;
656 }
657
658
659 /* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */
660
661 static _Float128
662 deval (_Float128 x, const _Float128 *p, int n)
663 {
664 _Float128 y;
665
666 p += n;
667 y = x + *p--;
668 do
669 {
670 y = y * x + *p--;
671 }
672 while (--n > 0);
673 return y;
674 }
675
676
677 /* Bessel function of the first kind, order zero. */
678
679 _Float128
680 __ieee754_j0l (_Float128 x)
681 {
682 _Float128 xx, xinv, z, p, q, c, s, cc, ss;
683
684 if (! isfinite (x))
685 {
686 if (x != x)
687 return x + x;
688 else
689 return 0;
690 }
691 if (x == 0)
692 return 1;
693
694 xx = fabsl (x);
695 if (xx <= 2)
696 {
697 if (xx < L(0x1p-57))
698 return 1;
699 /* 0 <= x <= 2 */
700 z = xx * xx;
701 p = z * z * neval (z, J0_2N, NJ0_2N) / deval (z, J0_2D, NJ0_2D);
702 p -= L(0.25) * z;
703 p += 1;
704 return p;
705 }
706
707 /* X = x - pi/4
708 cos(X) = cos(x) cos(pi/4) + sin(x) sin(pi/4)
709 = 1/sqrt(2) * (cos(x) + sin(x))
710 sin(X) = sin(x) cos(pi/4) - cos(x) sin(pi/4)
711 = 1/sqrt(2) * (sin(x) - cos(x))
712 sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
713 cf. Fdlibm. */
714 __sincosl (xx, &s, &c);
715 ss = s - c;
716 cc = s + c;
717 if (xx <= LDBL_MAX / 2)
718 {
719 z = -__cosl (xx + xx);
720 if ((s * c) < 0)
721 cc = z / ss;
722 else
723 ss = z / cc;
724 }
725
726 if (xx > L(0x1p256))
727 return ONEOSQPI * cc / sqrtl (xx);
728
729 xinv = 1 / xx;
730 z = xinv * xinv;
731 if (xinv <= 0.25)
732 {
733 if (xinv <= 0.125)
734 {
735 if (xinv <= 0.0625)
736 {
737 p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID);
738 q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID);
739 }
740 else
741 {
742 p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D);
743 q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D);
744 }
745 }
746 else if (xinv <= 0.1875)
747 {
748 p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D);
749 q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D);
750 }
751 else
752 {
753 p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D);
754 q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D);
755 }
756 } /* .25 */
757 else /* if (xinv <= 0.5) */
758 {
759 if (xinv <= 0.375)
760 {
761 if (xinv <= 0.3125)
762 {
763 p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D);
764 q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D);
765 }
766 else
767 {
768 p = neval (z, P2r7_3r2N, NP2r7_3r2N)
769 / deval (z, P2r7_3r2D, NP2r7_3r2D);
770 q = neval (z, Q2r7_3r2N, NQ2r7_3r2N)
771 / deval (z, Q2r7_3r2D, NQ2r7_3r2D);
772 }
773 }
774 else if (xinv <= 0.4375)
775 {
776 p = neval (z, P2r3_2r7N, NP2r3_2r7N)
777 / deval (z, P2r3_2r7D, NP2r3_2r7D);
778 q = neval (z, Q2r3_2r7N, NQ2r3_2r7N)
779 / deval (z, Q2r3_2r7D, NQ2r3_2r7D);
780 }
781 else
782 {
783 p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
784 q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
785 }
786 }
787 p = 1 + z * p;
788 q = z * xinv * q;
789 q = q - L(0.125) * xinv;
790 z = ONEOSQPI * (p * cc - q * ss) / sqrtl (xx);
791 return z;
792 }
793 libm_alias_finite (__ieee754_j0l, __j0l)
794
795
796 /* Y0(x) = 2/pi * log(x) * J0(x) + R(x^2)
797 Peak absolute error 1.7e-36 (relative where Y0 > 1)
798 0 <= x <= 2 */
799 #define NY0_2N 7
800 static const _Float128 Y0_2N[NY0_2N + 1] = {
801 L(-1.062023609591350692692296993537002558155E19),
802 L(2.542000883190248639104127452714966858866E19),
803 L(-1.984190771278515324281415820316054696545E18),
804 L(4.982586044371592942465373274440222033891E16),
805 L(-5.529326354780295177243773419090123407550E14),
806 L(3.013431465522152289279088265336861140391E12),
807 L(-7.959436160727126750732203098982718347785E9),
808 L(8.230845651379566339707130644134372793322E6),
809 };
810 #define NY0_2D 7
811 static const _Float128 Y0_2D[NY0_2D + 1] = {
812 L(1.438972634353286978700329883122253752192E20),
813 L(1.856409101981569254247700169486907405500E18),
814 L(1.219693352678218589553725579802986255614E16),
815 L(5.389428943282838648918475915779958097958E13),
816 L(1.774125762108874864433872173544743051653E11),
817 L(4.522104832545149534808218252434693007036E8),
818 L(8.872187401232943927082914504125234454930E5),
819 L(1.251945613186787532055610876304669413955E3),
820 /* 1.000000000000000000000000000000000000000E0 */
821 };
822
823 static const _Float128 U0 = L(-7.3804295108687225274343927948483016310862e-02);
824
825 /* Bessel function of the second kind, order zero. */
826
827 _Float128
828 __ieee754_y0l(_Float128 x)
829 {
830 _Float128 xx, xinv, z, p, q, c, s, cc, ss;
831
832 if (! isfinite (x))
833 return 1 / (x + x * x);
834 if (x <= 0)
835 {
836 if (x < 0)
837 return (zero / (zero * x));
838 return -1 / zero; /* -inf and divide by zero exception. */
839 }
840 xx = fabsl (x);
841 if (xx <= 0x1p-57)
842 return U0 + TWOOPI * __ieee754_logl (x);
843 if (xx <= 2)
844 {
845 /* 0 <= x <= 2 */
846 z = xx * xx;
847 p = neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D);
848 p = TWOOPI * __ieee754_logl (x) * __ieee754_j0l (x) + p;
849 return p;
850 }
851
852 /* X = x - pi/4
853 cos(X) = cos(x) cos(pi/4) + sin(x) sin(pi/4)
854 = 1/sqrt(2) * (cos(x) + sin(x))
855 sin(X) = sin(x) cos(pi/4) - cos(x) sin(pi/4)
856 = 1/sqrt(2) * (sin(x) - cos(x))
857 sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
858 cf. Fdlibm. */
859 __sincosl (x, &s, &c);
860 ss = s - c;
861 cc = s + c;
862 if (xx <= LDBL_MAX / 2)
863 {
864 z = -__cosl (x + x);
865 if ((s * c) < 0)
866 cc = z / ss;
867 else
868 ss = z / cc;
869 }
870
871 if (xx > L(0x1p256))
872 return ONEOSQPI * ss / sqrtl (x);
873
874 xinv = 1 / xx;
875 z = xinv * xinv;
876 if (xinv <= 0.25)
877 {
878 if (xinv <= 0.125)
879 {
880 if (xinv <= 0.0625)
881 {
882 p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID);
883 q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID);
884 }
885 else
886 {
887 p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D);
888 q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D);
889 }
890 }
891 else if (xinv <= 0.1875)
892 {
893 p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D);
894 q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D);
895 }
896 else
897 {
898 p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D);
899 q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D);
900 }
901 } /* .25 */
902 else /* if (xinv <= 0.5) */
903 {
904 if (xinv <= 0.375)
905 {
906 if (xinv <= 0.3125)
907 {
908 p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D);
909 q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D);
910 }
911 else
912 {
913 p = neval (z, P2r7_3r2N, NP2r7_3r2N)
914 / deval (z, P2r7_3r2D, NP2r7_3r2D);
915 q = neval (z, Q2r7_3r2N, NQ2r7_3r2N)
916 / deval (z, Q2r7_3r2D, NQ2r7_3r2D);
917 }
918 }
919 else if (xinv <= 0.4375)
920 {
921 p = neval (z, P2r3_2r7N, NP2r3_2r7N)
922 / deval (z, P2r3_2r7D, NP2r3_2r7D);
923 q = neval (z, Q2r3_2r7N, NQ2r3_2r7N)
924 / deval (z, Q2r3_2r7D, NQ2r3_2r7D);
925 }
926 else
927 {
928 p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
929 q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
930 }
931 }
932 p = 1 + z * p;
933 q = z * xinv * q;
934 q = q - L(0.125) * xinv;
935 z = ONEOSQPI * (p * ss + q * cc) / sqrtl (x);
936 return z;
937 }
938 libm_alias_finite (__ieee754_y0l, __y0l)