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