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