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