1/* j1l.c
2 *
3 * Bessel function of order one
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * __float128 x, y, j1q();
10 *
11 * y = j1q( x );
12 *
13 *
14 *
15 * DESCRIPTION:
16 *
17 * Returns Bessel function of first kind, order one 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 is
21 * J1(x) = .5x + x x^2 R(x^2)
22 *
23 * The second interval is further partitioned into eight equal segments
24 * of 1/x.
25 * J1(x) = sqrt(2/(pi x)) (P1(x) cos(X) - Q1(x) sin(X)),
26 * X = x - 3 pi / 4,
27 *
28 * and the auxiliary functions are given by
29 *
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)
32 *
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)).
35 *
36 *
37 *
38 * ACCURACY:
39 *
40 * Absolute error:
41 * arithmetic domain # trials peak rms
42 * IEEE 0, 30 100000 2.8e-34 2.7e-35
43 *
44 *
45 */
46
47/* y1l.c
48 *
49 * Bessel function of the second kind, order one
50 *
51 *
52 *
53 * SYNOPSIS:
54 *
55 * __float128, y, y1q();
56 *
57 * y = y1q( x );
58 *
59 *
60 *
61 * DESCRIPTION:
62 *
63 * Returns Bessel function of the second kind, of order
64 * one, of the argument.
65 *
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)),
71 * X = x - 3 pi / 4.
72 *
73 * ACCURACY:
74 *
75 * Absolute error, when y0(x) < 1; else relative error:
76 *
77 * arithmetic domain # trials peak rms
78 * IEEE 0, 30 100000 2.7e-34 2.9e-35
79 *
80 */
81
82/* Copyright 2001 by Stephen L. Moshier (moshier@na-net.onrl.gov).
83
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.
88
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.
93
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 */
97
98#include <errno.h>
99#include "quadmath-imp.h"
100
101/* 1 / sqrt(pi) */
102static const __float128 ONEOSQPI = 5.6418958354775628694807945156077258584405E-1Q;
103/* 2 / pi */
104static const __float128 TWOOPI = 6.3661977236758134307553505349005744813784E-1Q;
105static const __float128 zero = 0.0Q;
106
107/* J1(x) = .5x + x x^2 R(x^2)
108 Peak relative error 1.9e-35
109 0 <= x <= 2 */
110#define NJ0_2N 6
111static const __float128 J0_2N[NJ0_2N + 1] = {
112 -5.943799577386942855938508697619735179660E16Q,
113 1.812087021305009192259946997014044074711E15Q,
114 -2.761698314264509665075127515729146460895E13Q,
115 2.091089497823600978949389109350658815972E11Q,
116 -8.546413231387036372945453565654130054307E8Q,
117 1.797229225249742247475464052741320612261E6Q,
118 -1.559552840946694171346552770008812083969E3Q
119};
120#define NJ0_2D 6
121static const __float128 J0_2D[NJ0_2D + 1] = {
122 9.510079323819108569501613916191477479397E17Q,
123 1.063193817503280529676423936545854693915E16Q,
124 5.934143516050192600795972192791775226920E13Q,
125 2.168000911950620999091479265214368352883E11Q,
126 5.673775894803172808323058205986256928794E8Q,
127 1.080329960080981204840966206372671147224E6Q,
128 1.411951256636576283942477881535283304912E3Q,
129 /* 1.000000000000000000000000000000000000000E0Q */
130};
131
132/* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
133 0 <= 1/x <= .0625
134 Peak relative error 3.6e-36 */
135#define NP16_IN 9
136static const __float128 P16_IN[NP16_IN + 1] = {
137 5.143674369359646114999545149085139822905E-16Q,
138 4.836645664124562546056389268546233577376E-13Q,
139 1.730945562285804805325011561498453013673E-10Q,
140 3.047976856147077889834905908605310585810E-8Q,
141 2.855227609107969710407464739188141162386E-6Q,
142 1.439362407936705484122143713643023998457E-4Q,
143 3.774489768532936551500999699815873422073E-3Q,
144 4.723962172984642566142399678920790598426E-2Q,
145 2.359289678988743939925017240478818248735E-1Q,
146 3.032580002220628812728954785118117124520E-1Q,
147};
148#define NP16_ID 9
149static const __float128 P16_ID[NP16_ID + 1] = {
150 4.389268795186898018132945193912677177553E-15Q,
151 4.132671824807454334388868363256830961655E-12Q,
152 1.482133328179508835835963635130894413136E-9Q,
153 2.618941412861122118906353737117067376236E-7Q,
154 2.467854246740858470815714426201888034270E-5Q,
155 1.257192927368839847825938545925340230490E-3Q,
156 3.362739031941574274949719324644120720341E-2Q,
157 4.384458231338934105875343439265370178858E-1Q,
158 2.412830809841095249170909628197264854651E0Q,
159 4.176078204111348059102962617368214856874E0Q,
160 /* 1.000000000000000000000000000000000000000E0 */
161};
162
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 */
166#define NP8_16N 11
167static const __float128 P8_16N[NP8_16N + 1] = {
168 2.984612480763362345647303274082071598135E-16Q,
169 1.923651877544126103941232173085475682334E-13Q,
170 4.881258879388869396043760693256024307743E-11Q,
171 6.368866572475045408480898921866869811889E-9Q,
172 4.684818344104910450523906967821090796737E-7Q,
173 2.005177298271593587095982211091300382796E-5Q,
174 4.979808067163957634120681477207147536182E-4Q,
175 6.946005761642579085284689047091173581127E-3Q,
176 5.074601112955765012750207555985299026204E-2Q,
177 1.698599455896180893191766195194231825379E-1Q,
178 1.957536905259237627737222775573623779638E-1Q,
179 2.991314703282528370270179989044994319374E-2Q,
180};
181#define NP8_16D 10
182static const __float128 P8_16D[NP8_16D + 1] = {
183 2.546869316918069202079580939942463010937E-15Q,
184 1.644650111942455804019788382157745229955E-12Q,
185 4.185430770291694079925607420808011147173E-10Q,
186 5.485331966975218025368698195861074143153E-8Q,
187 4.062884421686912042335466327098932678905E-6Q,
188 1.758139661060905948870523641319556816772E-4Q,
189 4.445143889306356207566032244985607493096E-3Q,
190 6.391901016293512632765621532571159071158E-2Q,
191 4.933040207519900471177016015718145795434E-1Q,
192 1.839144086168947712971630337250761842976E0Q,
193 2.715120873995490920415616716916149586579E0Q,
194 /* 1.000000000000000000000000000000000000000E0 */
195};
196
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 */
200#define NP5_8N 10
201static const __float128 P5_8N[NP5_8N + 1] = {
202 2.837678373978003452653763806968237227234E-12Q,
203 9.726641165590364928442128579282742354806E-10Q,
204 1.284408003604131382028112171490633956539E-7Q,
205 8.524624695868291291250573339272194285008E-6Q,
206 3.111516908953172249853673787748841282846E-4Q,
207 6.423175156126364104172801983096596409176E-3Q,
208 7.430220589989104581004416356260692450652E-2Q,
209 4.608315409833682489016656279567605536619E-1Q,
210 1.396870223510964882676225042258855977512E0Q,
211 1.718500293904122365894630460672081526236E0Q,
212 5.465927698800862172307352821870223855365E-1Q
213};
214#define NP5_8D 10
215static const __float128 P5_8D[NP5_8D + 1] = {
216 2.421485545794616609951168511612060482715E-11Q,
217 8.329862750896452929030058039752327232310E-9Q,
218 1.106137992233383429630592081375289010720E-6Q,
219 7.405786153760681090127497796448503306939E-5Q,
220 2.740364785433195322492093333127633465227E-3Q,
221 5.781246470403095224872243564165254652198E-2Q,
222 6.927711353039742469918754111511109983546E-1Q,
223 4.558679283460430281188304515922826156690E0Q,
224 1.534468499844879487013168065728837900009E1Q,
225 2.313927430889218597919624843161569422745E1Q,
226 1.194506341319498844336768473218382828637E1Q,
227 /* 1.000000000000000000000000000000000000000E0 */
228};
229
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 */
233#define NP4_5N 10
234static const __float128 P4_5N[NP4_5N + 1] = {
235 1.846029078268368685834261260420933914621E-10Q,
236 3.916295939611376119377869680335444207768E-8Q,
237 3.122158792018920627984597530935323997312E-6Q,
238 1.218073444893078303994045653603392272450E-4Q,
239 2.536420827983485448140477159977981844883E-3Q,
240 2.883011322006690823959367922241169171315E-2Q,
241 1.755255190734902907438042414495469810830E-1Q,
242 5.379317079922628599870898285488723736599E-1Q,
243 7.284904050194300773890303361501726561938E-1Q,
244 3.270110346613085348094396323925000362813E-1Q,
245 1.804473805689725610052078464951722064757E-2Q,
246};
247#define NP4_5D 9
248static const __float128 P4_5D[NP4_5D + 1] = {
249 1.575278146806816970152174364308980863569E-9Q,
250 3.361289173657099516191331123405675054321E-7Q,
251 2.704692281550877810424745289838790693708E-5Q,
252 1.070854930483999749316546199273521063543E-3Q,
253 2.282373093495295842598097265627962125411E-2Q,
254 2.692025460665354148328762368240343249830E-1Q,
255 1.739892942593664447220951225734811133759E0Q,
256 5.890727576752230385342377570386657229324E0Q,
257 9.517442287057841500750256954117735128153E0Q,
258 6.100616353935338240775363403030137736013E0Q,
259 /* 1.000000000000000000000000000000000000000E0 */
260};
261
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 */
265#define NP3r2_4N 9
266static const __float128 P3r2_4N[NP3r2_4N + 1] = {
267 8.240803130988044478595580300846665863782E-8Q,
268 1.179418958381961224222969866406483744580E-5Q,
269 6.179787320956386624336959112503824397755E-4Q,
270 1.540270833608687596420595830747166658383E-2Q,
271 1.983904219491512618376375619598837355076E-1Q,
272 1.341465722692038870390470651608301155565E0Q,
273 4.617865326696612898792238245990854646057E0Q,
274 7.435574801812346424460233180412308000587E0Q,
275 4.671327027414635292514599201278557680420E0Q,
276 7.299530852495776936690976966995187714739E-1Q,
277};
278#define NP3r2_4D 9
279static const __float128 P3r2_4D[NP3r2_4D + 1] = {
280 7.032152009675729604487575753279187576521E-7Q,
281 1.015090352324577615777511269928856742848E-4Q,
282 5.394262184808448484302067955186308730620E-3Q,
283 1.375291438480256110455809354836988584325E-1Q,
284 1.836247144461106304788160919310404376670E0Q,
285 1.314378564254376655001094503090935880349E1Q,
286 4.957184590465712006934452500894672343488E1Q,
287 9.287394244300647738855415178790263465398E1Q,
288 7.652563275535900609085229286020552768399E1Q,
289 2.147042473003074533150718117770093209096E1Q,
290 /* 1.000000000000000000000000000000000000000E0 */
291};
292
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 */
296#define NP2r7_3r2N 9
297static const __float128 P2r7_3r2N[NP2r7_3r2N + 1] = {
298 4.599033469240421554219816935160627085991E-7Q,
299 4.665724440345003914596647144630893997284E-5Q,
300 1.684348845667764271596142716944374892756E-3Q,
301 2.802446446884455707845985913454440176223E-2Q,
302 2.321937586453963310008279956042545173930E-1Q,
303 9.640277413988055668692438709376437553804E-1Q,
304 1.911021064710270904508663334033003246028E0Q,
305 1.600811610164341450262992138893970224971E0Q,
306 4.266299218652587901171386591543457861138E-1Q,
307 1.316470424456061252962568223251247207325E-2Q,
308};
309#define NP2r7_3r2D 8
310static const __float128 P2r7_3r2D[NP2r7_3r2D + 1] = {
311 3.924508608545520758883457108453520099610E-6Q,
312 4.029707889408829273226495756222078039823E-4Q,
313 1.484629715787703260797886463307469600219E-2Q,
314 2.553136379967180865331706538897231588685E-1Q,
315 2.229457223891676394409880026887106228740E0Q,
316 1.005708903856384091956550845198392117318E1Q,
317 2.277082659664386953166629360352385889558E1Q,
318 2.384726835193630788249826630376533988245E1Q,
319 9.700989749041320895890113781610939632410E0Q,
320 /* 1.000000000000000000000000000000000000000E0 */
321};
322
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 */
326#define NP2r3_2r7N 9
327static const __float128 P2r3_2r7N[NP2r3_2r7N + 1] = {
328 3.916766777108274628543759603786857387402E-6Q,
329 3.212176636756546217390661984304645137013E-4Q,
330 9.255768488524816445220126081207248947118E-3Q,
331 1.214853146369078277453080641911700735354E-1Q,
332 7.855163309847214136198449861311404633665E-1Q,
333 2.520058073282978403655488662066019816540E0Q,
334 3.825136484837545257209234285382183711466E0Q,
335 2.432569427554248006229715163865569506873E0Q,
336 4.877934835018231178495030117729800489743E-1Q,
337 1.109902737860249670981355149101343427885E-2Q,
338};
339#define NP2r3_2r7D 8
340static const __float128 P2r3_2r7D[NP2r3_2r7D + 1] = {
341 3.342307880794065640312646341190547184461E-5Q,
342 2.782182891138893201544978009012096558265E-3Q,
343 8.221304931614200702142049236141249929207E-2Q,
344 1.123728246291165812392918571987858010949E0Q,
345 7.740482453652715577233858317133423434590E0Q,
346 2.737624677567945952953322566311201919139E1Q,
347 4.837181477096062403118304137851260715475E1Q,
348 3.941098643468580791437772701093795299274E1Q,
349 1.245821247166544627558323920382547533630E1Q,
350 /* 1.000000000000000000000000000000000000000E0 */
351};
352
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 */
356#define NP2_2r3N 8
357static const __float128 P2_2r3N[NP2_2r3N + 1] = {
358 3.397930802851248553545191160608731940751E-4Q,
359 2.104020902735482418784312825637833698217E-2Q,
360 4.442291771608095963935342749477836181939E-1Q,
361 4.131797328716583282869183304291833754967E0Q,
362 1.819920169779026500146134832455189917589E1Q,
363 3.781779616522937565300309684282401791291E1Q,
364 3.459605449728864218972931220783543410347E1Q,
365 1.173594248397603882049066603238568316561E1Q,
366 9.455702270242780642835086549285560316461E-1Q,
367};
368#define NP2_2r3D 8
369static const __float128 P2_2r3D[NP2_2r3D + 1] = {
370 2.899568897241432883079888249845707400614E-3Q,
371 1.831107138190848460767699919531132426356E-1Q,
372 3.999350044057883839080258832758908825165E0Q,
373 3.929041535867957938340569419874195303712E1Q,
374 1.884245613422523323068802689915538908291E2Q,
375 4.461469948819229734353852978424629815929E2Q,
376 5.004998753999796821224085972610636347903E2Q,
377 2.386342520092608513170837883757163414100E2Q,
378 3.791322528149347975999851588922424189957E1Q,
379 /* 1.000000000000000000000000000000000000000E0 */
380};
381
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
385 0 <= 1/x <= .0625 */
386#define NQ16_IN 10
387static const __float128 Q16_IN[NQ16_IN + 1] = {
388 -3.917420835712508001321875734030357393421E-18Q,
389 -4.440311387483014485304387406538069930457E-15Q,
390 -1.951635424076926487780929645954007139616E-12Q,
391 -4.318256438421012555040546775651612810513E-10Q,
392 -5.231244131926180765270446557146989238020E-8Q,
393 -3.540072702902043752460711989234732357653E-6Q,
394 -1.311017536555269966928228052917534882984E-4Q,
395 -2.495184669674631806622008769674827575088E-3Q,
396 -2.141868222987209028118086708697998506716E-2Q,
397 -6.184031415202148901863605871197272650090E-2Q,
398 -1.922298704033332356899546792898156493887E-2Q,
399};
400#define NQ16_ID 9
401static const __float128 Q16_ID[NQ16_ID + 1] = {
402 3.820418034066293517479619763498400162314E-17Q,
403 4.340702810799239909648911373329149354911E-14Q,
404 1.914985356383416140706179933075303538524E-11Q,
405 4.262333682610888819476498617261895474330E-9Q,
406 5.213481314722233980346462747902942182792E-7Q,
407 3.585741697694069399299005316809954590558E-5Q,
408 1.366513429642842006385029778105539457546E-3Q,
409 2.745282599850704662726337474371355160594E-2Q,
410 2.637644521611867647651200098449903330074E-1Q,
411 1.006953426110765984590782655598680488746E0Q,
412 /* 1.000000000000000000000000000000000000000E0 */
413 };
414
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 */
419#define NQ8_16N 11
420static const __float128 Q8_16N[NQ8_16N + 1] = {
421 -2.028630366670228670781362543615221542291E-17Q,
422 -1.519634620380959966438130374006858864624E-14Q,
423 -4.540596528116104986388796594639405114524E-12Q,
424 -7.085151756671466559280490913558388648274E-10Q,
425 -6.351062671323970823761883833531546885452E-8Q,
426 -3.390817171111032905297982523519503522491E-6Q,
427 -1.082340897018886970282138836861233213972E-4Q,
428 -2.020120801187226444822977006648252379508E-3Q,
429 -2.093169910981725694937457070649605557555E-2Q,
430 -1.092176538874275712359269481414448063393E-1Q,
431 -2.374790947854765809203590474789108718733E-1Q,
432 -1.365364204556573800719985118029601401323E-1Q,
433};
434#define NQ8_16D 11
435static const __float128 Q8_16D[NQ8_16D + 1] = {
436 1.978397614733632533581207058069628242280E-16Q,
437 1.487361156806202736877009608336766720560E-13Q,
438 4.468041406888412086042576067133365913456E-11Q,
439 7.027822074821007443672290507210594648877E-9Q,
440 6.375740580686101224127290062867976007374E-7Q,
441 3.466887658320002225888644977076410421940E-5Q,
442 1.138625640905289601186353909213719596986E-3Q,
443 2.224470799470414663443449818235008486439E-2Q,
444 2.487052928527244907490589787691478482358E-1Q,
445 1.483927406564349124649083853892380899217E0Q,
446 4.182773513276056975777258788903489507705E0Q,
447 4.419665392573449746043880892524360870944E0Q,
448 /* 1.000000000000000000000000000000000000000E0 */
449};
450
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 */
455#define NQ5_8N 10
456static const __float128 Q5_8N[NQ5_8N + 1] = {
457 -3.656082407740970534915918390488336879763E-13Q,
458 -1.344660308497244804752334556734121771023E-10Q,
459 -1.909765035234071738548629788698150760791E-8Q,
460 -1.366668038160120210269389551283666716453E-6Q,
461 -5.392327355984269366895210704976314135683E-5Q,
462 -1.206268245713024564674432357634540343884E-3Q,
463 -1.515456784370354374066417703736088291287E-2Q,
464 -1.022454301137286306933217746545237098518E-1Q,
465 -3.373438906472495080504907858424251082240E-1Q,
466 -4.510782522110845697262323973549178453405E-1Q,
467 -1.549000892545288676809660828213589804884E-1Q,
468};
469#define NQ5_8D 10
470static const __float128 Q5_8D[NQ5_8D + 1] = {
471 3.565550843359501079050699598913828460036E-12Q,
472 1.321016015556560621591847454285330528045E-9Q,
473 1.897542728662346479999969679234270605975E-7Q,
474 1.381720283068706710298734234287456219474E-5Q,
475 5.599248147286524662305325795203422873725E-4Q,
476 1.305442352653121436697064782499122164843E-2Q,
477 1.750234079626943298160445750078631894985E-1Q,
478 1.311420542073436520965439883806946678491E0Q,
479 5.162757689856842406744504211089724926650E0Q,
480 9.527760296384704425618556332087850581308E0Q,
481 6.604648207463236667912921642545100248584E0Q,
482 /* 1.000000000000000000000000000000000000000E0 */
483};
484
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 */
489#define NQ4_5N 10
490static const __float128 Q4_5N[NQ4_5N + 1] = {
491 -4.079513568708891749424783046520200903755E-11Q,
492 -9.326548104106791766891812583019664893311E-9Q,
493 -8.016795121318423066292906123815687003356E-7Q,
494 -3.372350544043594415609295225664186750995E-5Q,
495 -7.566238665947967882207277686375417983917E-4Q,
496 -9.248861580055565402130441618521591282617E-3Q,
497 -6.033106131055851432267702948850231270338E-2Q,
498 -1.966908754799996793730369265431584303447E-1Q,
499 -2.791062741179964150755788226623462207560E-1Q,
500 -1.255478605849190549914610121863534191666E-1Q,
501 -4.320429862021265463213168186061696944062E-3Q,
502};
503#define NQ4_5D 9
504static const __float128 Q4_5D[NQ4_5D + 1] = {
505 3.978497042580921479003851216297330701056E-10Q,
506 9.203304163828145809278568906420772246666E-8Q,
507 8.059685467088175644915010485174545743798E-6Q,
508 3.490187375993956409171098277561669167446E-4Q,
509 8.189109654456872150100501732073810028829E-3Q,
510 1.072572867311023640958725265762483033769E-1Q,
511 7.790606862409960053675717185714576937994E-1Q,
512 3.016049768232011196434185423512777656328E0Q,
513 5.722963851442769787733717162314477949360E0Q,
514 4.510527838428473279647251350931380867663E0Q,
515 /* 1.000000000000000000000000000000000000000E0 */
516};
517
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 */
522#define NQ3r2_4N 9
523static const __float128 Q3r2_4N[NQ3r2_4N + 1] = {
524 -1.087480809271383885936921889040388133627E-8Q,
525 -1.690067828697463740906962973479310170932E-6Q,
526 -9.608064416995105532790745641974762550982E-5Q,
527 -2.594198839156517191858208513873961837410E-3Q,
528 -3.610954144421543968160459863048062977822E-2Q,
529 -2.629866798251843212210482269563961685666E-1Q,
530 -9.709186825881775885917984975685752956660E-1Q,
531 -1.667521829918185121727268867619982417317E0Q,
532 -1.109255082925540057138766105229900943501E0Q,
533 -1.812932453006641348145049323713469043328E-1Q,
534};
535#define NQ3r2_4D 9
536static const __float128 Q3r2_4D[NQ3r2_4D + 1] = {
537 1.060552717496912381388763753841473407026E-7Q,
538 1.676928002024920520786883649102388708024E-5Q,
539 9.803481712245420839301400601140812255737E-4Q,
540 2.765559874262309494758505158089249012930E-2Q,
541 4.117921827792571791298862613287549140706E-1Q,
542 3.323769515244751267093378361930279161413E0Q,
543 1.436602494405814164724810151689705353670E1Q,
544 3.163087869617098638064881410646782408297E1Q,
545 3.198181264977021649489103980298349589419E1Q,
546 1.203649258862068431199471076202897823272E1Q,
547 /* 1.000000000000000000000000000000000000000E0 */
548};
549
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 */
554#define NQ2r7_3r2N 9
555static const __float128 Q2r7_3r2N[NQ2r7_3r2N + 1] = {
556 -1.723405393982209853244278760171643219530E-7Q,
557 -2.090508758514655456365709712333460087442E-5Q,
558 -9.140104013370974823232873472192719263019E-4Q,
559 -1.871349499990714843332742160292474780128E-2Q,
560 -1.948930738119938669637865956162512983416E-1Q,
561 -1.048764684978978127908439526343174139788E0Q,
562 -2.827714929925679500237476105843643064698E0Q,
563 -3.508761569156476114276988181329773987314E0Q,
564 -1.669332202790211090973255098624488308989E0Q,
565 -1.930796319299022954013840684651016077770E-1Q,
566};
567#define NQ2r7_3r2D 9
568static const __float128 Q2r7_3r2D[NQ2r7_3r2D + 1] = {
569 1.680730662300831976234547482334347983474E-6Q,
570 2.084241442440551016475972218719621841120E-4Q,
571 9.445316642108367479043541702688736295579E-3Q,
572 2.044637889456631896650179477133252184672E-1Q,
573 2.316091982244297350829522534435350078205E0Q,
574 1.412031891783015085196708811890448488865E1Q,
575 4.583830154673223384837091077279595496149E1Q,
576 7.549520609270909439885998474045974122261E1Q,
577 5.697605832808113367197494052388203310638E1Q,
578 1.601496240876192444526383314589371686234E1Q,
579 /* 1.000000000000000000000000000000000000000E0 */
580};
581
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 */
586#define NQ2r3_2r7N 9
587static const __float128 Q2r3_2r7N[NQ2r3_2r7N + 1] = {
588 -8.603042076329122085722385914954878953775E-7Q,
589 -7.701746260451647874214968882605186675720E-5Q,
590 -2.407932004380727587382493696877569654271E-3Q,
591 -3.403434217607634279028110636919987224188E-2Q,
592 -2.348707332185238159192422084985713102877E-1Q,
593 -7.957498841538254916147095255700637463207E-1Q,
594 -1.258469078442635106431098063707934348577E0Q,
595 -8.162415474676345812459353639449971369890E-1Q,
596 -1.581783890269379690141513949609572806898E-1Q,
597 -1.890595651683552228232308756569450822905E-3Q,
598};
599#define NQ2r3_2r7D 8
600static const __float128 Q2r3_2r7D[NQ2r3_2r7D + 1] = {
601 8.390017524798316921170710533381568175665E-6Q,
602 7.738148683730826286477254659973968763659E-4Q,
603 2.541480810958665794368759558791634341779E-2Q,
604 3.878879789711276799058486068562386244873E-1Q,
605 3.003783779325811292142957336802456109333E0Q,
606 1.206480374773322029883039064575464497400E1Q,
607 2.458414064785315978408974662900438351782E1Q,
608 2.367237826273668567199042088835448715228E1Q,
609 9.231451197519171090875569102116321676763E0Q,
610 /* 1.000000000000000000000000000000000000000E0 */
611};
612
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 */
617#define NQ2_2r3N 9
618static const __float128 Q2_2r3N[NQ2_2r3N + 1] = {
619 -5.552507516089087822166822364590806076174E-6Q,
620 -4.135067659799500521040944087433752970297E-4Q,
621 -1.059928728869218962607068840646564457980E-2Q,
622 -1.212070036005832342565792241385459023801E-1Q,
623 -6.688350110633603958684302153362735625156E-1Q,
624 -1.793587878197360221340277951304429821582E0Q,
625 -2.225407682237197485644647380483725045326E0Q,
626 -1.123402135458940189438898496348239744403E0Q,
627 -1.679187241566347077204805190763597299805E-1Q,
628 -1.458550613639093752909985189067233504148E-3Q,
629};
630#define NQ2_2r3D 8
631static const __float128 Q2_2r3D[NQ2_2r3D + 1] = {
632 5.415024336507980465169023996403597916115E-5Q,
633 4.179246497380453022046357404266022870788E-3Q,
634 1.136306384261959483095442402929502368598E-1Q,
635 1.422640343719842213484515445393284072830E0Q,
636 8.968786703393158374728850922289204805764E0Q,
637 2.914542473339246127533384118781216495934E1Q,
638 4.781605421020380669870197378210457054685E1Q,
639 3.693865837171883152382820584714795072937E1Q,
640 1.153220502744204904763115556224395893076E1Q,
641 /* 1.000000000000000000000000000000000000000E0 */
642};
643
644
645/* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
646
647static __float128
648neval (__float128 x, const __float128 *p, int n)
649{
650 __float128 y;
651
652 p += n;
653 y = *p--;
654 do
655 {
656 y = y * x + *p--;
657 }
658 while (--n > 0);
659 return y;
660}
661
662
663/* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */
664
665static __float128
666deval (__float128 x, const __float128 *p, int n)
667{
668 __float128 y;
669
670 p += n;
671 y = x + *p--;
672 do
673 {
674 y = y * x + *p--;
675 }
676 while (--n > 0);
677 return y;
678}
679
680
681/* Bessel function of the first kind, order one. */
682
683__float128
684j1q (__float128 x)
685{
686 __float128 xx, xinv, z, p, q, c, s, cc, ss;
687
688 if (! finiteq (x))
689 {
690 if (x != x)
691 return x + x;
692 else
693 return 0.0Q;
694 }
695 if (x == 0.0Q)
696 return x;
697 xx = fabsq (x);
698 if (xx <= 0x1p-58Q)
699 {
700 __float128 ret = x * 0.5Q;
701 math_check_force_underflow (ret);
702 if (ret == 0)
703 errno = ERANGE;
704 return ret;
705 }
706 if (xx <= 2.0Q)
707 {
708 /* 0 <= x <= 2 */
709 z = xx * xx;
710 p = xx * z * neval (z, J0_2N, NJ0_2N) / deval (z, J0_2D, NJ0_2D);
711 p += 0.5Q * xx;
712 if (x < 0)
713 p = -p;
714 return p;
715 }
716
717 /* X = x - 3 pi/4
718 cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
719 = 1/sqrt(2) * (-cos(x) + sin(x))
720 sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
721 = -1/sqrt(2) * (sin(x) + cos(x))
722 cf. Fdlibm. */
723 sincosq (xx, &s, &c);
724 ss = -s - c;
725 cc = s - c;
726 if (xx <= FLT128_MAX / 2.0Q)
727 {
728 z = cosq (xx + xx);
729 if ((s * c) > 0)
730 cc = z / ss;
731 else
732 ss = z / cc;
733 }
734
735 if (xx > 0x1p256Q)
736 {
737 z = ONEOSQPI * cc / sqrtq (xx);
738 if (x < 0)
739 z = -z;
740 return z;
741 }
742
743 xinv = 1.0Q / xx;
744 z = xinv * xinv;
745 if (xinv <= 0.25)
746 {
747 if (xinv <= 0.125)
748 {
749 if (xinv <= 0.0625)
750 {
751 p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID);
752 q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID);
753 }
754 else
755 {
756 p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D);
757 q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D);
758 }
759 }
760 else if (xinv <= 0.1875)
761 {
762 p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D);
763 q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D);
764 }
765 else
766 {
767 p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D);
768 q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D);
769 }
770 } /* .25 */
771 else /* if (xinv <= 0.5) */
772 {
773 if (xinv <= 0.375)
774 {
775 if (xinv <= 0.3125)
776 {
777 p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D);
778 q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D);
779 }
780 else
781 {
782 p = neval (z, P2r7_3r2N, NP2r7_3r2N)
783 / deval (z, P2r7_3r2D, NP2r7_3r2D);
784 q = neval (z, Q2r7_3r2N, NQ2r7_3r2N)
785 / deval (z, Q2r7_3r2D, NQ2r7_3r2D);
786 }
787 }
788 else if (xinv <= 0.4375)
789 {
790 p = neval (z, P2r3_2r7N, NP2r3_2r7N)
791 / deval (z, P2r3_2r7D, NP2r3_2r7D);
792 q = neval (z, Q2r3_2r7N, NQ2r3_2r7N)
793 / deval (z, Q2r3_2r7D, NQ2r3_2r7D);
794 }
795 else
796 {
797 p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
798 q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
799 }
800 }
801 p = 1.0Q + z * p;
802 q = z * q;
803 q = q * xinv + 0.375Q * xinv;
804 z = ONEOSQPI * (p * cc - q * ss) / sqrtq (xx);
805 if (x < 0)
806 z = -z;
807 return z;
808}
809
810
811/* Y1(x) = 2/pi * (log(x) * J1(x) - 1/x) + x R(x^2)
812 Peak relative error 6.2e-38
813 0 <= x <= 2 */
814#define NY0_2N 7
815static __float128 Y0_2N[NY0_2N + 1] = {
816 -6.804415404830253804408698161694720833249E19Q,
817 1.805450517967019908027153056150465849237E19Q,
818 -8.065747497063694098810419456383006737312E17Q,
819 1.401336667383028259295830955439028236299E16Q,
820 -1.171654432898137585000399489686629680230E14Q,
821 5.061267920943853732895341125243428129150E11Q,
822 -1.096677850566094204586208610960870217970E9Q,
823 9.541172044989995856117187515882879304461E5Q,
824};
825#define NY0_2D 7
826static __float128 Y0_2D[NY0_2D + 1] = {
827 3.470629591820267059538637461549677594549E20Q,
828 4.120796439009916326855848107545425217219E18Q,
829 2.477653371652018249749350657387030814542E16Q,
830 9.954678543353888958177169349272167762797E13Q,
831 2.957927997613630118216218290262851197754E11Q,
832 6.748421382188864486018861197614025972118E8Q,
833 1.173453425218010888004562071020305709319E6Q,
834 1.450335662961034949894009554536003377187E3Q,
835 /* 1.000000000000000000000000000000000000000E0 */
836};
837
838
839/* Bessel function of the second kind, order one. */
840
841__float128
842y1q (__float128 x)
843{
844 __float128 xx, xinv, z, p, q, c, s, cc, ss;
845
846 if (! finiteq (x))
847 return 1 / (x + x * x);
848 if (x <= 0.0Q)
849 {
850 if (x < 0.0Q)
851 return (zero / (zero * x));
852 return -1 / zero; /* -inf and divide by zero exception. */
853 }
854 xx = fabsq (x);
855 if (xx <= 0x1p-114)
856 {
857 z = -TWOOPI / x;
858 if (isinfq (z))
859 errno = ERANGE;
860 return z;
861 }
862 if (xx <= 2.0Q)
863 {
864 /* 0 <= x <= 2 */
865 /* FIXME: SET_RESTORE_ROUNDL (FE_TONEAREST); */
866 z = xx * xx;
867 p = xx * neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D);
868 p = -TWOOPI / xx + p;
869 p = TWOOPI * logq (x) * j1q (x) + p;
870 return p;
871 }
872
873 /* X = x - 3 pi/4
874 cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
875 = 1/sqrt(2) * (-cos(x) + sin(x))
876 sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
877 = -1/sqrt(2) * (sin(x) + cos(x))
878 cf. Fdlibm. */
879 sincosq (xx, &s, &c);
880 ss = -s - c;
881 cc = s - c;
882 if (xx <= FLT128_MAX / 2.0Q)
883 {
884 z = cosq (xx + xx);
885 if ((s * c) > 0)
886 cc = z / ss;
887 else
888 ss = z / cc;
889 }
890
891 if (xx > 0x1p256Q)
892 return ONEOSQPI * ss / sqrtq (xx);
893
894 xinv = 1.0Q / xx;
895 z = xinv * xinv;
896 if (xinv <= 0.25)
897 {
898 if (xinv <= 0.125)
899 {
900 if (xinv <= 0.0625)
901 {
902 p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID);
903 q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID);
904 }
905 else
906 {
907 p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D);
908 q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D);
909 }
910 }
911 else if (xinv <= 0.1875)
912 {
913 p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D);
914 q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D);
915 }
916 else
917 {
918 p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D);
919 q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D);
920 }
921 } /* .25 */
922 else /* if (xinv <= 0.5) */
923 {
924 if (xinv <= 0.375)
925 {
926 if (xinv <= 0.3125)
927 {
928 p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D);
929 q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D);
930 }
931 else
932 {
933 p = neval (z, P2r7_3r2N, NP2r7_3r2N)
934 / deval (z, P2r7_3r2D, NP2r7_3r2D);
935 q = neval (z, Q2r7_3r2N, NQ2r7_3r2N)
936 / deval (z, Q2r7_3r2D, NQ2r7_3r2D);
937 }
938 }
939 else if (xinv <= 0.4375)
940 {
941 p = neval (z, P2r3_2r7N, NP2r3_2r7N)
942 / deval (z, P2r3_2r7D, NP2r3_2r7D);
943 q = neval (z, Q2r3_2r7N, NQ2r3_2r7N)
944 / deval (z, Q2r3_2r7D, NQ2r3_2r7D);
945 }
946 else
947 {
948 p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
949 q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
950 }
951 }
952 p = 1.0Q + z * p;
953 q = z * q;
954 q = q * xinv + 0.375Q * xinv;
955 z = ONEOSQPI * (p * ss + q * cc) / sqrtq (xx);
956 return z;
957}
958