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