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, write to the Free Software
96 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
99 #include "math_private.h"
102 static const long double ONEOSQPI
= 5.6418958354775628694807945156077258584405E-1L;
104 static const long double TWOOPI
= 6.3661977236758134307553505349005744813784E-1L;
105 static const long double zero
= 0.0L;
107 /* J1(x) = .5x + x x^2 R(x^2)
108 Peak relative error 1.9e-35
111 static const long double J0_2N
[NJ0_2N
+ 1] = {
112 -5.943799577386942855938508697619735179660E16L
,
113 1.812087021305009192259946997014044074711E15L
,
114 -2.761698314264509665075127515729146460895E13L
,
115 2.091089497823600978949389109350658815972E11L
,
116 -8.546413231387036372945453565654130054307E8L
,
117 1.797229225249742247475464052741320612261E6L
,
118 -1.559552840946694171346552770008812083969E3L
121 static const long double J0_2D
[NJ0_2D
+ 1] = {
122 9.510079323819108569501613916191477479397E17L
,
123 1.063193817503280529676423936545854693915E16L
,
124 5.934143516050192600795972192791775226920E13L
,
125 2.168000911950620999091479265214368352883E11L
,
126 5.673775894803172808323058205986256928794E8L
,
127 1.080329960080981204840966206372671147224E6L
,
128 1.411951256636576283942477881535283304912E3L
,
129 /* 1.000000000000000000000000000000000000000E0L */
132 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
134 Peak relative error 3.6e-36 */
136 static const long double P16_IN
[NP16_IN
+ 1] = {
137 5.143674369359646114999545149085139822905E-16L,
138 4.836645664124562546056389268546233577376E-13L,
139 1.730945562285804805325011561498453013673E-10L,
140 3.047976856147077889834905908605310585810E-8L,
141 2.855227609107969710407464739188141162386E-6L,
142 1.439362407936705484122143713643023998457E-4L,
143 3.774489768532936551500999699815873422073E-3L,
144 4.723962172984642566142399678920790598426E-2L,
145 2.359289678988743939925017240478818248735E-1L,
146 3.032580002220628812728954785118117124520E-1L,
149 static const long double P16_ID
[NP16_ID
+ 1] = {
150 4.389268795186898018132945193912677177553E-15L,
151 4.132671824807454334388868363256830961655E-12L,
152 1.482133328179508835835963635130894413136E-9L,
153 2.618941412861122118906353737117067376236E-7L,
154 2.467854246740858470815714426201888034270E-5L,
155 1.257192927368839847825938545925340230490E-3L,
156 3.362739031941574274949719324644120720341E-2L,
157 4.384458231338934105875343439265370178858E-1L,
158 2.412830809841095249170909628197264854651E0L
,
159 4.176078204111348059102962617368214856874E0L
,
160 /* 1.000000000000000000000000000000000000000E0 */
163 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
164 0.0625 <= 1/x <= 0.125
165 Peak relative error 1.9e-36 */
167 static const long double P8_16N
[NP8_16N
+ 1] = {
168 2.984612480763362345647303274082071598135E-16L,
169 1.923651877544126103941232173085475682334E-13L,
170 4.881258879388869396043760693256024307743E-11L,
171 6.368866572475045408480898921866869811889E-9L,
172 4.684818344104910450523906967821090796737E-7L,
173 2.005177298271593587095982211091300382796E-5L,
174 4.979808067163957634120681477207147536182E-4L,
175 6.946005761642579085284689047091173581127E-3L,
176 5.074601112955765012750207555985299026204E-2L,
177 1.698599455896180893191766195194231825379E-1L,
178 1.957536905259237627737222775573623779638E-1L,
179 2.991314703282528370270179989044994319374E-2L,
182 static const long double P8_16D
[NP8_16D
+ 1] = {
183 2.546869316918069202079580939942463010937E-15L,
184 1.644650111942455804019788382157745229955E-12L,
185 4.185430770291694079925607420808011147173E-10L,
186 5.485331966975218025368698195861074143153E-8L,
187 4.062884421686912042335466327098932678905E-6L,
188 1.758139661060905948870523641319556816772E-4L,
189 4.445143889306356207566032244985607493096E-3L,
190 6.391901016293512632765621532571159071158E-2L,
191 4.933040207519900471177016015718145795434E-1L,
192 1.839144086168947712971630337250761842976E0L
,
193 2.715120873995490920415616716916149586579E0L
,
194 /* 1.000000000000000000000000000000000000000E0 */
197 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
198 0.125 <= 1/x <= 0.1875
199 Peak relative error 1.3e-36 */
201 static const long double P5_8N
[NP5_8N
+ 1] = {
202 2.837678373978003452653763806968237227234E-12L,
203 9.726641165590364928442128579282742354806E-10L,
204 1.284408003604131382028112171490633956539E-7L,
205 8.524624695868291291250573339272194285008E-6L,
206 3.111516908953172249853673787748841282846E-4L,
207 6.423175156126364104172801983096596409176E-3L,
208 7.430220589989104581004416356260692450652E-2L,
209 4.608315409833682489016656279567605536619E-1L,
210 1.396870223510964882676225042258855977512E0L
,
211 1.718500293904122365894630460672081526236E0L
,
212 5.465927698800862172307352821870223855365E-1L
215 static const long double P5_8D
[NP5_8D
+ 1] = {
216 2.421485545794616609951168511612060482715E-11L,
217 8.329862750896452929030058039752327232310E-9L,
218 1.106137992233383429630592081375289010720E-6L,
219 7.405786153760681090127497796448503306939E-5L,
220 2.740364785433195322492093333127633465227E-3L,
221 5.781246470403095224872243564165254652198E-2L,
222 6.927711353039742469918754111511109983546E-1L,
223 4.558679283460430281188304515922826156690E0L
,
224 1.534468499844879487013168065728837900009E1L
,
225 2.313927430889218597919624843161569422745E1L
,
226 1.194506341319498844336768473218382828637E1L
,
227 /* 1.000000000000000000000000000000000000000E0 */
230 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
231 Peak relative error 1.4e-36
232 0.1875 <= 1/x <= 0.25 */
234 static const long double P4_5N
[NP4_5N
+ 1] = {
235 1.846029078268368685834261260420933914621E-10L,
236 3.916295939611376119377869680335444207768E-8L,
237 3.122158792018920627984597530935323997312E-6L,
238 1.218073444893078303994045653603392272450E-4L,
239 2.536420827983485448140477159977981844883E-3L,
240 2.883011322006690823959367922241169171315E-2L,
241 1.755255190734902907438042414495469810830E-1L,
242 5.379317079922628599870898285488723736599E-1L,
243 7.284904050194300773890303361501726561938E-1L,
244 3.270110346613085348094396323925000362813E-1L,
245 1.804473805689725610052078464951722064757E-2L,
248 static const long double P4_5D
[NP4_5D
+ 1] = {
249 1.575278146806816970152174364308980863569E-9L,
250 3.361289173657099516191331123405675054321E-7L,
251 2.704692281550877810424745289838790693708E-5L,
252 1.070854930483999749316546199273521063543E-3L,
253 2.282373093495295842598097265627962125411E-2L,
254 2.692025460665354148328762368240343249830E-1L,
255 1.739892942593664447220951225734811133759E0L
,
256 5.890727576752230385342377570386657229324E0L
,
257 9.517442287057841500750256954117735128153E0L
,
258 6.100616353935338240775363403030137736013E0L
,
259 /* 1.000000000000000000000000000000000000000E0 */
262 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
263 Peak relative error 3.0e-36
264 0.25 <= 1/x <= 0.3125 */
266 static const long double P3r2_4N
[NP3r2_4N
+ 1] = {
267 8.240803130988044478595580300846665863782E-8L,
268 1.179418958381961224222969866406483744580E-5L,
269 6.179787320956386624336959112503824397755E-4L,
270 1.540270833608687596420595830747166658383E-2L,
271 1.983904219491512618376375619598837355076E-1L,
272 1.341465722692038870390470651608301155565E0L
,
273 4.617865326696612898792238245990854646057E0L
,
274 7.435574801812346424460233180412308000587E0L
,
275 4.671327027414635292514599201278557680420E0L
,
276 7.299530852495776936690976966995187714739E-1L,
279 static const long double P3r2_4D
[NP3r2_4D
+ 1] = {
280 7.032152009675729604487575753279187576521E-7L,
281 1.015090352324577615777511269928856742848E-4L,
282 5.394262184808448484302067955186308730620E-3L,
283 1.375291438480256110455809354836988584325E-1L,
284 1.836247144461106304788160919310404376670E0L
,
285 1.314378564254376655001094503090935880349E1L
,
286 4.957184590465712006934452500894672343488E1L
,
287 9.287394244300647738855415178790263465398E1L
,
288 7.652563275535900609085229286020552768399E1L
,
289 2.147042473003074533150718117770093209096E1L
,
290 /* 1.000000000000000000000000000000000000000E0 */
293 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
294 Peak relative error 1.0e-35
295 0.3125 <= 1/x <= 0.375 */
297 static const long double P2r7_3r2N
[NP2r7_3r2N
+ 1] = {
298 4.599033469240421554219816935160627085991E-7L,
299 4.665724440345003914596647144630893997284E-5L,
300 1.684348845667764271596142716944374892756E-3L,
301 2.802446446884455707845985913454440176223E-2L,
302 2.321937586453963310008279956042545173930E-1L,
303 9.640277413988055668692438709376437553804E-1L,
304 1.911021064710270904508663334033003246028E0L
,
305 1.600811610164341450262992138893970224971E0L
,
306 4.266299218652587901171386591543457861138E-1L,
307 1.316470424456061252962568223251247207325E-2L,
310 static const long double P2r7_3r2D
[NP2r7_3r2D
+ 1] = {
311 3.924508608545520758883457108453520099610E-6L,
312 4.029707889408829273226495756222078039823E-4L,
313 1.484629715787703260797886463307469600219E-2L,
314 2.553136379967180865331706538897231588685E-1L,
315 2.229457223891676394409880026887106228740E0L
,
316 1.005708903856384091956550845198392117318E1L
,
317 2.277082659664386953166629360352385889558E1L
,
318 2.384726835193630788249826630376533988245E1L
,
319 9.700989749041320895890113781610939632410E0L
,
320 /* 1.000000000000000000000000000000000000000E0 */
323 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
324 Peak relative error 1.7e-36
325 0.3125 <= 1/x <= 0.4375 */
327 static const long double P2r3_2r7N
[NP2r3_2r7N
+ 1] = {
328 3.916766777108274628543759603786857387402E-6L,
329 3.212176636756546217390661984304645137013E-4L,
330 9.255768488524816445220126081207248947118E-3L,
331 1.214853146369078277453080641911700735354E-1L,
332 7.855163309847214136198449861311404633665E-1L,
333 2.520058073282978403655488662066019816540E0L
,
334 3.825136484837545257209234285382183711466E0L
,
335 2.432569427554248006229715163865569506873E0L
,
336 4.877934835018231178495030117729800489743E-1L,
337 1.109902737860249670981355149101343427885E-2L,
340 static const long double P2r3_2r7D
[NP2r3_2r7D
+ 1] = {
341 3.342307880794065640312646341190547184461E-5L,
342 2.782182891138893201544978009012096558265E-3L,
343 8.221304931614200702142049236141249929207E-2L,
344 1.123728246291165812392918571987858010949E0L
,
345 7.740482453652715577233858317133423434590E0L
,
346 2.737624677567945952953322566311201919139E1L
,
347 4.837181477096062403118304137851260715475E1L
,
348 3.941098643468580791437772701093795299274E1L
,
349 1.245821247166544627558323920382547533630E1L
,
350 /* 1.000000000000000000000000000000000000000E0 */
353 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
354 Peak relative error 1.7e-35
355 0.4375 <= 1/x <= 0.5 */
357 static const long double P2_2r3N
[NP2_2r3N
+ 1] = {
358 3.397930802851248553545191160608731940751E-4L,
359 2.104020902735482418784312825637833698217E-2L,
360 4.442291771608095963935342749477836181939E-1L,
361 4.131797328716583282869183304291833754967E0L
,
362 1.819920169779026500146134832455189917589E1L
,
363 3.781779616522937565300309684282401791291E1L
,
364 3.459605449728864218972931220783543410347E1L
,
365 1.173594248397603882049066603238568316561E1L
,
366 9.455702270242780642835086549285560316461E-1L,
369 static const long double P2_2r3D
[NP2_2r3D
+ 1] = {
370 2.899568897241432883079888249845707400614E-3L,
371 1.831107138190848460767699919531132426356E-1L,
372 3.999350044057883839080258832758908825165E0L
,
373 3.929041535867957938340569419874195303712E1L
,
374 1.884245613422523323068802689915538908291E2L
,
375 4.461469948819229734353852978424629815929E2L
,
376 5.004998753999796821224085972610636347903E2L
,
377 2.386342520092608513170837883757163414100E2L
,
378 3.791322528149347975999851588922424189957E1L
,
379 /* 1.000000000000000000000000000000000000000E0 */
382 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
383 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
384 Peak relative error 8.0e-36
387 static const long double Q16_IN
[NQ16_IN
+ 1] = {
388 -3.917420835712508001321875734030357393421E-18L,
389 -4.440311387483014485304387406538069930457E-15L,
390 -1.951635424076926487780929645954007139616E-12L,
391 -4.318256438421012555040546775651612810513E-10L,
392 -5.231244131926180765270446557146989238020E-8L,
393 -3.540072702902043752460711989234732357653E-6L,
394 -1.311017536555269966928228052917534882984E-4L,
395 -2.495184669674631806622008769674827575088E-3L,
396 -2.141868222987209028118086708697998506716E-2L,
397 -6.184031415202148901863605871197272650090E-2L,
398 -1.922298704033332356899546792898156493887E-2L,
401 static const long double Q16_ID
[NQ16_ID
+ 1] = {
402 3.820418034066293517479619763498400162314E-17L,
403 4.340702810799239909648911373329149354911E-14L,
404 1.914985356383416140706179933075303538524E-11L,
405 4.262333682610888819476498617261895474330E-9L,
406 5.213481314722233980346462747902942182792E-7L,
407 3.585741697694069399299005316809954590558E-5L,
408 1.366513429642842006385029778105539457546E-3L,
409 2.745282599850704662726337474371355160594E-2L,
410 2.637644521611867647651200098449903330074E-1L,
411 1.006953426110765984590782655598680488746E0L
,
412 /* 1.000000000000000000000000000000000000000E0 */
415 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
416 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
417 Peak relative error 1.9e-36
418 0.0625 <= 1/x <= 0.125 */
420 static const long double Q8_16N
[NQ8_16N
+ 1] = {
421 -2.028630366670228670781362543615221542291E-17L,
422 -1.519634620380959966438130374006858864624E-14L,
423 -4.540596528116104986388796594639405114524E-12L,
424 -7.085151756671466559280490913558388648274E-10L,
425 -6.351062671323970823761883833531546885452E-8L,
426 -3.390817171111032905297982523519503522491E-6L,
427 -1.082340897018886970282138836861233213972E-4L,
428 -2.020120801187226444822977006648252379508E-3L,
429 -2.093169910981725694937457070649605557555E-2L,
430 -1.092176538874275712359269481414448063393E-1L,
431 -2.374790947854765809203590474789108718733E-1L,
432 -1.365364204556573800719985118029601401323E-1L,
435 static const long double Q8_16D
[NQ8_16D
+ 1] = {
436 1.978397614733632533581207058069628242280E-16L,
437 1.487361156806202736877009608336766720560E-13L,
438 4.468041406888412086042576067133365913456E-11L,
439 7.027822074821007443672290507210594648877E-9L,
440 6.375740580686101224127290062867976007374E-7L,
441 3.466887658320002225888644977076410421940E-5L,
442 1.138625640905289601186353909213719596986E-3L,
443 2.224470799470414663443449818235008486439E-2L,
444 2.487052928527244907490589787691478482358E-1L,
445 1.483927406564349124649083853892380899217E0L
,
446 4.182773513276056975777258788903489507705E0L
,
447 4.419665392573449746043880892524360870944E0L
,
448 /* 1.000000000000000000000000000000000000000E0 */
451 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
452 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
453 Peak relative error 1.5e-35
454 0.125 <= 1/x <= 0.1875 */
456 static const long double Q5_8N
[NQ5_8N
+ 1] = {
457 -3.656082407740970534915918390488336879763E-13L,
458 -1.344660308497244804752334556734121771023E-10L,
459 -1.909765035234071738548629788698150760791E-8L,
460 -1.366668038160120210269389551283666716453E-6L,
461 -5.392327355984269366895210704976314135683E-5L,
462 -1.206268245713024564674432357634540343884E-3L,
463 -1.515456784370354374066417703736088291287E-2L,
464 -1.022454301137286306933217746545237098518E-1L,
465 -3.373438906472495080504907858424251082240E-1L,
466 -4.510782522110845697262323973549178453405E-1L,
467 -1.549000892545288676809660828213589804884E-1L,
470 static const long double Q5_8D
[NQ5_8D
+ 1] = {
471 3.565550843359501079050699598913828460036E-12L,
472 1.321016015556560621591847454285330528045E-9L,
473 1.897542728662346479999969679234270605975E-7L,
474 1.381720283068706710298734234287456219474E-5L,
475 5.599248147286524662305325795203422873725E-4L,
476 1.305442352653121436697064782499122164843E-2L,
477 1.750234079626943298160445750078631894985E-1L,
478 1.311420542073436520965439883806946678491E0L
,
479 5.162757689856842406744504211089724926650E0L
,
480 9.527760296384704425618556332087850581308E0L
,
481 6.604648207463236667912921642545100248584E0L
,
482 /* 1.000000000000000000000000000000000000000E0 */
485 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
486 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
487 Peak relative error 1.3e-35
488 0.1875 <= 1/x <= 0.25 */
490 static const long double Q4_5N
[NQ4_5N
+ 1] = {
491 -4.079513568708891749424783046520200903755E-11L,
492 -9.326548104106791766891812583019664893311E-9L,
493 -8.016795121318423066292906123815687003356E-7L,
494 -3.372350544043594415609295225664186750995E-5L,
495 -7.566238665947967882207277686375417983917E-4L,
496 -9.248861580055565402130441618521591282617E-3L,
497 -6.033106131055851432267702948850231270338E-2L,
498 -1.966908754799996793730369265431584303447E-1L,
499 -2.791062741179964150755788226623462207560E-1L,
500 -1.255478605849190549914610121863534191666E-1L,
501 -4.320429862021265463213168186061696944062E-3L,
504 static const long double Q4_5D
[NQ4_5D
+ 1] = {
505 3.978497042580921479003851216297330701056E-10L,
506 9.203304163828145809278568906420772246666E-8L,
507 8.059685467088175644915010485174545743798E-6L,
508 3.490187375993956409171098277561669167446E-4L,
509 8.189109654456872150100501732073810028829E-3L,
510 1.072572867311023640958725265762483033769E-1L,
511 7.790606862409960053675717185714576937994E-1L,
512 3.016049768232011196434185423512777656328E0L
,
513 5.722963851442769787733717162314477949360E0L
,
514 4.510527838428473279647251350931380867663E0L
,
515 /* 1.000000000000000000000000000000000000000E0 */
518 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
519 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
520 Peak relative error 2.1e-35
521 0.25 <= 1/x <= 0.3125 */
523 static const long double Q3r2_4N
[NQ3r2_4N
+ 1] = {
524 -1.087480809271383885936921889040388133627E-8L,
525 -1.690067828697463740906962973479310170932E-6L,
526 -9.608064416995105532790745641974762550982E-5L,
527 -2.594198839156517191858208513873961837410E-3L,
528 -3.610954144421543968160459863048062977822E-2L,
529 -2.629866798251843212210482269563961685666E-1L,
530 -9.709186825881775885917984975685752956660E-1L,
531 -1.667521829918185121727268867619982417317E0L
,
532 -1.109255082925540057138766105229900943501E0L
,
533 -1.812932453006641348145049323713469043328E-1L,
536 static const long double Q3r2_4D
[NQ3r2_4D
+ 1] = {
537 1.060552717496912381388763753841473407026E-7L,
538 1.676928002024920520786883649102388708024E-5L,
539 9.803481712245420839301400601140812255737E-4L,
540 2.765559874262309494758505158089249012930E-2L,
541 4.117921827792571791298862613287549140706E-1L,
542 3.323769515244751267093378361930279161413E0L
,
543 1.436602494405814164724810151689705353670E1L
,
544 3.163087869617098638064881410646782408297E1L
,
545 3.198181264977021649489103980298349589419E1L
,
546 1.203649258862068431199471076202897823272E1L
,
547 /* 1.000000000000000000000000000000000000000E0 */
550 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
551 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
552 Peak relative error 1.6e-36
553 0.3125 <= 1/x <= 0.375 */
555 static const long double Q2r7_3r2N
[NQ2r7_3r2N
+ 1] = {
556 -1.723405393982209853244278760171643219530E-7L,
557 -2.090508758514655456365709712333460087442E-5L,
558 -9.140104013370974823232873472192719263019E-4L,
559 -1.871349499990714843332742160292474780128E-2L,
560 -1.948930738119938669637865956162512983416E-1L,
561 -1.048764684978978127908439526343174139788E0L
,
562 -2.827714929925679500237476105843643064698E0L
,
563 -3.508761569156476114276988181329773987314E0L
,
564 -1.669332202790211090973255098624488308989E0L
,
565 -1.930796319299022954013840684651016077770E-1L,
568 static const long double Q2r7_3r2D
[NQ2r7_3r2D
+ 1] = {
569 1.680730662300831976234547482334347983474E-6L,
570 2.084241442440551016475972218719621841120E-4L,
571 9.445316642108367479043541702688736295579E-3L,
572 2.044637889456631896650179477133252184672E-1L,
573 2.316091982244297350829522534435350078205E0L
,
574 1.412031891783015085196708811890448488865E1L
,
575 4.583830154673223384837091077279595496149E1L
,
576 7.549520609270909439885998474045974122261E1L
,
577 5.697605832808113367197494052388203310638E1L
,
578 1.601496240876192444526383314589371686234E1L
,
579 /* 1.000000000000000000000000000000000000000E0 */
582 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
583 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
584 Peak relative error 9.5e-36
585 0.375 <= 1/x <= 0.4375 */
587 static const long double Q2r3_2r7N
[NQ2r3_2r7N
+ 1] = {
588 -8.603042076329122085722385914954878953775E-7L,
589 -7.701746260451647874214968882605186675720E-5L,
590 -2.407932004380727587382493696877569654271E-3L,
591 -3.403434217607634279028110636919987224188E-2L,
592 -2.348707332185238159192422084985713102877E-1L,
593 -7.957498841538254916147095255700637463207E-1L,
594 -1.258469078442635106431098063707934348577E0L
,
595 -8.162415474676345812459353639449971369890E-1L,
596 -1.581783890269379690141513949609572806898E-1L,
597 -1.890595651683552228232308756569450822905E-3L,
600 static const long double Q2r3_2r7D
[NQ2r3_2r7D
+ 1] = {
601 8.390017524798316921170710533381568175665E-6L,
602 7.738148683730826286477254659973968763659E-4L,
603 2.541480810958665794368759558791634341779E-2L,
604 3.878879789711276799058486068562386244873E-1L,
605 3.003783779325811292142957336802456109333E0L
,
606 1.206480374773322029883039064575464497400E1L
,
607 2.458414064785315978408974662900438351782E1L
,
608 2.367237826273668567199042088835448715228E1L
,
609 9.231451197519171090875569102116321676763E0L
,
610 /* 1.000000000000000000000000000000000000000E0 */
613 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
614 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
615 Peak relative error 1.4e-36
616 0.4375 <= 1/x <= 0.5 */
618 static const long double Q2_2r3N
[NQ2_2r3N
+ 1] = {
619 -5.552507516089087822166822364590806076174E-6L,
620 -4.135067659799500521040944087433752970297E-4L,
621 -1.059928728869218962607068840646564457980E-2L,
622 -1.212070036005832342565792241385459023801E-1L,
623 -6.688350110633603958684302153362735625156E-1L,
624 -1.793587878197360221340277951304429821582E0L
,
625 -2.225407682237197485644647380483725045326E0L
,
626 -1.123402135458940189438898496348239744403E0L
,
627 -1.679187241566347077204805190763597299805E-1L,
628 -1.458550613639093752909985189067233504148E-3L,
631 static const long double Q2_2r3D
[NQ2_2r3D
+ 1] = {
632 5.415024336507980465169023996403597916115E-5L,
633 4.179246497380453022046357404266022870788E-3L,
634 1.136306384261959483095442402929502368598E-1L,
635 1.422640343719842213484515445393284072830E0L
,
636 8.968786703393158374728850922289204805764E0L
,
637 2.914542473339246127533384118781216495934E1L
,
638 4.781605421020380669870197378210457054685E1L
,
639 3.693865837171883152382820584714795072937E1L
,
640 1.153220502744204904763115556224395893076E1L
,
641 /* 1.000000000000000000000000000000000000000E0 */
645 /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
648 neval (long double x
, const long double *p
, int n
)
663 /* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */
666 deval (long double x
, const long double *p
, int n
)
681 /* Bessel function of the first kind, order one. */
684 __ieee754_j1l (long double x
)
686 long double xx
, xinv
, z
, p
, q
, c
, s
, cc
, ss
;
702 p
= xx
* z
* neval (z
, J0_2N
, NJ0_2N
) / deval (z
, J0_2D
, NJ0_2D
);
717 p
= neval (z
, P16_IN
, NP16_IN
) / deval (z
, P16_ID
, NP16_ID
);
718 q
= neval (z
, Q16_IN
, NQ16_IN
) / deval (z
, Q16_ID
, NQ16_ID
);
722 p
= neval (z
, P8_16N
, NP8_16N
) / deval (z
, P8_16D
, NP8_16D
);
723 q
= neval (z
, Q8_16N
, NQ8_16N
) / deval (z
, Q8_16D
, NQ8_16D
);
726 else if (xinv
<= 0.1875)
728 p
= neval (z
, P5_8N
, NP5_8N
) / deval (z
, P5_8D
, NP5_8D
);
729 q
= neval (z
, Q5_8N
, NQ5_8N
) / deval (z
, Q5_8D
, NQ5_8D
);
733 p
= neval (z
, P4_5N
, NP4_5N
) / deval (z
, P4_5D
, NP4_5D
);
734 q
= neval (z
, Q4_5N
, NQ4_5N
) / deval (z
, Q4_5D
, NQ4_5D
);
737 else /* if (xinv <= 0.5) */
743 p
= neval (z
, P3r2_4N
, NP3r2_4N
) / deval (z
, P3r2_4D
, NP3r2_4D
);
744 q
= neval (z
, Q3r2_4N
, NQ3r2_4N
) / deval (z
, Q3r2_4D
, NQ3r2_4D
);
748 p
= neval (z
, P2r7_3r2N
, NP2r7_3r2N
)
749 / deval (z
, P2r7_3r2D
, NP2r7_3r2D
);
750 q
= neval (z
, Q2r7_3r2N
, NQ2r7_3r2N
)
751 / deval (z
, Q2r7_3r2D
, NQ2r7_3r2D
);
754 else if (xinv
<= 0.4375)
756 p
= neval (z
, P2r3_2r7N
, NP2r3_2r7N
)
757 / deval (z
, P2r3_2r7D
, NP2r3_2r7D
);
758 q
= neval (z
, Q2r3_2r7N
, NQ2r3_2r7N
)
759 / deval (z
, Q2r3_2r7D
, NQ2r3_2r7D
);
763 p
= neval (z
, P2_2r3N
, NP2_2r3N
) / deval (z
, P2_2r3D
, NP2_2r3D
);
764 q
= neval (z
, Q2_2r3N
, NQ2_2r3N
) / deval (z
, Q2_2r3D
, NQ2_2r3D
);
769 q
= q
* xinv
+ 0.375L * xinv
;
771 cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
772 = 1/sqrt(2) * (-cos(x) + sin(x))
773 sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
774 = -1/sqrt(2) * (sin(x) + cos(x))
776 __sincosl (xx
, &s
, &c
);
779 z
= __cosl (xx
+ xx
);
784 z
= ONEOSQPI
* (p
* cc
- q
* ss
) / __ieee754_sqrtl (xx
);
791 /* Y1(x) = 2/pi * (log(x) * J1(x) - 1/x) + x R(x^2)
792 Peak relative error 6.2e-38
795 static long double Y0_2N
[NY0_2N
+ 1] = {
796 -6.804415404830253804408698161694720833249E19L
,
797 1.805450517967019908027153056150465849237E19L
,
798 -8.065747497063694098810419456383006737312E17L
,
799 1.401336667383028259295830955439028236299E16L
,
800 -1.171654432898137585000399489686629680230E14L
,
801 5.061267920943853732895341125243428129150E11L
,
802 -1.096677850566094204586208610960870217970E9L
,
803 9.541172044989995856117187515882879304461E5L
,
806 static long double Y0_2D
[NY0_2D
+ 1] = {
807 3.470629591820267059538637461549677594549E20L
,
808 4.120796439009916326855848107545425217219E18L
,
809 2.477653371652018249749350657387030814542E16L
,
810 9.954678543353888958177169349272167762797E13L
,
811 2.957927997613630118216218290262851197754E11L
,
812 6.748421382188864486018861197614025972118E8L
,
813 1.173453425218010888004562071020305709319E6L
,
814 1.450335662961034949894009554536003377187E3L
,
815 /* 1.000000000000000000000000000000000000000E0 */
819 /* Bessel function of the second kind, order one. */
822 __ieee754_y1l (long double x
)
824 long double xx
, xinv
, z
, p
, q
, c
, s
, cc
, ss
;
836 return (zero
/ (zero
* x
));
837 return -HUGE_VALL
+ x
;
844 p
= xx
* neval (z
, Y0_2N
, NY0_2N
) / deval (z
, Y0_2D
, NY0_2D
);
845 p
= -TWOOPI
/ xx
+ p
;
846 p
= TWOOPI
* __ieee754_logl (x
) * __ieee754_j1l (x
) + p
;
858 p
= neval (z
, P16_IN
, NP16_IN
) / deval (z
, P16_ID
, NP16_ID
);
859 q
= neval (z
, Q16_IN
, NQ16_IN
) / deval (z
, Q16_ID
, NQ16_ID
);
863 p
= neval (z
, P8_16N
, NP8_16N
) / deval (z
, P8_16D
, NP8_16D
);
864 q
= neval (z
, Q8_16N
, NQ8_16N
) / deval (z
, Q8_16D
, NQ8_16D
);
867 else if (xinv
<= 0.1875)
869 p
= neval (z
, P5_8N
, NP5_8N
) / deval (z
, P5_8D
, NP5_8D
);
870 q
= neval (z
, Q5_8N
, NQ5_8N
) / deval (z
, Q5_8D
, NQ5_8D
);
874 p
= neval (z
, P4_5N
, NP4_5N
) / deval (z
, P4_5D
, NP4_5D
);
875 q
= neval (z
, Q4_5N
, NQ4_5N
) / deval (z
, Q4_5D
, NQ4_5D
);
878 else /* if (xinv <= 0.5) */
884 p
= neval (z
, P3r2_4N
, NP3r2_4N
) / deval (z
, P3r2_4D
, NP3r2_4D
);
885 q
= neval (z
, Q3r2_4N
, NQ3r2_4N
) / deval (z
, Q3r2_4D
, NQ3r2_4D
);
889 p
= neval (z
, P2r7_3r2N
, NP2r7_3r2N
)
890 / deval (z
, P2r7_3r2D
, NP2r7_3r2D
);
891 q
= neval (z
, Q2r7_3r2N
, NQ2r7_3r2N
)
892 / deval (z
, Q2r7_3r2D
, NQ2r7_3r2D
);
895 else if (xinv
<= 0.4375)
897 p
= neval (z
, P2r3_2r7N
, NP2r3_2r7N
)
898 / deval (z
, P2r3_2r7D
, NP2r3_2r7D
);
899 q
= neval (z
, Q2r3_2r7N
, NQ2r3_2r7N
)
900 / deval (z
, Q2r3_2r7D
, NQ2r3_2r7D
);
904 p
= neval (z
, P2_2r3N
, NP2_2r3N
) / deval (z
, P2_2r3D
, NP2_2r3D
);
905 q
= neval (z
, Q2_2r3N
, NQ2_2r3N
) / deval (z
, Q2_2r3D
, NQ2_2r3D
);
910 q
= q
* xinv
+ 0.375L * xinv
;
912 cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
913 = 1/sqrt(2) * (-cos(x) + sin(x))
914 sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
915 = -1/sqrt(2) * (sin(x) + cos(x))
917 __sincosl (xx
, &s
, &c
);
920 z
= __cosl (xx
+ xx
);
925 z
= ONEOSQPI
* (p
* ss
+ q
* cc
) / __ieee754_sqrtl (xx
);