l4/pkg/libquadmath/lib/contrib/math/j0q.c

   1 /*                                                      j0l.c
   2  *
   3  *      Bessel function of order zero
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * __float128 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  * __float128 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 "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.0Q;
 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;
 685       else
 686         return 0.0Q;
 687     }
 688   if (x == 0.0Q)
 689     return 1.0Q;
 690
 691   xx = fabsq (x);
 692   if (xx <= 2.0Q)
 693     {
 694       /* 0 <= x <= 2 */
 695       z = xx * xx;
 696       p = z * z * neval (z, J0_2N, NJ0_2N) / deval (z, J0_2D, NJ0_2D);
 697       p -= 0.25Q * z;
 698       p += 1.0Q;
 699       return p;
 700     }
 701
 702   xinv = 1.0Q / xx;
 703   z = xinv * xinv;
 704   if (xinv <= 0.25)
 705     {
 706       if (xinv <= 0.125)
 707         {
 708           if (xinv <= 0.0625)
 709             {
 710               p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID);
 711               q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID);
 712             }
 713           else
 714             {
 715               p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D);
 716               q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D);
 717             }
 718         }
 719       else if (xinv <= 0.1875)
 720         {
 721           p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D);
 722           q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D);
 723         }
 724       else
 725         {
 726           p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D);
 727           q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D);
 728         }
 729     }                           /* .25 */
 730   else /* if (xinv <= 0.5) */
 731     {
 732       if (xinv <= 0.375)
 733         {
 734           if (xinv <= 0.3125)
 735             {
 736               p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D);
 737               q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D);
 738             }
 739           else
 740             {
 741               p = neval (z, P2r7_3r2N, NP2r7_3r2N)
 742                   / deval (z, P2r7_3r2D, NP2r7_3r2D);
 743               q = neval (z, Q2r7_3r2N, NQ2r7_3r2N)
 744                   / deval (z, Q2r7_3r2D, NQ2r7_3r2D);
 745             }
 746         }
 747       else if (xinv <= 0.4375)
 748         {
 749           p = neval (z, P2r3_2r7N, NP2r3_2r7N)
 750               / deval (z, P2r3_2r7D, NP2r3_2r7D);
 751           q = neval (z, Q2r3_2r7N, NQ2r3_2r7N)
 752               / deval (z, Q2r3_2r7D, NQ2r3_2r7D);
 753         }
 754       else
 755         {
 756           p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
 757           q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
 758         }
 759     }
 760   p = 1.0Q + z * p;
 761   q = z * xinv * q;
 762   q = q - 0.125Q * xinv;
 763   /* X = x - pi/4
 764      cos(X) = cos(x) cos(pi/4) + sin(x) sin(pi/4)
 765      = 1/sqrt(2) * (cos(x) + sin(x))
 766      sin(X) = sin(x) cos(pi/4) - cos(x) sin(pi/4)
 767      = 1/sqrt(2) * (sin(x) - cos(x))
 768      sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
 769      cf. Fdlibm.  */
 770   sincosq (xx, &s, &c);
 771   ss = s - c;
 772   cc = s + c;
 773   z = - cosq (xx + xx);
 774   if ((s * c) < 0)
 775     cc = z / ss;
 776   else
 777     ss = z / cc;
 778   z = ONEOSQPI * (p * cc - q * ss) / sqrtq (xx);
 779   return z;
 780 }
 781
 782
 783 /* Y0(x) = 2/pi * log(x) * J0(x) + R(x^2)
 784    Peak absolute error 1.7e-36 (relative where Y0 > 1)
 785    0 <= x <= 2   */
 786 #define NY0_2N 7
 787 static __float128 Y0_2N[NY0_2N + 1] = {
 788  -1.062023609591350692692296993537002558155E19Q,
 789   2.542000883190248639104127452714966858866E19Q,
 790  -1.984190771278515324281415820316054696545E18Q,
 791   4.982586044371592942465373274440222033891E16Q,
 792  -5.529326354780295177243773419090123407550E14Q,
 793   3.013431465522152289279088265336861140391E12Q,
 794  -7.959436160727126750732203098982718347785E9Q,
 795   8.230845651379566339707130644134372793322E6Q,
 796 };
 797 #define NY0_2D 7
 798 static __float128 Y0_2D[NY0_2D + 1] = {
 799   1.438972634353286978700329883122253752192E20Q,
 800   1.856409101981569254247700169486907405500E18Q,
 801   1.219693352678218589553725579802986255614E16Q,
 802   5.389428943282838648918475915779958097958E13Q,
 803   1.774125762108874864433872173544743051653E11Q,
 804   4.522104832545149534808218252434693007036E8Q,
 805   8.872187401232943927082914504125234454930E5Q,
 806   1.251945613186787532055610876304669413955E3Q,
 807  /* 1.000000000000000000000000000000000000000E0 */
 808 };
 809
 810 static const long double U0 = -7.3804295108687225274343927948483016310862e-02Q;
 811
 812 /* Bessel function of the second kind, order zero.  */
 813
 814 __float128
 815 y0q (__float128 x)
 816 {
 817   __float128 xx, xinv, z, p, q, c, s, cc, ss;
 818
 819   if (! finiteq (x))
 820     {
 821       if (x != x)
 822         return x;
 823       else
 824         return 0.0Q;
 825     }
 826   if (x <= 0.0Q)
 827     {
 828       if (x < 0.0Q)
 829         return (zero / (zero * x));
 830       return -HUGE_VALQ + x;
 831     }
 832   xx = fabsq (x);
 833   if (xx <= 0x1p-57)
 834     return U0 + TWOOPI * logq (x);
 835   if (xx <= 2.0Q)
 836     {
 837       /* 0 <= x <= 2 */
 838       z = xx * xx;
 839       p = neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D);
 840       p = TWOOPI * logq (x) * j0q (x) + p;
 841       return p;
 842     }
 843
 844   xinv = 1.0Q / xx;
 845   z = xinv * xinv;
 846   if (xinv <= 0.25)
 847     {
 848       if (xinv <= 0.125)
 849         {
 850           if (xinv <= 0.0625)
 851             {
 852               p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID);
 853               q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID);
 854             }
 855           else
 856             {
 857               p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D);
 858               q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D);
 859             }
 860         }
 861       else if (xinv <= 0.1875)
 862         {
 863           p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D);
 864           q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D);
 865         }
 866       else
 867         {
 868           p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D);
 869           q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D);
 870         }
 871     }                           /* .25 */
 872   else /* if (xinv <= 0.5) */
 873     {
 874       if (xinv <= 0.375)
 875         {
 876           if (xinv <= 0.3125)
 877             {
 878               p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D);
 879               q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D);
 880             }
 881           else
 882             {
 883               p = neval (z, P2r7_3r2N, NP2r7_3r2N)
 884                   / deval (z, P2r7_3r2D, NP2r7_3r2D);
 885               q = neval (z, Q2r7_3r2N, NQ2r7_3r2N)
 886                   / deval (z, Q2r7_3r2D, NQ2r7_3r2D);
 887             }
 888         }
 889       else if (xinv <= 0.4375)
 890         {
 891           p = neval (z, P2r3_2r7N, NP2r3_2r7N)
 892               / deval (z, P2r3_2r7D, NP2r3_2r7D);
 893           q = neval (z, Q2r3_2r7N, NQ2r3_2r7N)
 894               / deval (z, Q2r3_2r7D, NQ2r3_2r7D);
 895         }
 896       else
 897         {
 898           p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
 899           q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
 900         }
 901     }
 902   p = 1.0Q + z * p;
 903   q = z * xinv * q;
 904   q = q - 0.125Q * xinv;
 905   /* X = x - pi/4
 906      cos(X) = cos(x) cos(pi/4) + sin(x) sin(pi/4)
 907      = 1/sqrt(2) * (cos(x) + sin(x))
 908      sin(X) = sin(x) cos(pi/4) - cos(x) sin(pi/4)
 909      = 1/sqrt(2) * (sin(x) - cos(x))
 910      sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
 911      cf. Fdlibm.  */
 912   sincosq (x, &s, &c);
 913   ss = s - c;
 914   cc = s + c;
 915   z = - cosq (x + x);
 916   if ((s * c) < 0)
 917     cc = z / ss;
 918   else
 919     ss = z / cc;
 920   z = ONEOSQPI * (p * ss + q * cc) / sqrtq (x);
 921   return z;
 922 }