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