Mercurial > hg > CbC > CbC_gcc

comparison libquadmath/math/j1q.c @ 68:561a7518be6b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
update gcc-4.6
author Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
date Sun, 21 Aug 2011 07:07:55 +0900
parents
children 04ced10e8804
comparison
equal deleted inserted replaced
67:f6334be47118 68:561a7518be6b
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, write to the Free Software
96 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
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.0Q;
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.000000000000000000000000000000000000000E0Q */
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
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
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
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;
691 else
692 return 0.0Q;
693 }
694 if (x == 0.0Q)
695 return x;
696 xx = fabsq (x);
697 if (xx <= 2.0Q)
698 {
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;
706 }
707
708 xinv = 1.0Q / xx;
709 z = xinv * xinv;
710 if (xinv <= 0.25)
711 {
712 if (xinv <= 0.125)
713 {
714 if (xinv <= 0.0625)
715 {
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);
718 }
719 else
720 {
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);
723 }
724 }
725 else if (xinv <= 0.1875)
726 {
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);
729 }
730 else
731 {
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);
734 }
735 } /* .25 */
736 else /* if (xinv <= 0.5) */
737 {
738 if (xinv <= 0.375)
739 {
740 if (xinv <= 0.3125)
741 {
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);
744 }
745 else
746 {
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);
751 }
752 }
753 else if (xinv <= 0.4375)
754 {
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);
759 }
760 else
761 {
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);
764 }
765 }
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;
787 }
788
789
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,
803 };
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 */
815 };
816
817
818 /* Bessel function of the second kind, order one. */
819
820 __float128
821 y1q (__float128 x)
822 {
823 __float128 xx, xinv, z, p, q, c, s, cc, ss;
824
825 if (! finiteq (x))
826 {
827 if (x != x)
828 return x;
829 else
830 return 0.0Q;
831 }
832 if (x <= 0.0Q)
833 {
834 if (x < 0.0Q)
835 return (zero / (zero * x));
836 return -HUGE_VALQ + x;
837 }
838 xx = fabsq (x);
839 if (xx <= 2.0Q)
840 {
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;
847 }
848
849 xinv = 1.0Q / xx;
850 z = xinv * xinv;
851 if (xinv <= 0.25)
852 {
853 if (xinv <= 0.125)
854 {
855 if (xinv <= 0.0625)
856 {
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);
859 }
860 else
861 {
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);
864 }
865 }
866 else if (xinv <= 0.1875)
867 {
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);
870 }
871 else
872 {
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);
875 }
876 } /* .25 */
877 else /* if (xinv <= 0.5) */
878 {
879 if (xinv <= 0.375)
880 {
881 if (xinv <= 0.3125)
882 {
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);
885 }
886 else
887 {
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);
892 }
893 }
894 else if (xinv <= 0.4375)
895 {
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);
900 }
901 else
902 {
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);
905 }
906 }
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;
926 }