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     <http://www.gnu.org/licenses/>.  */
  93
  94 #include <math.h>
  95 #include <math_private.h>
  96 #include <float.h>
  97
  98 /* 1 / sqrt(pi) */
  99 static const _Float128 ONEOSQPI = L(5.6418958354775628694807945156077258584405E-1);
 100 /* 2 / pi */
 101 static const _Float128 TWOOPI = L(6.3661977236758134307553505349005744813784E-1);
 102 static 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
 108 static 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
 118 static 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
 133 static 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
 146 static 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
 164 static 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
 178 static 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
 197 static 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
 211 static 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
 229 static 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
 242 static 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
 260 static 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
 273 static 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
 291 static 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
 304 static 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
 321 static 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
 334 static 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
 351 static 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
 363 static 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
 381 static 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
 395 static 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
 414 static 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
 429 static 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
 450 static 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
 464 static 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
 484 static 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
 498 static 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
 517 static 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
 531 static 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
 550 static 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
 563 static 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
 582 static 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
 595 static 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
 613 static 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
 626 static 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
 642 static _Float128
 643 neval (_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
 660 static _Float128
 661 deval (_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 }
 792 strong_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
 799 static const _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
 810 static const _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
 822 static 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 }
 937 strong_alias (__ieee754_y0l, __y0l_finite)