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