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