1 |
740 |
jeremybenn |
/* 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 |
|
|
|
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 |
|
|
|
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 |
|
|
|
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 |
|
|
|
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 |
|
|
}
|