split checksum into cfg checksum and line checksum
[official-gcc.git] / libquadmath / math / j1q.c
blobf6bf2a256ce97e4d641ab012f6023ab802e27654
1 /* j1l.c
3 * Bessel function of order one
7 * SYNOPSIS:
9 * long double x, y, j1l();
11 * y = j1l( x );
15 * DESCRIPTION:
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
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,
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)).
38 * ACCURACY:
40 * Absolute error:
41 * arithmetic domain # trials peak rms
42 * IEEE 0, 30 100000 2.8e-34 2.7e-35
47 /* y1l.c
49 * Bessel function of the second kind, order one
53 * SYNOPSIS:
55 * double x, y, y1l();
57 * y = y1l( x );
61 * DESCRIPTION:
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)),
71 * X = x - 3 pi / 4.
73 * ACCURACY:
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 */
98 #include "quadmath-imp.h"
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.0Q;
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
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.000000000000000000000000000000000000000E0Q */
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,
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 */
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,
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 */
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
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 */
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,
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 */
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,
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 */
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,
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 */
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,
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 */
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,
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 */
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,
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 */
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,
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 */
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,
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 */
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,
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 */
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,
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 */
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,
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 */
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,
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 */
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,
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 */
644 /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
646 static __float128
647 neval (__float128 x, const __float128 *p, int n)
649 __float128 y;
651 p += n;
652 y = *p--;
655 y = y * x + *p--;
657 while (--n > 0);
658 return y;
662 /* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */
664 static __float128
665 deval (__float128 x, const __float128 *p, int n)
667 __float128 y;
669 p += n;
670 y = x + *p--;
673 y = y * x + *p--;
675 while (--n > 0);
676 return y;
680 /* Bessel function of the first kind, order one. */
682 __float128
683 j1q (__float128 x)
685 __float128 xx, xinv, z, p, q, c, s, cc, ss;
687 if (! finiteq (x))
689 if (x != x)
690 return x;
691 else
692 return 0.0Q;
694 if (x == 0.0Q)
695 return x;
696 xx = fabsq (x);
697 if (xx <= 2.0Q)
699 /* 0 <= x <= 2 */
700 z = xx * xx;
701 p = xx * z * neval (z, J0_2N, NJ0_2N) / deval (z, J0_2D, NJ0_2D);
702 p += 0.5Q * xx;
703 if (x < 0)
704 p = -p;
705 return p;
708 xinv = 1.0Q / xx;
709 z = xinv * xinv;
710 if (xinv <= 0.25)
712 if (xinv <= 0.125)
714 if (xinv <= 0.0625)
716 p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID);
717 q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID);
719 else
721 p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D);
722 q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D);
725 else if (xinv <= 0.1875)
727 p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D);
728 q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D);
730 else
732 p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D);
733 q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D);
735 } /* .25 */
736 else /* if (xinv <= 0.5) */
738 if (xinv <= 0.375)
740 if (xinv <= 0.3125)
742 p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D);
743 q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D);
745 else
747 p = neval (z, P2r7_3r2N, NP2r7_3r2N)
748 / deval (z, P2r7_3r2D, NP2r7_3r2D);
749 q = neval (z, Q2r7_3r2N, NQ2r7_3r2N)
750 / deval (z, Q2r7_3r2D, NQ2r7_3r2D);
753 else if (xinv <= 0.4375)
755 p = neval (z, P2r3_2r7N, NP2r3_2r7N)
756 / deval (z, P2r3_2r7D, NP2r3_2r7D);
757 q = neval (z, Q2r3_2r7N, NQ2r3_2r7N)
758 / deval (z, Q2r3_2r7D, NQ2r3_2r7D);
760 else
762 p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
763 q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
766 p = 1.0Q + z * p;
767 q = z * q;
768 q = q * xinv + 0.375Q * xinv;
769 /* X = x - 3 pi/4
770 cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
771 = 1/sqrt(2) * (-cos(x) + sin(x))
772 sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
773 = -1/sqrt(2) * (sin(x) + cos(x))
774 cf. Fdlibm. */
775 sincosq (xx, &s, &c);
776 ss = -s - c;
777 cc = s - c;
778 z = cosq (xx + xx);
779 if ((s * c) > 0)
780 cc = z / ss;
781 else
782 ss = z / cc;
783 z = ONEOSQPI * (p * cc - q * ss) / sqrtq (xx);
784 if (x < 0)
785 z = -z;
786 return z;
790 /* Y1(x) = 2/pi * (log(x) * J1(x) - 1/x) + x R(x^2)
791 Peak relative error 6.2e-38
792 0 <= x <= 2 */
793 #define NY0_2N 7
794 static __float128 Y0_2N[NY0_2N + 1] = {
795 -6.804415404830253804408698161694720833249E19Q,
796 1.805450517967019908027153056150465849237E19Q,
797 -8.065747497063694098810419456383006737312E17Q,
798 1.401336667383028259295830955439028236299E16Q,
799 -1.171654432898137585000399489686629680230E14Q,
800 5.061267920943853732895341125243428129150E11Q,
801 -1.096677850566094204586208610960870217970E9Q,
802 9.541172044989995856117187515882879304461E5Q,
804 #define NY0_2D 7
805 static __float128 Y0_2D[NY0_2D + 1] = {
806 3.470629591820267059538637461549677594549E20Q,
807 4.120796439009916326855848107545425217219E18Q,
808 2.477653371652018249749350657387030814542E16Q,
809 9.954678543353888958177169349272167762797E13Q,
810 2.957927997613630118216218290262851197754E11Q,
811 6.748421382188864486018861197614025972118E8Q,
812 1.173453425218010888004562071020305709319E6Q,
813 1.450335662961034949894009554536003377187E3Q,
814 /* 1.000000000000000000000000000000000000000E0 */
818 /* Bessel function of the second kind, order one. */
820 __float128
821 y1q (__float128 x)
823 __float128 xx, xinv, z, p, q, c, s, cc, ss;
825 if (! finiteq (x))
827 if (x != x)
828 return x;
829 else
830 return 0.0Q;
832 if (x <= 0.0Q)
834 if (x < 0.0Q)
835 return (zero / (zero * x));
836 return -HUGE_VALQ + x;
838 xx = fabsq (x);
839 if (xx <= 2.0Q)
841 /* 0 <= x <= 2 */
842 z = xx * xx;
843 p = xx * neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D);
844 p = -TWOOPI / xx + p;
845 p = TWOOPI * logq (x) * j1q (x) + p;
846 return p;
849 xinv = 1.0Q / xx;
850 z = xinv * xinv;
851 if (xinv <= 0.25)
853 if (xinv <= 0.125)
855 if (xinv <= 0.0625)
857 p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID);
858 q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID);
860 else
862 p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D);
863 q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D);
866 else if (xinv <= 0.1875)
868 p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D);
869 q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D);
871 else
873 p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D);
874 q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D);
876 } /* .25 */
877 else /* if (xinv <= 0.5) */
879 if (xinv <= 0.375)
881 if (xinv <= 0.3125)
883 p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D);
884 q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D);
886 else
888 p = neval (z, P2r7_3r2N, NP2r7_3r2N)
889 / deval (z, P2r7_3r2D, NP2r7_3r2D);
890 q = neval (z, Q2r7_3r2N, NQ2r7_3r2N)
891 / deval (z, Q2r7_3r2D, NQ2r7_3r2D);
894 else if (xinv <= 0.4375)
896 p = neval (z, P2r3_2r7N, NP2r3_2r7N)
897 / deval (z, P2r3_2r7D, NP2r3_2r7D);
898 q = neval (z, Q2r3_2r7N, NQ2r3_2r7N)
899 / deval (z, Q2r3_2r7D, NQ2r3_2r7D);
901 else
903 p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
904 q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
907 p = 1.0Q + z * p;
908 q = z * q;
909 q = q * xinv + 0.375Q * xinv;
910 /* X = x - 3 pi/4
911 cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
912 = 1/sqrt(2) * (-cos(x) + sin(x))
913 sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
914 = -1/sqrt(2) * (sin(x) + cos(x))
915 cf. Fdlibm. */
916 sincosq (xx, &s, &c);
917 ss = -s - c;
918 cc = s - c;
919 z = cosq (xx + xx);
920 if ((s * c) > 0)
921 cc = z / ss;
922 else
923 ss = z / cc;
924 z = ONEOSQPI * (p * ss + q * cc) / sqrtq (xx);
925 return z;