]>
Commit | Line | Data |
---|---|---|
f524e4a2 AJ |
1 | /* j0l.c |
2 | * | |
3 | * Bessel function of order zero | |
4 | * | |
5 | * | |
6 | * | |
7 | * SYNOPSIS: | |
8 | * | |
9 | * long double x, y, j0l(); | |
10 | * | |
11 | * y = j0l( x ); | |
12 | * | |
13 | * | |
14 | * | |
15 | * DESCRIPTION: | |
16 | * | |
17 | * Returns Bessel function of first kind, order zero 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 | |
21 | * is J0(x) = 1 - x^2 / 4 + x^4 R(x^2) | |
22 | * The second interval is further partitioned into eight equal segments | |
23 | * of 1/x. | |
24 | * | |
25 | * J0(x) = sqrt(2/(pi x)) (P0(x) cos(X) - Q0(x) sin(X)), | |
26 | * X = x - pi/4, | |
27 | * | |
28 | * and the auxiliary functions are given by | |
29 | * | |
30 | * J0(x)cos(X) + Y0(x)sin(X) = sqrt( 2/(pi x)) P0(x), | |
31 | * P0(x) = 1 + 1/x^2 R(1/x^2) | |
32 | * | |
33 | * Y0(x)cos(X) - J0(x)sin(X) = sqrt( 2/(pi x)) Q0(x), | |
34 | * Q0(x) = 1/x (-.125 + 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 1.7e-34 2.4e-35 | |
43 | * | |
44 | * | |
45 | */ | |
46 | ||
47 | /* y0l.c | |
48 | * | |
49 | * Bessel function of the second kind, order zero | |
50 | * | |
51 | * | |
52 | * | |
53 | * SYNOPSIS: | |
54 | * | |
55 | * double x, y, y0l(); | |
56 | * | |
57 | * y = y0l( x ); | |
58 | * | |
59 | * | |
60 | * | |
61 | * DESCRIPTION: | |
62 | * | |
63 | * Returns Bessel function of the second kind, of order | |
64 | * zero, of the argument. | |
65 | * | |
66 | * The approximation is the same as for J0(x), and | |
67 | * Y0(x) = sqrt(2/(pi x)) (P0(x) sin(X) + Q0(x) cos(X)). | |
68 | * | |
69 | * ACCURACY: | |
70 | * | |
71 | * Absolute error, when y0(x) < 1; else relative error: | |
72 | * | |
73 | * arithmetic domain # trials peak rms | |
74 | * IEEE 0, 30 100000 3.0e-34 2.7e-35 | |
75 | * | |
76 | */ | |
77 | ||
cc7375ce RM |
78 | /* Copyright 2001 by Stephen L. Moshier (moshier@na-net.ornl.gov). |
79 | ||
80 | This library is free software; you can redistribute it and/or | |
81 | modify it under the terms of the GNU Lesser General Public | |
82 | License as published by the Free Software Foundation; either | |
83 | version 2.1 of the License, or (at your option) any later version. | |
84 | ||
85 | This library is distributed in the hope that it will be useful, | |
86 | but WITHOUT ANY WARRANTY; without even the implied warranty of | |
87 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
88 | Lesser General Public License for more details. | |
89 | ||
90 | You should have received a copy of the GNU Lesser General Public | |
59ba27a6 PE |
91 | License along with this library; if not, see |
92 | <http://www.gnu.org/licenses/>. */ | |
f524e4a2 | 93 | |
1ed0291c RH |
94 | #include <math.h> |
95 | #include <math_private.h> | |
98c48fe5 | 96 | #include <float.h> |
f524e4a2 AJ |
97 | |
98 | /* 1 / sqrt(pi) */ | |
1f5649f8 | 99 | static const long double ONEOSQPI = 5.6418958354775628694807945156077258584405E-1L; |
f524e4a2 | 100 | /* 2 / pi */ |
1f5649f8 UD |
101 | static const long double TWOOPI = 6.3661977236758134307553505349005744813784E-1L; |
102 | static const long double zero = 0.0L; | |
f524e4a2 AJ |
103 | |
104 | /* J0(x) = 1 - x^2/4 + x^2 x^2 R(x^2) | |
105 | Peak relative error 3.4e-37 | |
106 | 0 <= x <= 2 */ | |
107 | #define NJ0_2N 6 | |
1f5649f8 | 108 | static const long double J0_2N[NJ0_2N + 1] = { |
f524e4a2 AJ |
109 | 3.133239376997663645548490085151484674892E16L, |
110 | -5.479944965767990821079467311839107722107E14L, | |
111 | 6.290828903904724265980249871997551894090E12L, | |
112 | -3.633750176832769659849028554429106299915E10L, | |
113 | 1.207743757532429576399485415069244807022E8L, | |
114 | -2.107485999925074577174305650549367415465E5L, | |
115 | 1.562826808020631846245296572935547005859E2L, | |
116 | }; | |
117 | #define NJ0_2D 6 | |
1f5649f8 | 118 | static const long double J0_2D[NJ0_2D + 1] = { |
f524e4a2 AJ |
119 | 2.005273201278504733151033654496928968261E18L, |
120 | 2.063038558793221244373123294054149790864E16L, | |
121 | 1.053350447931127971406896594022010524994E14L, | |
122 | 3.496556557558702583143527876385508882310E11L, | |
123 | 8.249114511878616075860654484367133976306E8L, | |
124 | 1.402965782449571800199759247964242790589E6L, | |
125 | 1.619910762853439600957801751815074787351E3L, | |
126 | /* 1.000000000000000000000000000000000000000E0 */ | |
127 | }; | |
128 | ||
129 | /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2), | |
130 | 0 <= 1/x <= .0625 | |
131 | Peak relative error 3.3e-36 */ | |
132 | #define NP16_IN 9 | |
1f5649f8 | 133 | static const long double P16_IN[NP16_IN + 1] = { |
f524e4a2 AJ |
134 | -1.901689868258117463979611259731176301065E-16L, |
135 | -1.798743043824071514483008340803573980931E-13L, | |
136 | -6.481746687115262291873324132944647438959E-11L, | |
137 | -1.150651553745409037257197798528294248012E-8L, | |
138 | -1.088408467297401082271185599507222695995E-6L, | |
139 | -5.551996725183495852661022587879817546508E-5L, | |
140 | -1.477286941214245433866838787454880214736E-3L, | |
141 | -1.882877976157714592017345347609200402472E-2L, | |
142 | -9.620983176855405325086530374317855880515E-2L, | |
143 | -1.271468546258855781530458854476627766233E-1L, | |
144 | }; | |
145 | #define NP16_ID 9 | |
1f5649f8 | 146 | static const long double P16_ID[NP16_ID + 1] = { |
f524e4a2 AJ |
147 | 2.704625590411544837659891569420764475007E-15L, |
148 | 2.562526347676857624104306349421985403573E-12L, | |
149 | 9.259137589952741054108665570122085036246E-10L, | |
150 | 1.651044705794378365237454962653430805272E-7L, | |
151 | 1.573561544138733044977714063100859136660E-5L, | |
152 | 8.134482112334882274688298469629884804056E-4L, | |
153 | 2.219259239404080863919375103673593571689E-2L, | |
154 | 2.976990606226596289580242451096393862792E-1L, | |
155 | 1.713895630454693931742734911930937246254E0L, | |
156 | 3.231552290717904041465898249160757368855E0L, | |
157 | /* 1.000000000000000000000000000000000000000E0 */ | |
158 | }; | |
159 | ||
160 | /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2) | |
161 | 0.0625 <= 1/x <= 0.125 | |
162 | Peak relative error 2.4e-35 */ | |
163 | #define NP8_16N 10 | |
1f5649f8 | 164 | static const long double P8_16N[NP8_16N + 1] = { |
f524e4a2 AJ |
165 | -2.335166846111159458466553806683579003632E-15L, |
166 | -1.382763674252402720401020004169367089975E-12L, | |
167 | -3.192160804534716696058987967592784857907E-10L, | |
168 | -3.744199606283752333686144670572632116899E-8L, | |
169 | -2.439161236879511162078619292571922772224E-6L, | |
170 | -9.068436986859420951664151060267045346549E-5L, | |
171 | -1.905407090637058116299757292660002697359E-3L, | |
172 | -2.164456143936718388053842376884252978872E-2L, | |
173 | -1.212178415116411222341491717748696499966E-1L, | |
174 | -2.782433626588541494473277445959593334494E-1L, | |
175 | -1.670703190068873186016102289227646035035E-1L, | |
176 | }; | |
177 | #define NP8_16D 10 | |
1f5649f8 | 178 | static const long double P8_16D[NP8_16D + 1] = { |
f524e4a2 AJ |
179 | 3.321126181135871232648331450082662856743E-14L, |
180 | 1.971894594837650840586859228510007703641E-11L, | |
181 | 4.571144364787008285981633719513897281690E-9L, | |
182 | 5.396419143536287457142904742849052402103E-7L, | |
183 | 3.551548222385845912370226756036899901549E-5L, | |
184 | 1.342353874566932014705609788054598013516E-3L, | |
185 | 2.899133293006771317589357444614157734385E-2L, | |
186 | 3.455374978185770197704507681491574261545E-1L, | |
187 | 2.116616964297512311314454834712634820514E0L, | |
188 | 5.850768316827915470087758636881584174432E0L, | |
189 | 5.655273858938766830855753983631132928968E0L, | |
190 | /* 1.000000000000000000000000000000000000000E0 */ | |
191 | }; | |
192 | ||
193 | /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2) | |
194 | 0.125 <= 1/x <= 0.1875 | |
195 | Peak relative error 2.7e-35 */ | |
196 | #define NP5_8N 10 | |
1f5649f8 | 197 | static const long double P5_8N[NP5_8N + 1] = { |
f524e4a2 AJ |
198 | -1.270478335089770355749591358934012019596E-12L, |
199 | -4.007588712145412921057254992155810347245E-10L, | |
200 | -4.815187822989597568124520080486652009281E-8L, | |
201 | -2.867070063972764880024598300408284868021E-6L, | |
202 | -9.218742195161302204046454768106063638006E-5L, | |
203 | -1.635746821447052827526320629828043529997E-3L, | |
204 | -1.570376886640308408247709616497261011707E-2L, | |
205 | -7.656484795303305596941813361786219477807E-2L, | |
206 | -1.659371030767513274944805479908858628053E-1L, | |
207 | -1.185340550030955660015841796219919804915E-1L, | |
208 | -8.920026499909994671248893388013790366712E-3L, | |
209 | }; | |
210 | #define NP5_8D 9 | |
1f5649f8 | 211 | static const long double P5_8D[NP5_8D + 1] = { |
f524e4a2 AJ |
212 | 1.806902521016705225778045904631543990314E-11L, |
213 | 5.728502760243502431663549179135868966031E-9L, | |
214 | 6.938168504826004255287618819550667978450E-7L, | |
215 | 4.183769964807453250763325026573037785902E-5L, | |
216 | 1.372660678476925468014882230851637878587E-3L, | |
217 | 2.516452105242920335873286419212708961771E-2L, | |
218 | 2.550502712902647803796267951846557316182E-1L, | |
219 | 1.365861559418983216913629123778747617072E0L, | |
220 | 3.523825618308783966723472468855042541407E0L, | |
221 | 3.656365803506136165615111349150536282434E0L, | |
222 | /* 1.000000000000000000000000000000000000000E0 */ | |
223 | }; | |
224 | ||
225 | /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2) | |
226 | Peak relative error 3.5e-35 | |
227 | 0.1875 <= 1/x <= 0.25 */ | |
228 | #define NP4_5N 9 | |
1f5649f8 | 229 | static const long double P4_5N[NP4_5N + 1] = { |
f524e4a2 AJ |
230 | -9.791405771694098960254468859195175708252E-10L, |
231 | -1.917193059944531970421626610188102836352E-7L, | |
232 | -1.393597539508855262243816152893982002084E-5L, | |
233 | -4.881863490846771259880606911667479860077E-4L, | |
234 | -8.946571245022470127331892085881699269853E-3L, | |
235 | -8.707474232568097513415336886103899434251E-2L, | |
236 | -4.362042697474650737898551272505525973766E-1L, | |
237 | -1.032712171267523975431451359962375617386E0L, | |
238 | -9.630502683169895107062182070514713702346E-1L, | |
239 | -2.251804386252969656586810309252357233320E-1L, | |
240 | }; | |
241 | #define NP4_5D 9 | |
1f5649f8 | 242 | static const long double P4_5D[NP4_5D + 1] = { |
f524e4a2 AJ |
243 | 1.392555487577717669739688337895791213139E-8L, |
244 | 2.748886559120659027172816051276451376854E-6L, | |
245 | 2.024717710644378047477189849678576659290E-4L, | |
246 | 7.244868609350416002930624752604670292469E-3L, | |
247 | 1.373631762292244371102989739300382152416E-1L, | |
248 | 1.412298581400224267910294815260613240668E0L, | |
249 | 7.742495637843445079276397723849017617210E0L, | |
250 | 2.138429269198406512028307045259503811861E1L, | |
251 | 2.651547684548423476506826951831712762610E1L, | |
252 | 1.167499382465291931571685222882909166935E1L, | |
253 | /* 1.000000000000000000000000000000000000000E0 */ | |
254 | }; | |
255 | ||
256 | /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2) | |
257 | Peak relative error 2.3e-36 | |
258 | 0.25 <= 1/x <= 0.3125 */ | |
259 | #define NP3r2_4N 9 | |
1f5649f8 | 260 | static const long double P3r2_4N[NP3r2_4N + 1] = { |
f524e4a2 AJ |
261 | -2.589155123706348361249809342508270121788E-8L, |
262 | -3.746254369796115441118148490849195516593E-6L, | |
263 | -1.985595497390808544622893738135529701062E-4L, | |
264 | -5.008253705202932091290132760394976551426E-3L, | |
265 | -6.529469780539591572179155511840853077232E-2L, | |
266 | -4.468736064761814602927408833818990271514E-1L, | |
267 | -1.556391252586395038089729428444444823380E0L, | |
268 | -2.533135309840530224072920725976994981638E0L, | |
269 | -1.605509621731068453869408718565392869560E0L, | |
270 | -2.518966692256192789269859830255724429375E-1L, | |
271 | }; | |
272 | #define NP3r2_4D 9 | |
1f5649f8 | 273 | static const long double P3r2_4D[NP3r2_4D + 1] = { |
f524e4a2 AJ |
274 | 3.682353957237979993646169732962573930237E-7L, |
275 | 5.386741661883067824698973455566332102029E-5L, | |
276 | 2.906881154171822780345134853794241037053E-3L, | |
277 | 7.545832595801289519475806339863492074126E-2L, | |
278 | 1.029405357245594877344360389469584526654E0L, | |
279 | 7.565706120589873131187989560509757626725E0L, | |
280 | 2.951172890699569545357692207898667665796E1L, | |
281 | 5.785723537170311456298467310529815457536E1L, | |
282 | 5.095621464598267889126015412522773474467E1L, | |
283 | 1.602958484169953109437547474953308401442E1L, | |
284 | /* 1.000000000000000000000000000000000000000E0 */ | |
285 | }; | |
286 | ||
287 | /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2) | |
288 | Peak relative error 1.0e-35 | |
289 | 0.3125 <= 1/x <= 0.375 */ | |
290 | #define NP2r7_3r2N 9 | |
1f5649f8 | 291 | static const long double P2r7_3r2N[NP2r7_3r2N + 1] = { |
f524e4a2 AJ |
292 | -1.917322340814391131073820537027234322550E-7L, |
293 | -1.966595744473227183846019639723259011906E-5L, | |
294 | -7.177081163619679403212623526632690465290E-4L, | |
295 | -1.206467373860974695661544653741899755695E-2L, | |
296 | -1.008656452188539812154551482286328107316E-1L, | |
297 | -4.216016116408810856620947307438823892707E-1L, | |
298 | -8.378631013025721741744285026537009814161E-1L, | |
299 | -6.973895635309960850033762745957946272579E-1L, | |
300 | -1.797864718878320770670740413285763554812E-1L, | |
301 | -4.098025357743657347681137871388402849581E-3L, | |
302 | }; | |
303 | #define NP2r7_3r2D 8 | |
1f5649f8 | 304 | static const long double P2r7_3r2D[NP2r7_3r2D + 1] = { |
f524e4a2 AJ |
305 | 2.726858489303036441686496086962545034018E-6L, |
306 | 2.840430827557109238386808968234848081424E-4L, | |
307 | 1.063826772041781947891481054529454088832E-2L, | |
308 | 1.864775537138364773178044431045514405468E-1L, | |
309 | 1.665660052857205170440952607701728254211E0L, | |
310 | 7.723745889544331153080842168958348568395E0L, | |
311 | 1.810726427571829798856428548102077799835E1L, | |
312 | 1.986460672157794440666187503833545388527E1L, | |
313 | 8.645503204552282306364296517220055815488E0L, | |
314 | /* 1.000000000000000000000000000000000000000E0 */ | |
315 | }; | |
316 | ||
317 | /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2) | |
318 | Peak relative error 1.3e-36 | |
319 | 0.3125 <= 1/x <= 0.4375 */ | |
320 | #define NP2r3_2r7N 9 | |
1f5649f8 | 321 | static const long double P2r3_2r7N[NP2r3_2r7N + 1] = { |
f524e4a2 AJ |
322 | -1.594642785584856746358609622003310312622E-6L, |
323 | -1.323238196302221554194031733595194539794E-4L, | |
324 | -3.856087818696874802689922536987100372345E-3L, | |
325 | -5.113241710697777193011470733601522047399E-2L, | |
326 | -3.334229537209911914449990372942022350558E-1L, | |
327 | -1.075703518198127096179198549659283422832E0L, | |
328 | -1.634174803414062725476343124267110981807E0L, | |
329 | -1.030133247434119595616826842367268304880E0L, | |
330 | -1.989811539080358501229347481000707289391E-1L, | |
331 | -3.246859189246653459359775001466924610236E-3L, | |
332 | }; | |
333 | #define NP2r3_2r7D 8 | |
1f5649f8 | 334 | static const long double P2r3_2r7D[NP2r3_2r7D + 1] = { |
f524e4a2 AJ |
335 | 2.267936634217251403663034189684284173018E-5L, |
336 | 1.918112982168673386858072491437971732237E-3L, | |
337 | 5.771704085468423159125856786653868219522E-2L, | |
338 | 8.056124451167969333717642810661498890507E-1L, | |
339 | 5.687897967531010276788680634413789328776E0L, | |
340 | 2.072596760717695491085444438270778394421E1L, | |
341 | 3.801722099819929988585197088613160496684E1L, | |
342 | 3.254620235902912339534998592085115836829E1L, | |
343 | 1.104847772130720331801884344645060675036E1L, | |
344 | /* 1.000000000000000000000000000000000000000E0 */ | |
345 | }; | |
346 | ||
347 | /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2) | |
348 | Peak relative error 1.2e-35 | |
349 | 0.4375 <= 1/x <= 0.5 */ | |
350 | #define NP2_2r3N 8 | |
1f5649f8 | 351 | static const long double P2_2r3N[NP2_2r3N + 1] = { |
f524e4a2 AJ |
352 | -1.001042324337684297465071506097365389123E-4L, |
353 | -6.289034524673365824853547252689991418981E-3L, | |
354 | -1.346527918018624234373664526930736205806E-1L, | |
355 | -1.268808313614288355444506172560463315102E0L, | |
356 | -5.654126123607146048354132115649177406163E0L, | |
357 | -1.186649511267312652171775803270911971693E1L, | |
358 | -1.094032424931998612551588246779200724257E1L, | |
359 | -3.728792136814520055025256353193674625267E0L, | |
360 | -3.000348318524471807839934764596331810608E-1L, | |
361 | }; | |
362 | #define NP2_2r3D 8 | |
1f5649f8 | 363 | static const long double P2_2r3D[NP2_2r3D + 1] = { |
f524e4a2 AJ |
364 | 1.423705538269770974803901422532055612980E-3L, |
365 | 9.171476630091439978533535167485230575894E-2L, | |
366 | 2.049776318166637248868444600215942828537E0L, | |
367 | 2.068970329743769804547326701946144899583E1L, | |
368 | 1.025103500560831035592731539565060347709E2L, | |
369 | 2.528088049697570728252145557167066708284E2L, | |
370 | 2.992160327587558573740271294804830114205E2L, | |
371 | 1.540193761146551025832707739468679973036E2L, | |
372 | 2.779516701986912132637672140709452502650E1L, | |
373 | /* 1.000000000000000000000000000000000000000E0 */ | |
374 | }; | |
375 | ||
376 | /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x), | |
377 | Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2)) | |
378 | Peak relative error 2.2e-35 | |
379 | 0 <= 1/x <= .0625 */ | |
380 | #define NQ16_IN 10 | |
1f5649f8 | 381 | static const long double Q16_IN[NQ16_IN + 1] = { |
f524e4a2 AJ |
382 | 2.343640834407975740545326632205999437469E-18L, |
383 | 2.667978112927811452221176781536278257448E-15L, | |
384 | 1.178415018484555397390098879501969116536E-12L, | |
385 | 2.622049767502719728905924701288614016597E-10L, | |
386 | 3.196908059607618864801313380896308968673E-8L, | |
387 | 2.179466154171673958770030655199434798494E-6L, | |
388 | 8.139959091628545225221976413795645177291E-5L, | |
389 | 1.563900725721039825236927137885747138654E-3L, | |
390 | 1.355172364265825167113562519307194840307E-2L, | |
391 | 3.928058355906967977269780046844768588532E-2L, | |
392 | 1.107891967702173292405380993183694932208E-2L, | |
393 | }; | |
394 | #define NQ16_ID 9 | |
1f5649f8 | 395 | static const long double Q16_ID[NQ16_ID + 1] = { |
f524e4a2 AJ |
396 | 3.199850952578356211091219295199301766718E-17L, |
397 | 3.652601488020654842194486058637953363918E-14L, | |
398 | 1.620179741394865258354608590461839031281E-11L, | |
399 | 3.629359209474609630056463248923684371426E-9L, | |
400 | 4.473680923894354600193264347733477363305E-7L, | |
401 | 3.106368086644715743265603656011050476736E-5L, | |
402 | 1.198239259946770604954664925153424252622E-3L, | |
403 | 2.446041004004283102372887804475767568272E-2L, | |
404 | 2.403235525011860603014707768815113698768E-1L, | |
405 | 9.491006790682158612266270665136910927149E-1L, | |
406 | /* 1.000000000000000000000000000000000000000E0 */ | |
407 | }; | |
408 | ||
409 | /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x), | |
410 | Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2)) | |
411 | Peak relative error 5.1e-36 | |
412 | 0.0625 <= 1/x <= 0.125 */ | |
413 | #define NQ8_16N 11 | |
1f5649f8 | 414 | static const long double Q8_16N[NQ8_16N + 1] = { |
f524e4a2 AJ |
415 | 1.001954266485599464105669390693597125904E-17L, |
416 | 7.545499865295034556206475956620160007849E-15L, | |
417 | 2.267838684785673931024792538193202559922E-12L, | |
418 | 3.561909705814420373609574999542459912419E-10L, | |
419 | 3.216201422768092505214730633842924944671E-8L, | |
420 | 1.731194793857907454569364622452058554314E-6L, | |
421 | 5.576944613034537050396518509871004586039E-5L, | |
422 | 1.051787760316848982655967052985391418146E-3L, | |
423 | 1.102852974036687441600678598019883746959E-2L, | |
424 | 5.834647019292460494254225988766702933571E-2L, | |
425 | 1.290281921604364618912425380717127576529E-1L, | |
426 | 7.598886310387075708640370806458926458301E-2L, | |
427 | }; | |
428 | #define NQ8_16D 11 | |
1f5649f8 | 429 | static const long double Q8_16D[NQ8_16D + 1] = { |
f524e4a2 AJ |
430 | 1.368001558508338469503329967729951830843E-16L, |
431 | 1.034454121857542147020549303317348297289E-13L, | |
432 | 3.128109209247090744354764050629381674436E-11L, | |
433 | 4.957795214328501986562102573522064468671E-9L, | |
434 | 4.537872468606711261992676606899273588899E-7L, | |
435 | 2.493639207101727713192687060517509774182E-5L, | |
436 | 8.294957278145328349785532236663051405805E-4L, | |
437 | 1.646471258966713577374948205279380115839E-2L, | |
438 | 1.878910092770966718491814497982191447073E-1L, | |
439 | 1.152641605706170353727903052525652504075E0L, | |
440 | 3.383550240669773485412333679367792932235E0L, | |
441 | 3.823875252882035706910024716609908473970E0L, | |
442 | /* 1.000000000000000000000000000000000000000E0 */ | |
443 | }; | |
444 | ||
445 | /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x), | |
446 | Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2)) | |
447 | Peak relative error 3.9e-35 | |
448 | 0.125 <= 1/x <= 0.1875 */ | |
449 | #define NQ5_8N 10 | |
1f5649f8 | 450 | static const long double Q5_8N[NQ5_8N + 1] = { |
f524e4a2 AJ |
451 | 1.750399094021293722243426623211733898747E-13L, |
452 | 6.483426211748008735242909236490115050294E-11L, | |
453 | 9.279430665656575457141747875716899958373E-9L, | |
454 | 6.696634968526907231258534757736576340266E-7L, | |
455 | 2.666560823798895649685231292142838188061E-5L, | |
456 | 6.025087697259436271271562769707550594540E-4L, | |
457 | 7.652807734168613251901945778921336353485E-3L, | |
458 | 5.226269002589406461622551452343519078905E-2L, | |
459 | 1.748390159751117658969324896330142895079E-1L, | |
460 | 2.378188719097006494782174902213083589660E-1L, | |
461 | 8.383984859679804095463699702165659216831E-2L, | |
462 | }; | |
463 | #define NQ5_8D 10 | |
1f5649f8 | 464 | static const long double Q5_8D[NQ5_8D + 1] = { |
f524e4a2 AJ |
465 | 2.389878229704327939008104855942987615715E-12L, |
466 | 8.926142817142546018703814194987786425099E-10L, | |
467 | 1.294065862406745901206588525833274399038E-7L, | |
468 | 9.524139899457666250828752185212769682191E-6L, | |
469 | 3.908332488377770886091936221573123353489E-4L, | |
470 | 9.250427033957236609624199884089916836748E-3L, | |
471 | 1.263420066165922645975830877751588421451E-1L, | |
472 | 9.692527053860420229711317379861733180654E-1L, | |
473 | 3.937813834630430172221329298841520707954E0L, | |
474 | 7.603126427436356534498908111445191312181E0L, | |
475 | 5.670677653334105479259958485084550934305E0L, | |
476 | /* 1.000000000000000000000000000000000000000E0 */ | |
477 | }; | |
478 | ||
479 | /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x), | |
480 | Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2)) | |
481 | Peak relative error 3.2e-35 | |
482 | 0.1875 <= 1/x <= 0.25 */ | |
483 | #define NQ4_5N 10 | |
1f5649f8 | 484 | static const long double Q4_5N[NQ4_5N + 1] = { |
f524e4a2 AJ |
485 | 2.233870042925895644234072357400122854086E-11L, |
486 | 5.146223225761993222808463878999151699792E-9L, | |
487 | 4.459114531468296461688753521109797474523E-7L, | |
488 | 1.891397692931537975547242165291668056276E-5L, | |
489 | 4.279519145911541776938964806470674565504E-4L, | |
490 | 5.275239415656560634702073291768904783989E-3L, | |
491 | 3.468698403240744801278238473898432608887E-2L, | |
492 | 1.138773146337708415188856882915457888274E-1L, | |
493 | 1.622717518946443013587108598334636458955E-1L, | |
494 | 7.249040006390586123760992346453034628227E-2L, | |
495 | 1.941595365256460232175236758506411486667E-3L, | |
496 | }; | |
497 | #define NQ4_5D 9 | |
1f5649f8 | 498 | static const long double Q4_5D[NQ4_5D + 1] = { |
f524e4a2 AJ |
499 | 3.049977232266999249626430127217988047453E-10L, |
500 | 7.120883230531035857746096928889676144099E-8L, | |
501 | 6.301786064753734446784637919554359588859E-6L, | |
502 | 2.762010530095069598480766869426308077192E-4L, | |
503 | 6.572163250572867859316828886203406361251E-3L, | |
504 | 8.752566114841221958200215255461843397776E-2L, | |
505 | 6.487654992874805093499285311075289932664E-1L, | |
506 | 2.576550017826654579451615283022812801435E0L, | |
507 | 5.056392229924022835364779562707348096036E0L, | |
508 | 4.179770081068251464907531367859072157773E0L, | |
509 | /* 1.000000000000000000000000000000000000000E0 */ | |
510 | }; | |
511 | ||
512 | /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x), | |
513 | Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2)) | |
514 | Peak relative error 1.4e-36 | |
515 | 0.25 <= 1/x <= 0.3125 */ | |
516 | #define NQ3r2_4N 10 | |
1f5649f8 | 517 | static const long double Q3r2_4N[NQ3r2_4N + 1] = { |
f524e4a2 AJ |
518 | 6.126167301024815034423262653066023684411E-10L, |
519 | 1.043969327113173261820028225053598975128E-7L, | |
520 | 6.592927270288697027757438170153763220190E-6L, | |
521 | 2.009103660938497963095652951912071336730E-4L, | |
522 | 3.220543385492643525985862356352195896964E-3L, | |
523 | 2.774405975730545157543417650436941650990E-2L, | |
524 | 1.258114008023826384487378016636555041129E-1L, | |
525 | 2.811724258266902502344701449984698323860E-1L, | |
526 | 2.691837665193548059322831687432415014067E-1L, | |
527 | 7.949087384900985370683770525312735605034E-2L, | |
528 | 1.229509543620976530030153018986910810747E-3L, | |
529 | }; | |
530 | #define NQ3r2_4D 9 | |
1f5649f8 | 531 | static const long double Q3r2_4D[NQ3r2_4D + 1] = { |
f524e4a2 AJ |
532 | 8.364260446128475461539941389210166156568E-9L, |
533 | 1.451301850638956578622154585560759862764E-6L, | |
534 | 9.431830010924603664244578867057141839463E-5L, | |
535 | 3.004105101667433434196388593004526182741E-3L, | |
536 | 5.148157397848271739710011717102773780221E-2L, | |
537 | 4.901089301726939576055285374953887874895E-1L, | |
538 | 2.581760991981709901216967665934142240346E0L, | |
539 | 7.257105880775059281391729708630912791847E0L, | |
540 | 1.006014717326362868007913423810737369312E1L, | |
541 | 5.879416600465399514404064187445293212470E0L, | |
542 | /* 1.000000000000000000000000000000000000000E0*/ | |
543 | }; | |
544 | ||
545 | /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x), | |
546 | Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2)) | |
547 | Peak relative error 3.8e-36 | |
548 | 0.3125 <= 1/x <= 0.375 */ | |
549 | #define NQ2r7_3r2N 9 | |
1f5649f8 | 550 | static const long double Q2r7_3r2N[NQ2r7_3r2N + 1] = { |
f524e4a2 AJ |
551 | 7.584861620402450302063691901886141875454E-8L, |
552 | 9.300939338814216296064659459966041794591E-6L, | |
553 | 4.112108906197521696032158235392604947895E-4L, | |
554 | 8.515168851578898791897038357239630654431E-3L, | |
555 | 8.971286321017307400142720556749573229058E-2L, | |
556 | 4.885856732902956303343015636331874194498E-1L, | |
557 | 1.334506268733103291656253500506406045846E0L, | |
558 | 1.681207956863028164179042145803851824654E0L, | |
559 | 8.165042692571721959157677701625853772271E-1L, | |
560 | 9.805848115375053300608712721986235900715E-2L, | |
561 | }; | |
562 | #define NQ2r7_3r2D 9 | |
1f5649f8 | 563 | static const long double Q2r7_3r2D[NQ2r7_3r2D + 1] = { |
f524e4a2 AJ |
564 | 1.035586492113036586458163971239438078160E-6L, |
565 | 1.301999337731768381683593636500979713689E-4L, | |
566 | 5.993695702564527062553071126719088859654E-3L, | |
567 | 1.321184892887881883489141186815457808785E-1L, | |
568 | 1.528766555485015021144963194165165083312E0L, | |
569 | 9.561463309176490874525827051566494939295E0L, | |
570 | 3.203719484883967351729513662089163356911E1L, | |
571 | 5.497294687660930446641539152123568668447E1L, | |
572 | 4.391158169390578768508675452986948391118E1L, | |
573 | 1.347836630730048077907818943625789418378E1L, | |
574 | /* 1.000000000000000000000000000000000000000E0 */ | |
575 | }; | |
576 | ||
577 | /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x), | |
578 | Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2)) | |
579 | Peak relative error 2.2e-35 | |
580 | 0.375 <= 1/x <= 0.4375 */ | |
581 | #define NQ2r3_2r7N 9 | |
1f5649f8 | 582 | static const long double Q2r3_2r7N[NQ2r3_2r7N + 1] = { |
f524e4a2 AJ |
583 | 4.455027774980750211349941766420190722088E-7L, |
584 | 4.031998274578520170631601850866780366466E-5L, | |
585 | 1.273987274325947007856695677491340636339E-3L, | |
586 | 1.818754543377448509897226554179659122873E-2L, | |
587 | 1.266748858326568264126353051352269875352E-1L, | |
588 | 4.327578594728723821137731555139472880414E-1L, | |
589 | 6.892532471436503074928194969154192615359E-1L, | |
590 | 4.490775818438716873422163588640262036506E-1L, | |
591 | 8.649615949297322440032000346117031581572E-2L, | |
592 | 7.261345286655345047417257611469066147561E-4L, | |
593 | }; | |
594 | #define NQ2r3_2r7D 8 | |
1f5649f8 | 595 | static const long double Q2r3_2r7D[NQ2r3_2r7D + 1] = { |
f524e4a2 AJ |
596 | 6.082600739680555266312417978064954793142E-6L, |
597 | 5.693622538165494742945717226571441747567E-4L, | |
598 | 1.901625907009092204458328768129666975975E-2L, | |
599 | 2.958689532697857335456896889409923371570E-1L, | |
600 | 2.343124711045660081603809437993368799568E0L, | |
601 | 9.665894032187458293568704885528192804376E0L, | |
602 | 2.035273104990617136065743426322454881353E1L, | |
603 | 2.044102010478792896815088858740075165531E1L, | |
604 | 8.445937177863155827844146643468706599304E0L, | |
605 | /* 1.000000000000000000000000000000000000000E0 */ | |
606 | }; | |
607 | ||
608 | /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x), | |
609 | Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2)) | |
610 | Peak relative error 3.1e-36 | |
611 | 0.4375 <= 1/x <= 0.5 */ | |
612 | #define NQ2_2r3N 9 | |
1f5649f8 | 613 | static const long double Q2_2r3N[NQ2_2r3N + 1] = { |
f524e4a2 AJ |
614 | 2.817566786579768804844367382809101929314E-6L, |
615 | 2.122772176396691634147024348373539744935E-4L, | |
616 | 5.501378031780457828919593905395747517585E-3L, | |
617 | 6.355374424341762686099147452020466524659E-2L, | |
618 | 3.539652320122661637429658698954748337223E-1L, | |
619 | 9.571721066119617436343740541777014319695E-1L, | |
620 | 1.196258777828426399432550698612171955305E0L, | |
621 | 6.069388659458926158392384709893753793967E-1L, | |
622 | 9.026746127269713176512359976978248763621E-2L, | |
623 | 5.317668723070450235320878117210807236375E-4L, | |
624 | }; | |
625 | #define NQ2_2r3D 8 | |
1f5649f8 | 626 | static const long double Q2_2r3D[NQ2_2r3D + 1] = { |
f524e4a2 AJ |
627 | 3.846924354014260866793741072933159380158E-5L, |
628 | 3.017562820057704325510067178327449946763E-3L, | |
629 | 8.356305620686867949798885808540444210935E-2L, | |
630 | 1.068314930499906838814019619594424586273E0L, | |
631 | 6.900279623894821067017966573640732685233E0L, | |
632 | 2.307667390886377924509090271780839563141E1L, | |
633 | 3.921043465412723970791036825401273528513E1L, | |
634 | 3.167569478939719383241775717095729233436E1L, | |
635 | 1.051023841699200920276198346301543665909E1L, | |
636 | /* 1.000000000000000000000000000000000000000E0*/ | |
637 | }; | |
638 | ||
639 | ||
640 | /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */ | |
641 | ||
642 | static long double | |
60a06b7c | 643 | neval (long double x, const long double *p, int n) |
f524e4a2 AJ |
644 | { |
645 | long double y; | |
646 | ||
647 | p += n; | |
648 | y = *p--; | |
649 | do | |
650 | { | |
651 | y = y * x + *p--; | |
652 | } | |
653 | while (--n > 0); | |
654 | return y; | |
655 | } | |
656 | ||
657 | ||
658 | /* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */ | |
659 | ||
660 | static long double | |
60a06b7c | 661 | deval (long double x, const long double *p, int n) |
f524e4a2 AJ |
662 | { |
663 | long double y; | |
664 | ||
665 | p += n; | |
666 | y = x + *p--; | |
667 | do | |
668 | { | |
669 | y = y * x + *p--; | |
670 | } | |
671 | while (--n > 0); | |
672 | return y; | |
673 | } | |
674 | ||
675 | ||
676 | /* Bessel function of the first kind, order zero. */ | |
677 | ||
678 | long double | |
679 | __ieee754_j0l (long double x) | |
680 | { | |
681 | long double xx, xinv, z, p, q, c, s, cc, ss; | |
682 | ||
d81f90cc | 683 | if (! isfinite (x)) |
f524e4a2 AJ |
684 | { |
685 | if (x != x) | |
d73e7bdb | 686 | return x + x; |
f524e4a2 AJ |
687 | else |
688 | return 0.0L; | |
689 | } | |
690 | if (x == 0.0L) | |
691 | return 1.0L; | |
692 | ||
693 | xx = fabsl (x); | |
694 | if (xx <= 2.0L) | |
695 | { | |
0a90a8f2 JM |
696 | if (xx < 0x1p-57L) |
697 | return 1.0L; | |
f524e4a2 AJ |
698 | /* 0 <= x <= 2 */ |
699 | z = xx * xx; | |
700 | p = z * z * neval (z, J0_2N, NJ0_2N) / deval (z, J0_2D, NJ0_2D); | |
701 | p -= 0.25L * z; | |
702 | p += 1.0L; | |
703 | return p; | |
704 | } | |
705 | ||
2a185d32 JM |
706 | /* X = x - pi/4 |
707 | cos(X) = cos(x) cos(pi/4) + sin(x) sin(pi/4) | |
708 | = 1/sqrt(2) * (cos(x) + sin(x)) | |
709 | sin(X) = sin(x) cos(pi/4) - cos(x) sin(pi/4) | |
710 | = 1/sqrt(2) * (sin(x) - cos(x)) | |
711 | sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) | |
712 | cf. Fdlibm. */ | |
713 | __sincosl (xx, &s, &c); | |
714 | ss = s - c; | |
715 | cc = s + c; | |
98c48fe5 JM |
716 | if (xx <= LDBL_MAX / 2.0L) |
717 | { | |
718 | z = -__cosl (xx + xx); | |
719 | if ((s * c) < 0) | |
720 | cc = z / ss; | |
721 | else | |
722 | ss = z / cc; | |
723 | } | |
2a185d32 JM |
724 | |
725 | if (xx > 0x1p256L) | |
726 | return ONEOSQPI * cc / __ieee754_sqrtl (xx); | |
727 | ||
f524e4a2 AJ |
728 | xinv = 1.0L / xx; |
729 | z = xinv * xinv; | |
730 | if (xinv <= 0.25) | |
731 | { | |
732 | if (xinv <= 0.125) | |
733 | { | |
734 | if (xinv <= 0.0625) | |
735 | { | |
736 | p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID); | |
737 | q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID); | |
738 | } | |
739 | else | |
740 | { | |
741 | p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D); | |
742 | q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D); | |
743 | } | |
744 | } | |
745 | else if (xinv <= 0.1875) | |
746 | { | |
747 | p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D); | |
748 | q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D); | |
749 | } | |
750 | else | |
751 | { | |
752 | p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D); | |
753 | q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D); | |
754 | } | |
755 | } /* .25 */ | |
756 | else /* if (xinv <= 0.5) */ | |
757 | { | |
758 | if (xinv <= 0.375) | |
759 | { | |
760 | if (xinv <= 0.3125) | |
761 | { | |
762 | p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D); | |
763 | q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D); | |
764 | } | |
765 | else | |
766 | { | |
767 | p = neval (z, P2r7_3r2N, NP2r7_3r2N) | |
768 | / deval (z, P2r7_3r2D, NP2r7_3r2D); | |
769 | q = neval (z, Q2r7_3r2N, NQ2r7_3r2N) | |
770 | / deval (z, Q2r7_3r2D, NQ2r7_3r2D); | |
771 | } | |
772 | } | |
773 | else if (xinv <= 0.4375) | |
774 | { | |
775 | p = neval (z, P2r3_2r7N, NP2r3_2r7N) | |
776 | / deval (z, P2r3_2r7D, NP2r3_2r7D); | |
777 | q = neval (z, Q2r3_2r7N, NQ2r3_2r7N) | |
778 | / deval (z, Q2r3_2r7D, NQ2r3_2r7D); | |
779 | } | |
780 | else | |
781 | { | |
782 | p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D); | |
783 | q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D); | |
784 | } | |
785 | } | |
786 | p = 1.0L + z * p; | |
787 | q = z * xinv * q; | |
788 | q = q - 0.125L * xinv; | |
246ec411 | 789 | z = ONEOSQPI * (p * cc - q * ss) / __ieee754_sqrtl (xx); |
f524e4a2 AJ |
790 | return z; |
791 | } | |
0ac5ae23 | 792 | strong_alias (__ieee754_j0l, __j0l_finite) |
f524e4a2 AJ |
793 | |
794 | ||
795 | /* Y0(x) = 2/pi * log(x) * J0(x) + R(x^2) | |
796 | Peak absolute error 1.7e-36 (relative where Y0 > 1) | |
797 | 0 <= x <= 2 */ | |
798 | #define NY0_2N 7 | |
799 | static long double Y0_2N[NY0_2N + 1] = { | |
800 | -1.062023609591350692692296993537002558155E19L, | |
801 | 2.542000883190248639104127452714966858866E19L, | |
802 | -1.984190771278515324281415820316054696545E18L, | |
803 | 4.982586044371592942465373274440222033891E16L, | |
804 | -5.529326354780295177243773419090123407550E14L, | |
805 | 3.013431465522152289279088265336861140391E12L, | |
806 | -7.959436160727126750732203098982718347785E9L, | |
807 | 8.230845651379566339707130644134372793322E6L, | |
808 | }; | |
809 | #define NY0_2D 7 | |
810 | static long double Y0_2D[NY0_2D + 1] = { | |
811 | 1.438972634353286978700329883122253752192E20L, | |
812 | 1.856409101981569254247700169486907405500E18L, | |
813 | 1.219693352678218589553725579802986255614E16L, | |
814 | 5.389428943282838648918475915779958097958E13L, | |
815 | 1.774125762108874864433872173544743051653E11L, | |
816 | 4.522104832545149534808218252434693007036E8L, | |
817 | 8.872187401232943927082914504125234454930E5L, | |
818 | 1.251945613186787532055610876304669413955E3L, | |
819 | /* 1.000000000000000000000000000000000000000E0 */ | |
820 | }; | |
821 | ||
05b227bd | 822 | static const long double U0 = -7.3804295108687225274343927948483016310862e-02L; |
f524e4a2 AJ |
823 | |
824 | /* Bessel function of the second kind, order zero. */ | |
825 | ||
826 | long double | |
827 | __ieee754_y0l(long double x) | |
828 | { | |
829 | long double xx, xinv, z, p, q, c, s, cc, ss; | |
830 | ||
d81f90cc | 831 | if (! isfinite (x)) |
f524e4a2 AJ |
832 | { |
833 | if (x != x) | |
d73e7bdb | 834 | return x + x; |
f524e4a2 AJ |
835 | else |
836 | return 0.0L; | |
837 | } | |
52e1b618 UD |
838 | if (x <= 0.0L) |
839 | { | |
840 | if (x < 0.0L) | |
1f510b3f AJ |
841 | return (zero / (zero * x)); |
842 | return -HUGE_VALL + x; | |
52e1b618 | 843 | } |
f524e4a2 | 844 | xx = fabsl (x); |
05b227bd DM |
845 | if (xx <= 0x1p-57) |
846 | return U0 + TWOOPI * __ieee754_logl (x); | |
f524e4a2 AJ |
847 | if (xx <= 2.0L) |
848 | { | |
849 | /* 0 <= x <= 2 */ | |
850 | z = xx * xx; | |
851 | p = neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D); | |
246ec411 | 852 | p = TWOOPI * __ieee754_logl (x) * __ieee754_j0l (x) + p; |
f524e4a2 AJ |
853 | return p; |
854 | } | |
855 | ||
2a185d32 JM |
856 | /* X = x - pi/4 |
857 | cos(X) = cos(x) cos(pi/4) + sin(x) sin(pi/4) | |
858 | = 1/sqrt(2) * (cos(x) + sin(x)) | |
859 | sin(X) = sin(x) cos(pi/4) - cos(x) sin(pi/4) | |
860 | = 1/sqrt(2) * (sin(x) - cos(x)) | |
861 | sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) | |
862 | cf. Fdlibm. */ | |
863 | __sincosl (x, &s, &c); | |
864 | ss = s - c; | |
865 | cc = s + c; | |
98c48fe5 JM |
866 | if (xx <= LDBL_MAX / 2.0L) |
867 | { | |
868 | z = -__cosl (x + x); | |
869 | if ((s * c) < 0) | |
870 | cc = z / ss; | |
871 | else | |
872 | ss = z / cc; | |
873 | } | |
2a185d32 JM |
874 | |
875 | if (xx > 0x1p256L) | |
876 | return ONEOSQPI * ss / __ieee754_sqrtl (x); | |
877 | ||
f524e4a2 AJ |
878 | xinv = 1.0L / xx; |
879 | z = xinv * xinv; | |
880 | if (xinv <= 0.25) | |
881 | { | |
882 | if (xinv <= 0.125) | |
883 | { | |
884 | if (xinv <= 0.0625) | |
885 | { | |
886 | p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID); | |
887 | q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID); | |
888 | } | |
889 | else | |
890 | { | |
891 | p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D); | |
892 | q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D); | |
893 | } | |
894 | } | |
895 | else if (xinv <= 0.1875) | |
896 | { | |
897 | p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D); | |
898 | q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D); | |
899 | } | |
900 | else | |
901 | { | |
902 | p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D); | |
903 | q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D); | |
904 | } | |
905 | } /* .25 */ | |
906 | else /* if (xinv <= 0.5) */ | |
907 | { | |
908 | if (xinv <= 0.375) | |
909 | { | |
910 | if (xinv <= 0.3125) | |
911 | { | |
912 | p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D); | |
913 | q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D); | |
914 | } | |
915 | else | |
916 | { | |
917 | p = neval (z, P2r7_3r2N, NP2r7_3r2N) | |
918 | / deval (z, P2r7_3r2D, NP2r7_3r2D); | |
919 | q = neval (z, Q2r7_3r2N, NQ2r7_3r2N) | |
920 | / deval (z, Q2r7_3r2D, NQ2r7_3r2D); | |
921 | } | |
922 | } | |
923 | else if (xinv <= 0.4375) | |
924 | { | |
925 | p = neval (z, P2r3_2r7N, NP2r3_2r7N) | |
926 | / deval (z, P2r3_2r7D, NP2r3_2r7D); | |
927 | q = neval (z, Q2r3_2r7N, NQ2r3_2r7N) | |
928 | / deval (z, Q2r3_2r7D, NQ2r3_2r7D); | |
929 | } | |
930 | else | |
931 | { | |
932 | p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D); | |
933 | q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D); | |
934 | } | |
935 | } | |
936 | p = 1.0L + z * p; | |
937 | q = z * xinv * q; | |
938 | q = q - 0.125L * xinv; | |
246ec411 | 939 | z = ONEOSQPI * (p * ss + q * cc) / __ieee754_sqrtl (x); |
f524e4a2 AJ |
940 | return z; |
941 | } | |
0ac5ae23 | 942 | strong_alias (__ieee754_y0l, __y0l_finite) |