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