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 "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.0Q;
 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.000000000000000000000000000000000000000E0Q */
 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;
 691       else
 692         return 0.0Q;
 693     }
 694   if (x == 0.0Q)
 695     return x;
 696   xx = fabsq (x);
 697   if (xx <= 2.0Q)
 698     {
 699       /* 0 <= x <= 2 */
 700       z = xx * xx;
 701       p = xx * z * neval (z, J0_2N, NJ0_2N) / deval (z, J0_2D, NJ0_2D);
 702       p += 0.5Q * xx;
 703       if (x < 0)
 704         p = -p;
 705       return p;
 706     }
 707
 708   xinv = 1.0Q / xx;
 709   z = xinv * xinv;
 710   if (xinv <= 0.25)
 711     {
 712       if (xinv <= 0.125)
 713         {
 714           if (xinv <= 0.0625)
 715             {
 716               p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID);
 717               q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID);
 718             }
 719           else
 720             {
 721               p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D);
 722               q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D);
 723             }
 724         }
 725       else if (xinv <= 0.1875)
 726         {
 727           p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D);
 728           q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D);
 729         }
 730       else
 731         {
 732           p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D);
 733           q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D);
 734         }
 735     }                           /* .25 */
 736   else /* if (xinv <= 0.5) */
 737     {
 738       if (xinv <= 0.375)
 739         {
 740           if (xinv <= 0.3125)
 741             {
 742               p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D);
 743               q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D);
 744             }
 745           else
 746             {
 747               p = neval (z, P2r7_3r2N, NP2r7_3r2N)
 748                   / deval (z, P2r7_3r2D, NP2r7_3r2D);
 749               q = neval (z, Q2r7_3r2N, NQ2r7_3r2N)
 750                   / deval (z, Q2r7_3r2D, NQ2r7_3r2D);
 751             }
 752         }
 753       else if (xinv <= 0.4375)
 754         {
 755           p = neval (z, P2r3_2r7N, NP2r3_2r7N)
 756               / deval (z, P2r3_2r7D, NP2r3_2r7D);
 757           q = neval (z, Q2r3_2r7N, NQ2r3_2r7N)
 758               / deval (z, Q2r3_2r7D, NQ2r3_2r7D);
 759         }
 760       else
 761         {
 762           p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
 763           q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
 764         }
 765     }
 766   p = 1.0Q + z * p;
 767   q = z * q;
 768   q = q * xinv + 0.375Q * xinv;
 769   /* X = x - 3 pi/4
 770      cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
 771      = 1/sqrt(2) * (-cos(x) + sin(x))
 772      sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
 773      = -1/sqrt(2) * (sin(x) + cos(x))
 774      cf. Fdlibm.  */
 775   sincosq (xx, &s, &c);
 776   ss = -s - c;
 777   cc = s - c;
 778   z = cosq (xx + xx);
 779   if ((s * c) > 0)
 780     cc = z / ss;
 781   else
 782     ss = z / cc;
 783   z = ONEOSQPI * (p * cc - q * ss) / sqrtq (xx);
 784   if (x < 0)
 785     z = -z;
 786   return z;
 787 }
 788
 789
 790 /* Y1(x) = 2/pi * (log(x) * J1(x) - 1/x) + x R(x^2)
 791    Peak relative error 6.2e-38
 792    0 <= x <= 2   */
 793 #define NY0_2N 7
 794 static __float128 Y0_2N[NY0_2N + 1] = {
 795   -6.804415404830253804408698161694720833249E19Q,
 796   1.805450517967019908027153056150465849237E19Q,
 797   -8.065747497063694098810419456383006737312E17Q,
 798   1.401336667383028259295830955439028236299E16Q,
 799   -1.171654432898137585000399489686629680230E14Q,
 800   5.061267920943853732895341125243428129150E11Q,
 801   -1.096677850566094204586208610960870217970E9Q,
 802   9.541172044989995856117187515882879304461E5Q,
 803 };
 804 #define NY0_2D 7
 805 static __float128 Y0_2D[NY0_2D + 1] = {
 806   3.470629591820267059538637461549677594549E20Q,
 807   4.120796439009916326855848107545425217219E18Q,
 808   2.477653371652018249749350657387030814542E16Q,
 809   9.954678543353888958177169349272167762797E13Q,
 810   2.957927997613630118216218290262851197754E11Q,
 811   6.748421382188864486018861197614025972118E8Q,
 812   1.173453425218010888004562071020305709319E6Q,
 813   1.450335662961034949894009554536003377187E3Q,
 814   /* 1.000000000000000000000000000000000000000E0 */
 815 };
 816
 817
 818 /* Bessel function of the second kind, order one.  */
 819
 820 __float128
 821 y1q (__float128 x)
 822 {
 823   __float128 xx, xinv, z, p, q, c, s, cc, ss;
 824
 825   if (! finiteq (x))
 826     {
 827       if (x != x)
 828         return x;
 829       else
 830         return 0.0Q;
 831     }
 832   if (x <= 0.0Q)
 833     {
 834       if (x < 0.0Q)
 835         return (zero / (zero * x));
 836       return -HUGE_VALQ + x;
 837     }
 838   xx = fabsq (x);
 839   if (xx <= 0x1p-114)
 840     return -TWOOPI / x;
 841   if (xx <= 2.0Q)
 842    {
 843       /* 0 <= x <= 2 */
 844       z = xx * xx;
 845       p = xx * neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D);
 846       p = -TWOOPI / xx + p;
 847       p = TWOOPI * logq (x) * j1q (x) + p;
 848       return p;
 849     }
 850
 851   xinv = 1.0Q / xx;
 852   z = xinv * xinv;
 853   if (xinv <= 0.25)
 854     {
 855       if (xinv <= 0.125)
 856         {
 857           if (xinv <= 0.0625)
 858             {
 859               p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID);
 860               q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID);
 861             }
 862           else
 863             {
 864               p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D);
 865               q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D);
 866             }
 867         }
 868       else if (xinv <= 0.1875)
 869         {
 870           p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D);
 871           q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D);
 872         }
 873       else
 874         {
 875           p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D);
 876           q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D);
 877         }
 878     }                           /* .25 */
 879   else /* if (xinv <= 0.5) */
 880     {
 881       if (xinv <= 0.375)
 882         {
 883           if (xinv <= 0.3125)
 884             {
 885               p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D);
 886               q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D);
 887             }
 888           else
 889             {
 890               p = neval (z, P2r7_3r2N, NP2r7_3r2N)
 891                   / deval (z, P2r7_3r2D, NP2r7_3r2D);
 892               q = neval (z, Q2r7_3r2N, NQ2r7_3r2N)
 893                   / deval (z, Q2r7_3r2D, NQ2r7_3r2D);
 894             }
 895         }
 896       else if (xinv <= 0.4375)
 897         {
 898           p = neval (z, P2r3_2r7N, NP2r3_2r7N)
 899               / deval (z, P2r3_2r7D, NP2r3_2r7D);
 900           q = neval (z, Q2r3_2r7N, NQ2r3_2r7N)
 901               / deval (z, Q2r3_2r7D, NQ2r3_2r7D);
 902         }
 903       else
 904         {
 905           p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
 906           q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
 907         }
 908     }
 909   p = 1.0Q + z * p;
 910   q = z * q;
 911   q = q * xinv + 0.375Q * xinv;
 912   /* X = x - 3 pi/4
 913      cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
 914      = 1/sqrt(2) * (-cos(x) + sin(x))
 915      sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
 916      = -1/sqrt(2) * (sin(x) + cos(x))
 917      cf. Fdlibm.  */
 918   sincosq (xx, &s, &c);
 919   ss = -s - c;
 920   cc = s - c;
 921   z = cosq (xx + xx);
 922   if ((s * c) > 0)
 923     cc = z / ss;
 924   else
 925     ss = z / cc;
 926   z = ONEOSQPI * (p * ss + q * cc) / sqrtq (xx);
 927   return z;
 928 }