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) */ 101static const __float128 ONEOSQPI = 5.6418958354775628694807945156077258584405E-1Q; 102/* 2 / pi */ 103static const __float128 TWOOPI = 6.3661977236758134307553505349005744813784E-1Q; 104static 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 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.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 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; 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 794static __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 805static __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 821y1q (__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} 929