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     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)