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) */ 101static const __float128 ONEOSQPI = 5.6418958354775628694807945156077258584405E-1Q; 102/* 2 / pi */ 103static const __float128 TWOOPI = 6.3661977236758134307553505349005744813784E-1Q; 104static 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 110static 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 120static 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 135static 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 148static 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 166static 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 181static 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 200static 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 214static 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 233static 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 247static 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 265static 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 278static 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 296static 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 309static 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 326static 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 339static 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 356static 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 368static 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 386static 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 400static 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 419static 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 434static 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 455static 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 469static 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 489static 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 503static 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 522static 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 535static 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 554static 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 567static 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 586static 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 599static 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 617static 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 630static 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 646static __float128 647neval (__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 664static __float128 665deval (__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 683j1q (__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 815static 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 826static 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 842y1q (__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} 958