1 /* j1l.c
2 *
3 * Bessel function of order one
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * long double x, y, j1l();
10 *
11 * y = j1l( 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 * double x, y, y1l();
56 *
57 * y = y1l( 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, see
96 <http://www.gnu.org/licenses/>. */
97
98 #include "quadmath-imp.h"
99
100 /* 1 / sqrt(pi) */
101 static const __float128 ONEOSQPI = 5.6418958354775628694807945156077258584405E-1Q;
102 /* 2 / pi */
103 static const __float128 TWOOPI = 6.3661977236758134307553505349005744813784E-1Q;
104 static const __float128 zero = 0;
105
106 /* J1(x) = .5x + x x^2 R(x^2)
107 Peak relative error 1.9e-35
108 0 <= x <= 2 */
109 #define NJ0_2N 6
110 static const __float128 J0_2N[NJ0_2N + 1] = {
111 -5.943799577386942855938508697619735179660E16Q,
112 1.812087021305009192259946997014044074711E15Q,
113 -2.761698314264509665075127515729146460895E13Q,
114 2.091089497823600978949389109350658815972E11Q,
115 -8.546413231387036372945453565654130054307E8Q,
116 1.797229225249742247475464052741320612261E6Q,
117 -1.559552840946694171346552770008812083969E3Q
118 };
119 #define NJ0_2D 6
120 static const __float128 J0_2D[NJ0_2D + 1] = {
121 9.510079323819108569501613916191477479397E17Q,
122 1.063193817503280529676423936545854693915E16Q,
123 5.934143516050192600795972192791775226920E13Q,
124 2.168000911950620999091479265214368352883E11Q,
125 5.673775894803172808323058205986256928794E8Q,
126 1.080329960080981204840966206372671147224E6Q,
127 1.411951256636576283942477881535283304912E3Q,
128 /* 1.000000000000000000000000000000000000000E0L */
129 };
130
131 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
132 0 <= 1/x <= .0625
133 Peak relative error 3.6e-36 */
134 #define NP16_IN 9
135 static const __float128 P16_IN[NP16_IN + 1] = {
136 5.143674369359646114999545149085139822905E-16Q,
137 4.836645664124562546056389268546233577376E-13Q,
138 1.730945562285804805325011561498453013673E-10Q,
139 3.047976856147077889834905908605310585810E-8Q,
140 2.855227609107969710407464739188141162386E-6Q,
141 1.439362407936705484122143713643023998457E-4Q,
142 3.774489768532936551500999699815873422073E-3Q,
143 4.723962172984642566142399678920790598426E-2Q,
144 2.359289678988743939925017240478818248735E-1Q,
145 3.032580002220628812728954785118117124520E-1Q,
146 };
147 #define NP16_ID 9
148 static const __float128 P16_ID[NP16_ID + 1] = {
149 4.389268795186898018132945193912677177553E-15Q,
150 4.132671824807454334388868363256830961655E-12Q,
151 1.482133328179508835835963635130894413136E-9Q,
152 2.618941412861122118906353737117067376236E-7Q,
153 2.467854246740858470815714426201888034270E-5Q,
154 1.257192927368839847825938545925340230490E-3Q,
155 3.362739031941574274949719324644120720341E-2Q,
156 4.384458231338934105875343439265370178858E-1Q,
157 2.412830809841095249170909628197264854651E0Q,
158 4.176078204111348059102962617368214856874E0Q,
159 /* 1.000000000000000000000000000000000000000E0 */
160 };
161
162 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
163 0.0625 <= 1/x <= 0.125
164 Peak relative error 1.9e-36 */
165 #define NP8_16N 11
166 static const __float128 P8_16N[NP8_16N + 1] = {
167 2.984612480763362345647303274082071598135E-16Q,
168 1.923651877544126103941232173085475682334E-13Q,
169 4.881258879388869396043760693256024307743E-11Q,
170 6.368866572475045408480898921866869811889E-9Q,
171 4.684818344104910450523906967821090796737E-7Q,
172 2.005177298271593587095982211091300382796E-5Q,
173 4.979808067163957634120681477207147536182E-4Q,
174 6.946005761642579085284689047091173581127E-3Q,
175 5.074601112955765012750207555985299026204E-2Q,
176 1.698599455896180893191766195194231825379E-1Q,
177 1.957536905259237627737222775573623779638E-1Q,
178 2.991314703282528370270179989044994319374E-2Q,
179 };
180 #define NP8_16D 10
181 static const __float128 P8_16D[NP8_16D + 1] = {
182 2.546869316918069202079580939942463010937E-15Q,
183 1.644650111942455804019788382157745229955E-12Q,
184 4.185430770291694079925607420808011147173E-10Q,
185 5.485331966975218025368698195861074143153E-8Q,
186 4.062884421686912042335466327098932678905E-6Q,
187 1.758139661060905948870523641319556816772E-4Q,
188 4.445143889306356207566032244985607493096E-3Q,
189 6.391901016293512632765621532571159071158E-2Q,
190 4.933040207519900471177016015718145795434E-1Q,
191 1.839144086168947712971630337250761842976E0Q,
192 2.715120873995490920415616716916149586579E0Q,
193 /* 1.000000000000000000000000000000000000000E0 */
194 };
195
196 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
197 0.125 <= 1/x <= 0.1875
198 Peak relative error 1.3e-36 */
199 #define NP5_8N 10
200 static const __float128 P5_8N[NP5_8N + 1] = {
201 2.837678373978003452653763806968237227234E-12Q,
202 9.726641165590364928442128579282742354806E-10Q,
203 1.284408003604131382028112171490633956539E-7Q,
204 8.524624695868291291250573339272194285008E-6Q,
205 3.111516908953172249853673787748841282846E-4Q,
206 6.423175156126364104172801983096596409176E-3Q,
207 7.430220589989104581004416356260692450652E-2Q,
208 4.608315409833682489016656279567605536619E-1Q,
209 1.396870223510964882676225042258855977512E0Q,
210 1.718500293904122365894630460672081526236E0Q,
211 5.465927698800862172307352821870223855365E-1Q
212 };
213 #define NP5_8D 10
214 static const __float128 P5_8D[NP5_8D + 1] = {
215 2.421485545794616609951168511612060482715E-11Q,
216 8.329862750896452929030058039752327232310E-9Q,
217 1.106137992233383429630592081375289010720E-6Q,
218 7.405786153760681090127497796448503306939E-5Q,
219 2.740364785433195322492093333127633465227E-3Q,
220 5.781246470403095224872243564165254652198E-2Q,
221 6.927711353039742469918754111511109983546E-1Q,
222 4.558679283460430281188304515922826156690E0Q,
223 1.534468499844879487013168065728837900009E1Q,
224 2.313927430889218597919624843161569422745E1Q,
225 1.194506341319498844336768473218382828637E1Q,
226 /* 1.000000000000000000000000000000000000000E0 */
227 };
228
229 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
230 Peak relative error 1.4e-36
231 0.1875 <= 1/x <= 0.25 */
232 #define NP4_5N 10
233 static const __float128 P4_5N[NP4_5N + 1] = {
234 1.846029078268368685834261260420933914621E-10Q,
235 3.916295939611376119377869680335444207768E-8Q,
236 3.122158792018920627984597530935323997312E-6Q,
237 1.218073444893078303994045653603392272450E-4Q,
238 2.536420827983485448140477159977981844883E-3Q,
239 2.883011322006690823959367922241169171315E-2Q,
240 1.755255190734902907438042414495469810830E-1Q,
241 5.379317079922628599870898285488723736599E-1Q,
242 7.284904050194300773890303361501726561938E-1Q,
243 3.270110346613085348094396323925000362813E-1Q,
244 1.804473805689725610052078464951722064757E-2Q,
245 };
246 #define NP4_5D 9
247 static const __float128 P4_5D[NP4_5D + 1] = {
248 1.575278146806816970152174364308980863569E-9Q,
249 3.361289173657099516191331123405675054321E-7Q,
250 2.704692281550877810424745289838790693708E-5Q,
251 1.070854930483999749316546199273521063543E-3Q,
252 2.282373093495295842598097265627962125411E-2Q,
253 2.692025460665354148328762368240343249830E-1Q,
254 1.739892942593664447220951225734811133759E0Q,
255 5.890727576752230385342377570386657229324E0Q,
256 9.517442287057841500750256954117735128153E0Q,
257 6.100616353935338240775363403030137736013E0Q,
258 /* 1.000000000000000000000000000000000000000E0 */
259 };
260
261 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
262 Peak relative error 3.0e-36
263 0.25 <= 1/x <= 0.3125 */
264 #define NP3r2_4N 9
265 static const __float128 P3r2_4N[NP3r2_4N + 1] = {
266 8.240803130988044478595580300846665863782E-8Q,
267 1.179418958381961224222969866406483744580E-5Q,
268 6.179787320956386624336959112503824397755E-4Q,
269 1.540270833608687596420595830747166658383E-2Q,
270 1.983904219491512618376375619598837355076E-1Q,
271 1.341465722692038870390470651608301155565E0Q,
272 4.617865326696612898792238245990854646057E0Q,
273 7.435574801812346424460233180412308000587E0Q,
274 4.671327027414635292514599201278557680420E0Q,
275 7.299530852495776936690976966995187714739E-1Q,
276 };
277 #define NP3r2_4D 9
278 static const __float128 P3r2_4D[NP3r2_4D + 1] = {
279 7.032152009675729604487575753279187576521E-7Q,
280 1.015090352324577615777511269928856742848E-4Q,
281 5.394262184808448484302067955186308730620E-3Q,
282 1.375291438480256110455809354836988584325E-1Q,
283 1.836247144461106304788160919310404376670E0Q,
284 1.314378564254376655001094503090935880349E1Q,
285 4.957184590465712006934452500894672343488E1Q,
286 9.287394244300647738855415178790263465398E1Q,
287 7.652563275535900609085229286020552768399E1Q,
288 2.147042473003074533150718117770093209096E1Q,
289 /* 1.000000000000000000000000000000000000000E0 */
290 };
291
292 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
293 Peak relative error 1.0e-35
294 0.3125 <= 1/x <= 0.375 */
295 #define NP2r7_3r2N 9
296 static const __float128 P2r7_3r2N[NP2r7_3r2N + 1] = {
297 4.599033469240421554219816935160627085991E-7Q,
298 4.665724440345003914596647144630893997284E-5Q,
299 1.684348845667764271596142716944374892756E-3Q,
300 2.802446446884455707845985913454440176223E-2Q,
301 2.321937586453963310008279956042545173930E-1Q,
302 9.640277413988055668692438709376437553804E-1Q,
303 1.911021064710270904508663334033003246028E0Q,
304 1.600811610164341450262992138893970224971E0Q,
305 4.266299218652587901171386591543457861138E-1Q,
306 1.316470424456061252962568223251247207325E-2Q,
307 };
308 #define NP2r7_3r2D 8
309 static const __float128 P2r7_3r2D[NP2r7_3r2D + 1] = {
310 3.924508608545520758883457108453520099610E-6Q,
311 4.029707889408829273226495756222078039823E-4Q,
312 1.484629715787703260797886463307469600219E-2Q,
313 2.553136379967180865331706538897231588685E-1Q,
314 2.229457223891676394409880026887106228740E0Q,
315 1.005708903856384091956550845198392117318E1Q,
316 2.277082659664386953166629360352385889558E1Q,
317 2.384726835193630788249826630376533988245E1Q,
318 9.700989749041320895890113781610939632410E0Q,
319 /* 1.000000000000000000000000000000000000000E0 */
320 };
321
322 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
323 Peak relative error 1.7e-36
324 0.3125 <= 1/x <= 0.4375 */
325 #define NP2r3_2r7N 9
326 static const __float128 P2r3_2r7N[NP2r3_2r7N + 1] = {
327 3.916766777108274628543759603786857387402E-6Q,
328 3.212176636756546217390661984304645137013E-4Q,
329 9.255768488524816445220126081207248947118E-3Q,
330 1.214853146369078277453080641911700735354E-1Q,
331 7.855163309847214136198449861311404633665E-1Q,
332 2.520058073282978403655488662066019816540E0Q,
333 3.825136484837545257209234285382183711466E0Q,
334 2.432569427554248006229715163865569506873E0Q,
335 4.877934835018231178495030117729800489743E-1Q,
336 1.109902737860249670981355149101343427885E-2Q,
337 };
338 #define NP2r3_2r7D 8
339 static const __float128 P2r3_2r7D[NP2r3_2r7D + 1] = {
340 3.342307880794065640312646341190547184461E-5Q,
341 2.782182891138893201544978009012096558265E-3Q,
342 8.221304931614200702142049236141249929207E-2Q,
343 1.123728246291165812392918571987858010949E0Q,
344 7.740482453652715577233858317133423434590E0Q,
345 2.737624677567945952953322566311201919139E1Q,
346 4.837181477096062403118304137851260715475E1Q,
347 3.941098643468580791437772701093795299274E1Q,
348 1.245821247166544627558323920382547533630E1Q,
349 /* 1.000000000000000000000000000000000000000E0 */
350 };
351
352 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
353 Peak relative error 1.7e-35
354 0.4375 <= 1/x <= 0.5 */
355 #define NP2_2r3N 8
356 static const __float128 P2_2r3N[NP2_2r3N + 1] = {
357 3.397930802851248553545191160608731940751E-4Q,
358 2.104020902735482418784312825637833698217E-2Q,
359 4.442291771608095963935342749477836181939E-1Q,
360 4.131797328716583282869183304291833754967E0Q,
361 1.819920169779026500146134832455189917589E1Q,
362 3.781779616522937565300309684282401791291E1Q,
363 3.459605449728864218972931220783543410347E1Q,
364 1.173594248397603882049066603238568316561E1Q,
365 9.455702270242780642835086549285560316461E-1Q,
366 };
367 #define NP2_2r3D 8
368 static const __float128 P2_2r3D[NP2_2r3D + 1] = {
369 2.899568897241432883079888249845707400614E-3Q,
370 1.831107138190848460767699919531132426356E-1Q,
371 3.999350044057883839080258832758908825165E0Q,
372 3.929041535867957938340569419874195303712E1Q,
373 1.884245613422523323068802689915538908291E2Q,
374 4.461469948819229734353852978424629815929E2Q,
375 5.004998753999796821224085972610636347903E2Q,
376 2.386342520092608513170837883757163414100E2Q,
377 3.791322528149347975999851588922424189957E1Q,
378 /* 1.000000000000000000000000000000000000000E0 */
379 };
380
381 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
382 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
383 Peak relative error 8.0e-36
384 0 <= 1/x <= .0625 */
385 #define NQ16_IN 10
386 static const __float128 Q16_IN[NQ16_IN + 1] = {
387 -3.917420835712508001321875734030357393421E-18Q,
388 -4.440311387483014485304387406538069930457E-15Q,
389 -1.951635424076926487780929645954007139616E-12Q,
390 -4.318256438421012555040546775651612810513E-10Q,
391 -5.231244131926180765270446557146989238020E-8Q,
392 -3.540072702902043752460711989234732357653E-6Q,
393 -1.311017536555269966928228052917534882984E-4Q,
394 -2.495184669674631806622008769674827575088E-3Q,
395 -2.141868222987209028118086708697998506716E-2Q,
396 -6.184031415202148901863605871197272650090E-2Q,
397 -1.922298704033332356899546792898156493887E-2Q,
398 };
399 #define NQ16_ID 9
400 static const __float128 Q16_ID[NQ16_ID + 1] = {
401 3.820418034066293517479619763498400162314E-17Q,
402 4.340702810799239909648911373329149354911E-14Q,
403 1.914985356383416140706179933075303538524E-11Q,
404 4.262333682610888819476498617261895474330E-9Q,
405 5.213481314722233980346462747902942182792E-7Q,
406 3.585741697694069399299005316809954590558E-5Q,
407 1.366513429642842006385029778105539457546E-3Q,
408 2.745282599850704662726337474371355160594E-2Q,
409 2.637644521611867647651200098449903330074E-1Q,
410 1.006953426110765984590782655598680488746E0Q,
411 /* 1.000000000000000000000000000000000000000E0 */
412 };
413
414 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
415 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
416 Peak relative error 1.9e-36
417 0.0625 <= 1/x <= 0.125 */
418 #define NQ8_16N 11
419 static const __float128 Q8_16N[NQ8_16N + 1] = {
420 -2.028630366670228670781362543615221542291E-17Q,
421 -1.519634620380959966438130374006858864624E-14Q,
422 -4.540596528116104986388796594639405114524E-12Q,
423 -7.085151756671466559280490913558388648274E-10Q,
424 -6.351062671323970823761883833531546885452E-8Q,
425 -3.390817171111032905297982523519503522491E-6Q,
426 -1.082340897018886970282138836861233213972E-4Q,
427 -2.020120801187226444822977006648252379508E-3Q,
428 -2.093169910981725694937457070649605557555E-2Q,
429 -1.092176538874275712359269481414448063393E-1Q,
430 -2.374790947854765809203590474789108718733E-1Q,
431 -1.365364204556573800719985118029601401323E-1Q,
432 };
433 #define NQ8_16D 11
434 static const __float128 Q8_16D[NQ8_16D + 1] = {
435 1.978397614733632533581207058069628242280E-16Q,
436 1.487361156806202736877009608336766720560E-13Q,
437 4.468041406888412086042576067133365913456E-11Q,
438 7.027822074821007443672290507210594648877E-9Q,
439 6.375740580686101224127290062867976007374E-7Q,
440 3.466887658320002225888644977076410421940E-5Q,
441 1.138625640905289601186353909213719596986E-3Q,
442 2.224470799470414663443449818235008486439E-2Q,
443 2.487052928527244907490589787691478482358E-1Q,
444 1.483927406564349124649083853892380899217E0Q,
445 4.182773513276056975777258788903489507705E0Q,
446 4.419665392573449746043880892524360870944E0Q,
447 /* 1.000000000000000000000000000000000000000E0 */
448 };
449
450 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
451 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
452 Peak relative error 1.5e-35
453 0.125 <= 1/x <= 0.1875 */
454 #define NQ5_8N 10
455 static const __float128 Q5_8N[NQ5_8N + 1] = {
456 -3.656082407740970534915918390488336879763E-13Q,
457 -1.344660308497244804752334556734121771023E-10Q,
458 -1.909765035234071738548629788698150760791E-8Q,
459 -1.366668038160120210269389551283666716453E-6Q,
460 -5.392327355984269366895210704976314135683E-5Q,
461 -1.206268245713024564674432357634540343884E-3Q,
462 -1.515456784370354374066417703736088291287E-2Q,
463 -1.022454301137286306933217746545237098518E-1Q,
464 -3.373438906472495080504907858424251082240E-1Q,
465 -4.510782522110845697262323973549178453405E-1Q,
466 -1.549000892545288676809660828213589804884E-1Q,
467 };
468 #define NQ5_8D 10
469 static const __float128 Q5_8D[NQ5_8D + 1] = {
470 3.565550843359501079050699598913828460036E-12Q,
471 1.321016015556560621591847454285330528045E-9Q,
472 1.897542728662346479999969679234270605975E-7Q,
473 1.381720283068706710298734234287456219474E-5Q,
474 5.599248147286524662305325795203422873725E-4Q,
475 1.305442352653121436697064782499122164843E-2Q,
476 1.750234079626943298160445750078631894985E-1Q,
477 1.311420542073436520965439883806946678491E0Q,
478 5.162757689856842406744504211089724926650E0Q,
479 9.527760296384704425618556332087850581308E0Q,
480 6.604648207463236667912921642545100248584E0Q,
481 /* 1.000000000000000000000000000000000000000E0 */
482 };
483
484 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
485 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
486 Peak relative error 1.3e-35
487 0.1875 <= 1/x <= 0.25 */
488 #define NQ4_5N 10
489 static const __float128 Q4_5N[NQ4_5N + 1] = {
490 -4.079513568708891749424783046520200903755E-11Q,
491 -9.326548104106791766891812583019664893311E-9Q,
492 -8.016795121318423066292906123815687003356E-7Q,
493 -3.372350544043594415609295225664186750995E-5Q,
494 -7.566238665947967882207277686375417983917E-4Q,
495 -9.248861580055565402130441618521591282617E-3Q,
496 -6.033106131055851432267702948850231270338E-2Q,
497 -1.966908754799996793730369265431584303447E-1Q,
498 -2.791062741179964150755788226623462207560E-1Q,
499 -1.255478605849190549914610121863534191666E-1Q,
500 -4.320429862021265463213168186061696944062E-3Q,
501 };
502 #define NQ4_5D 9
503 static const __float128 Q4_5D[NQ4_5D + 1] = {
504 3.978497042580921479003851216297330701056E-10Q,
505 9.203304163828145809278568906420772246666E-8Q,
506 8.059685467088175644915010485174545743798E-6Q,
507 3.490187375993956409171098277561669167446E-4Q,
508 8.189109654456872150100501732073810028829E-3Q,
509 1.072572867311023640958725265762483033769E-1Q,
510 7.790606862409960053675717185714576937994E-1Q,
511 3.016049768232011196434185423512777656328E0Q,
512 5.722963851442769787733717162314477949360E0Q,
513 4.510527838428473279647251350931380867663E0Q,
514 /* 1.000000000000000000000000000000000000000E0 */
515 };
516
517 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
518 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
519 Peak relative error 2.1e-35
520 0.25 <= 1/x <= 0.3125 */
521 #define NQ3r2_4N 9
522 static const __float128 Q3r2_4N[NQ3r2_4N + 1] = {
523 -1.087480809271383885936921889040388133627E-8Q,
524 -1.690067828697463740906962973479310170932E-6Q,
525 -9.608064416995105532790745641974762550982E-5Q,
526 -2.594198839156517191858208513873961837410E-3Q,
527 -3.610954144421543968160459863048062977822E-2Q,
528 -2.629866798251843212210482269563961685666E-1Q,
529 -9.709186825881775885917984975685752956660E-1Q,
530 -1.667521829918185121727268867619982417317E0Q,
531 -1.109255082925540057138766105229900943501E0Q,
532 -1.812932453006641348145049323713469043328E-1Q,
533 };
534 #define NQ3r2_4D 9
535 static const __float128 Q3r2_4D[NQ3r2_4D + 1] = {
536 1.060552717496912381388763753841473407026E-7Q,
537 1.676928002024920520786883649102388708024E-5Q,
538 9.803481712245420839301400601140812255737E-4Q,
539 2.765559874262309494758505158089249012930E-2Q,
540 4.117921827792571791298862613287549140706E-1Q,
541 3.323769515244751267093378361930279161413E0Q,
542 1.436602494405814164724810151689705353670E1Q,
543 3.163087869617098638064881410646782408297E1Q,
544 3.198181264977021649489103980298349589419E1Q,
545 1.203649258862068431199471076202897823272E1Q,
546 /* 1.000000000000000000000000000000000000000E0 */
547 };
548
549 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
550 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
551 Peak relative error 1.6e-36
552 0.3125 <= 1/x <= 0.375 */
553 #define NQ2r7_3r2N 9
554 static const __float128 Q2r7_3r2N[NQ2r7_3r2N + 1] = {
555 -1.723405393982209853244278760171643219530E-7Q,
556 -2.090508758514655456365709712333460087442E-5Q,
557 -9.140104013370974823232873472192719263019E-4Q,
558 -1.871349499990714843332742160292474780128E-2Q,
559 -1.948930738119938669637865956162512983416E-1Q,
560 -1.048764684978978127908439526343174139788E0Q,
561 -2.827714929925679500237476105843643064698E0Q,
562 -3.508761569156476114276988181329773987314E0Q,
563 -1.669332202790211090973255098624488308989E0Q,
564 -1.930796319299022954013840684651016077770E-1Q,
565 };
566 #define NQ2r7_3r2D 9
567 static const __float128 Q2r7_3r2D[NQ2r7_3r2D + 1] = {
568 1.680730662300831976234547482334347983474E-6Q,
569 2.084241442440551016475972218719621841120E-4Q,
570 9.445316642108367479043541702688736295579E-3Q,
571 2.044637889456631896650179477133252184672E-1Q,
572 2.316091982244297350829522534435350078205E0Q,
573 1.412031891783015085196708811890448488865E1Q,
574 4.583830154673223384837091077279595496149E1Q,
575 7.549520609270909439885998474045974122261E1Q,
576 5.697605832808113367197494052388203310638E1Q,
577 1.601496240876192444526383314589371686234E1Q,
578 /* 1.000000000000000000000000000000000000000E0 */
579 };
580
581 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
582 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
583 Peak relative error 9.5e-36
584 0.375 <= 1/x <= 0.4375 */
585 #define NQ2r3_2r7N 9
586 static const __float128 Q2r3_2r7N[NQ2r3_2r7N + 1] = {
587 -8.603042076329122085722385914954878953775E-7Q,
588 -7.701746260451647874214968882605186675720E-5Q,
589 -2.407932004380727587382493696877569654271E-3Q,
590 -3.403434217607634279028110636919987224188E-2Q,
591 -2.348707332185238159192422084985713102877E-1Q,
592 -7.957498841538254916147095255700637463207E-1Q,
593 -1.258469078442635106431098063707934348577E0Q,
594 -8.162415474676345812459353639449971369890E-1Q,
595 -1.581783890269379690141513949609572806898E-1Q,
596 -1.890595651683552228232308756569450822905E-3Q,
597 };
598 #define NQ2r3_2r7D 8
599 static const __float128 Q2r3_2r7D[NQ2r3_2r7D + 1] = {
600 8.390017524798316921170710533381568175665E-6Q,
601 7.738148683730826286477254659973968763659E-4Q,
602 2.541480810958665794368759558791634341779E-2Q,
603 3.878879789711276799058486068562386244873E-1Q,
604 3.003783779325811292142957336802456109333E0Q,
605 1.206480374773322029883039064575464497400E1Q,
606 2.458414064785315978408974662900438351782E1Q,
607 2.367237826273668567199042088835448715228E1Q,
608 9.231451197519171090875569102116321676763E0Q,
609 /* 1.000000000000000000000000000000000000000E0 */
610 };
611
612 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
613 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
614 Peak relative error 1.4e-36
615 0.4375 <= 1/x <= 0.5 */
616 #define NQ2_2r3N 9
617 static const __float128 Q2_2r3N[NQ2_2r3N + 1] = {
618 -5.552507516089087822166822364590806076174E-6Q,
619 -4.135067659799500521040944087433752970297E-4Q,
620 -1.059928728869218962607068840646564457980E-2Q,
621 -1.212070036005832342565792241385459023801E-1Q,
622 -6.688350110633603958684302153362735625156E-1Q,
623 -1.793587878197360221340277951304429821582E0Q,
624 -2.225407682237197485644647380483725045326E0Q,
625 -1.123402135458940189438898496348239744403E0Q,
626 -1.679187241566347077204805190763597299805E-1Q,
627 -1.458550613639093752909985189067233504148E-3Q,
628 };
629 #define NQ2_2r3D 8
630 static const __float128 Q2_2r3D[NQ2_2r3D + 1] = {
631 5.415024336507980465169023996403597916115E-5Q,
632 4.179246497380453022046357404266022870788E-3Q,
633 1.136306384261959483095442402929502368598E-1Q,
634 1.422640343719842213484515445393284072830E0Q,
635 8.968786703393158374728850922289204805764E0Q,
636 2.914542473339246127533384118781216495934E1Q,
637 4.781605421020380669870197378210457054685E1Q,
638 3.693865837171883152382820584714795072937E1Q,
639 1.153220502744204904763115556224395893076E1Q,
640 /* 1.000000000000000000000000000000000000000E0 */
641 };
642
643
644 /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
645
646 static __float128
neval(__float128 x,const __float128 * p,int n)647 neval (__float128 x, const __float128 *p, int n)
648 {
649 __float128 y;
650
651 p += n;
652 y = *p--;
653 do
654 {
655 y = y * x + *p--;
656 }
657 while (--n > 0);
658 return y;
659 }
660
661
662 /* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */
663
664 static __float128
deval(__float128 x,const __float128 * p,int n)665 deval (__float128 x, const __float128 *p, int n)
666 {
667 __float128 y;
668
669 p += n;
670 y = x + *p--;
671 do
672 {
673 y = y * x + *p--;
674 }
675 while (--n > 0);
676 return y;
677 }
678
679
680 /* Bessel function of the first kind, order one. */
681
682 __float128
j1q(__float128 x)683 j1q (__float128 x)
684 {
685 __float128 xx, xinv, z, p, q, c, s, cc, ss;
686
687 if (! finiteq (x))
688 {
689 if (x != x)
690 return x + x;
691 else
692 return 0;
693 }
694 if (x == 0)
695 return x;
696 xx = fabsq (x);
697 if (xx <= 0x1p-58Q)
698 {
699 __float128 ret = x * 0.5Q;
700 math_check_force_underflow (ret);
701 if (ret == 0)
702 errno = ERANGE;
703 return ret;
704 }
705 if (xx <= 2)
706 {
707 /* 0 <= x <= 2 */
708 z = xx * xx;
709 p = xx * z * neval (z, J0_2N, NJ0_2N) / deval (z, J0_2D, NJ0_2D);
710 p += 0.5Q * xx;
711 if (x < 0)
712 p = -p;
713 return p;
714 }
715
716 /* X = x - 3 pi/4
717 cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
718 = 1/sqrt(2) * (-cos(x) + sin(x))
719 sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
720 = -1/sqrt(2) * (sin(x) + cos(x))
721 cf. Fdlibm. */
722 sincosq (xx, &s, &c);
723 ss = -s - c;
724 cc = s - c;
725 if (xx <= FLT128_MAX / 2)
726 {
727 z = cosq (xx + xx);
728 if ((s * c) > 0)
729 cc = z / ss;
730 else
731 ss = z / cc;
732 }
733
734 if (xx > 0x1p256Q)
735 {
736 z = ONEOSQPI * cc / sqrtq (xx);
737 if (x < 0)
738 z = -z;
739 return z;
740 }
741
742 xinv = 1 / xx;
743 z = xinv * xinv;
744 if (xinv <= 0.25)
745 {
746 if (xinv <= 0.125)
747 {
748 if (xinv <= 0.0625)
749 {
750 p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID);
751 q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID);
752 }
753 else
754 {
755 p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D);
756 q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D);
757 }
758 }
759 else if (xinv <= 0.1875)
760 {
761 p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D);
762 q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D);
763 }
764 else
765 {
766 p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D);
767 q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D);
768 }
769 } /* .25 */
770 else /* if (xinv <= 0.5) */
771 {
772 if (xinv <= 0.375)
773 {
774 if (xinv <= 0.3125)
775 {
776 p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D);
777 q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D);
778 }
779 else
780 {
781 p = neval (z, P2r7_3r2N, NP2r7_3r2N)
782 / deval (z, P2r7_3r2D, NP2r7_3r2D);
783 q = neval (z, Q2r7_3r2N, NQ2r7_3r2N)
784 / deval (z, Q2r7_3r2D, NQ2r7_3r2D);
785 }
786 }
787 else if (xinv <= 0.4375)
788 {
789 p = neval (z, P2r3_2r7N, NP2r3_2r7N)
790 / deval (z, P2r3_2r7D, NP2r3_2r7D);
791 q = neval (z, Q2r3_2r7N, NQ2r3_2r7N)
792 / deval (z, Q2r3_2r7D, NQ2r3_2r7D);
793 }
794 else
795 {
796 p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
797 q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
798 }
799 }
800 p = 1 + z * p;
801 q = z * q;
802 q = q * xinv + 0.375Q * xinv;
803 z = ONEOSQPI * (p * cc - q * ss) / sqrtq (xx);
804 if (x < 0)
805 z = -z;
806 return z;
807 }
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
815 static const __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
826 static const __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
y1q(__float128 x)842 y1q (__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)
849 {
850 if (x < 0)
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)
863 {
864 /* 0 <= x <= 2 */
865 SET_RESTORE_ROUNDF128 (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)
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 / 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 + 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