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