libquadmath/math/j1q.c

   1 /*                                                      j1l.c
   2  *
   3  *      Bessel function of order one
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * __float128 x, y, j1q();
  10  *
  11  * y = j1q( 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  * __float128, y, y1q();
  56  *
  57  * y = y1q( 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 <errno.h>
  99 #include "quadmath-imp.h"
 100
 101 /* 1 / sqrt(pi) */
 102 static const __float128 ONEOSQPI = 5.6418958354775628694807945156077258584405E-1Q;
 103 /* 2 / pi */
 104 static const __float128 TWOOPI = 6.3661977236758134307553505349005744813784E-1Q;
 105 static const __float128 zero = 0.0Q;
 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 __float128 J0_2N[NJ0_2N + 1] = {
 112  -5.943799577386942855938508697619735179660E16Q,
 113   1.812087021305009192259946997014044074711E15Q,
 114  -2.761698314264509665075127515729146460895E13Q,
 115   2.091089497823600978949389109350658815972E11Q,
 116  -8.546413231387036372945453565654130054307E8Q,
 117   1.797229225249742247475464052741320612261E6Q,
 118  -1.559552840946694171346552770008812083969E3Q
 119 };
 120 #define NJ0_2D 6
 121 static const __float128 J0_2D[NJ0_2D + 1] = {
 122   9.510079323819108569501613916191477479397E17Q,
 123   1.063193817503280529676423936545854693915E16Q,
 124   5.934143516050192600795972192791775226920E13Q,
 125   2.168000911950620999091479265214368352883E11Q,
 126   5.673775894803172808323058205986256928794E8Q,
 127   1.080329960080981204840966206372671147224E6Q,
 128   1.411951256636576283942477881535283304912E3Q,
 129  /* 1.000000000000000000000000000000000000000E0Q */
 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 __float128 P16_IN[NP16_IN + 1] = {
 137   5.143674369359646114999545149085139822905E-16Q,
 138   4.836645664124562546056389268546233577376E-13Q,
 139   1.730945562285804805325011561498453013673E-10Q,
 140   3.047976856147077889834905908605310585810E-8Q,
 141   2.855227609107969710407464739188141162386E-6Q,
 142   1.439362407936705484122143713643023998457E-4Q,
 143   3.774489768532936551500999699815873422073E-3Q,
 144   4.723962172984642566142399678920790598426E-2Q,
 145   2.359289678988743939925017240478818248735E-1Q,
 146   3.032580002220628812728954785118117124520E-1Q,
 147 };
 148 #define NP16_ID 9
 149 static const __float128 P16_ID[NP16_ID + 1] = {
 150   4.389268795186898018132945193912677177553E-15Q,
 151   4.132671824807454334388868363256830961655E-12Q,
 152   1.482133328179508835835963635130894413136E-9Q,
 153   2.618941412861122118906353737117067376236E-7Q,
 154   2.467854246740858470815714426201888034270E-5Q,
 155   1.257192927368839847825938545925340230490E-3Q,
 156   3.362739031941574274949719324644120720341E-2Q,
 157   4.384458231338934105875343439265370178858E-1Q,
 158   2.412830809841095249170909628197264854651E0Q,
 159   4.176078204111348059102962617368214856874E0Q,
 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 __float128 P8_16N[NP8_16N + 1] = {
 168   2.984612480763362345647303274082071598135E-16Q,
 169   1.923651877544126103941232173085475682334E-13Q,
 170   4.881258879388869396043760693256024307743E-11Q,
 171   6.368866572475045408480898921866869811889E-9Q,
 172   4.684818344104910450523906967821090796737E-7Q,
 173   2.005177298271593587095982211091300382796E-5Q,
 174   4.979808067163957634120681477207147536182E-4Q,
 175   6.946005761642579085284689047091173581127E-3Q,
 176   5.074601112955765012750207555985299026204E-2Q,
 177   1.698599455896180893191766195194231825379E-1Q,
 178   1.957536905259237627737222775573623779638E-1Q,
 179   2.991314703282528370270179989044994319374E-2Q,
 180 };
 181 #define NP8_16D 10
 182 static const __float128 P8_16D[NP8_16D + 1] = {
 183   2.546869316918069202079580939942463010937E-15Q,
 184   1.644650111942455804019788382157745229955E-12Q,
 185   4.185430770291694079925607420808011147173E-10Q,
 186   5.485331966975218025368698195861074143153E-8Q,
 187   4.062884421686912042335466327098932678905E-6Q,
 188   1.758139661060905948870523641319556816772E-4Q,
 189   4.445143889306356207566032244985607493096E-3Q,
 190   6.391901016293512632765621532571159071158E-2Q,
 191   4.933040207519900471177016015718145795434E-1Q,
 192   1.839144086168947712971630337250761842976E0Q,
 193   2.715120873995490920415616716916149586579E0Q,
 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 __float128 P5_8N[NP5_8N + 1] = {
 202   2.837678373978003452653763806968237227234E-12Q,
 203   9.726641165590364928442128579282742354806E-10Q,
 204   1.284408003604131382028112171490633956539E-7Q,
 205   8.524624695868291291250573339272194285008E-6Q,
 206   3.111516908953172249853673787748841282846E-4Q,
 207   6.423175156126364104172801983096596409176E-3Q,
 208   7.430220589989104581004416356260692450652E-2Q,
 209   4.608315409833682489016656279567605536619E-1Q,
 210   1.396870223510964882676225042258855977512E0Q,
 211   1.718500293904122365894630460672081526236E0Q,
 212   5.465927698800862172307352821870223855365E-1Q
 213 };
 214 #define NP5_8D 10
 215 static const __float128 P5_8D[NP5_8D + 1] = {
 216   2.421485545794616609951168511612060482715E-11Q,
 217   8.329862750896452929030058039752327232310E-9Q,
 218   1.106137992233383429630592081375289010720E-6Q,
 219   7.405786153760681090127497796448503306939E-5Q,
 220   2.740364785433195322492093333127633465227E-3Q,
 221   5.781246470403095224872243564165254652198E-2Q,
 222   6.927711353039742469918754111511109983546E-1Q,
 223   4.558679283460430281188304515922826156690E0Q,
 224   1.534468499844879487013168065728837900009E1Q,
 225   2.313927430889218597919624843161569422745E1Q,
 226   1.194506341319498844336768473218382828637E1Q,
 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 __float128 P4_5N[NP4_5N + 1] = {
 235   1.846029078268368685834261260420933914621E-10Q,
 236   3.916295939611376119377869680335444207768E-8Q,
 237   3.122158792018920627984597530935323997312E-6Q,
 238   1.218073444893078303994045653603392272450E-4Q,
 239   2.536420827983485448140477159977981844883E-3Q,
 240   2.883011322006690823959367922241169171315E-2Q,
 241   1.755255190734902907438042414495469810830E-1Q,
 242   5.379317079922628599870898285488723736599E-1Q,
 243   7.284904050194300773890303361501726561938E-1Q,
 244   3.270110346613085348094396323925000362813E-1Q,
 245   1.804473805689725610052078464951722064757E-2Q,
 246 };
 247 #define NP4_5D 9
 248 static const __float128 P4_5D[NP4_5D + 1] = {
 249   1.575278146806816970152174364308980863569E-9Q,
 250   3.361289173657099516191331123405675054321E-7Q,
 251   2.704692281550877810424745289838790693708E-5Q,
 252   1.070854930483999749316546199273521063543E-3Q,
 253   2.282373093495295842598097265627962125411E-2Q,
 254   2.692025460665354148328762368240343249830E-1Q,
 255   1.739892942593664447220951225734811133759E0Q,
 256   5.890727576752230385342377570386657229324E0Q,
 257   9.517442287057841500750256954117735128153E0Q,
 258   6.100616353935338240775363403030137736013E0Q,
 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 __float128 P3r2_4N[NP3r2_4N + 1] = {
 267   8.240803130988044478595580300846665863782E-8Q,
 268   1.179418958381961224222969866406483744580E-5Q,
 269   6.179787320956386624336959112503824397755E-4Q,
 270   1.540270833608687596420595830747166658383E-2Q,
 271   1.983904219491512618376375619598837355076E-1Q,
 272   1.341465722692038870390470651608301155565E0Q,
 273   4.617865326696612898792238245990854646057E0Q,
 274   7.435574801812346424460233180412308000587E0Q,
 275   4.671327027414635292514599201278557680420E0Q,
 276   7.299530852495776936690976966995187714739E-1Q,
 277 };
 278 #define NP3r2_4D 9
 279 static const __float128 P3r2_4D[NP3r2_4D + 1] = {
 280   7.032152009675729604487575753279187576521E-7Q,
 281   1.015090352324577615777511269928856742848E-4Q,
 282   5.394262184808448484302067955186308730620E-3Q,
 283   1.375291438480256110455809354836988584325E-1Q,
 284   1.836247144461106304788160919310404376670E0Q,
 285   1.314378564254376655001094503090935880349E1Q,
 286   4.957184590465712006934452500894672343488E1Q,
 287   9.287394244300647738855415178790263465398E1Q,
 288   7.652563275535900609085229286020552768399E1Q,
 289   2.147042473003074533150718117770093209096E1Q,
 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 __float128 P2r7_3r2N[NP2r7_3r2N + 1] = {
 298   4.599033469240421554219816935160627085991E-7Q,
 299   4.665724440345003914596647144630893997284E-5Q,
 300   1.684348845667764271596142716944374892756E-3Q,
 301   2.802446446884455707845985913454440176223E-2Q,
 302   2.321937586453963310008279956042545173930E-1Q,
 303   9.640277413988055668692438709376437553804E-1Q,
 304   1.911021064710270904508663334033003246028E0Q,
 305   1.600811610164341450262992138893970224971E0Q,
 306   4.266299218652587901171386591543457861138E-1Q,
 307   1.316470424456061252962568223251247207325E-2Q,
 308 };
 309 #define NP2r7_3r2D 8
 310 static const __float128 P2r7_3r2D[NP2r7_3r2D + 1] = {
 311   3.924508608545520758883457108453520099610E-6Q,
 312   4.029707889408829273226495756222078039823E-4Q,
 313   1.484629715787703260797886463307469600219E-2Q,
 314   2.553136379967180865331706538897231588685E-1Q,
 315   2.229457223891676394409880026887106228740E0Q,
 316   1.005708903856384091956550845198392117318E1Q,
 317   2.277082659664386953166629360352385889558E1Q,
 318   2.384726835193630788249826630376533988245E1Q,
 319   9.700989749041320895890113781610939632410E0Q,
 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 __float128 P2r3_2r7N[NP2r3_2r7N + 1] = {
 328   3.916766777108274628543759603786857387402E-6Q,
 329   3.212176636756546217390661984304645137013E-4Q,
 330   9.255768488524816445220126081207248947118E-3Q,
 331   1.214853146369078277453080641911700735354E-1Q,
 332   7.855163309847214136198449861311404633665E-1Q,
 333   2.520058073282978403655488662066019816540E0Q,
 334   3.825136484837545257209234285382183711466E0Q,
 335   2.432569427554248006229715163865569506873E0Q,
 336   4.877934835018231178495030117729800489743E-1Q,
 337   1.109902737860249670981355149101343427885E-2Q,
 338 };
 339 #define NP2r3_2r7D 8
 340 static const __float128 P2r3_2r7D[NP2r3_2r7D + 1] = {
 341   3.342307880794065640312646341190547184461E-5Q,
 342   2.782182891138893201544978009012096558265E-3Q,
 343   8.221304931614200702142049236141249929207E-2Q,
 344   1.123728246291165812392918571987858010949E0Q,
 345   7.740482453652715577233858317133423434590E0Q,
 346   2.737624677567945952953322566311201919139E1Q,
 347   4.837181477096062403118304137851260715475E1Q,
 348   3.941098643468580791437772701093795299274E1Q,
 349   1.245821247166544627558323920382547533630E1Q,
 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 __float128 P2_2r3N[NP2_2r3N + 1] = {
 358   3.397930802851248553545191160608731940751E-4Q,
 359   2.104020902735482418784312825637833698217E-2Q,
 360   4.442291771608095963935342749477836181939E-1Q,
 361   4.131797328716583282869183304291833754967E0Q,
 362   1.819920169779026500146134832455189917589E1Q,
 363   3.781779616522937565300309684282401791291E1Q,
 364   3.459605449728864218972931220783543410347E1Q,
 365   1.173594248397603882049066603238568316561E1Q,
 366   9.455702270242780642835086549285560316461E-1Q,
 367 };
 368 #define NP2_2r3D 8
 369 static const __float128 P2_2r3D[NP2_2r3D + 1] = {
 370   2.899568897241432883079888249845707400614E-3Q,
 371   1.831107138190848460767699919531132426356E-1Q,
 372   3.999350044057883839080258832758908825165E0Q,
 373   3.929041535867957938340569419874195303712E1Q,
 374   1.884245613422523323068802689915538908291E2Q,
 375   4.461469948819229734353852978424629815929E2Q,
 376   5.004998753999796821224085972610636347903E2Q,
 377   2.386342520092608513170837883757163414100E2Q,
 378   3.791322528149347975999851588922424189957E1Q,
 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 __float128 Q16_IN[NQ16_IN + 1] = {
 388   -3.917420835712508001321875734030357393421E-18Q,
 389   -4.440311387483014485304387406538069930457E-15Q,
 390   -1.951635424076926487780929645954007139616E-12Q,
 391   -4.318256438421012555040546775651612810513E-10Q,
 392   -5.231244131926180765270446557146989238020E-8Q,
 393   -3.540072702902043752460711989234732357653E-6Q,
 394   -1.311017536555269966928228052917534882984E-4Q,
 395   -2.495184669674631806622008769674827575088E-3Q,
 396   -2.141868222987209028118086708697998506716E-2Q,
 397   -6.184031415202148901863605871197272650090E-2Q,
 398   -1.922298704033332356899546792898156493887E-2Q,
 399 };
 400 #define NQ16_ID 9
 401 static const __float128 Q16_ID[NQ16_ID + 1] = {
 402   3.820418034066293517479619763498400162314E-17Q,
 403   4.340702810799239909648911373329149354911E-14Q,
 404   1.914985356383416140706179933075303538524E-11Q,
 405   4.262333682610888819476498617261895474330E-9Q,
 406   5.213481314722233980346462747902942182792E-7Q,
 407   3.585741697694069399299005316809954590558E-5Q,
 408   1.366513429642842006385029778105539457546E-3Q,
 409   2.745282599850704662726337474371355160594E-2Q,
 410   2.637644521611867647651200098449903330074E-1Q,
 411   1.006953426110765984590782655598680488746E0Q,
 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 __float128 Q8_16N[NQ8_16N + 1] = {
 421   -2.028630366670228670781362543615221542291E-17Q,
 422   -1.519634620380959966438130374006858864624E-14Q,
 423   -4.540596528116104986388796594639405114524E-12Q,
 424   -7.085151756671466559280490913558388648274E-10Q,
 425   -6.351062671323970823761883833531546885452E-8Q,
 426   -3.390817171111032905297982523519503522491E-6Q,
 427   -1.082340897018886970282138836861233213972E-4Q,
 428   -2.020120801187226444822977006648252379508E-3Q,
 429   -2.093169910981725694937457070649605557555E-2Q,
 430   -1.092176538874275712359269481414448063393E-1Q,
 431   -2.374790947854765809203590474789108718733E-1Q,
 432   -1.365364204556573800719985118029601401323E-1Q,
 433 };
 434 #define NQ8_16D 11
 435 static const __float128 Q8_16D[NQ8_16D + 1] = {
 436   1.978397614733632533581207058069628242280E-16Q,
 437   1.487361156806202736877009608336766720560E-13Q,
 438   4.468041406888412086042576067133365913456E-11Q,
 439   7.027822074821007443672290507210594648877E-9Q,
 440   6.375740580686101224127290062867976007374E-7Q,
 441   3.466887658320002225888644977076410421940E-5Q,
 442   1.138625640905289601186353909213719596986E-3Q,
 443   2.224470799470414663443449818235008486439E-2Q,
 444   2.487052928527244907490589787691478482358E-1Q,
 445   1.483927406564349124649083853892380899217E0Q,
 446   4.182773513276056975777258788903489507705E0Q,
 447   4.419665392573449746043880892524360870944E0Q,
 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 __float128 Q5_8N[NQ5_8N + 1] = {
 457   -3.656082407740970534915918390488336879763E-13Q,
 458   -1.344660308497244804752334556734121771023E-10Q,
 459   -1.909765035234071738548629788698150760791E-8Q,
 460   -1.366668038160120210269389551283666716453E-6Q,
 461   -5.392327355984269366895210704976314135683E-5Q,
 462   -1.206268245713024564674432357634540343884E-3Q,
 463   -1.515456784370354374066417703736088291287E-2Q,
 464   -1.022454301137286306933217746545237098518E-1Q,
 465   -3.373438906472495080504907858424251082240E-1Q,
 466   -4.510782522110845697262323973549178453405E-1Q,
 467   -1.549000892545288676809660828213589804884E-1Q,
 468 };
 469 #define NQ5_8D 10
 470 static const __float128 Q5_8D[NQ5_8D + 1] = {
 471   3.565550843359501079050699598913828460036E-12Q,
 472   1.321016015556560621591847454285330528045E-9Q,
 473   1.897542728662346479999969679234270605975E-7Q,
 474   1.381720283068706710298734234287456219474E-5Q,
 475   5.599248147286524662305325795203422873725E-4Q,
 476   1.305442352653121436697064782499122164843E-2Q,
 477   1.750234079626943298160445750078631894985E-1Q,
 478   1.311420542073436520965439883806946678491E0Q,
 479   5.162757689856842406744504211089724926650E0Q,
 480   9.527760296384704425618556332087850581308E0Q,
 481   6.604648207463236667912921642545100248584E0Q,
 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 __float128 Q4_5N[NQ4_5N + 1] = {
 491   -4.079513568708891749424783046520200903755E-11Q,
 492   -9.326548104106791766891812583019664893311E-9Q,
 493   -8.016795121318423066292906123815687003356E-7Q,
 494   -3.372350544043594415609295225664186750995E-5Q,
 495   -7.566238665947967882207277686375417983917E-4Q,
 496   -9.248861580055565402130441618521591282617E-3Q,
 497   -6.033106131055851432267702948850231270338E-2Q,
 498   -1.966908754799996793730369265431584303447E-1Q,
 499   -2.791062741179964150755788226623462207560E-1Q,
 500   -1.255478605849190549914610121863534191666E-1Q,
 501   -4.320429862021265463213168186061696944062E-3Q,
 502 };
 503 #define NQ4_5D 9
 504 static const __float128 Q4_5D[NQ4_5D + 1] = {
 505   3.978497042580921479003851216297330701056E-10Q,
 506   9.203304163828145809278568906420772246666E-8Q,
 507   8.059685467088175644915010485174545743798E-6Q,
 508   3.490187375993956409171098277561669167446E-4Q,
 509   8.189109654456872150100501732073810028829E-3Q,
 510   1.072572867311023640958725265762483033769E-1Q,
 511   7.790606862409960053675717185714576937994E-1Q,
 512   3.016049768232011196434185423512777656328E0Q,
 513   5.722963851442769787733717162314477949360E0Q,
 514   4.510527838428473279647251350931380867663E0Q,
 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 __float128 Q3r2_4N[NQ3r2_4N + 1] = {
 524   -1.087480809271383885936921889040388133627E-8Q,
 525   -1.690067828697463740906962973479310170932E-6Q,
 526   -9.608064416995105532790745641974762550982E-5Q,
 527   -2.594198839156517191858208513873961837410E-3Q,
 528   -3.610954144421543968160459863048062977822E-2Q,
 529   -2.629866798251843212210482269563961685666E-1Q,
 530   -9.709186825881775885917984975685752956660E-1Q,
 531   -1.667521829918185121727268867619982417317E0Q,
 532   -1.109255082925540057138766105229900943501E0Q,
 533   -1.812932453006641348145049323713469043328E-1Q,
 534 };
 535 #define NQ3r2_4D 9
 536 static const __float128 Q3r2_4D[NQ3r2_4D + 1] = {
 537   1.060552717496912381388763753841473407026E-7Q,
 538   1.676928002024920520786883649102388708024E-5Q,
 539   9.803481712245420839301400601140812255737E-4Q,
 540   2.765559874262309494758505158089249012930E-2Q,
 541   4.117921827792571791298862613287549140706E-1Q,
 542   3.323769515244751267093378361930279161413E0Q,
 543   1.436602494405814164724810151689705353670E1Q,
 544   3.163087869617098638064881410646782408297E1Q,
 545   3.198181264977021649489103980298349589419E1Q,
 546   1.203649258862068431199471076202897823272E1Q,
 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 __float128 Q2r7_3r2N[NQ2r7_3r2N + 1] = {
 556   -1.723405393982209853244278760171643219530E-7Q,
 557   -2.090508758514655456365709712333460087442E-5Q,
 558   -9.140104013370974823232873472192719263019E-4Q,
 559   -1.871349499990714843332742160292474780128E-2Q,
 560   -1.948930738119938669637865956162512983416E-1Q,
 561   -1.048764684978978127908439526343174139788E0Q,
 562   -2.827714929925679500237476105843643064698E0Q,
 563   -3.508761569156476114276988181329773987314E0Q,
 564   -1.669332202790211090973255098624488308989E0Q,
 565   -1.930796319299022954013840684651016077770E-1Q,
 566 };
 567 #define NQ2r7_3r2D 9
 568 static const __float128 Q2r7_3r2D[NQ2r7_3r2D + 1] = {
 569   1.680730662300831976234547482334347983474E-6Q,
 570   2.084241442440551016475972218719621841120E-4Q,
 571   9.445316642108367479043541702688736295579E-3Q,
 572   2.044637889456631896650179477133252184672E-1Q,
 573   2.316091982244297350829522534435350078205E0Q,
 574   1.412031891783015085196708811890448488865E1Q,
 575   4.583830154673223384837091077279595496149E1Q,
 576   7.549520609270909439885998474045974122261E1Q,
 577   5.697605832808113367197494052388203310638E1Q,
 578   1.601496240876192444526383314589371686234E1Q,
 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 __float128 Q2r3_2r7N[NQ2r3_2r7N + 1] = {
 588   -8.603042076329122085722385914954878953775E-7Q,
 589   -7.701746260451647874214968882605186675720E-5Q,
 590   -2.407932004380727587382493696877569654271E-3Q,
 591   -3.403434217607634279028110636919987224188E-2Q,
 592   -2.348707332185238159192422084985713102877E-1Q,
 593   -7.957498841538254916147095255700637463207E-1Q,
 594   -1.258469078442635106431098063707934348577E0Q,
 595   -8.162415474676345812459353639449971369890E-1Q,
 596   -1.581783890269379690141513949609572806898E-1Q,
 597   -1.890595651683552228232308756569450822905E-3Q,
 598 };
 599 #define NQ2r3_2r7D 8
 600 static const __float128 Q2r3_2r7D[NQ2r3_2r7D + 1] = {
 601   8.390017524798316921170710533381568175665E-6Q,
 602   7.738148683730826286477254659973968763659E-4Q,
 603   2.541480810958665794368759558791634341779E-2Q,
 604   3.878879789711276799058486068562386244873E-1Q,
 605   3.003783779325811292142957336802456109333E0Q,
 606   1.206480374773322029883039064575464497400E1Q,
 607   2.458414064785315978408974662900438351782E1Q,
 608   2.367237826273668567199042088835448715228E1Q,
 609   9.231451197519171090875569102116321676763E0Q,
 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 __float128 Q2_2r3N[NQ2_2r3N + 1] = {
 619   -5.552507516089087822166822364590806076174E-6Q,
 620   -4.135067659799500521040944087433752970297E-4Q,
 621   -1.059928728869218962607068840646564457980E-2Q,
 622   -1.212070036005832342565792241385459023801E-1Q,
 623   -6.688350110633603958684302153362735625156E-1Q,
 624   -1.793587878197360221340277951304429821582E0Q,
 625   -2.225407682237197485644647380483725045326E0Q,
 626   -1.123402135458940189438898496348239744403E0Q,
 627   -1.679187241566347077204805190763597299805E-1Q,
 628   -1.458550613639093752909985189067233504148E-3Q,
 629 };
 630 #define NQ2_2r3D 8
 631 static const __float128 Q2_2r3D[NQ2_2r3D + 1] = {
 632   5.415024336507980465169023996403597916115E-5Q,
 633   4.179246497380453022046357404266022870788E-3Q,
 634   1.136306384261959483095442402929502368598E-1Q,
 635   1.422640343719842213484515445393284072830E0Q,
 636   8.968786703393158374728850922289204805764E0Q,
 637   2.914542473339246127533384118781216495934E1Q,
 638   4.781605421020380669870197378210457054685E1Q,
 639   3.693865837171883152382820584714795072937E1Q,
 640   1.153220502744204904763115556224395893076E1Q,
 641   /* 1.000000000000000000000000000000000000000E0 */
 642 };
 643
 644
 645 /* Evaluate P[n] x^n  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
 646
 647 static __float128
 648 neval (__float128 x, const __float128 *p, int n)
 649 {
 650   __float128 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 __float128
 666 deval (__float128 x, const __float128 *p, int n)
 667 {
 668   __float128 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 __float128
 684 j1q (__float128 x)
 685 {
 686   __float128 xx, xinv, z, p, q, c, s, cc, ss;
 687
 688   if (! finiteq (x))
 689     {
 690       if (x != x)
 691         return x + x;
 692       else
 693         return 0.0Q;
 694     }
 695   if (x == 0.0Q)
 696     return x;
 697   xx = fabsq (x);
 698   if (xx <= 0x1p-58Q)
 699     {
 700       __float128 ret = x * 0.5Q;
 701       math_check_force_underflow (ret);
 702       if (ret == 0)
 703         errno = ERANGE;
 704       return ret;
 705     }
 706   if (xx <= 2.0Q)
 707     {
 708       /* 0 <= x <= 2 */
 709       z = xx * xx;
 710       p = xx * z * neval (z, J0_2N, NJ0_2N) / deval (z, J0_2D, NJ0_2D);
 711       p += 0.5Q * xx;
 712       if (x < 0)
 713         p = -p;
 714       return p;
 715     }
 716
 717   /* X = x - 3 pi/4
 718      cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
 719      = 1/sqrt(2) * (-cos(x) + sin(x))
 720      sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
 721      = -1/sqrt(2) * (sin(x) + cos(x))
 722      cf. Fdlibm.  */
 723   sincosq (xx, &s, &c);
 724   ss = -s - c;
 725   cc = s - c;
 726   if (xx <= FLT128_MAX / 2.0Q)
 727     {
 728       z = cosq (xx + xx);
 729       if ((s * c) > 0)
 730         cc = z / ss;
 731       else
 732         ss = z / cc;
 733     }
 734
 735   if (xx > 0x1p256Q)
 736     {
 737       z = ONEOSQPI * cc / sqrtq (xx);
 738       if (x < 0)
 739         z = -z;
 740       return z;
 741     }
 742
 743   xinv = 1.0Q / xx;
 744   z = xinv * xinv;
 745   if (xinv <= 0.25)
 746     {
 747       if (xinv <= 0.125)
 748         {
 749           if (xinv <= 0.0625)
 750             {
 751               p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID);
 752               q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID);
 753             }
 754           else
 755             {
 756               p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D);
 757               q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D);
 758             }
 759         }
 760       else if (xinv <= 0.1875)
 761         {
 762           p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D);
 763           q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D);
 764         }
 765       else
 766         {
 767           p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D);
 768           q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D);
 769         }
 770     }                           /* .25 */
 771   else /* if (xinv <= 0.5) */
 772     {
 773       if (xinv <= 0.375)
 774         {
 775           if (xinv <= 0.3125)
 776             {
 777               p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D);
 778               q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D);
 779             }
 780           else
 781             {
 782               p = neval (z, P2r7_3r2N, NP2r7_3r2N)
 783                   / deval (z, P2r7_3r2D, NP2r7_3r2D);
 784               q = neval (z, Q2r7_3r2N, NQ2r7_3r2N)
 785                   / deval (z, Q2r7_3r2D, NQ2r7_3r2D);
 786             }
 787         }
 788       else if (xinv <= 0.4375)
 789         {
 790           p = neval (z, P2r3_2r7N, NP2r3_2r7N)
 791               / deval (z, P2r3_2r7D, NP2r3_2r7D);
 792           q = neval (z, Q2r3_2r7N, NQ2r3_2r7N)
 793               / deval (z, Q2r3_2r7D, NQ2r3_2r7D);
 794         }
 795       else
 796         {
 797           p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
 798           q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
 799         }
 800     }
 801   p = 1.0Q + z * p;
 802   q = z * q;
 803   q = q * xinv + 0.375Q * xinv;
 804   z = ONEOSQPI * (p * cc - q * ss) / sqrtq (xx);
 805   if (x < 0)
 806     z = -z;
 807   return z;
 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 __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 __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.0Q)
 849     {
 850       if (x < 0.0Q)
 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.0Q)
 863    {
 864       /* 0 <= x <= 2 */
 865       /* FIXME:  SET_RESTORE_ROUNDL (FE_TONEAREST);  */
 866       z = xx * xx;
 867       p = xx * neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D);
 868       p = -TWOOPI / xx + p;
 869       p = TWOOPI * 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.0Q)
 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.0Q / 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.0Q + 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 }