]>
git.ipfire.org Git - thirdparty/glibc.git/blob - sysdeps/ieee754/ldbl-128/e_j1l.c
174f35df0762adaaca55384cd3e0e713c519c737
3 * Bessel function of order one
9 * long double x, y, j1l();
17 * Returns Bessel function of first kind, order one 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 is
21 * J1(x) = .5x + x x^2 R(x^2)
23 * The second interval is further partitioned into eight equal segments
25 * J1(x) = sqrt(2/(pi x)) (P1(x) cos(X) - Q1(x) sin(X)),
28 * and the auxiliary functions are given by
30 * J1(x)cos(X) + Y1(x)sin(X) = sqrt( 2/(pi x)) P1(x),
31 * P1(x) = 1 + 1/x^2 R(1/x^2)
33 * Y1(x)cos(X) - J1(x)sin(X) = sqrt( 2/(pi x)) Q1(x),
34 * Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)).
41 * arithmetic domain # trials peak rms
42 * IEEE 0, 30 100000 2.8e-34 2.7e-35
49 * Bessel function of the second kind, order one
63 * Returns Bessel function of the second kind, of order
64 * one, of the argument.
66 * The domain is divided into two major intervals [0, 2] and
67 * (2, infinity). In the first interval the rational approximation is
68 * Y1(x) = 2/pi * (log(x) * J1(x) - 1/x) + x R(x^2) .
69 * In the second interval the approximation is the same as for J1(x), and
70 * Y1(x) = sqrt(2/(pi x)) (P1(x) sin(X) + Q1(x) cos(X)),
75 * Absolute error, when y0(x) < 1; else relative error:
77 * arithmetic domain # trials peak rms
78 * IEEE 0, 30 100000 2.7e-34 2.9e-35
82 /* Copyright 2001 by Stephen L. Moshier (moshier@na-net.onrl.gov).
84 This library is free software; you can redistribute it and/or
85 modify it under the terms of the GNU Lesser General Public
86 License as published by the Free Software Foundation; either
87 version 2.1 of the License, or (at your option) any later version.
89 This library is distributed in the hope that it will be useful,
90 but WITHOUT ANY WARRANTY; without even the implied warranty of
91 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
92 Lesser General Public License for more details.
94 You should have received a copy of the GNU Lesser General Public
95 License along with this library; if not, see
96 <https://www.gnu.org/licenses/>. */
100 #include <math_private.h>
101 #include <fenv_private.h>
102 #include <math-underflow.h>
106 static const _Float128 ONEOSQPI
= L(5.6418958354775628694807945156077258584405E-1);
108 static const _Float128 TWOOPI
= L(6.3661977236758134307553505349005744813784E-1);
109 static const _Float128 zero
= 0;
111 /* J1(x) = .5x + x x^2 R(x^2)
112 Peak relative error 1.9e-35
115 static const _Float128 J0_2N
[NJ0_2N
+ 1] = {
116 L(-5.943799577386942855938508697619735179660E16
),
117 L(1.812087021305009192259946997014044074711E15
),
118 L(-2.761698314264509665075127515729146460895E13
),
119 L(2.091089497823600978949389109350658815972E11
),
120 L(-8.546413231387036372945453565654130054307E8
),
121 L(1.797229225249742247475464052741320612261E6
),
122 L(-1.559552840946694171346552770008812083969E3
)
125 static const _Float128 J0_2D
[NJ0_2D
+ 1] = {
126 L(9.510079323819108569501613916191477479397E17
),
127 L(1.063193817503280529676423936545854693915E16
),
128 L(5.934143516050192600795972192791775226920E13
),
129 L(2.168000911950620999091479265214368352883E11
),
130 L(5.673775894803172808323058205986256928794E8
),
131 L(1.080329960080981204840966206372671147224E6
),
132 L(1.411951256636576283942477881535283304912E3
),
133 /* 1.000000000000000000000000000000000000000E0L */
136 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
138 Peak relative error 3.6e-36 */
140 static const _Float128 P16_IN
[NP16_IN
+ 1] = {
141 L(5.143674369359646114999545149085139822905E-16),
142 L(4.836645664124562546056389268546233577376E-13),
143 L(1.730945562285804805325011561498453013673E-10),
144 L(3.047976856147077889834905908605310585810E-8),
145 L(2.855227609107969710407464739188141162386E-6),
146 L(1.439362407936705484122143713643023998457E-4),
147 L(3.774489768532936551500999699815873422073E-3),
148 L(4.723962172984642566142399678920790598426E-2),
149 L(2.359289678988743939925017240478818248735E-1),
150 L(3.032580002220628812728954785118117124520E-1),
153 static const _Float128 P16_ID
[NP16_ID
+ 1] = {
154 L(4.389268795186898018132945193912677177553E-15),
155 L(4.132671824807454334388868363256830961655E-12),
156 L(1.482133328179508835835963635130894413136E-9),
157 L(2.618941412861122118906353737117067376236E-7),
158 L(2.467854246740858470815714426201888034270E-5),
159 L(1.257192927368839847825938545925340230490E-3),
160 L(3.362739031941574274949719324644120720341E-2),
161 L(4.384458231338934105875343439265370178858E-1),
162 L(2.412830809841095249170909628197264854651E0
),
163 L(4.176078204111348059102962617368214856874E0
),
164 /* 1.000000000000000000000000000000000000000E0 */
167 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
168 0.0625 <= 1/x <= 0.125
169 Peak relative error 1.9e-36 */
171 static const _Float128 P8_16N
[NP8_16N
+ 1] = {
172 L(2.984612480763362345647303274082071598135E-16),
173 L(1.923651877544126103941232173085475682334E-13),
174 L(4.881258879388869396043760693256024307743E-11),
175 L(6.368866572475045408480898921866869811889E-9),
176 L(4.684818344104910450523906967821090796737E-7),
177 L(2.005177298271593587095982211091300382796E-5),
178 L(4.979808067163957634120681477207147536182E-4),
179 L(6.946005761642579085284689047091173581127E-3),
180 L(5.074601112955765012750207555985299026204E-2),
181 L(1.698599455896180893191766195194231825379E-1),
182 L(1.957536905259237627737222775573623779638E-1),
183 L(2.991314703282528370270179989044994319374E-2),
186 static const _Float128 P8_16D
[NP8_16D
+ 1] = {
187 L(2.546869316918069202079580939942463010937E-15),
188 L(1.644650111942455804019788382157745229955E-12),
189 L(4.185430770291694079925607420808011147173E-10),
190 L(5.485331966975218025368698195861074143153E-8),
191 L(4.062884421686912042335466327098932678905E-6),
192 L(1.758139661060905948870523641319556816772E-4),
193 L(4.445143889306356207566032244985607493096E-3),
194 L(6.391901016293512632765621532571159071158E-2),
195 L(4.933040207519900471177016015718145795434E-1),
196 L(1.839144086168947712971630337250761842976E0
),
197 L(2.715120873995490920415616716916149586579E0
),
198 /* 1.000000000000000000000000000000000000000E0 */
201 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
202 0.125 <= 1/x <= 0.1875
203 Peak relative error 1.3e-36 */
205 static const _Float128 P5_8N
[NP5_8N
+ 1] = {
206 L(2.837678373978003452653763806968237227234E-12),
207 L(9.726641165590364928442128579282742354806E-10),
208 L(1.284408003604131382028112171490633956539E-7),
209 L(8.524624695868291291250573339272194285008E-6),
210 L(3.111516908953172249853673787748841282846E-4),
211 L(6.423175156126364104172801983096596409176E-3),
212 L(7.430220589989104581004416356260692450652E-2),
213 L(4.608315409833682489016656279567605536619E-1),
214 L(1.396870223510964882676225042258855977512E0
),
215 L(1.718500293904122365894630460672081526236E0
),
216 L(5.465927698800862172307352821870223855365E-1)
219 static const _Float128 P5_8D
[NP5_8D
+ 1] = {
220 L(2.421485545794616609951168511612060482715E-11),
221 L(8.329862750896452929030058039752327232310E-9),
222 L(1.106137992233383429630592081375289010720E-6),
223 L(7.405786153760681090127497796448503306939E-5),
224 L(2.740364785433195322492093333127633465227E-3),
225 L(5.781246470403095224872243564165254652198E-2),
226 L(6.927711353039742469918754111511109983546E-1),
227 L(4.558679283460430281188304515922826156690E0
),
228 L(1.534468499844879487013168065728837900009E1
),
229 L(2.313927430889218597919624843161569422745E1
),
230 L(1.194506341319498844336768473218382828637E1
),
231 /* 1.000000000000000000000000000000000000000E0 */
234 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
235 Peak relative error 1.4e-36
236 0.1875 <= 1/x <= 0.25 */
238 static const _Float128 P4_5N
[NP4_5N
+ 1] = {
239 L(1.846029078268368685834261260420933914621E-10),
240 L(3.916295939611376119377869680335444207768E-8),
241 L(3.122158792018920627984597530935323997312E-6),
242 L(1.218073444893078303994045653603392272450E-4),
243 L(2.536420827983485448140477159977981844883E-3),
244 L(2.883011322006690823959367922241169171315E-2),
245 L(1.755255190734902907438042414495469810830E-1),
246 L(5.379317079922628599870898285488723736599E-1),
247 L(7.284904050194300773890303361501726561938E-1),
248 L(3.270110346613085348094396323925000362813E-1),
249 L(1.804473805689725610052078464951722064757E-2),
252 static const _Float128 P4_5D
[NP4_5D
+ 1] = {
253 L(1.575278146806816970152174364308980863569E-9),
254 L(3.361289173657099516191331123405675054321E-7),
255 L(2.704692281550877810424745289838790693708E-5),
256 L(1.070854930483999749316546199273521063543E-3),
257 L(2.282373093495295842598097265627962125411E-2),
258 L(2.692025460665354148328762368240343249830E-1),
259 L(1.739892942593664447220951225734811133759E0
),
260 L(5.890727576752230385342377570386657229324E0
),
261 L(9.517442287057841500750256954117735128153E0
),
262 L(6.100616353935338240775363403030137736013E0
),
263 /* 1.000000000000000000000000000000000000000E0 */
266 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
267 Peak relative error 3.0e-36
268 0.25 <= 1/x <= 0.3125 */
270 static const _Float128 P3r2_4N
[NP3r2_4N
+ 1] = {
271 L(8.240803130988044478595580300846665863782E-8),
272 L(1.179418958381961224222969866406483744580E-5),
273 L(6.179787320956386624336959112503824397755E-4),
274 L(1.540270833608687596420595830747166658383E-2),
275 L(1.983904219491512618376375619598837355076E-1),
276 L(1.341465722692038870390470651608301155565E0
),
277 L(4.617865326696612898792238245990854646057E0
),
278 L(7.435574801812346424460233180412308000587E0
),
279 L(4.671327027414635292514599201278557680420E0
),
280 L(7.299530852495776936690976966995187714739E-1),
283 static const _Float128 P3r2_4D
[NP3r2_4D
+ 1] = {
284 L(7.032152009675729604487575753279187576521E-7),
285 L(1.015090352324577615777511269928856742848E-4),
286 L(5.394262184808448484302067955186308730620E-3),
287 L(1.375291438480256110455809354836988584325E-1),
288 L(1.836247144461106304788160919310404376670E0
),
289 L(1.314378564254376655001094503090935880349E1
),
290 L(4.957184590465712006934452500894672343488E1
),
291 L(9.287394244300647738855415178790263465398E1
),
292 L(7.652563275535900609085229286020552768399E1
),
293 L(2.147042473003074533150718117770093209096E1
),
294 /* 1.000000000000000000000000000000000000000E0 */
297 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
298 Peak relative error 1.0e-35
299 0.3125 <= 1/x <= 0.375 */
301 static const _Float128 P2r7_3r2N
[NP2r7_3r2N
+ 1] = {
302 L(4.599033469240421554219816935160627085991E-7),
303 L(4.665724440345003914596647144630893997284E-5),
304 L(1.684348845667764271596142716944374892756E-3),
305 L(2.802446446884455707845985913454440176223E-2),
306 L(2.321937586453963310008279956042545173930E-1),
307 L(9.640277413988055668692438709376437553804E-1),
308 L(1.911021064710270904508663334033003246028E0
),
309 L(1.600811610164341450262992138893970224971E0
),
310 L(4.266299218652587901171386591543457861138E-1),
311 L(1.316470424456061252962568223251247207325E-2),
314 static const _Float128 P2r7_3r2D
[NP2r7_3r2D
+ 1] = {
315 L(3.924508608545520758883457108453520099610E-6),
316 L(4.029707889408829273226495756222078039823E-4),
317 L(1.484629715787703260797886463307469600219E-2),
318 L(2.553136379967180865331706538897231588685E-1),
319 L(2.229457223891676394409880026887106228740E0
),
320 L(1.005708903856384091956550845198392117318E1
),
321 L(2.277082659664386953166629360352385889558E1
),
322 L(2.384726835193630788249826630376533988245E1
),
323 L(9.700989749041320895890113781610939632410E0
),
324 /* 1.000000000000000000000000000000000000000E0 */
327 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
328 Peak relative error 1.7e-36
329 0.3125 <= 1/x <= 0.4375 */
331 static const _Float128 P2r3_2r7N
[NP2r3_2r7N
+ 1] = {
332 L(3.916766777108274628543759603786857387402E-6),
333 L(3.212176636756546217390661984304645137013E-4),
334 L(9.255768488524816445220126081207248947118E-3),
335 L(1.214853146369078277453080641911700735354E-1),
336 L(7.855163309847214136198449861311404633665E-1),
337 L(2.520058073282978403655488662066019816540E0
),
338 L(3.825136484837545257209234285382183711466E0
),
339 L(2.432569427554248006229715163865569506873E0
),
340 L(4.877934835018231178495030117729800489743E-1),
341 L(1.109902737860249670981355149101343427885E-2),
344 static const _Float128 P2r3_2r7D
[NP2r3_2r7D
+ 1] = {
345 L(3.342307880794065640312646341190547184461E-5),
346 L(2.782182891138893201544978009012096558265E-3),
347 L(8.221304931614200702142049236141249929207E-2),
348 L(1.123728246291165812392918571987858010949E0
),
349 L(7.740482453652715577233858317133423434590E0
),
350 L(2.737624677567945952953322566311201919139E1
),
351 L(4.837181477096062403118304137851260715475E1
),
352 L(3.941098643468580791437772701093795299274E1
),
353 L(1.245821247166544627558323920382547533630E1
),
354 /* 1.000000000000000000000000000000000000000E0 */
357 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
358 Peak relative error 1.7e-35
359 0.4375 <= 1/x <= 0.5 */
361 static const _Float128 P2_2r3N
[NP2_2r3N
+ 1] = {
362 L(3.397930802851248553545191160608731940751E-4),
363 L(2.104020902735482418784312825637833698217E-2),
364 L(4.442291771608095963935342749477836181939E-1),
365 L(4.131797328716583282869183304291833754967E0
),
366 L(1.819920169779026500146134832455189917589E1
),
367 L(3.781779616522937565300309684282401791291E1
),
368 L(3.459605449728864218972931220783543410347E1
),
369 L(1.173594248397603882049066603238568316561E1
),
370 L(9.455702270242780642835086549285560316461E-1),
373 static const _Float128 P2_2r3D
[NP2_2r3D
+ 1] = {
374 L(2.899568897241432883079888249845707400614E-3),
375 L(1.831107138190848460767699919531132426356E-1),
376 L(3.999350044057883839080258832758908825165E0
),
377 L(3.929041535867957938340569419874195303712E1
),
378 L(1.884245613422523323068802689915538908291E2
),
379 L(4.461469948819229734353852978424629815929E2
),
380 L(5.004998753999796821224085972610636347903E2
),
381 L(2.386342520092608513170837883757163414100E2
),
382 L(3.791322528149347975999851588922424189957E1
),
383 /* 1.000000000000000000000000000000000000000E0 */
386 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
387 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
388 Peak relative error 8.0e-36
391 static const _Float128 Q16_IN
[NQ16_IN
+ 1] = {
392 L(-3.917420835712508001321875734030357393421E-18),
393 L(-4.440311387483014485304387406538069930457E-15),
394 L(-1.951635424076926487780929645954007139616E-12),
395 L(-4.318256438421012555040546775651612810513E-10),
396 L(-5.231244131926180765270446557146989238020E-8),
397 L(-3.540072702902043752460711989234732357653E-6),
398 L(-1.311017536555269966928228052917534882984E-4),
399 L(-2.495184669674631806622008769674827575088E-3),
400 L(-2.141868222987209028118086708697998506716E-2),
401 L(-6.184031415202148901863605871197272650090E-2),
402 L(-1.922298704033332356899546792898156493887E-2),
405 static const _Float128 Q16_ID
[NQ16_ID
+ 1] = {
406 L(3.820418034066293517479619763498400162314E-17),
407 L(4.340702810799239909648911373329149354911E-14),
408 L(1.914985356383416140706179933075303538524E-11),
409 L(4.262333682610888819476498617261895474330E-9),
410 L(5.213481314722233980346462747902942182792E-7),
411 L(3.585741697694069399299005316809954590558E-5),
412 L(1.366513429642842006385029778105539457546E-3),
413 L(2.745282599850704662726337474371355160594E-2),
414 L(2.637644521611867647651200098449903330074E-1),
415 L(1.006953426110765984590782655598680488746E0
),
416 /* 1.000000000000000000000000000000000000000E0 */
419 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
420 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
421 Peak relative error 1.9e-36
422 0.0625 <= 1/x <= 0.125 */
424 static const _Float128 Q8_16N
[NQ8_16N
+ 1] = {
425 L(-2.028630366670228670781362543615221542291E-17),
426 L(-1.519634620380959966438130374006858864624E-14),
427 L(-4.540596528116104986388796594639405114524E-12),
428 L(-7.085151756671466559280490913558388648274E-10),
429 L(-6.351062671323970823761883833531546885452E-8),
430 L(-3.390817171111032905297982523519503522491E-6),
431 L(-1.082340897018886970282138836861233213972E-4),
432 L(-2.020120801187226444822977006648252379508E-3),
433 L(-2.093169910981725694937457070649605557555E-2),
434 L(-1.092176538874275712359269481414448063393E-1),
435 L(-2.374790947854765809203590474789108718733E-1),
436 L(-1.365364204556573800719985118029601401323E-1),
439 static const _Float128 Q8_16D
[NQ8_16D
+ 1] = {
440 L(1.978397614733632533581207058069628242280E-16),
441 L(1.487361156806202736877009608336766720560E-13),
442 L(4.468041406888412086042576067133365913456E-11),
443 L(7.027822074821007443672290507210594648877E-9),
444 L(6.375740580686101224127290062867976007374E-7),
445 L(3.466887658320002225888644977076410421940E-5),
446 L(1.138625640905289601186353909213719596986E-3),
447 L(2.224470799470414663443449818235008486439E-2),
448 L(2.487052928527244907490589787691478482358E-1),
449 L(1.483927406564349124649083853892380899217E0
),
450 L(4.182773513276056975777258788903489507705E0
),
451 L(4.419665392573449746043880892524360870944E0
),
452 /* 1.000000000000000000000000000000000000000E0 */
455 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
456 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
457 Peak relative error 1.5e-35
458 0.125 <= 1/x <= 0.1875 */
460 static const _Float128 Q5_8N
[NQ5_8N
+ 1] = {
461 L(-3.656082407740970534915918390488336879763E-13),
462 L(-1.344660308497244804752334556734121771023E-10),
463 L(-1.909765035234071738548629788698150760791E-8),
464 L(-1.366668038160120210269389551283666716453E-6),
465 L(-5.392327355984269366895210704976314135683E-5),
466 L(-1.206268245713024564674432357634540343884E-3),
467 L(-1.515456784370354374066417703736088291287E-2),
468 L(-1.022454301137286306933217746545237098518E-1),
469 L(-3.373438906472495080504907858424251082240E-1),
470 L(-4.510782522110845697262323973549178453405E-1),
471 L(-1.549000892545288676809660828213589804884E-1),
474 static const _Float128 Q5_8D
[NQ5_8D
+ 1] = {
475 L(3.565550843359501079050699598913828460036E-12),
476 L(1.321016015556560621591847454285330528045E-9),
477 L(1.897542728662346479999969679234270605975E-7),
478 L(1.381720283068706710298734234287456219474E-5),
479 L(5.599248147286524662305325795203422873725E-4),
480 L(1.305442352653121436697064782499122164843E-2),
481 L(1.750234079626943298160445750078631894985E-1),
482 L(1.311420542073436520965439883806946678491E0
),
483 L(5.162757689856842406744504211089724926650E0
),
484 L(9.527760296384704425618556332087850581308E0
),
485 L(6.604648207463236667912921642545100248584E0
),
486 /* 1.000000000000000000000000000000000000000E0 */
489 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
490 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
491 Peak relative error 1.3e-35
492 0.1875 <= 1/x <= 0.25 */
494 static const _Float128 Q4_5N
[NQ4_5N
+ 1] = {
495 L(-4.079513568708891749424783046520200903755E-11),
496 L(-9.326548104106791766891812583019664893311E-9),
497 L(-8.016795121318423066292906123815687003356E-7),
498 L(-3.372350544043594415609295225664186750995E-5),
499 L(-7.566238665947967882207277686375417983917E-4),
500 L(-9.248861580055565402130441618521591282617E-3),
501 L(-6.033106131055851432267702948850231270338E-2),
502 L(-1.966908754799996793730369265431584303447E-1),
503 L(-2.791062741179964150755788226623462207560E-1),
504 L(-1.255478605849190549914610121863534191666E-1),
505 L(-4.320429862021265463213168186061696944062E-3),
508 static const _Float128 Q4_5D
[NQ4_5D
+ 1] = {
509 L(3.978497042580921479003851216297330701056E-10),
510 L(9.203304163828145809278568906420772246666E-8),
511 L(8.059685467088175644915010485174545743798E-6),
512 L(3.490187375993956409171098277561669167446E-4),
513 L(8.189109654456872150100501732073810028829E-3),
514 L(1.072572867311023640958725265762483033769E-1),
515 L(7.790606862409960053675717185714576937994E-1),
516 L(3.016049768232011196434185423512777656328E0
),
517 L(5.722963851442769787733717162314477949360E0
),
518 L(4.510527838428473279647251350931380867663E0
),
519 /* 1.000000000000000000000000000000000000000E0 */
522 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
523 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
524 Peak relative error 2.1e-35
525 0.25 <= 1/x <= 0.3125 */
527 static const _Float128 Q3r2_4N
[NQ3r2_4N
+ 1] = {
528 L(-1.087480809271383885936921889040388133627E-8),
529 L(-1.690067828697463740906962973479310170932E-6),
530 L(-9.608064416995105532790745641974762550982E-5),
531 L(-2.594198839156517191858208513873961837410E-3),
532 L(-3.610954144421543968160459863048062977822E-2),
533 L(-2.629866798251843212210482269563961685666E-1),
534 L(-9.709186825881775885917984975685752956660E-1),
535 L(-1.667521829918185121727268867619982417317E0
),
536 L(-1.109255082925540057138766105229900943501E0
),
537 L(-1.812932453006641348145049323713469043328E-1),
540 static const _Float128 Q3r2_4D
[NQ3r2_4D
+ 1] = {
541 L(1.060552717496912381388763753841473407026E-7),
542 L(1.676928002024920520786883649102388708024E-5),
543 L(9.803481712245420839301400601140812255737E-4),
544 L(2.765559874262309494758505158089249012930E-2),
545 L(4.117921827792571791298862613287549140706E-1),
546 L(3.323769515244751267093378361930279161413E0
),
547 L(1.436602494405814164724810151689705353670E1
),
548 L(3.163087869617098638064881410646782408297E1
),
549 L(3.198181264977021649489103980298349589419E1
),
550 L(1.203649258862068431199471076202897823272E1
),
551 /* 1.000000000000000000000000000000000000000E0 */
554 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
555 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
556 Peak relative error 1.6e-36
557 0.3125 <= 1/x <= 0.375 */
559 static const _Float128 Q2r7_3r2N
[NQ2r7_3r2N
+ 1] = {
560 L(-1.723405393982209853244278760171643219530E-7),
561 L(-2.090508758514655456365709712333460087442E-5),
562 L(-9.140104013370974823232873472192719263019E-4),
563 L(-1.871349499990714843332742160292474780128E-2),
564 L(-1.948930738119938669637865956162512983416E-1),
565 L(-1.048764684978978127908439526343174139788E0
),
566 L(-2.827714929925679500237476105843643064698E0
),
567 L(-3.508761569156476114276988181329773987314E0
),
568 L(-1.669332202790211090973255098624488308989E0
),
569 L(-1.930796319299022954013840684651016077770E-1),
572 static const _Float128 Q2r7_3r2D
[NQ2r7_3r2D
+ 1] = {
573 L(1.680730662300831976234547482334347983474E-6),
574 L(2.084241442440551016475972218719621841120E-4),
575 L(9.445316642108367479043541702688736295579E-3),
576 L(2.044637889456631896650179477133252184672E-1),
577 L(2.316091982244297350829522534435350078205E0
),
578 L(1.412031891783015085196708811890448488865E1
),
579 L(4.583830154673223384837091077279595496149E1
),
580 L(7.549520609270909439885998474045974122261E1
),
581 L(5.697605832808113367197494052388203310638E1
),
582 L(1.601496240876192444526383314589371686234E1
),
583 /* 1.000000000000000000000000000000000000000E0 */
586 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
587 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
588 Peak relative error 9.5e-36
589 0.375 <= 1/x <= 0.4375 */
591 static const _Float128 Q2r3_2r7N
[NQ2r3_2r7N
+ 1] = {
592 L(-8.603042076329122085722385914954878953775E-7),
593 L(-7.701746260451647874214968882605186675720E-5),
594 L(-2.407932004380727587382493696877569654271E-3),
595 L(-3.403434217607634279028110636919987224188E-2),
596 L(-2.348707332185238159192422084985713102877E-1),
597 L(-7.957498841538254916147095255700637463207E-1),
598 L(-1.258469078442635106431098063707934348577E0
),
599 L(-8.162415474676345812459353639449971369890E-1),
600 L(-1.581783890269379690141513949609572806898E-1),
601 L(-1.890595651683552228232308756569450822905E-3),
604 static const _Float128 Q2r3_2r7D
[NQ2r3_2r7D
+ 1] = {
605 L(8.390017524798316921170710533381568175665E-6),
606 L(7.738148683730826286477254659973968763659E-4),
607 L(2.541480810958665794368759558791634341779E-2),
608 L(3.878879789711276799058486068562386244873E-1),
609 L(3.003783779325811292142957336802456109333E0
),
610 L(1.206480374773322029883039064575464497400E1
),
611 L(2.458414064785315978408974662900438351782E1
),
612 L(2.367237826273668567199042088835448715228E1
),
613 L(9.231451197519171090875569102116321676763E0
),
614 /* 1.000000000000000000000000000000000000000E0 */
617 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
618 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
619 Peak relative error 1.4e-36
620 0.4375 <= 1/x <= 0.5 */
622 static const _Float128 Q2_2r3N
[NQ2_2r3N
+ 1] = {
623 L(-5.552507516089087822166822364590806076174E-6),
624 L(-4.135067659799500521040944087433752970297E-4),
625 L(-1.059928728869218962607068840646564457980E-2),
626 L(-1.212070036005832342565792241385459023801E-1),
627 L(-6.688350110633603958684302153362735625156E-1),
628 L(-1.793587878197360221340277951304429821582E0
),
629 L(-2.225407682237197485644647380483725045326E0
),
630 L(-1.123402135458940189438898496348239744403E0
),
631 L(-1.679187241566347077204805190763597299805E-1),
632 L(-1.458550613639093752909985189067233504148E-3),
635 static const _Float128 Q2_2r3D
[NQ2_2r3D
+ 1] = {
636 L(5.415024336507980465169023996403597916115E-5),
637 L(4.179246497380453022046357404266022870788E-3),
638 L(1.136306384261959483095442402929502368598E-1),
639 L(1.422640343719842213484515445393284072830E0
),
640 L(8.968786703393158374728850922289204805764E0
),
641 L(2.914542473339246127533384118781216495934E1
),
642 L(4.781605421020380669870197378210457054685E1
),
643 L(3.693865837171883152382820584714795072937E1
),
644 L(1.153220502744204904763115556224395893076E1
),
645 /* 1.000000000000000000000000000000000000000E0 */
649 /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
652 neval (_Float128 x
, const _Float128
*p
, int n
)
667 /* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */
670 deval (_Float128 x
, const _Float128
*p
, int n
)
685 /* Bessel function of the first kind, order one. */
688 __ieee754_j1l (_Float128 x
)
690 _Float128 xx
, xinv
, z
, p
, q
, c
, s
, cc
, ss
;
702 if (xx
<= L(0x1p
-58))
704 _Float128 ret
= x
* L(0.5);
705 math_check_force_underflow (ret
);
707 __set_errno (ERANGE
);
714 p
= xx
* z
* neval (z
, J0_2N
, NJ0_2N
) / deval (z
, J0_2D
, NJ0_2D
);
722 cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
723 = 1/sqrt(2) * (-cos(x) + sin(x))
724 sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
725 = -1/sqrt(2) * (sin(x) + cos(x))
727 __sincosl (xx
, &s
, &c
);
730 if (xx
<= LDBL_MAX
/ 2)
732 z
= __cosl (xx
+ xx
);
741 z
= ONEOSQPI
* cc
/ sqrtl (xx
);
755 p
= neval (z
, P16_IN
, NP16_IN
) / deval (z
, P16_ID
, NP16_ID
);
756 q
= neval (z
, Q16_IN
, NQ16_IN
) / deval (z
, Q16_ID
, NQ16_ID
);
760 p
= neval (z
, P8_16N
, NP8_16N
) / deval (z
, P8_16D
, NP8_16D
);
761 q
= neval (z
, Q8_16N
, NQ8_16N
) / deval (z
, Q8_16D
, NQ8_16D
);
764 else if (xinv
<= 0.1875)
766 p
= neval (z
, P5_8N
, NP5_8N
) / deval (z
, P5_8D
, NP5_8D
);
767 q
= neval (z
, Q5_8N
, NQ5_8N
) / deval (z
, Q5_8D
, NQ5_8D
);
771 p
= neval (z
, P4_5N
, NP4_5N
) / deval (z
, P4_5D
, NP4_5D
);
772 q
= neval (z
, Q4_5N
, NQ4_5N
) / deval (z
, Q4_5D
, NQ4_5D
);
775 else /* if (xinv <= 0.5) */
781 p
= neval (z
, P3r2_4N
, NP3r2_4N
) / deval (z
, P3r2_4D
, NP3r2_4D
);
782 q
= neval (z
, Q3r2_4N
, NQ3r2_4N
) / deval (z
, Q3r2_4D
, NQ3r2_4D
);
786 p
= neval (z
, P2r7_3r2N
, NP2r7_3r2N
)
787 / deval (z
, P2r7_3r2D
, NP2r7_3r2D
);
788 q
= neval (z
, Q2r7_3r2N
, NQ2r7_3r2N
)
789 / deval (z
, Q2r7_3r2D
, NQ2r7_3r2D
);
792 else if (xinv
<= 0.4375)
794 p
= neval (z
, P2r3_2r7N
, NP2r3_2r7N
)
795 / deval (z
, P2r3_2r7D
, NP2r3_2r7D
);
796 q
= neval (z
, Q2r3_2r7N
, NQ2r3_2r7N
)
797 / deval (z
, Q2r3_2r7D
, NQ2r3_2r7D
);
801 p
= neval (z
, P2_2r3N
, NP2_2r3N
) / deval (z
, P2_2r3D
, NP2_2r3D
);
802 q
= neval (z
, Q2_2r3N
, NQ2_2r3N
) / deval (z
, Q2_2r3D
, NQ2_2r3D
);
807 q
= q
* xinv
+ L(0.375) * xinv
;
808 z
= ONEOSQPI
* (p
* cc
- q
* ss
) / sqrtl (xx
);
813 strong_alias (__ieee754_j1l
, __j1l_finite
)
816 /* Y1(x) = 2/pi * (log(x) * J1(x) - 1/x) + x R(x^2)
817 Peak relative error 6.2e-38
820 static const _Float128 Y0_2N
[NY0_2N
+ 1] = {
821 L(-6.804415404830253804408698161694720833249E19
),
822 L(1.805450517967019908027153056150465849237E19
),
823 L(-8.065747497063694098810419456383006737312E17
),
824 L(1.401336667383028259295830955439028236299E16
),
825 L(-1.171654432898137585000399489686629680230E14
),
826 L(5.061267920943853732895341125243428129150E11
),
827 L(-1.096677850566094204586208610960870217970E9
),
828 L(9.541172044989995856117187515882879304461E5
),
831 static const _Float128 Y0_2D
[NY0_2D
+ 1] = {
832 L(3.470629591820267059538637461549677594549E20
),
833 L(4.120796439009916326855848107545425217219E18
),
834 L(2.477653371652018249749350657387030814542E16
),
835 L(9.954678543353888958177169349272167762797E13
),
836 L(2.957927997613630118216218290262851197754E11
),
837 L(6.748421382188864486018861197614025972118E8
),
838 L(1.173453425218010888004562071020305709319E6
),
839 L(1.450335662961034949894009554536003377187E3
),
840 /* 1.000000000000000000000000000000000000000E0 */
844 /* Bessel function of the second kind, order one. */
847 __ieee754_y1l (_Float128 x
)
849 _Float128 xx
, xinv
, z
, p
, q
, c
, s
, cc
, ss
;
852 return 1 / (x
+ x
* x
);
856 return (zero
/ (zero
* x
));
857 return -1 / zero
; /* -inf and divide by zero exception. */
864 __set_errno (ERANGE
);
870 SET_RESTORE_ROUNDL (FE_TONEAREST
);
872 p
= xx
* neval (z
, Y0_2N
, NY0_2N
) / deval (z
, Y0_2D
, NY0_2D
);
873 p
= -TWOOPI
/ xx
+ p
;
874 p
= TWOOPI
* __ieee754_logl (x
) * __ieee754_j1l (x
) + p
;
879 cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
880 = 1/sqrt(2) * (-cos(x) + sin(x))
881 sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
882 = -1/sqrt(2) * (sin(x) + cos(x))
884 __sincosl (xx
, &s
, &c
);
887 if (xx
<= LDBL_MAX
/ 2)
889 z
= __cosl (xx
+ xx
);
897 return ONEOSQPI
* ss
/ sqrtl (xx
);
907 p
= neval (z
, P16_IN
, NP16_IN
) / deval (z
, P16_ID
, NP16_ID
);
908 q
= neval (z
, Q16_IN
, NQ16_IN
) / deval (z
, Q16_ID
, NQ16_ID
);
912 p
= neval (z
, P8_16N
, NP8_16N
) / deval (z
, P8_16D
, NP8_16D
);
913 q
= neval (z
, Q8_16N
, NQ8_16N
) / deval (z
, Q8_16D
, NQ8_16D
);
916 else if (xinv
<= 0.1875)
918 p
= neval (z
, P5_8N
, NP5_8N
) / deval (z
, P5_8D
, NP5_8D
);
919 q
= neval (z
, Q5_8N
, NQ5_8N
) / deval (z
, Q5_8D
, NQ5_8D
);
923 p
= neval (z
, P4_5N
, NP4_5N
) / deval (z
, P4_5D
, NP4_5D
);
924 q
= neval (z
, Q4_5N
, NQ4_5N
) / deval (z
, Q4_5D
, NQ4_5D
);
927 else /* if (xinv <= 0.5) */
933 p
= neval (z
, P3r2_4N
, NP3r2_4N
) / deval (z
, P3r2_4D
, NP3r2_4D
);
934 q
= neval (z
, Q3r2_4N
, NQ3r2_4N
) / deval (z
, Q3r2_4D
, NQ3r2_4D
);
938 p
= neval (z
, P2r7_3r2N
, NP2r7_3r2N
)
939 / deval (z
, P2r7_3r2D
, NP2r7_3r2D
);
940 q
= neval (z
, Q2r7_3r2N
, NQ2r7_3r2N
)
941 / deval (z
, Q2r7_3r2D
, NQ2r7_3r2D
);
944 else if (xinv
<= 0.4375)
946 p
= neval (z
, P2r3_2r7N
, NP2r3_2r7N
)
947 / deval (z
, P2r3_2r7D
, NP2r3_2r7D
);
948 q
= neval (z
, Q2r3_2r7N
, NQ2r3_2r7N
)
949 / deval (z
, Q2r3_2r7D
, NQ2r3_2r7D
);
953 p
= neval (z
, P2_2r3N
, NP2_2r3N
) / deval (z
, P2_2r3D
, NP2_2r3D
);
954 q
= neval (z
, Q2_2r3N
, NQ2_2r3N
) / deval (z
, Q2_2r3D
, NQ2_2r3D
);
959 q
= q
* xinv
+ L(0.375) * xinv
;
960 z
= ONEOSQPI
* (p
* ss
+ q
* cc
) / sqrtl (xx
);
963 strong_alias (__ieee754_y1l
, __y1l_finite
)