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