1 /* j0l.c 2 * 3 * Bessel function of order zero 4 * 5 * 6 * 7 * SYNOPSIS: 8 * 9 * long double x, y, j0l(); 10 * 11 * y = j0l( x ); 12 * 13 * 14 * 15 * DESCRIPTION: 16 * 17 * Returns Bessel function of first kind, order zero 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 21 * is J0(x) = 1 - x^2 / 4 + x^4 R(x^2) 22 * The second interval is further partitioned into eight equal segments 23 * of 1/x. 24 * 25 * J0(x) = sqrt(2/(pi x)) (P0(x) cos(X) - Q0(x) sin(X)), 26 * X = x - pi/4, 27 * 28 * and the auxiliary functions are given by 29 * 30 * J0(x)cos(X) + Y0(x)sin(X) = sqrt( 2/(pi x)) P0(x), 31 * P0(x) = 1 + 1/x^2 R(1/x^2) 32 * 33 * Y0(x)cos(X) - J0(x)sin(X) = sqrt( 2/(pi x)) Q0(x), 34 * Q0(x) = 1/x (-.125 + 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 1.7e-34 2.4e-35 43 * 44 * 45 */ 46 47 /* y0l.c 48 * 49 * Bessel function of the second kind, order zero 50 * 51 * 52 * 53 * SYNOPSIS: 54 * 55 * double x, y, y0l(); 56 * 57 * y = y0l( x ); 58 * 59 * 60 * 61 * DESCRIPTION: 62 * 63 * Returns Bessel function of the second kind, of order 64 * zero, of the argument. 65 * 66 * The approximation is the same as for J0(x), and 67 * Y0(x) = sqrt(2/(pi x)) (P0(x) sin(X) + Q0(x) cos(X)). 68 * 69 * ACCURACY: 70 * 71 * Absolute error, when y0(x) < 1; else relative error: 72 * 73 * arithmetic domain # trials peak rms 74 * IEEE 0, 30 100000 3.0e-34 2.7e-35 75 * 76 */ 77 78 /* Copyright 2001 by Stephen L. Moshier (moshier@na-net.ornl.gov). 79 80 This library is free software; you can redistribute it and/or 81 modify it under the terms of the GNU Lesser General Public 82 License as published by the Free Software Foundation; either 83 version 2.1 of the License, or (at your option) any later version. 84 85 This library is distributed in the hope that it will be useful, 86 but WITHOUT ANY WARRANTY; without even the implied warranty of 87 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 88 Lesser General Public License for more details. 89 90 You should have received a copy of the GNU Lesser General Public 91 License along with this library; if not, see 92 <http://www.gnu.org/licenses/>. */ 93 94 #include "quadmath-imp.h" 95 96 /* 1 / sqrt(pi) */ 97 static const __float128 ONEOSQPI = 5.6418958354775628694807945156077258584405E-1Q; 98 /* 2 / pi */ 99 static const __float128 TWOOPI = 6.3661977236758134307553505349005744813784E-1Q; 100 static const __float128 zero = 0; 101 102 /* J0(x) = 1 - x^2/4 + x^2 x^2 R(x^2) 103 Peak relative error 3.4e-37 104 0 <= x <= 2 */ 105 #define NJ0_2N 6 106 static const __float128 J0_2N[NJ0_2N + 1] = { 107 3.133239376997663645548490085151484674892E16Q, 108 -5.479944965767990821079467311839107722107E14Q, 109 6.290828903904724265980249871997551894090E12Q, 110 -3.633750176832769659849028554429106299915E10Q, 111 1.207743757532429576399485415069244807022E8Q, 112 -2.107485999925074577174305650549367415465E5Q, 113 1.562826808020631846245296572935547005859E2Q, 114 }; 115 #define NJ0_2D 6 116 static const __float128 J0_2D[NJ0_2D + 1] = { 117 2.005273201278504733151033654496928968261E18Q, 118 2.063038558793221244373123294054149790864E16Q, 119 1.053350447931127971406896594022010524994E14Q, 120 3.496556557558702583143527876385508882310E11Q, 121 8.249114511878616075860654484367133976306E8Q, 122 1.402965782449571800199759247964242790589E6Q, 123 1.619910762853439600957801751815074787351E3Q, 124 /* 1.000000000000000000000000000000000000000E0 */ 125 }; 126 127 /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2), 128 0 <= 1/x <= .0625 129 Peak relative error 3.3e-36 */ 130 #define NP16_IN 9 131 static const __float128 P16_IN[NP16_IN + 1] = { 132 -1.901689868258117463979611259731176301065E-16Q, 133 -1.798743043824071514483008340803573980931E-13Q, 134 -6.481746687115262291873324132944647438959E-11Q, 135 -1.150651553745409037257197798528294248012E-8Q, 136 -1.088408467297401082271185599507222695995E-6Q, 137 -5.551996725183495852661022587879817546508E-5Q, 138 -1.477286941214245433866838787454880214736E-3Q, 139 -1.882877976157714592017345347609200402472E-2Q, 140 -9.620983176855405325086530374317855880515E-2Q, 141 -1.271468546258855781530458854476627766233E-1Q, 142 }; 143 #define NP16_ID 9 144 static const __float128 P16_ID[NP16_ID + 1] = { 145 2.704625590411544837659891569420764475007E-15Q, 146 2.562526347676857624104306349421985403573E-12Q, 147 9.259137589952741054108665570122085036246E-10Q, 148 1.651044705794378365237454962653430805272E-7Q, 149 1.573561544138733044977714063100859136660E-5Q, 150 8.134482112334882274688298469629884804056E-4Q, 151 2.219259239404080863919375103673593571689E-2Q, 152 2.976990606226596289580242451096393862792E-1Q, 153 1.713895630454693931742734911930937246254E0Q, 154 3.231552290717904041465898249160757368855E0Q, 155 /* 1.000000000000000000000000000000000000000E0 */ 156 }; 157 158 /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2) 159 0.0625 <= 1/x <= 0.125 160 Peak relative error 2.4e-35 */ 161 #define NP8_16N 10 162 static const __float128 P8_16N[NP8_16N + 1] = { 163 -2.335166846111159458466553806683579003632E-15Q, 164 -1.382763674252402720401020004169367089975E-12Q, 165 -3.192160804534716696058987967592784857907E-10Q, 166 -3.744199606283752333686144670572632116899E-8Q, 167 -2.439161236879511162078619292571922772224E-6Q, 168 -9.068436986859420951664151060267045346549E-5Q, 169 -1.905407090637058116299757292660002697359E-3Q, 170 -2.164456143936718388053842376884252978872E-2Q, 171 -1.212178415116411222341491717748696499966E-1Q, 172 -2.782433626588541494473277445959593334494E-1Q, 173 -1.670703190068873186016102289227646035035E-1Q, 174 }; 175 #define NP8_16D 10 176 static const __float128 P8_16D[NP8_16D + 1] = { 177 3.321126181135871232648331450082662856743E-14Q, 178 1.971894594837650840586859228510007703641E-11Q, 179 4.571144364787008285981633719513897281690E-9Q, 180 5.396419143536287457142904742849052402103E-7Q, 181 3.551548222385845912370226756036899901549E-5Q, 182 1.342353874566932014705609788054598013516E-3Q, 183 2.899133293006771317589357444614157734385E-2Q, 184 3.455374978185770197704507681491574261545E-1Q, 185 2.116616964297512311314454834712634820514E0Q, 186 5.850768316827915470087758636881584174432E0Q, 187 5.655273858938766830855753983631132928968E0Q, 188 /* 1.000000000000000000000000000000000000000E0 */ 189 }; 190 191 /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2) 192 0.125 <= 1/x <= 0.1875 193 Peak relative error 2.7e-35 */ 194 #define NP5_8N 10 195 static const __float128 P5_8N[NP5_8N + 1] = { 196 -1.270478335089770355749591358934012019596E-12Q, 197 -4.007588712145412921057254992155810347245E-10Q, 198 -4.815187822989597568124520080486652009281E-8Q, 199 -2.867070063972764880024598300408284868021E-6Q, 200 -9.218742195161302204046454768106063638006E-5Q, 201 -1.635746821447052827526320629828043529997E-3Q, 202 -1.570376886640308408247709616497261011707E-2Q, 203 -7.656484795303305596941813361786219477807E-2Q, 204 -1.659371030767513274944805479908858628053E-1Q, 205 -1.185340550030955660015841796219919804915E-1Q, 206 -8.920026499909994671248893388013790366712E-3Q, 207 }; 208 #define NP5_8D 9 209 static const __float128 P5_8D[NP5_8D + 1] = { 210 1.806902521016705225778045904631543990314E-11Q, 211 5.728502760243502431663549179135868966031E-9Q, 212 6.938168504826004255287618819550667978450E-7Q, 213 4.183769964807453250763325026573037785902E-5Q, 214 1.372660678476925468014882230851637878587E-3Q, 215 2.516452105242920335873286419212708961771E-2Q, 216 2.550502712902647803796267951846557316182E-1Q, 217 1.365861559418983216913629123778747617072E0Q, 218 3.523825618308783966723472468855042541407E0Q, 219 3.656365803506136165615111349150536282434E0Q, 220 /* 1.000000000000000000000000000000000000000E0 */ 221 }; 222 223 /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2) 224 Peak relative error 3.5e-35 225 0.1875 <= 1/x <= 0.25 */ 226 #define NP4_5N 9 227 static const __float128 P4_5N[NP4_5N + 1] = { 228 -9.791405771694098960254468859195175708252E-10Q, 229 -1.917193059944531970421626610188102836352E-7Q, 230 -1.393597539508855262243816152893982002084E-5Q, 231 -4.881863490846771259880606911667479860077E-4Q, 232 -8.946571245022470127331892085881699269853E-3Q, 233 -8.707474232568097513415336886103899434251E-2Q, 234 -4.362042697474650737898551272505525973766E-1Q, 235 -1.032712171267523975431451359962375617386E0Q, 236 -9.630502683169895107062182070514713702346E-1Q, 237 -2.251804386252969656586810309252357233320E-1Q, 238 }; 239 #define NP4_5D 9 240 static const __float128 P4_5D[NP4_5D + 1] = { 241 1.392555487577717669739688337895791213139E-8Q, 242 2.748886559120659027172816051276451376854E-6Q, 243 2.024717710644378047477189849678576659290E-4Q, 244 7.244868609350416002930624752604670292469E-3Q, 245 1.373631762292244371102989739300382152416E-1Q, 246 1.412298581400224267910294815260613240668E0Q, 247 7.742495637843445079276397723849017617210E0Q, 248 2.138429269198406512028307045259503811861E1Q, 249 2.651547684548423476506826951831712762610E1Q, 250 1.167499382465291931571685222882909166935E1Q, 251 /* 1.000000000000000000000000000000000000000E0 */ 252 }; 253 254 /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2) 255 Peak relative error 2.3e-36 256 0.25 <= 1/x <= 0.3125 */ 257 #define NP3r2_4N 9 258 static const __float128 P3r2_4N[NP3r2_4N + 1] = { 259 -2.589155123706348361249809342508270121788E-8Q, 260 -3.746254369796115441118148490849195516593E-6Q, 261 -1.985595497390808544622893738135529701062E-4Q, 262 -5.008253705202932091290132760394976551426E-3Q, 263 -6.529469780539591572179155511840853077232E-2Q, 264 -4.468736064761814602927408833818990271514E-1Q, 265 -1.556391252586395038089729428444444823380E0Q, 266 -2.533135309840530224072920725976994981638E0Q, 267 -1.605509621731068453869408718565392869560E0Q, 268 -2.518966692256192789269859830255724429375E-1Q, 269 }; 270 #define NP3r2_4D 9 271 static const __float128 P3r2_4D[NP3r2_4D + 1] = { 272 3.682353957237979993646169732962573930237E-7Q, 273 5.386741661883067824698973455566332102029E-5Q, 274 2.906881154171822780345134853794241037053E-3Q, 275 7.545832595801289519475806339863492074126E-2Q, 276 1.029405357245594877344360389469584526654E0Q, 277 7.565706120589873131187989560509757626725E0Q, 278 2.951172890699569545357692207898667665796E1Q, 279 5.785723537170311456298467310529815457536E1Q, 280 5.095621464598267889126015412522773474467E1Q, 281 1.602958484169953109437547474953308401442E1Q, 282 /* 1.000000000000000000000000000000000000000E0 */ 283 }; 284 285 /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2) 286 Peak relative error 1.0e-35 287 0.3125 <= 1/x <= 0.375 */ 288 #define NP2r7_3r2N 9 289 static const __float128 P2r7_3r2N[NP2r7_3r2N + 1] = { 290 -1.917322340814391131073820537027234322550E-7Q, 291 -1.966595744473227183846019639723259011906E-5Q, 292 -7.177081163619679403212623526632690465290E-4Q, 293 -1.206467373860974695661544653741899755695E-2Q, 294 -1.008656452188539812154551482286328107316E-1Q, 295 -4.216016116408810856620947307438823892707E-1Q, 296 -8.378631013025721741744285026537009814161E-1Q, 297 -6.973895635309960850033762745957946272579E-1Q, 298 -1.797864718878320770670740413285763554812E-1Q, 299 -4.098025357743657347681137871388402849581E-3Q, 300 }; 301 #define NP2r7_3r2D 8 302 static const __float128 P2r7_3r2D[NP2r7_3r2D + 1] = { 303 2.726858489303036441686496086962545034018E-6Q, 304 2.840430827557109238386808968234848081424E-4Q, 305 1.063826772041781947891481054529454088832E-2Q, 306 1.864775537138364773178044431045514405468E-1Q, 307 1.665660052857205170440952607701728254211E0Q, 308 7.723745889544331153080842168958348568395E0Q, 309 1.810726427571829798856428548102077799835E1Q, 310 1.986460672157794440666187503833545388527E1Q, 311 8.645503204552282306364296517220055815488E0Q, 312 /* 1.000000000000000000000000000000000000000E0 */ 313 }; 314 315 /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2) 316 Peak relative error 1.3e-36 317 0.3125 <= 1/x <= 0.4375 */ 318 #define NP2r3_2r7N 9 319 static const __float128 P2r3_2r7N[NP2r3_2r7N + 1] = { 320 -1.594642785584856746358609622003310312622E-6Q, 321 -1.323238196302221554194031733595194539794E-4Q, 322 -3.856087818696874802689922536987100372345E-3Q, 323 -5.113241710697777193011470733601522047399E-2Q, 324 -3.334229537209911914449990372942022350558E-1Q, 325 -1.075703518198127096179198549659283422832E0Q, 326 -1.634174803414062725476343124267110981807E0Q, 327 -1.030133247434119595616826842367268304880E0Q, 328 -1.989811539080358501229347481000707289391E-1Q, 329 -3.246859189246653459359775001466924610236E-3Q, 330 }; 331 #define NP2r3_2r7D 8 332 static const __float128 P2r3_2r7D[NP2r3_2r7D + 1] = { 333 2.267936634217251403663034189684284173018E-5Q, 334 1.918112982168673386858072491437971732237E-3Q, 335 5.771704085468423159125856786653868219522E-2Q, 336 8.056124451167969333717642810661498890507E-1Q, 337 5.687897967531010276788680634413789328776E0Q, 338 2.072596760717695491085444438270778394421E1Q, 339 3.801722099819929988585197088613160496684E1Q, 340 3.254620235902912339534998592085115836829E1Q, 341 1.104847772130720331801884344645060675036E1Q, 342 /* 1.000000000000000000000000000000000000000E0 */ 343 }; 344 345 /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2) 346 Peak relative error 1.2e-35 347 0.4375 <= 1/x <= 0.5 */ 348 #define NP2_2r3N 8 349 static const __float128 P2_2r3N[NP2_2r3N + 1] = { 350 -1.001042324337684297465071506097365389123E-4Q, 351 -6.289034524673365824853547252689991418981E-3Q, 352 -1.346527918018624234373664526930736205806E-1Q, 353 -1.268808313614288355444506172560463315102E0Q, 354 -5.654126123607146048354132115649177406163E0Q, 355 -1.186649511267312652171775803270911971693E1Q, 356 -1.094032424931998612551588246779200724257E1Q, 357 -3.728792136814520055025256353193674625267E0Q, 358 -3.000348318524471807839934764596331810608E-1Q, 359 }; 360 #define NP2_2r3D 8 361 static const __float128 P2_2r3D[NP2_2r3D + 1] = { 362 1.423705538269770974803901422532055612980E-3Q, 363 9.171476630091439978533535167485230575894E-2Q, 364 2.049776318166637248868444600215942828537E0Q, 365 2.068970329743769804547326701946144899583E1Q, 366 1.025103500560831035592731539565060347709E2Q, 367 2.528088049697570728252145557167066708284E2Q, 368 2.992160327587558573740271294804830114205E2Q, 369 1.540193761146551025832707739468679973036E2Q, 370 2.779516701986912132637672140709452502650E1Q, 371 /* 1.000000000000000000000000000000000000000E0 */ 372 }; 373 374 /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x), 375 Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2)) 376 Peak relative error 2.2e-35 377 0 <= 1/x <= .0625 */ 378 #define NQ16_IN 10 379 static const __float128 Q16_IN[NQ16_IN + 1] = { 380 2.343640834407975740545326632205999437469E-18Q, 381 2.667978112927811452221176781536278257448E-15Q, 382 1.178415018484555397390098879501969116536E-12Q, 383 2.622049767502719728905924701288614016597E-10Q, 384 3.196908059607618864801313380896308968673E-8Q, 385 2.179466154171673958770030655199434798494E-6Q, 386 8.139959091628545225221976413795645177291E-5Q, 387 1.563900725721039825236927137885747138654E-3Q, 388 1.355172364265825167113562519307194840307E-2Q, 389 3.928058355906967977269780046844768588532E-2Q, 390 1.107891967702173292405380993183694932208E-2Q, 391 }; 392 #define NQ16_ID 9 393 static const __float128 Q16_ID[NQ16_ID + 1] = { 394 3.199850952578356211091219295199301766718E-17Q, 395 3.652601488020654842194486058637953363918E-14Q, 396 1.620179741394865258354608590461839031281E-11Q, 397 3.629359209474609630056463248923684371426E-9Q, 398 4.473680923894354600193264347733477363305E-7Q, 399 3.106368086644715743265603656011050476736E-5Q, 400 1.198239259946770604954664925153424252622E-3Q, 401 2.446041004004283102372887804475767568272E-2Q, 402 2.403235525011860603014707768815113698768E-1Q, 403 9.491006790682158612266270665136910927149E-1Q, 404 /* 1.000000000000000000000000000000000000000E0 */ 405 }; 406 407 /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x), 408 Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2)) 409 Peak relative error 5.1e-36 410 0.0625 <= 1/x <= 0.125 */ 411 #define NQ8_16N 11 412 static const __float128 Q8_16N[NQ8_16N + 1] = { 413 1.001954266485599464105669390693597125904E-17Q, 414 7.545499865295034556206475956620160007849E-15Q, 415 2.267838684785673931024792538193202559922E-12Q, 416 3.561909705814420373609574999542459912419E-10Q, 417 3.216201422768092505214730633842924944671E-8Q, 418 1.731194793857907454569364622452058554314E-6Q, 419 5.576944613034537050396518509871004586039E-5Q, 420 1.051787760316848982655967052985391418146E-3Q, 421 1.102852974036687441600678598019883746959E-2Q, 422 5.834647019292460494254225988766702933571E-2Q, 423 1.290281921604364618912425380717127576529E-1Q, 424 7.598886310387075708640370806458926458301E-2Q, 425 }; 426 #define NQ8_16D 11 427 static const __float128 Q8_16D[NQ8_16D + 1] = { 428 1.368001558508338469503329967729951830843E-16Q, 429 1.034454121857542147020549303317348297289E-13Q, 430 3.128109209247090744354764050629381674436E-11Q, 431 4.957795214328501986562102573522064468671E-9Q, 432 4.537872468606711261992676606899273588899E-7Q, 433 2.493639207101727713192687060517509774182E-5Q, 434 8.294957278145328349785532236663051405805E-4Q, 435 1.646471258966713577374948205279380115839E-2Q, 436 1.878910092770966718491814497982191447073E-1Q, 437 1.152641605706170353727903052525652504075E0Q, 438 3.383550240669773485412333679367792932235E0Q, 439 3.823875252882035706910024716609908473970E0Q, 440 /* 1.000000000000000000000000000000000000000E0 */ 441 }; 442 443 /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x), 444 Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2)) 445 Peak relative error 3.9e-35 446 0.125 <= 1/x <= 0.1875 */ 447 #define NQ5_8N 10 448 static const __float128 Q5_8N[NQ5_8N + 1] = { 449 1.750399094021293722243426623211733898747E-13Q, 450 6.483426211748008735242909236490115050294E-11Q, 451 9.279430665656575457141747875716899958373E-9Q, 452 6.696634968526907231258534757736576340266E-7Q, 453 2.666560823798895649685231292142838188061E-5Q, 454 6.025087697259436271271562769707550594540E-4Q, 455 7.652807734168613251901945778921336353485E-3Q, 456 5.226269002589406461622551452343519078905E-2Q, 457 1.748390159751117658969324896330142895079E-1Q, 458 2.378188719097006494782174902213083589660E-1Q, 459 8.383984859679804095463699702165659216831E-2Q, 460 }; 461 #define NQ5_8D 10 462 static const __float128 Q5_8D[NQ5_8D + 1] = { 463 2.389878229704327939008104855942987615715E-12Q, 464 8.926142817142546018703814194987786425099E-10Q, 465 1.294065862406745901206588525833274399038E-7Q, 466 9.524139899457666250828752185212769682191E-6Q, 467 3.908332488377770886091936221573123353489E-4Q, 468 9.250427033957236609624199884089916836748E-3Q, 469 1.263420066165922645975830877751588421451E-1Q, 470 9.692527053860420229711317379861733180654E-1Q, 471 3.937813834630430172221329298841520707954E0Q, 472 7.603126427436356534498908111445191312181E0Q, 473 5.670677653334105479259958485084550934305E0Q, 474 /* 1.000000000000000000000000000000000000000E0 */ 475 }; 476 477 /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x), 478 Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2)) 479 Peak relative error 3.2e-35 480 0.1875 <= 1/x <= 0.25 */ 481 #define NQ4_5N 10 482 static const __float128 Q4_5N[NQ4_5N + 1] = { 483 2.233870042925895644234072357400122854086E-11Q, 484 5.146223225761993222808463878999151699792E-9Q, 485 4.459114531468296461688753521109797474523E-7Q, 486 1.891397692931537975547242165291668056276E-5Q, 487 4.279519145911541776938964806470674565504E-4Q, 488 5.275239415656560634702073291768904783989E-3Q, 489 3.468698403240744801278238473898432608887E-2Q, 490 1.138773146337708415188856882915457888274E-1Q, 491 1.622717518946443013587108598334636458955E-1Q, 492 7.249040006390586123760992346453034628227E-2Q, 493 1.941595365256460232175236758506411486667E-3Q, 494 }; 495 #define NQ4_5D 9 496 static const __float128 Q4_5D[NQ4_5D + 1] = { 497 3.049977232266999249626430127217988047453E-10Q, 498 7.120883230531035857746096928889676144099E-8Q, 499 6.301786064753734446784637919554359588859E-6Q, 500 2.762010530095069598480766869426308077192E-4Q, 501 6.572163250572867859316828886203406361251E-3Q, 502 8.752566114841221958200215255461843397776E-2Q, 503 6.487654992874805093499285311075289932664E-1Q, 504 2.576550017826654579451615283022812801435E0Q, 505 5.056392229924022835364779562707348096036E0Q, 506 4.179770081068251464907531367859072157773E0Q, 507 /* 1.000000000000000000000000000000000000000E0 */ 508 }; 509 510 /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x), 511 Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2)) 512 Peak relative error 1.4e-36 513 0.25 <= 1/x <= 0.3125 */ 514 #define NQ3r2_4N 10 515 static const __float128 Q3r2_4N[NQ3r2_4N + 1] = { 516 6.126167301024815034423262653066023684411E-10Q, 517 1.043969327113173261820028225053598975128E-7Q, 518 6.592927270288697027757438170153763220190E-6Q, 519 2.009103660938497963095652951912071336730E-4Q, 520 3.220543385492643525985862356352195896964E-3Q, 521 2.774405975730545157543417650436941650990E-2Q, 522 1.258114008023826384487378016636555041129E-1Q, 523 2.811724258266902502344701449984698323860E-1Q, 524 2.691837665193548059322831687432415014067E-1Q, 525 7.949087384900985370683770525312735605034E-2Q, 526 1.229509543620976530030153018986910810747E-3Q, 527 }; 528 #define NQ3r2_4D 9 529 static const __float128 Q3r2_4D[NQ3r2_4D + 1] = { 530 8.364260446128475461539941389210166156568E-9Q, 531 1.451301850638956578622154585560759862764E-6Q, 532 9.431830010924603664244578867057141839463E-5Q, 533 3.004105101667433434196388593004526182741E-3Q, 534 5.148157397848271739710011717102773780221E-2Q, 535 4.901089301726939576055285374953887874895E-1Q, 536 2.581760991981709901216967665934142240346E0Q, 537 7.257105880775059281391729708630912791847E0Q, 538 1.006014717326362868007913423810737369312E1Q, 539 5.879416600465399514404064187445293212470E0Q, 540 /* 1.000000000000000000000000000000000000000E0*/ 541 }; 542 543 /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x), 544 Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2)) 545 Peak relative error 3.8e-36 546 0.3125 <= 1/x <= 0.375 */ 547 #define NQ2r7_3r2N 9 548 static const __float128 Q2r7_3r2N[NQ2r7_3r2N + 1] = { 549 7.584861620402450302063691901886141875454E-8Q, 550 9.300939338814216296064659459966041794591E-6Q, 551 4.112108906197521696032158235392604947895E-4Q, 552 8.515168851578898791897038357239630654431E-3Q, 553 8.971286321017307400142720556749573229058E-2Q, 554 4.885856732902956303343015636331874194498E-1Q, 555 1.334506268733103291656253500506406045846E0Q, 556 1.681207956863028164179042145803851824654E0Q, 557 8.165042692571721959157677701625853772271E-1Q, 558 9.805848115375053300608712721986235900715E-2Q, 559 }; 560 #define NQ2r7_3r2D 9 561 static const __float128 Q2r7_3r2D[NQ2r7_3r2D + 1] = { 562 1.035586492113036586458163971239438078160E-6Q, 563 1.301999337731768381683593636500979713689E-4Q, 564 5.993695702564527062553071126719088859654E-3Q, 565 1.321184892887881883489141186815457808785E-1Q, 566 1.528766555485015021144963194165165083312E0Q, 567 9.561463309176490874525827051566494939295E0Q, 568 3.203719484883967351729513662089163356911E1Q, 569 5.497294687660930446641539152123568668447E1Q, 570 4.391158169390578768508675452986948391118E1Q, 571 1.347836630730048077907818943625789418378E1Q, 572 /* 1.000000000000000000000000000000000000000E0 */ 573 }; 574 575 /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x), 576 Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2)) 577 Peak relative error 2.2e-35 578 0.375 <= 1/x <= 0.4375 */ 579 #define NQ2r3_2r7N 9 580 static const __float128 Q2r3_2r7N[NQ2r3_2r7N + 1] = { 581 4.455027774980750211349941766420190722088E-7Q, 582 4.031998274578520170631601850866780366466E-5Q, 583 1.273987274325947007856695677491340636339E-3Q, 584 1.818754543377448509897226554179659122873E-2Q, 585 1.266748858326568264126353051352269875352E-1Q, 586 4.327578594728723821137731555139472880414E-1Q, 587 6.892532471436503074928194969154192615359E-1Q, 588 4.490775818438716873422163588640262036506E-1Q, 589 8.649615949297322440032000346117031581572E-2Q, 590 7.261345286655345047417257611469066147561E-4Q, 591 }; 592 #define NQ2r3_2r7D 8 593 static const __float128 Q2r3_2r7D[NQ2r3_2r7D + 1] = { 594 6.082600739680555266312417978064954793142E-6Q, 595 5.693622538165494742945717226571441747567E-4Q, 596 1.901625907009092204458328768129666975975E-2Q, 597 2.958689532697857335456896889409923371570E-1Q, 598 2.343124711045660081603809437993368799568E0Q, 599 9.665894032187458293568704885528192804376E0Q, 600 2.035273104990617136065743426322454881353E1Q, 601 2.044102010478792896815088858740075165531E1Q, 602 8.445937177863155827844146643468706599304E0Q, 603 /* 1.000000000000000000000000000000000000000E0 */ 604 }; 605 606 /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x), 607 Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2)) 608 Peak relative error 3.1e-36 609 0.4375 <= 1/x <= 0.5 */ 610 #define NQ2_2r3N 9 611 static const __float128 Q2_2r3N[NQ2_2r3N + 1] = { 612 2.817566786579768804844367382809101929314E-6Q, 613 2.122772176396691634147024348373539744935E-4Q, 614 5.501378031780457828919593905395747517585E-3Q, 615 6.355374424341762686099147452020466524659E-2Q, 616 3.539652320122661637429658698954748337223E-1Q, 617 9.571721066119617436343740541777014319695E-1Q, 618 1.196258777828426399432550698612171955305E0Q, 619 6.069388659458926158392384709893753793967E-1Q, 620 9.026746127269713176512359976978248763621E-2Q, 621 5.317668723070450235320878117210807236375E-4Q, 622 }; 623 #define NQ2_2r3D 8 624 static const __float128 Q2_2r3D[NQ2_2r3D + 1] = { 625 3.846924354014260866793741072933159380158E-5Q, 626 3.017562820057704325510067178327449946763E-3Q, 627 8.356305620686867949798885808540444210935E-2Q, 628 1.068314930499906838814019619594424586273E0Q, 629 6.900279623894821067017966573640732685233E0Q, 630 2.307667390886377924509090271780839563141E1Q, 631 3.921043465412723970791036825401273528513E1Q, 632 3.167569478939719383241775717095729233436E1Q, 633 1.051023841699200920276198346301543665909E1Q, 634 /* 1.000000000000000000000000000000000000000E0*/ 635 }; 636 637 638 /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */ 639 640 static __float128 641 neval (__float128 x, const __float128 *p, int n) 642 { 643 __float128 y; 644 645 p += n; 646 y = *p--; 647 do 648 { 649 y = y * x + *p--; 650 } 651 while (--n > 0); 652 return y; 653 } 654 655 656 /* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */ 657 658 static __float128 659 deval (__float128 x, const __float128 *p, int n) 660 { 661 __float128 y; 662 663 p += n; 664 y = x + *p--; 665 do 666 { 667 y = y * x + *p--; 668 } 669 while (--n > 0); 670 return y; 671 } 672 673 674 /* Bessel function of the first kind, order zero. */ 675 676 __float128 677 j0q (__float128 x) 678 { 679 __float128 xx, xinv, z, p, q, c, s, cc, ss; 680 681 if (! finiteq (x)) 682 { 683 if (x != x) 684 return x + x; 685 else 686 return 0; 687 } 688 if (x == 0) 689 return 1; 690 691 xx = fabsq (x); 692 if (xx <= 2) 693 { 694 if (xx < 0x1p-57Q) 695 return 1; 696 /* 0 <= x <= 2 */ 697 z = xx * xx; 698 p = z * z * neval (z, J0_2N, NJ0_2N) / deval (z, J0_2D, NJ0_2D); 699 p -= 0.25Q * z; 700 p += 1; 701 return p; 702 } 703 704 /* X = x - pi/4 705 cos(X) = cos(x) cos(pi/4) + sin(x) sin(pi/4) 706 = 1/sqrt(2) * (cos(x) + sin(x)) 707 sin(X) = sin(x) cos(pi/4) - cos(x) sin(pi/4) 708 = 1/sqrt(2) * (sin(x) - cos(x)) 709 sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) 710 cf. Fdlibm. */ 711 sincosq (xx, &s, &c); 712 ss = s - c; 713 cc = s + c; 714 if (xx <= FLT128_MAX / 2) 715 { 716 z = -cosq (xx + xx); 717 if ((s * c) < 0) 718 cc = z / ss; 719 else 720 ss = z / cc; 721 } 722 723 if (xx > 0x1p256Q) 724 return ONEOSQPI * cc / sqrtq (xx); 725 726 xinv = 1 / xx; 727 z = xinv * xinv; 728 if (xinv <= 0.25) 729 { 730 if (xinv <= 0.125) 731 { 732 if (xinv <= 0.0625) 733 { 734 p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID); 735 q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID); 736 } 737 else 738 { 739 p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D); 740 q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D); 741 } 742 } 743 else if (xinv <= 0.1875) 744 { 745 p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D); 746 q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D); 747 } 748 else 749 { 750 p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D); 751 q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D); 752 } 753 } /* .25 */ 754 else /* if (xinv <= 0.5) */ 755 { 756 if (xinv <= 0.375) 757 { 758 if (xinv <= 0.3125) 759 { 760 p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D); 761 q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D); 762 } 763 else 764 { 765 p = neval (z, P2r7_3r2N, NP2r7_3r2N) 766 / deval (z, P2r7_3r2D, NP2r7_3r2D); 767 q = neval (z, Q2r7_3r2N, NQ2r7_3r2N) 768 / deval (z, Q2r7_3r2D, NQ2r7_3r2D); 769 } 770 } 771 else if (xinv <= 0.4375) 772 { 773 p = neval (z, P2r3_2r7N, NP2r3_2r7N) 774 / deval (z, P2r3_2r7D, NP2r3_2r7D); 775 q = neval (z, Q2r3_2r7N, NQ2r3_2r7N) 776 / deval (z, Q2r3_2r7D, NQ2r3_2r7D); 777 } 778 else 779 { 780 p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D); 781 q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D); 782 } 783 } 784 p = 1 + z * p; 785 q = z * xinv * q; 786 q = q - 0.125Q * xinv; 787 z = ONEOSQPI * (p * cc - q * ss) / sqrtq (xx); 788 return z; 789 } 790 791 792 793 /* Y0(x) = 2/pi * log(x) * J0(x) + R(x^2) 794 Peak absolute error 1.7e-36 (relative where Y0 > 1) 795 0 <= x <= 2 */ 796 #define NY0_2N 7 797 static const __float128 Y0_2N[NY0_2N + 1] = { 798 -1.062023609591350692692296993537002558155E19Q, 799 2.542000883190248639104127452714966858866E19Q, 800 -1.984190771278515324281415820316054696545E18Q, 801 4.982586044371592942465373274440222033891E16Q, 802 -5.529326354780295177243773419090123407550E14Q, 803 3.013431465522152289279088265336861140391E12Q, 804 -7.959436160727126750732203098982718347785E9Q, 805 8.230845651379566339707130644134372793322E6Q, 806 }; 807 #define NY0_2D 7 808 static const __float128 Y0_2D[NY0_2D + 1] = { 809 1.438972634353286978700329883122253752192E20Q, 810 1.856409101981569254247700169486907405500E18Q, 811 1.219693352678218589553725579802986255614E16Q, 812 5.389428943282838648918475915779958097958E13Q, 813 1.774125762108874864433872173544743051653E11Q, 814 4.522104832545149534808218252434693007036E8Q, 815 8.872187401232943927082914504125234454930E5Q, 816 1.251945613186787532055610876304669413955E3Q, 817 /* 1.000000000000000000000000000000000000000E0 */ 818 }; 819 820 static const __float128 U0 = -7.3804295108687225274343927948483016310862e-02Q; 821 822 /* Bessel function of the second kind, order zero. */ 823 824 __float128 825 y0q(__float128 x) 826 { 827 __float128 xx, xinv, z, p, q, c, s, cc, ss; 828 829 if (! finiteq (x)) 830 return 1 / (x + x * x); 831 if (x <= 0) 832 { 833 if (x < 0) 834 return (zero / (zero * x)); 835 return -1 / zero; /* -inf and divide by zero exception. */ 836 } 837 xx = fabsq (x); 838 if (xx <= 0x1p-57) 839 return U0 + TWOOPI * logq (x); 840 if (xx <= 2) 841 { 842 /* 0 <= x <= 2 */ 843 z = xx * xx; 844 p = neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D); 845 p = TWOOPI * logq (x) * j0q (x) + p; 846 return p; 847 } 848 849 /* X = x - pi/4 850 cos(X) = cos(x) cos(pi/4) + sin(x) sin(pi/4) 851 = 1/sqrt(2) * (cos(x) + sin(x)) 852 sin(X) = sin(x) cos(pi/4) - cos(x) sin(pi/4) 853 = 1/sqrt(2) * (sin(x) - cos(x)) 854 sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) 855 cf. Fdlibm. */ 856 sincosq (x, &s, &c); 857 ss = s - c; 858 cc = s + c; 859 if (xx <= FLT128_MAX / 2) 860 { 861 z = -cosq (x + x); 862 if ((s * c) < 0) 863 cc = z / ss; 864 else 865 ss = z / cc; 866 } 867 868 if (xx > 0x1p256Q) 869 return ONEOSQPI * ss / sqrtq (x); 870 871 xinv = 1 / xx; 872 z = xinv * xinv; 873 if (xinv <= 0.25) 874 { 875 if (xinv <= 0.125) 876 { 877 if (xinv <= 0.0625) 878 { 879 p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID); 880 q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID); 881 } 882 else 883 { 884 p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D); 885 q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D); 886 } 887 } 888 else if (xinv <= 0.1875) 889 { 890 p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D); 891 q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D); 892 } 893 else 894 { 895 p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D); 896 q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D); 897 } 898 } /* .25 */ 899 else /* if (xinv <= 0.5) */ 900 { 901 if (xinv <= 0.375) 902 { 903 if (xinv <= 0.3125) 904 { 905 p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D); 906 q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D); 907 } 908 else 909 { 910 p = neval (z, P2r7_3r2N, NP2r7_3r2N) 911 / deval (z, P2r7_3r2D, NP2r7_3r2D); 912 q = neval (z, Q2r7_3r2N, NQ2r7_3r2N) 913 / deval (z, Q2r7_3r2D, NQ2r7_3r2D); 914 } 915 } 916 else if (xinv <= 0.4375) 917 { 918 p = neval (z, P2r3_2r7N, NP2r3_2r7N) 919 / deval (z, P2r3_2r7D, NP2r3_2r7D); 920 q = neval (z, Q2r3_2r7N, NQ2r3_2r7N) 921 / deval (z, Q2r3_2r7D, NQ2r3_2r7D); 922 } 923 else 924 { 925 p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D); 926 q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D); 927 } 928 } 929 p = 1 + z * p; 930 q = z * xinv * q; 931 q = q - 0.125Q * xinv; 932 z = ONEOSQPI * (p * ss + q * cc) / sqrtq (x); 933 return z; 934 } 935