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