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