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
 790
 791 /* Y1(x) = 2/pi * (log(x) * J1(x) - 1/x) + x R(x^2)
 792    Peak relative error 6.2e-38
 793    0 <= x <= 2   */
 794 #define NY0_2N 7
 795 static long double Y0_2N[NY0_2N + 1] = {
 796   -6.804415404830253804408698161694720833249E19L,
 797   1.805450517967019908027153056150465849237E19L,
 798   -8.065747497063694098810419456383006737312E17L,
 799   1.401336667383028259295830955439028236299E16L,
 800   -1.171654432898137585000399489686629680230E14L,
 801   5.061267920943853732895341125243428129150E11L,
 802   -1.096677850566094204586208610960870217970E9L,
 803   9.541172044989995856117187515882879304461E5L,
 804 };
 805 #define NY0_2D 7
 806 static long double Y0_2D[NY0_2D + 1] = {
 807   3.470629591820267059538637461549677594549E20L,
 808   4.120796439009916326855848107545425217219E18L,
 809   2.477653371652018249749350657387030814542E16L,
 810   9.954678543353888958177169349272167762797E13L,
 811   2.957927997613630118216218290262851197754E11L,
 812   6.748421382188864486018861197614025972118E8L,
 813   1.173453425218010888004562071020305709319E6L,
 814   1.450335662961034949894009554536003377187E3L,
 815   /* 1.000000000000000000000000000000000000000E0 */
 816 };
 817
 818
 819 /* Bessel function of the second kind, order one.  */
 820
 821 long double
 822 __ieee754_y1l (long double x)
 823 {
 824   long double xx, xinv, z, p, q, c, s, cc, ss;
 825
 826   if (! __finitel (x))
 827     {
 828       if (x != x)
 829         return x;
 830       else
 831         return 0.0L;
 832     }
 833   if (x <= 0.0L)
 834     {
 835       if (x < 0.0L)
 836         return (zero / (zero * x));
 837       return -HUGE_VALL + x;
 838     }
 839   xx = fabsl (x);
 840   if (xx <= 2.0L)
 841     {
 842       /* 0 <= x <= 2 */
 843       z = xx * xx;
 844       p = xx * neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D);
 845       p = -TWOOPI / xx + p;
 846       p = TWOOPI * __ieee754_logl (x) * __ieee754_j1l (x) + p;
 847       return p;
 848     }
 849
 850   xinv = 1.0L / xx;
 851   z = xinv * xinv;
 852   if (xinv <= 0.25)
 853     {
 854       if (xinv <= 0.125)
 855         {
 856           if (xinv <= 0.0625)
 857             {
 858               p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID);
 859               q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID);
 860             }
 861           else
 862             {
 863               p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D);
 864               q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D);
 865             }
 866         }
 867       else if (xinv <= 0.1875)
 868         {
 869           p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D);
 870           q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D);
 871         }
 872       else
 873         {
 874           p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D);
 875           q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D);
 876         }
 877     }                           /* .25 */
 878   else /* if (xinv <= 0.5) */
 879     {
 880       if (xinv <= 0.375)
 881         {
 882           if (xinv <= 0.3125)
 883             {
 884               p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D);
 885               q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D);
 886             }
 887           else
 888             {
 889               p = neval (z, P2r7_3r2N, NP2r7_3r2N)
 890                   / deval (z, P2r7_3r2D, NP2r7_3r2D);
 891               q = neval (z, Q2r7_3r2N, NQ2r7_3r2N)
 892                   / deval (z, Q2r7_3r2D, NQ2r7_3r2D);
 893             }
 894         }
 895       else if (xinv <= 0.4375)
 896         {
 897           p = neval (z, P2r3_2r7N, NP2r3_2r7N)
 898               / deval (z, P2r3_2r7D, NP2r3_2r7D);
 899           q = neval (z, Q2r3_2r7N, NQ2r3_2r7N)
 900               / deval (z, Q2r3_2r7D, NQ2r3_2r7D);
 901         }
 902       else
 903         {
 904           p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
 905           q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
 906         }
 907     }
 908   p = 1.0L + z * p;
 909   q = z * q;
 910   q = q * xinv + 0.375L * xinv;
 911   /* X = x - 3 pi/4
 912      cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
 913      = 1/sqrt(2) * (-cos(x) + sin(x))
 914      sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
 915      = -1/sqrt(2) * (sin(x) + cos(x))
 916      cf. Fdlibm.  */
 917   __sincosl (xx, &s, &c);
 918   ss = -s - c;
 919   cc = s - c;
 920   z = __cosl (xx + xx);
 921   if ((s * c) > 0)
 922     cc = z / ss;
 923   else
 924     ss = z / cc;
 925   z = ONEOSQPI * (p * ss + q * cc) / __ieee754_sqrtl (xx);
 926   return z;
 927 }