3 * Bessel function of order zero
9 * __float128 x, y, j0l();
17 * Returns Bessel function of first kind, order zero of the argument.
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
25 * J0(x) = sqrt(2/(pi x)) (P0(x) cos(X) - Q0(x) sin(X)),
28 * and the auxiliary functions are given by
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)
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))
41 * arithmetic domain # trials peak rms
42 * IEEE 0, 30 100000 1.7e-34 2.4e-35
49 * Bessel function of the second kind, order zero
55 * __float128 x, y, y0l();
63 * Returns Bessel function of the second kind, of order
64 * zero, of the argument.
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)).
71 * Absolute error, when y0(x) < 1; else relative error:
73 * arithmetic domain # trials peak rms
74 * IEEE 0, 30 100000 3.0e-34 2.7e-35
78 /* Copyright 2001 by Stephen L. Moshier (moshier@na-net.ornl.gov).
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.
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.
90 You should have received a copy of the GNU Lesser General Public
91 License along with this library; if not, write to the Free Software
92 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
94 #include "quadmath-imp.h"
97 static const __float128 ONEOSQPI
= 5.6418958354775628694807945156077258584405E-1Q
;
99 static const __float128 TWOOPI
= 6.3661977236758134307553505349005744813784E-1Q
;
100 static const __float128 zero
= 0.0Q
;
102 /* J0(x) = 1 - x^2/4 + x^2 x^2 R(x^2)
103 Peak relative error 3.4e-37
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
,
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 */
127 /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2),
129 Peak relative error 3.3e-36 */
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
,
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 */
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 */
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
,
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 */
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 */
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
,
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 */
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 */
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
,
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 */
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 */
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
,
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 */
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 */
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
,
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 */
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 */
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
,
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 */
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 */
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
,
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 */
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
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
,
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 */
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 */
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
,
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 */
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 */
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
,
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 */
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 */
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
,
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 */
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 */
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
,
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*/
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 */
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
,
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 */
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 */
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
,
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 */
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 */
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
,
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*/
638 /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
641 neval (__float128 x
, const __float128
*p
, int n
)
656 /* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */
659 deval (__float128 x
, const __float128
*p
, int n
)
674 /* Bessel function of the first kind, order zero. */
679 __float128 xx
, xinv
, z
, p
, q
, c
, s
, cc
, ss
;
698 p
= z
* z
* neval (z
, J0_2N
, NJ0_2N
) / deval (z
, J0_2D
, NJ0_2D
);
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))
711 sincosq (xx
, &s
, &c
);
714 if (xx
<= FLT128_MAX
/ 2.0Q
)
724 return ONEOSQPI
* cc
/ sqrtq (xx
);
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
);
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
);
743 else if (xinv
<= 0.1875)
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
);
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
);
754 else /* if (xinv <= 0.5) */
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
);
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
);
771 else if (xinv
<= 0.4375)
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
);
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
);
786 q
= q
- 0.125Q
* xinv
;
787 z
= ONEOSQPI
* (p
* cc
- q
* ss
) / sqrtq (xx
);
792 /* Y0(x) = 2/pi * log(x) * J0(x) + R(x^2)
793 Peak absolute error 1.7e-36 (relative where Y0 > 1)
796 static __float128 Y0_2N
[NY0_2N
+ 1] = {
797 -1.062023609591350692692296993537002558155E19Q
,
798 2.542000883190248639104127452714966858866E19Q
,
799 -1.984190771278515324281415820316054696545E18Q
,
800 4.982586044371592942465373274440222033891E16Q
,
801 -5.529326354780295177243773419090123407550E14Q
,
802 3.013431465522152289279088265336861140391E12Q
,
803 -7.959436160727126750732203098982718347785E9Q
,
804 8.230845651379566339707130644134372793322E6Q
,
807 static __float128 Y0_2D
[NY0_2D
+ 1] = {
808 1.438972634353286978700329883122253752192E20Q
,
809 1.856409101981569254247700169486907405500E18Q
,
810 1.219693352678218589553725579802986255614E16Q
,
811 5.389428943282838648918475915779958097958E13Q
,
812 1.774125762108874864433872173544743051653E11Q
,
813 4.522104832545149534808218252434693007036E8Q
,
814 8.872187401232943927082914504125234454930E5Q
,
815 1.251945613186787532055610876304669413955E3Q
,
816 /* 1.000000000000000000000000000000000000000E0 */
819 static const __float128 U0
= -7.3804295108687225274343927948483016310862e-02Q
;
821 /* Bessel function of the second kind, order zero. */
826 __float128 xx
, xinv
, z
, p
, q
, c
, s
, cc
, ss
;
829 return 1 / (x
+ x
* x
);
833 return (zero
/ (zero
* x
));
834 return -1 / zero
; /* -inf and divide by zero exception. */
838 return U0
+ TWOOPI
* logq (x
);
843 p
= neval (z
, Y0_2N
, NY0_2N
) / deval (z
, Y0_2D
, NY0_2D
);
844 p
= TWOOPI
* logq (x
) * j0q (x
) + p
;
849 cos(X) = cos(x) cos(pi/4) + sin(x) sin(pi/4)
850 = 1/sqrt(2) * (cos(x) + sin(x))
851 sin(X) = sin(x) cos(pi/4) - cos(x) sin(pi/4)
852 = 1/sqrt(2) * (sin(x) - cos(x))
853 sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
858 if (xx
<= FLT128_MAX
/ 2.0Q
)
868 return ONEOSQPI
* ss
/ sqrtq (x
);
878 p
= neval (z
, P16_IN
, NP16_IN
) / deval (z
, P16_ID
, NP16_ID
);
879 q
= neval (z
, Q16_IN
, NQ16_IN
) / deval (z
, Q16_ID
, NQ16_ID
);
883 p
= neval (z
, P8_16N
, NP8_16N
) / deval (z
, P8_16D
, NP8_16D
);
884 q
= neval (z
, Q8_16N
, NQ8_16N
) / deval (z
, Q8_16D
, NQ8_16D
);
887 else if (xinv
<= 0.1875)
889 p
= neval (z
, P5_8N
, NP5_8N
) / deval (z
, P5_8D
, NP5_8D
);
890 q
= neval (z
, Q5_8N
, NQ5_8N
) / deval (z
, Q5_8D
, NQ5_8D
);
894 p
= neval (z
, P4_5N
, NP4_5N
) / deval (z
, P4_5D
, NP4_5D
);
895 q
= neval (z
, Q4_5N
, NQ4_5N
) / deval (z
, Q4_5D
, NQ4_5D
);
898 else /* if (xinv <= 0.5) */
904 p
= neval (z
, P3r2_4N
, NP3r2_4N
) / deval (z
, P3r2_4D
, NP3r2_4D
);
905 q
= neval (z
, Q3r2_4N
, NQ3r2_4N
) / deval (z
, Q3r2_4D
, NQ3r2_4D
);
909 p
= neval (z
, P2r7_3r2N
, NP2r7_3r2N
)
910 / deval (z
, P2r7_3r2D
, NP2r7_3r2D
);
911 q
= neval (z
, Q2r7_3r2N
, NQ2r7_3r2N
)
912 / deval (z
, Q2r7_3r2D
, NQ2r7_3r2D
);
915 else if (xinv
<= 0.4375)
917 p
= neval (z
, P2r3_2r7N
, NP2r3_2r7N
)
918 / deval (z
, P2r3_2r7D
, NP2r3_2r7D
);
919 q
= neval (z
, Q2r3_2r7N
, NQ2r3_2r7N
)
920 / deval (z
, Q2r3_2r7D
, NQ2r3_2r7D
);
924 p
= neval (z
, P2_2r3N
, NP2_2r3N
) / deval (z
, P2_2r3D
, NP2_2r3D
);
925 q
= neval (z
, Q2_2r3N
, NQ2_2r3N
) / deval (z
, Q2_2r3D
, NQ2_2r3D
);
930 q
= q
- 0.125Q
* xinv
;
931 z
= ONEOSQPI
* (p
* ss
+ q
* cc
) / sqrtq (x
);