sysdeps/ieee754/ldbl-128/e_j1l.c

   1 /*                                                      j1l.c
   2  *
   3  *      Bessel function of order one
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * long double x, y, j1l();
  10  *
  11  * y = j1l( x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Returns Bessel function of first kind, order one 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 is
  21  * J1(x) = .5x + x x^2 R(x^2)
  22  *
  23  * The second interval is further partitioned into eight equal segments
  24  * of 1/x.
  25  * J1(x) = sqrt(2/(pi x)) (P1(x) cos(X) - Q1(x) sin(X)),
  26  * X = x - 3 pi / 4,
  27  *
  28  * and the auxiliary functions are given by
  29  *
  30  * J1(x)cos(X) + Y1(x)sin(X) = sqrt( 2/(pi x)) P1(x),
  31  * P1(x) = 1 + 1/x^2 R(1/x^2)
  32  *
  33  * Y1(x)cos(X) - J1(x)sin(X) = sqrt( 2/(pi x)) Q1(x),
  34  * Q1(x) = 1/x (.375 + 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      2.8e-34      2.7e-35
  43  *
  44  *
  45  */
  46
  47 /*                                                      y1l.c
  48  *
  49  *      Bessel function of the second kind, order one
  50  *
  51  *
  52  *
  53  * SYNOPSIS:
  54  *
  55  * double x, y, y1l();
  56  *
  57  * y = y1l( x );
  58  *
  59  *
  60  *
  61  * DESCRIPTION:
  62  *
  63  * Returns Bessel function of the second kind, of order
  64  * one, of the argument.
  65  *
  66  * The domain is divided into two major intervals [0, 2] and
  67  * (2, infinity). In the first interval the rational approximation is
  68  * Y1(x) = 2/pi * (log(x) * J1(x) - 1/x) + x R(x^2) .
  69  * In the second interval the approximation is the same as for J1(x), and
  70  * Y1(x) = sqrt(2/(pi x)) (P1(x) sin(X) + Q1(x) cos(X)),
  71  * X = x - 3 pi / 4.
  72  *
  73  * ACCURACY:
  74  *
  75  *  Absolute error, when y0(x) < 1; else relative error:
  76  *
  77  * arithmetic   domain     # trials      peak         rms
  78  *    IEEE      0, 30       100000      2.7e-34     2.9e-35
  79  *
  80  */
  81
  82 /* Copyright 2001 by Stephen L. Moshier (moshier@na-net.onrl.gov).
  83
  84     This library is free software; you can redistribute it and/or
  85     modify it under the terms of the GNU Lesser General Public
  86     License as published by the Free Software Foundation; either
  87     version 2.1 of the License, or (at your option) any later version.
  88
  89     This library is distributed in the hope that it will be useful,
  90     but WITHOUT ANY WARRANTY; without even the implied warranty of
  91     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  92     Lesser General Public License for more details.
  93
  94     You should have received a copy of the GNU Lesser General Public
  95     License along with this library; if not, see
  96     <http://www.gnu.org/licenses/>.  */
  97
  98 #include <errno.h>
  99 #include <math.h>
 100 #include <math_private.h>
 101 #include <float.h>
 102
 103 /* 1 / sqrt(pi) */
 104 static const _Float128 ONEOSQPI = L(5.6418958354775628694807945156077258584405E-1);
 105 /* 2 / pi */
 106 static const _Float128 TWOOPI = L(6.3661977236758134307553505349005744813784E-1);
 107 static const _Float128 zero = 0;
 108
 109 /* J1(x) = .5x + x x^2 R(x^2)
 110    Peak relative error 1.9e-35
 111    0 <= x <= 2  */
 112 #define NJ0_2N 6
 113 static const _Float128 J0_2N[NJ0_2N + 1] = {
 114  L(-5.943799577386942855938508697619735179660E16),
 115   L(1.812087021305009192259946997014044074711E15),
 116  L(-2.761698314264509665075127515729146460895E13),
 117   L(2.091089497823600978949389109350658815972E11),
 118  L(-8.546413231387036372945453565654130054307E8),
 119   L(1.797229225249742247475464052741320612261E6),
 120  L(-1.559552840946694171346552770008812083969E3)
 121 };
 122 #define NJ0_2D 6
 123 static const _Float128 J0_2D[NJ0_2D + 1] = {
 124   L(9.510079323819108569501613916191477479397E17),
 125   L(1.063193817503280529676423936545854693915E16),
 126   L(5.934143516050192600795972192791775226920E13),
 127   L(2.168000911950620999091479265214368352883E11),
 128   L(5.673775894803172808323058205986256928794E8),
 129   L(1.080329960080981204840966206372671147224E6),
 130   L(1.411951256636576283942477881535283304912E3),
 131  /* 1.000000000000000000000000000000000000000E0L */
 132 };
 133
 134 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
 135    0 <= 1/x <= .0625
 136    Peak relative error 3.6e-36  */
 137 #define NP16_IN 9
 138 static const _Float128 P16_IN[NP16_IN + 1] = {
 139   L(5.143674369359646114999545149085139822905E-16),
 140   L(4.836645664124562546056389268546233577376E-13),
 141   L(1.730945562285804805325011561498453013673E-10),
 142   L(3.047976856147077889834905908605310585810E-8),
 143   L(2.855227609107969710407464739188141162386E-6),
 144   L(1.439362407936705484122143713643023998457E-4),
 145   L(3.774489768532936551500999699815873422073E-3),
 146   L(4.723962172984642566142399678920790598426E-2),
 147   L(2.359289678988743939925017240478818248735E-1),
 148   L(3.032580002220628812728954785118117124520E-1),
 149 };
 150 #define NP16_ID 9
 151 static const _Float128 P16_ID[NP16_ID + 1] = {
 152   L(4.389268795186898018132945193912677177553E-15),
 153   L(4.132671824807454334388868363256830961655E-12),
 154   L(1.482133328179508835835963635130894413136E-9),
 155   L(2.618941412861122118906353737117067376236E-7),
 156   L(2.467854246740858470815714426201888034270E-5),
 157   L(1.257192927368839847825938545925340230490E-3),
 158   L(3.362739031941574274949719324644120720341E-2),
 159   L(4.384458231338934105875343439265370178858E-1),
 160   L(2.412830809841095249170909628197264854651E0),
 161   L(4.176078204111348059102962617368214856874E0),
 162  /* 1.000000000000000000000000000000000000000E0 */
 163 };
 164
 165 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
 166     0.0625 <= 1/x <= 0.125
 167     Peak relative error 1.9e-36  */
 168 #define NP8_16N 11
 169 static const _Float128 P8_16N[NP8_16N + 1] = {
 170   L(2.984612480763362345647303274082071598135E-16),
 171   L(1.923651877544126103941232173085475682334E-13),
 172   L(4.881258879388869396043760693256024307743E-11),
 173   L(6.368866572475045408480898921866869811889E-9),
 174   L(4.684818344104910450523906967821090796737E-7),
 175   L(2.005177298271593587095982211091300382796E-5),
 176   L(4.979808067163957634120681477207147536182E-4),
 177   L(6.946005761642579085284689047091173581127E-3),
 178   L(5.074601112955765012750207555985299026204E-2),
 179   L(1.698599455896180893191766195194231825379E-1),
 180   L(1.957536905259237627737222775573623779638E-1),
 181   L(2.991314703282528370270179989044994319374E-2),
 182 };
 183 #define NP8_16D 10
 184 static const _Float128 P8_16D[NP8_16D + 1] = {
 185   L(2.546869316918069202079580939942463010937E-15),
 186   L(1.644650111942455804019788382157745229955E-12),
 187   L(4.185430770291694079925607420808011147173E-10),
 188   L(5.485331966975218025368698195861074143153E-8),
 189   L(4.062884421686912042335466327098932678905E-6),
 190   L(1.758139661060905948870523641319556816772E-4),
 191   L(4.445143889306356207566032244985607493096E-3),
 192   L(6.391901016293512632765621532571159071158E-2),
 193   L(4.933040207519900471177016015718145795434E-1),
 194   L(1.839144086168947712971630337250761842976E0),
 195   L(2.715120873995490920415616716916149586579E0),
 196  /* 1.000000000000000000000000000000000000000E0 */
 197 };
 198
 199 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
 200   0.125 <= 1/x <= 0.1875
 201   Peak relative error 1.3e-36  */
 202 #define NP5_8N 10
 203 static const _Float128 P5_8N[NP5_8N + 1] = {
 204   L(2.837678373978003452653763806968237227234E-12),
 205   L(9.726641165590364928442128579282742354806E-10),
 206   L(1.284408003604131382028112171490633956539E-7),
 207   L(8.524624695868291291250573339272194285008E-6),
 208   L(3.111516908953172249853673787748841282846E-4),
 209   L(6.423175156126364104172801983096596409176E-3),
 210   L(7.430220589989104581004416356260692450652E-2),
 211   L(4.608315409833682489016656279567605536619E-1),
 212   L(1.396870223510964882676225042258855977512E0),
 213   L(1.718500293904122365894630460672081526236E0),
 214   L(5.465927698800862172307352821870223855365E-1)
 215 };
 216 #define NP5_8D 10
 217 static const _Float128 P5_8D[NP5_8D + 1] = {
 218   L(2.421485545794616609951168511612060482715E-11),
 219   L(8.329862750896452929030058039752327232310E-9),
 220   L(1.106137992233383429630592081375289010720E-6),
 221   L(7.405786153760681090127497796448503306939E-5),
 222   L(2.740364785433195322492093333127633465227E-3),
 223   L(5.781246470403095224872243564165254652198E-2),
 224   L(6.927711353039742469918754111511109983546E-1),
 225   L(4.558679283460430281188304515922826156690E0),
 226   L(1.534468499844879487013168065728837900009E1),
 227   L(2.313927430889218597919624843161569422745E1),
 228   L(1.194506341319498844336768473218382828637E1),
 229  /* 1.000000000000000000000000000000000000000E0 */
 230 };
 231
 232 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
 233    Peak relative error 1.4e-36
 234    0.1875 <= 1/x <= 0.25  */
 235 #define NP4_5N 10
 236 static const _Float128 P4_5N[NP4_5N + 1] = {
 237   L(1.846029078268368685834261260420933914621E-10),
 238   L(3.916295939611376119377869680335444207768E-8),
 239   L(3.122158792018920627984597530935323997312E-6),
 240   L(1.218073444893078303994045653603392272450E-4),
 241   L(2.536420827983485448140477159977981844883E-3),
 242   L(2.883011322006690823959367922241169171315E-2),
 243   L(1.755255190734902907438042414495469810830E-1),
 244   L(5.379317079922628599870898285488723736599E-1),
 245   L(7.284904050194300773890303361501726561938E-1),
 246   L(3.270110346613085348094396323925000362813E-1),
 247   L(1.804473805689725610052078464951722064757E-2),
 248 };
 249 #define NP4_5D 9
 250 static const _Float128 P4_5D[NP4_5D + 1] = {
 251   L(1.575278146806816970152174364308980863569E-9),
 252   L(3.361289173657099516191331123405675054321E-7),
 253   L(2.704692281550877810424745289838790693708E-5),
 254   L(1.070854930483999749316546199273521063543E-3),
 255   L(2.282373093495295842598097265627962125411E-2),
 256   L(2.692025460665354148328762368240343249830E-1),
 257   L(1.739892942593664447220951225734811133759E0),
 258   L(5.890727576752230385342377570386657229324E0),
 259   L(9.517442287057841500750256954117735128153E0),
 260   L(6.100616353935338240775363403030137736013E0),
 261  /* 1.000000000000000000000000000000000000000E0 */
 262 };
 263
 264 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
 265    Peak relative error 3.0e-36
 266    0.25 <= 1/x <= 0.3125  */
 267 #define NP3r2_4N 9
 268 static const _Float128 P3r2_4N[NP3r2_4N + 1] = {
 269   L(8.240803130988044478595580300846665863782E-8),
 270   L(1.179418958381961224222969866406483744580E-5),
 271   L(6.179787320956386624336959112503824397755E-4),
 272   L(1.540270833608687596420595830747166658383E-2),
 273   L(1.983904219491512618376375619598837355076E-1),
 274   L(1.341465722692038870390470651608301155565E0),
 275   L(4.617865326696612898792238245990854646057E0),
 276   L(7.435574801812346424460233180412308000587E0),
 277   L(4.671327027414635292514599201278557680420E0),
 278   L(7.299530852495776936690976966995187714739E-1),
 279 };
 280 #define NP3r2_4D 9
 281 static const _Float128 P3r2_4D[NP3r2_4D + 1] = {
 282   L(7.032152009675729604487575753279187576521E-7),
 283   L(1.015090352324577615777511269928856742848E-4),
 284   L(5.394262184808448484302067955186308730620E-3),
 285   L(1.375291438480256110455809354836988584325E-1),
 286   L(1.836247144461106304788160919310404376670E0),
 287   L(1.314378564254376655001094503090935880349E1),
 288   L(4.957184590465712006934452500894672343488E1),
 289   L(9.287394244300647738855415178790263465398E1),
 290   L(7.652563275535900609085229286020552768399E1),
 291   L(2.147042473003074533150718117770093209096E1),
 292  /* 1.000000000000000000000000000000000000000E0 */
 293 };
 294
 295 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
 296    Peak relative error 1.0e-35
 297    0.3125 <= 1/x <= 0.375  */
 298 #define NP2r7_3r2N 9
 299 static const _Float128 P2r7_3r2N[NP2r7_3r2N + 1] = {
 300   L(4.599033469240421554219816935160627085991E-7),
 301   L(4.665724440345003914596647144630893997284E-5),
 302   L(1.684348845667764271596142716944374892756E-3),
 303   L(2.802446446884455707845985913454440176223E-2),
 304   L(2.321937586453963310008279956042545173930E-1),
 305   L(9.640277413988055668692438709376437553804E-1),
 306   L(1.911021064710270904508663334033003246028E0),
 307   L(1.600811610164341450262992138893970224971E0),
 308   L(4.266299218652587901171386591543457861138E-1),
 309   L(1.316470424456061252962568223251247207325E-2),
 310 };
 311 #define NP2r7_3r2D 8
 312 static const _Float128 P2r7_3r2D[NP2r7_3r2D + 1] = {
 313   L(3.924508608545520758883457108453520099610E-6),
 314   L(4.029707889408829273226495756222078039823E-4),
 315   L(1.484629715787703260797886463307469600219E-2),
 316   L(2.553136379967180865331706538897231588685E-1),
 317   L(2.229457223891676394409880026887106228740E0),
 318   L(1.005708903856384091956550845198392117318E1),
 319   L(2.277082659664386953166629360352385889558E1),
 320   L(2.384726835193630788249826630376533988245E1),
 321   L(9.700989749041320895890113781610939632410E0),
 322  /* 1.000000000000000000000000000000000000000E0 */
 323 };
 324
 325 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
 326    Peak relative error 1.7e-36
 327    0.3125 <= 1/x <= 0.4375  */
 328 #define NP2r3_2r7N 9
 329 static const _Float128 P2r3_2r7N[NP2r3_2r7N + 1] = {
 330   L(3.916766777108274628543759603786857387402E-6),
 331   L(3.212176636756546217390661984304645137013E-4),
 332   L(9.255768488524816445220126081207248947118E-3),
 333   L(1.214853146369078277453080641911700735354E-1),
 334   L(7.855163309847214136198449861311404633665E-1),
 335   L(2.520058073282978403655488662066019816540E0),
 336   L(3.825136484837545257209234285382183711466E0),
 337   L(2.432569427554248006229715163865569506873E0),
 338   L(4.877934835018231178495030117729800489743E-1),
 339   L(1.109902737860249670981355149101343427885E-2),
 340 };
 341 #define NP2r3_2r7D 8
 342 static const _Float128 P2r3_2r7D[NP2r3_2r7D + 1] = {
 343   L(3.342307880794065640312646341190547184461E-5),
 344   L(2.782182891138893201544978009012096558265E-3),
 345   L(8.221304931614200702142049236141249929207E-2),
 346   L(1.123728246291165812392918571987858010949E0),
 347   L(7.740482453652715577233858317133423434590E0),
 348   L(2.737624677567945952953322566311201919139E1),
 349   L(4.837181477096062403118304137851260715475E1),
 350   L(3.941098643468580791437772701093795299274E1),
 351   L(1.245821247166544627558323920382547533630E1),
 352  /* 1.000000000000000000000000000000000000000E0 */
 353 };
 354
 355 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
 356    Peak relative error 1.7e-35
 357    0.4375 <= 1/x <= 0.5  */
 358 #define NP2_2r3N 8
 359 static const _Float128 P2_2r3N[NP2_2r3N + 1] = {
 360   L(3.397930802851248553545191160608731940751E-4),
 361   L(2.104020902735482418784312825637833698217E-2),
 362   L(4.442291771608095963935342749477836181939E-1),
 363   L(4.131797328716583282869183304291833754967E0),
 364   L(1.819920169779026500146134832455189917589E1),
 365   L(3.781779616522937565300309684282401791291E1),
 366   L(3.459605449728864218972931220783543410347E1),
 367   L(1.173594248397603882049066603238568316561E1),
 368   L(9.455702270242780642835086549285560316461E-1),
 369 };
 370 #define NP2_2r3D 8
 371 static const _Float128 P2_2r3D[NP2_2r3D + 1] = {
 372   L(2.899568897241432883079888249845707400614E-3),
 373   L(1.831107138190848460767699919531132426356E-1),
 374   L(3.999350044057883839080258832758908825165E0),
 375   L(3.929041535867957938340569419874195303712E1),
 376   L(1.884245613422523323068802689915538908291E2),
 377   L(4.461469948819229734353852978424629815929E2),
 378   L(5.004998753999796821224085972610636347903E2),
 379   L(2.386342520092608513170837883757163414100E2),
 380   L(3.791322528149347975999851588922424189957E1),
 381  /* 1.000000000000000000000000000000000000000E0 */
 382 };
 383
 384 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
 385    Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
 386    Peak relative error 8.0e-36
 387    0 <= 1/x <= .0625  */
 388 #define NQ16_IN 10
 389 static const _Float128 Q16_IN[NQ16_IN + 1] = {
 390   L(-3.917420835712508001321875734030357393421E-18),
 391   L(-4.440311387483014485304387406538069930457E-15),
 392   L(-1.951635424076926487780929645954007139616E-12),
 393   L(-4.318256438421012555040546775651612810513E-10),
 394   L(-5.231244131926180765270446557146989238020E-8),
 395   L(-3.540072702902043752460711989234732357653E-6),
 396   L(-1.311017536555269966928228052917534882984E-4),
 397   L(-2.495184669674631806622008769674827575088E-3),
 398   L(-2.141868222987209028118086708697998506716E-2),
 399   L(-6.184031415202148901863605871197272650090E-2),
 400   L(-1.922298704033332356899546792898156493887E-2),
 401 };
 402 #define NQ16_ID 9
 403 static const _Float128 Q16_ID[NQ16_ID + 1] = {
 404   L(3.820418034066293517479619763498400162314E-17),
 405   L(4.340702810799239909648911373329149354911E-14),
 406   L(1.914985356383416140706179933075303538524E-11),
 407   L(4.262333682610888819476498617261895474330E-9),
 408   L(5.213481314722233980346462747902942182792E-7),
 409   L(3.585741697694069399299005316809954590558E-5),
 410   L(1.366513429642842006385029778105539457546E-3),
 411   L(2.745282599850704662726337474371355160594E-2),
 412   L(2.637644521611867647651200098449903330074E-1),
 413   L(1.006953426110765984590782655598680488746E0),
 414  /* 1.000000000000000000000000000000000000000E0 */
 415  };
 416
 417 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
 418    Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
 419    Peak relative error 1.9e-36
 420    0.0625 <= 1/x <= 0.125  */
 421 #define NQ8_16N 11
 422 static const _Float128 Q8_16N[NQ8_16N + 1] = {
 423   L(-2.028630366670228670781362543615221542291E-17),
 424   L(-1.519634620380959966438130374006858864624E-14),
 425   L(-4.540596528116104986388796594639405114524E-12),
 426   L(-7.085151756671466559280490913558388648274E-10),
 427   L(-6.351062671323970823761883833531546885452E-8),
 428   L(-3.390817171111032905297982523519503522491E-6),
 429   L(-1.082340897018886970282138836861233213972E-4),
 430   L(-2.020120801187226444822977006648252379508E-3),
 431   L(-2.093169910981725694937457070649605557555E-2),
 432   L(-1.092176538874275712359269481414448063393E-1),
 433   L(-2.374790947854765809203590474789108718733E-1),
 434   L(-1.365364204556573800719985118029601401323E-1),
 435 };
 436 #define NQ8_16D 11
 437 static const _Float128 Q8_16D[NQ8_16D + 1] = {
 438   L(1.978397614733632533581207058069628242280E-16),
 439   L(1.487361156806202736877009608336766720560E-13),
 440   L(4.468041406888412086042576067133365913456E-11),
 441   L(7.027822074821007443672290507210594648877E-9),
 442   L(6.375740580686101224127290062867976007374E-7),
 443   L(3.466887658320002225888644977076410421940E-5),
 444   L(1.138625640905289601186353909213719596986E-3),
 445   L(2.224470799470414663443449818235008486439E-2),
 446   L(2.487052928527244907490589787691478482358E-1),
 447   L(1.483927406564349124649083853892380899217E0),
 448   L(4.182773513276056975777258788903489507705E0),
 449   L(4.419665392573449746043880892524360870944E0),
 450  /* 1.000000000000000000000000000000000000000E0 */
 451 };
 452
 453 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
 454    Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
 455    Peak relative error 1.5e-35
 456    0.125 <= 1/x <= 0.1875  */
 457 #define NQ5_8N 10
 458 static const _Float128 Q5_8N[NQ5_8N + 1] = {
 459   L(-3.656082407740970534915918390488336879763E-13),
 460   L(-1.344660308497244804752334556734121771023E-10),
 461   L(-1.909765035234071738548629788698150760791E-8),
 462   L(-1.366668038160120210269389551283666716453E-6),
 463   L(-5.392327355984269366895210704976314135683E-5),
 464   L(-1.206268245713024564674432357634540343884E-3),
 465   L(-1.515456784370354374066417703736088291287E-2),
 466   L(-1.022454301137286306933217746545237098518E-1),
 467   L(-3.373438906472495080504907858424251082240E-1),
 468   L(-4.510782522110845697262323973549178453405E-1),
 469   L(-1.549000892545288676809660828213589804884E-1),
 470 };
 471 #define NQ5_8D 10
 472 static const _Float128 Q5_8D[NQ5_8D + 1] = {
 473   L(3.565550843359501079050699598913828460036E-12),
 474   L(1.321016015556560621591847454285330528045E-9),
 475   L(1.897542728662346479999969679234270605975E-7),
 476   L(1.381720283068706710298734234287456219474E-5),
 477   L(5.599248147286524662305325795203422873725E-4),
 478   L(1.305442352653121436697064782499122164843E-2),
 479   L(1.750234079626943298160445750078631894985E-1),
 480   L(1.311420542073436520965439883806946678491E0),
 481   L(5.162757689856842406744504211089724926650E0),
 482   L(9.527760296384704425618556332087850581308E0),
 483   L(6.604648207463236667912921642545100248584E0),
 484  /* 1.000000000000000000000000000000000000000E0 */
 485 };
 486
 487 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
 488    Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
 489    Peak relative error 1.3e-35
 490    0.1875 <= 1/x <= 0.25  */
 491 #define NQ4_5N 10
 492 static const _Float128 Q4_5N[NQ4_5N + 1] = {
 493   L(-4.079513568708891749424783046520200903755E-11),
 494   L(-9.326548104106791766891812583019664893311E-9),
 495   L(-8.016795121318423066292906123815687003356E-7),
 496   L(-3.372350544043594415609295225664186750995E-5),
 497   L(-7.566238665947967882207277686375417983917E-4),
 498   L(-9.248861580055565402130441618521591282617E-3),
 499   L(-6.033106131055851432267702948850231270338E-2),
 500   L(-1.966908754799996793730369265431584303447E-1),
 501   L(-2.791062741179964150755788226623462207560E-1),
 502   L(-1.255478605849190549914610121863534191666E-1),
 503   L(-4.320429862021265463213168186061696944062E-3),
 504 };
 505 #define NQ4_5D 9
 506 static const _Float128 Q4_5D[NQ4_5D + 1] = {
 507   L(3.978497042580921479003851216297330701056E-10),
 508   L(9.203304163828145809278568906420772246666E-8),
 509   L(8.059685467088175644915010485174545743798E-6),
 510   L(3.490187375993956409171098277561669167446E-4),
 511   L(8.189109654456872150100501732073810028829E-3),
 512   L(1.072572867311023640958725265762483033769E-1),
 513   L(7.790606862409960053675717185714576937994E-1),
 514   L(3.016049768232011196434185423512777656328E0),
 515   L(5.722963851442769787733717162314477949360E0),
 516   L(4.510527838428473279647251350931380867663E0),
 517  /* 1.000000000000000000000000000000000000000E0 */
 518 };
 519
 520 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
 521    Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
 522    Peak relative error 2.1e-35
 523    0.25 <= 1/x <= 0.3125  */
 524 #define NQ3r2_4N 9
 525 static const _Float128 Q3r2_4N[NQ3r2_4N + 1] = {
 526   L(-1.087480809271383885936921889040388133627E-8),
 527   L(-1.690067828697463740906962973479310170932E-6),
 528   L(-9.608064416995105532790745641974762550982E-5),
 529   L(-2.594198839156517191858208513873961837410E-3),
 530   L(-3.610954144421543968160459863048062977822E-2),
 531   L(-2.629866798251843212210482269563961685666E-1),
 532   L(-9.709186825881775885917984975685752956660E-1),
 533   L(-1.667521829918185121727268867619982417317E0),
 534   L(-1.109255082925540057138766105229900943501E0),
 535   L(-1.812932453006641348145049323713469043328E-1),
 536 };
 537 #define NQ3r2_4D 9
 538 static const _Float128 Q3r2_4D[NQ3r2_4D + 1] = {
 539   L(1.060552717496912381388763753841473407026E-7),
 540   L(1.676928002024920520786883649102388708024E-5),
 541   L(9.803481712245420839301400601140812255737E-4),
 542   L(2.765559874262309494758505158089249012930E-2),
 543   L(4.117921827792571791298862613287549140706E-1),
 544   L(3.323769515244751267093378361930279161413E0),
 545   L(1.436602494405814164724810151689705353670E1),
 546   L(3.163087869617098638064881410646782408297E1),
 547   L(3.198181264977021649489103980298349589419E1),
 548   L(1.203649258862068431199471076202897823272E1),
 549  /* 1.000000000000000000000000000000000000000E0  */
 550 };
 551
 552 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
 553    Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
 554    Peak relative error 1.6e-36
 555    0.3125 <= 1/x <= 0.375  */
 556 #define NQ2r7_3r2N 9
 557 static const _Float128 Q2r7_3r2N[NQ2r7_3r2N + 1] = {
 558   L(-1.723405393982209853244278760171643219530E-7),
 559   L(-2.090508758514655456365709712333460087442E-5),
 560   L(-9.140104013370974823232873472192719263019E-4),
 561   L(-1.871349499990714843332742160292474780128E-2),
 562   L(-1.948930738119938669637865956162512983416E-1),
 563   L(-1.048764684978978127908439526343174139788E0),
 564   L(-2.827714929925679500237476105843643064698E0),
 565   L(-3.508761569156476114276988181329773987314E0),
 566   L(-1.669332202790211090973255098624488308989E0),
 567   L(-1.930796319299022954013840684651016077770E-1),
 568 };
 569 #define NQ2r7_3r2D 9
 570 static const _Float128 Q2r7_3r2D[NQ2r7_3r2D + 1] = {
 571   L(1.680730662300831976234547482334347983474E-6),
 572   L(2.084241442440551016475972218719621841120E-4),
 573   L(9.445316642108367479043541702688736295579E-3),
 574   L(2.044637889456631896650179477133252184672E-1),
 575   L(2.316091982244297350829522534435350078205E0),
 576   L(1.412031891783015085196708811890448488865E1),
 577   L(4.583830154673223384837091077279595496149E1),
 578   L(7.549520609270909439885998474045974122261E1),
 579   L(5.697605832808113367197494052388203310638E1),
 580   L(1.601496240876192444526383314589371686234E1),
 581   /* 1.000000000000000000000000000000000000000E0 */
 582 };
 583
 584 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
 585    Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
 586    Peak relative error 9.5e-36
 587    0.375 <= 1/x <= 0.4375  */
 588 #define NQ2r3_2r7N 9
 589 static const _Float128 Q2r3_2r7N[NQ2r3_2r7N + 1] = {
 590   L(-8.603042076329122085722385914954878953775E-7),
 591   L(-7.701746260451647874214968882605186675720E-5),
 592   L(-2.407932004380727587382493696877569654271E-3),
 593   L(-3.403434217607634279028110636919987224188E-2),
 594   L(-2.348707332185238159192422084985713102877E-1),
 595   L(-7.957498841538254916147095255700637463207E-1),
 596   L(-1.258469078442635106431098063707934348577E0),
 597   L(-8.162415474676345812459353639449971369890E-1),
 598   L(-1.581783890269379690141513949609572806898E-1),
 599   L(-1.890595651683552228232308756569450822905E-3),
 600 };
 601 #define NQ2r3_2r7D 8
 602 static const _Float128 Q2r3_2r7D[NQ2r3_2r7D + 1] = {
 603   L(8.390017524798316921170710533381568175665E-6),
 604   L(7.738148683730826286477254659973968763659E-4),
 605   L(2.541480810958665794368759558791634341779E-2),
 606   L(3.878879789711276799058486068562386244873E-1),
 607   L(3.003783779325811292142957336802456109333E0),
 608   L(1.206480374773322029883039064575464497400E1),
 609   L(2.458414064785315978408974662900438351782E1),
 610   L(2.367237826273668567199042088835448715228E1),
 611   L(9.231451197519171090875569102116321676763E0),
 612  /* 1.000000000000000000000000000000000000000E0 */
 613 };
 614
 615 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
 616    Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
 617    Peak relative error 1.4e-36
 618    0.4375 <= 1/x <= 0.5  */
 619 #define NQ2_2r3N 9
 620 static const _Float128 Q2_2r3N[NQ2_2r3N + 1] = {
 621   L(-5.552507516089087822166822364590806076174E-6),
 622   L(-4.135067659799500521040944087433752970297E-4),
 623   L(-1.059928728869218962607068840646564457980E-2),
 624   L(-1.212070036005832342565792241385459023801E-1),
 625   L(-6.688350110633603958684302153362735625156E-1),
 626   L(-1.793587878197360221340277951304429821582E0),
 627   L(-2.225407682237197485644647380483725045326E0),
 628   L(-1.123402135458940189438898496348239744403E0),
 629   L(-1.679187241566347077204805190763597299805E-1),
 630   L(-1.458550613639093752909985189067233504148E-3),
 631 };
 632 #define NQ2_2r3D 8
 633 static const _Float128 Q2_2r3D[NQ2_2r3D + 1] = {
 634   L(5.415024336507980465169023996403597916115E-5),
 635   L(4.179246497380453022046357404266022870788E-3),
 636   L(1.136306384261959483095442402929502368598E-1),
 637   L(1.422640343719842213484515445393284072830E0),
 638   L(8.968786703393158374728850922289204805764E0),
 639   L(2.914542473339246127533384118781216495934E1),
 640   L(4.781605421020380669870197378210457054685E1),
 641   L(3.693865837171883152382820584714795072937E1),
 642   L(1.153220502744204904763115556224395893076E1),
 643   /* 1.000000000000000000000000000000000000000E0 */
 644 };
 645
 646
 647 /* Evaluate P[n] x^n  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
 648
 649 static _Float128
 650 neval (_Float128 x, const _Float128 *p, int n)
 651 {
 652   _Float128 y;
 653
 654   p += n;
 655   y = *p--;
 656   do
 657     {
 658       y = y * x + *p--;
 659     }
 660   while (--n > 0);
 661   return y;
 662 }
 663
 664
 665 /* Evaluate x^n+1  +  P[n] x^(n)  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
 666
 667 static _Float128
 668 deval (_Float128 x, const _Float128 *p, int n)
 669 {
 670   _Float128 y;
 671
 672   p += n;
 673   y = x + *p--;
 674   do
 675     {
 676       y = y * x + *p--;
 677     }
 678   while (--n > 0);
 679   return y;
 680 }
 681
 682
 683 /* Bessel function of the first kind, order one.  */
 684
 685 _Float128
 686 __ieee754_j1l (_Float128 x)
 687 {
 688   _Float128 xx, xinv, z, p, q, c, s, cc, ss;
 689
 690   if (! isfinite (x))
 691     {
 692       if (x != x)
 693         return x + x;
 694       else
 695         return 0;
 696     }
 697   if (x == 0)
 698     return x;
 699   xx = fabsl (x);
 700   if (xx <= L(0x1p-58))
 701     {
 702       _Float128 ret = x * L(0.5);
 703       math_check_force_underflow (ret);
 704       if (ret == 0)
 705         __set_errno (ERANGE);
 706       return ret;
 707     }
 708   if (xx <= 2)
 709     {
 710       /* 0 <= x <= 2 */
 711       z = xx * xx;
 712       p = xx * z * neval (z, J0_2N, NJ0_2N) / deval (z, J0_2D, NJ0_2D);
 713       p += L(0.5) * xx;
 714       if (x < 0)
 715         p = -p;
 716       return p;
 717     }
 718
 719   /* X = x - 3 pi/4
 720      cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
 721      = 1/sqrt(2) * (-cos(x) + sin(x))
 722      sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
 723      = -1/sqrt(2) * (sin(x) + cos(x))
 724      cf. Fdlibm.  */
 725   __sincosl (xx, &s, &c);
 726   ss = -s - c;
 727   cc = s - c;
 728   if (xx <= LDBL_MAX / 2)
 729     {
 730       z = __cosl (xx + xx);
 731       if ((s * c) > 0)
 732         cc = z / ss;
 733       else
 734         ss = z / cc;
 735     }
 736
 737   if (xx > L(0x1p256))
 738     {
 739       z = ONEOSQPI * cc / __ieee754_sqrtl (xx);
 740       if (x < 0)
 741         z = -z;
 742       return z;
 743     }
 744
 745   xinv = 1 / xx;
 746   z = xinv * xinv;
 747   if (xinv <= 0.25)
 748     {
 749       if (xinv <= 0.125)
 750         {
 751           if (xinv <= 0.0625)
 752             {
 753               p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID);
 754               q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID);
 755             }
 756           else
 757             {
 758               p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D);
 759               q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D);
 760             }
 761         }
 762       else if (xinv <= 0.1875)
 763         {
 764           p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D);
 765           q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D);
 766         }
 767       else
 768         {
 769           p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D);
 770           q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D);
 771         }
 772     }                           /* .25 */
 773   else /* if (xinv <= 0.5) */
 774     {
 775       if (xinv <= 0.375)
 776         {
 777           if (xinv <= 0.3125)
 778             {
 779               p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D);
 780               q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D);
 781             }
 782           else
 783             {
 784               p = neval (z, P2r7_3r2N, NP2r7_3r2N)
 785                   / deval (z, P2r7_3r2D, NP2r7_3r2D);
 786               q = neval (z, Q2r7_3r2N, NQ2r7_3r2N)
 787                   / deval (z, Q2r7_3r2D, NQ2r7_3r2D);
 788             }
 789         }
 790       else if (xinv <= 0.4375)
 791         {
 792           p = neval (z, P2r3_2r7N, NP2r3_2r7N)
 793               / deval (z, P2r3_2r7D, NP2r3_2r7D);
 794           q = neval (z, Q2r3_2r7N, NQ2r3_2r7N)
 795               / deval (z, Q2r3_2r7D, NQ2r3_2r7D);
 796         }
 797       else
 798         {
 799           p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
 800           q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
 801         }
 802     }
 803   p = 1 + z * p;
 804   q = z * q;
 805   q = q * xinv + L(0.375) * xinv;
 806   z = ONEOSQPI * (p * cc - q * ss) / __ieee754_sqrtl (xx);
 807   if (x < 0)
 808     z = -z;
 809   return z;
 810 }
 811 strong_alias (__ieee754_j1l, __j1l_finite)
 812
 813
 814 /* Y1(x) = 2/pi * (log(x) * J1(x) - 1/x) + x R(x^2)
 815    Peak relative error 6.2e-38
 816    0 <= x <= 2   */
 817 #define NY0_2N 7
 818 static _Float128 Y0_2N[NY0_2N + 1] = {
 819   L(-6.804415404830253804408698161694720833249E19),
 820   L(1.805450517967019908027153056150465849237E19),
 821   L(-8.065747497063694098810419456383006737312E17),
 822   L(1.401336667383028259295830955439028236299E16),
 823   L(-1.171654432898137585000399489686629680230E14),
 824   L(5.061267920943853732895341125243428129150E11),
 825   L(-1.096677850566094204586208610960870217970E9),
 826   L(9.541172044989995856117187515882879304461E5),
 827 };
 828 #define NY0_2D 7
 829 static _Float128 Y0_2D[NY0_2D + 1] = {
 830   L(3.470629591820267059538637461549677594549E20),
 831   L(4.120796439009916326855848107545425217219E18),
 832   L(2.477653371652018249749350657387030814542E16),
 833   L(9.954678543353888958177169349272167762797E13),
 834   L(2.957927997613630118216218290262851197754E11),
 835   L(6.748421382188864486018861197614025972118E8),
 836   L(1.173453425218010888004562071020305709319E6),
 837   L(1.450335662961034949894009554536003377187E3),
 838   /* 1.000000000000000000000000000000000000000E0 */
 839 };
 840
 841
 842 /* Bessel function of the second kind, order one.  */
 843
 844 _Float128
 845 __ieee754_y1l (_Float128 x)
 846 {
 847   _Float128 xx, xinv, z, p, q, c, s, cc, ss;
 848
 849   if (! isfinite (x))
 850     return 1 / (x + x * x);
 851   if (x <= 0)
 852     {
 853       if (x < 0)
 854         return (zero / (zero * x));
 855       return -1 / zero; /* -inf and divide by zero exception.  */
 856     }
 857   xx = fabsl (x);
 858   if (xx <= 0x1p-114)
 859     {
 860       z = -TWOOPI / x;
 861       if (isinf (z))
 862         __set_errno (ERANGE);
 863       return z;
 864     }
 865   if (xx <= 2)
 866     {
 867       /* 0 <= x <= 2 */
 868       SET_RESTORE_ROUNDL (FE_TONEAREST);
 869       z = xx * xx;
 870       p = xx * neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D);
 871       p = -TWOOPI / xx + p;
 872       p = TWOOPI * __ieee754_logl (x) * __ieee754_j1l (x) + p;
 873       return p;
 874     }
 875
 876   /* X = x - 3 pi/4
 877      cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
 878      = 1/sqrt(2) * (-cos(x) + sin(x))
 879      sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
 880      = -1/sqrt(2) * (sin(x) + cos(x))
 881      cf. Fdlibm.  */
 882   __sincosl (xx, &s, &c);
 883   ss = -s - c;
 884   cc = s - c;
 885   if (xx <= LDBL_MAX / 2)
 886     {
 887       z = __cosl (xx + xx);
 888       if ((s * c) > 0)
 889         cc = z / ss;
 890       else
 891         ss = z / cc;
 892     }
 893
 894   if (xx > L(0x1p256))
 895     return ONEOSQPI * ss / __ieee754_sqrtl (xx);
 896
 897   xinv = 1 / xx;
 898   z = xinv * xinv;
 899   if (xinv <= 0.25)
 900     {
 901       if (xinv <= 0.125)
 902         {
 903           if (xinv <= 0.0625)
 904             {
 905               p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID);
 906               q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID);
 907             }
 908           else
 909             {
 910               p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D);
 911               q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D);
 912             }
 913         }
 914       else if (xinv <= 0.1875)
 915         {
 916           p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D);
 917           q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D);
 918         }
 919       else
 920         {
 921           p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D);
 922           q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D);
 923         }
 924     }                           /* .25 */
 925   else /* if (xinv <= 0.5) */
 926     {
 927       if (xinv <= 0.375)
 928         {
 929           if (xinv <= 0.3125)
 930             {
 931               p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D);
 932               q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D);
 933             }
 934           else
 935             {
 936               p = neval (z, P2r7_3r2N, NP2r7_3r2N)
 937                   / deval (z, P2r7_3r2D, NP2r7_3r2D);
 938               q = neval (z, Q2r7_3r2N, NQ2r7_3r2N)
 939                   / deval (z, Q2r7_3r2D, NQ2r7_3r2D);
 940             }
 941         }
 942       else if (xinv <= 0.4375)
 943         {
 944           p = neval (z, P2r3_2r7N, NP2r3_2r7N)
 945               / deval (z, P2r3_2r7D, NP2r3_2r7D);
 946           q = neval (z, Q2r3_2r7N, NQ2r3_2r7N)
 947               / deval (z, Q2r3_2r7D, NQ2r3_2r7D);
 948         }
 949       else
 950         {
 951           p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
 952           q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
 953         }
 954     }
 955   p = 1 + z * p;
 956   q = z * q;
 957   q = q * xinv + L(0.375) * xinv;
 958   z = ONEOSQPI * (p * ss + q * cc) / __ieee754_sqrtl (xx);
 959   return z;
 960 }
 961 strong_alias (__ieee754_y1l, __y1l_finite)