Don't use __glibcxx_assert to check class invariants
[official-gcc.git] / libquadmath / math / j1q.c
blob5eb705084e2722667003bc2fbfaf49817cf41857
1 /* j1l.c
3 * Bessel function of order one
7 * SYNOPSIS:
9 * __float128 x, y, j1q();
11 * y = j1q( x );
15 * DESCRIPTION:
17 * Returns Bessel function of first kind, order one of the argument.
19 * The domain is divided into two major intervals [0, 2] and
20 * (2, infinity). In the first interval the rational approximation is
21 * J1(x) = .5x + x x^2 R(x^2)
23 * The second interval is further partitioned into eight equal segments
24 * of 1/x.
25 * J1(x) = sqrt(2/(pi x)) (P1(x) cos(X) - Q1(x) sin(X)),
26 * X = x - 3 pi / 4,
28 * and the auxiliary functions are given by
30 * J1(x)cos(X) + Y1(x)sin(X) = sqrt( 2/(pi x)) P1(x),
31 * P1(x) = 1 + 1/x^2 R(1/x^2)
33 * Y1(x)cos(X) - J1(x)sin(X) = sqrt( 2/(pi x)) Q1(x),
34 * Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)).
38 * ACCURACY:
40 * Absolute error:
41 * arithmetic domain # trials peak rms
42 * IEEE 0, 30 100000 2.8e-34 2.7e-35
47 /* y1l.c
49 * Bessel function of the second kind, order one
53 * SYNOPSIS:
55 * __float128, y, y1q();
57 * y = y1q( x );
61 * DESCRIPTION:
63 * Returns Bessel function of the second kind, of order
64 * one, of the argument.
66 * The domain is divided into two major intervals [0, 2] and
67 * (2, infinity). In the first interval the rational approximation is
68 * Y1(x) = 2/pi * (log(x) * J1(x) - 1/x) + x R(x^2) .
69 * In the second interval the approximation is the same as for J1(x), and
70 * Y1(x) = sqrt(2/(pi x)) (P1(x) sin(X) + Q1(x) cos(X)),
71 * X = x - 3 pi / 4.
73 * ACCURACY:
75 * Absolute error, when y0(x) < 1; else relative error:
77 * arithmetic domain # trials peak rms
78 * IEEE 0, 30 100000 2.7e-34 2.9e-35
82 /* Copyright 2001 by Stephen L. Moshier (moshier@na-net.onrl.gov).
84 This library is free software; you can redistribute it and/or
85 modify it under the terms of the GNU Lesser General Public
86 License as published by the Free Software Foundation; either
87 version 2.1 of the License, or (at your option) any later version.
89 This library is distributed in the hope that it will be useful,
90 but WITHOUT ANY WARRANTY; without even the implied warranty of
91 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
92 Lesser General Public License for more details.
94 You should have received a copy of the GNU Lesser General Public
95 License along with this library; if not, write to the Free Software
96 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
98 #include <errno.h>
99 #include "quadmath-imp.h"
101 /* 1 / sqrt(pi) */
102 static const __float128 ONEOSQPI = 5.6418958354775628694807945156077258584405E-1Q;
103 /* 2 / pi */
104 static const __float128 TWOOPI = 6.3661977236758134307553505349005744813784E-1Q;
105 static const __float128 zero = 0.0Q;
107 /* J1(x) = .5x + x x^2 R(x^2)
108 Peak relative error 1.9e-35
109 0 <= x <= 2 */
110 #define NJ0_2N 6
111 static const __float128 J0_2N[NJ0_2N + 1] = {
112 -5.943799577386942855938508697619735179660E16Q,
113 1.812087021305009192259946997014044074711E15Q,
114 -2.761698314264509665075127515729146460895E13Q,
115 2.091089497823600978949389109350658815972E11Q,
116 -8.546413231387036372945453565654130054307E8Q,
117 1.797229225249742247475464052741320612261E6Q,
118 -1.559552840946694171346552770008812083969E3Q
120 #define NJ0_2D 6
121 static const __float128 J0_2D[NJ0_2D + 1] = {
122 9.510079323819108569501613916191477479397E17Q,
123 1.063193817503280529676423936545854693915E16Q,
124 5.934143516050192600795972192791775226920E13Q,
125 2.168000911950620999091479265214368352883E11Q,
126 5.673775894803172808323058205986256928794E8Q,
127 1.080329960080981204840966206372671147224E6Q,
128 1.411951256636576283942477881535283304912E3Q,
129 /* 1.000000000000000000000000000000000000000E0Q */
132 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
133 0 <= 1/x <= .0625
134 Peak relative error 3.6e-36 */
135 #define NP16_IN 9
136 static const __float128 P16_IN[NP16_IN + 1] = {
137 5.143674369359646114999545149085139822905E-16Q,
138 4.836645664124562546056389268546233577376E-13Q,
139 1.730945562285804805325011561498453013673E-10Q,
140 3.047976856147077889834905908605310585810E-8Q,
141 2.855227609107969710407464739188141162386E-6Q,
142 1.439362407936705484122143713643023998457E-4Q,
143 3.774489768532936551500999699815873422073E-3Q,
144 4.723962172984642566142399678920790598426E-2Q,
145 2.359289678988743939925017240478818248735E-1Q,
146 3.032580002220628812728954785118117124520E-1Q,
148 #define NP16_ID 9
149 static const __float128 P16_ID[NP16_ID + 1] = {
150 4.389268795186898018132945193912677177553E-15Q,
151 4.132671824807454334388868363256830961655E-12Q,
152 1.482133328179508835835963635130894413136E-9Q,
153 2.618941412861122118906353737117067376236E-7Q,
154 2.467854246740858470815714426201888034270E-5Q,
155 1.257192927368839847825938545925340230490E-3Q,
156 3.362739031941574274949719324644120720341E-2Q,
157 4.384458231338934105875343439265370178858E-1Q,
158 2.412830809841095249170909628197264854651E0Q,
159 4.176078204111348059102962617368214856874E0Q,
160 /* 1.000000000000000000000000000000000000000E0 */
163 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
164 0.0625 <= 1/x <= 0.125
165 Peak relative error 1.9e-36 */
166 #define NP8_16N 11
167 static const __float128 P8_16N[NP8_16N + 1] = {
168 2.984612480763362345647303274082071598135E-16Q,
169 1.923651877544126103941232173085475682334E-13Q,
170 4.881258879388869396043760693256024307743E-11Q,
171 6.368866572475045408480898921866869811889E-9Q,
172 4.684818344104910450523906967821090796737E-7Q,
173 2.005177298271593587095982211091300382796E-5Q,
174 4.979808067163957634120681477207147536182E-4Q,
175 6.946005761642579085284689047091173581127E-3Q,
176 5.074601112955765012750207555985299026204E-2Q,
177 1.698599455896180893191766195194231825379E-1Q,
178 1.957536905259237627737222775573623779638E-1Q,
179 2.991314703282528370270179989044994319374E-2Q,
181 #define NP8_16D 10
182 static const __float128 P8_16D[NP8_16D + 1] = {
183 2.546869316918069202079580939942463010937E-15Q,
184 1.644650111942455804019788382157745229955E-12Q,
185 4.185430770291694079925607420808011147173E-10Q,
186 5.485331966975218025368698195861074143153E-8Q,
187 4.062884421686912042335466327098932678905E-6Q,
188 1.758139661060905948870523641319556816772E-4Q,
189 4.445143889306356207566032244985607493096E-3Q,
190 6.391901016293512632765621532571159071158E-2Q,
191 4.933040207519900471177016015718145795434E-1Q,
192 1.839144086168947712971630337250761842976E0Q,
193 2.715120873995490920415616716916149586579E0Q,
194 /* 1.000000000000000000000000000000000000000E0 */
197 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
198 0.125 <= 1/x <= 0.1875
199 Peak relative error 1.3e-36 */
200 #define NP5_8N 10
201 static const __float128 P5_8N[NP5_8N + 1] = {
202 2.837678373978003452653763806968237227234E-12Q,
203 9.726641165590364928442128579282742354806E-10Q,
204 1.284408003604131382028112171490633956539E-7Q,
205 8.524624695868291291250573339272194285008E-6Q,
206 3.111516908953172249853673787748841282846E-4Q,
207 6.423175156126364104172801983096596409176E-3Q,
208 7.430220589989104581004416356260692450652E-2Q,
209 4.608315409833682489016656279567605536619E-1Q,
210 1.396870223510964882676225042258855977512E0Q,
211 1.718500293904122365894630460672081526236E0Q,
212 5.465927698800862172307352821870223855365E-1Q
214 #define NP5_8D 10
215 static const __float128 P5_8D[NP5_8D + 1] = {
216 2.421485545794616609951168511612060482715E-11Q,
217 8.329862750896452929030058039752327232310E-9Q,
218 1.106137992233383429630592081375289010720E-6Q,
219 7.405786153760681090127497796448503306939E-5Q,
220 2.740364785433195322492093333127633465227E-3Q,
221 5.781246470403095224872243564165254652198E-2Q,
222 6.927711353039742469918754111511109983546E-1Q,
223 4.558679283460430281188304515922826156690E0Q,
224 1.534468499844879487013168065728837900009E1Q,
225 2.313927430889218597919624843161569422745E1Q,
226 1.194506341319498844336768473218382828637E1Q,
227 /* 1.000000000000000000000000000000000000000E0 */
230 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
231 Peak relative error 1.4e-36
232 0.1875 <= 1/x <= 0.25 */
233 #define NP4_5N 10
234 static const __float128 P4_5N[NP4_5N + 1] = {
235 1.846029078268368685834261260420933914621E-10Q,
236 3.916295939611376119377869680335444207768E-8Q,
237 3.122158792018920627984597530935323997312E-6Q,
238 1.218073444893078303994045653603392272450E-4Q,
239 2.536420827983485448140477159977981844883E-3Q,
240 2.883011322006690823959367922241169171315E-2Q,
241 1.755255190734902907438042414495469810830E-1Q,
242 5.379317079922628599870898285488723736599E-1Q,
243 7.284904050194300773890303361501726561938E-1Q,
244 3.270110346613085348094396323925000362813E-1Q,
245 1.804473805689725610052078464951722064757E-2Q,
247 #define NP4_5D 9
248 static const __float128 P4_5D[NP4_5D + 1] = {
249 1.575278146806816970152174364308980863569E-9Q,
250 3.361289173657099516191331123405675054321E-7Q,
251 2.704692281550877810424745289838790693708E-5Q,
252 1.070854930483999749316546199273521063543E-3Q,
253 2.282373093495295842598097265627962125411E-2Q,
254 2.692025460665354148328762368240343249830E-1Q,
255 1.739892942593664447220951225734811133759E0Q,
256 5.890727576752230385342377570386657229324E0Q,
257 9.517442287057841500750256954117735128153E0Q,
258 6.100616353935338240775363403030137736013E0Q,
259 /* 1.000000000000000000000000000000000000000E0 */
262 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
263 Peak relative error 3.0e-36
264 0.25 <= 1/x <= 0.3125 */
265 #define NP3r2_4N 9
266 static const __float128 P3r2_4N[NP3r2_4N + 1] = {
267 8.240803130988044478595580300846665863782E-8Q,
268 1.179418958381961224222969866406483744580E-5Q,
269 6.179787320956386624336959112503824397755E-4Q,
270 1.540270833608687596420595830747166658383E-2Q,
271 1.983904219491512618376375619598837355076E-1Q,
272 1.341465722692038870390470651608301155565E0Q,
273 4.617865326696612898792238245990854646057E0Q,
274 7.435574801812346424460233180412308000587E0Q,
275 4.671327027414635292514599201278557680420E0Q,
276 7.299530852495776936690976966995187714739E-1Q,
278 #define NP3r2_4D 9
279 static const __float128 P3r2_4D[NP3r2_4D + 1] = {
280 7.032152009675729604487575753279187576521E-7Q,
281 1.015090352324577615777511269928856742848E-4Q,
282 5.394262184808448484302067955186308730620E-3Q,
283 1.375291438480256110455809354836988584325E-1Q,
284 1.836247144461106304788160919310404376670E0Q,
285 1.314378564254376655001094503090935880349E1Q,
286 4.957184590465712006934452500894672343488E1Q,
287 9.287394244300647738855415178790263465398E1Q,
288 7.652563275535900609085229286020552768399E1Q,
289 2.147042473003074533150718117770093209096E1Q,
290 /* 1.000000000000000000000000000000000000000E0 */
293 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
294 Peak relative error 1.0e-35
295 0.3125 <= 1/x <= 0.375 */
296 #define NP2r7_3r2N 9
297 static const __float128 P2r7_3r2N[NP2r7_3r2N + 1] = {
298 4.599033469240421554219816935160627085991E-7Q,
299 4.665724440345003914596647144630893997284E-5Q,
300 1.684348845667764271596142716944374892756E-3Q,
301 2.802446446884455707845985913454440176223E-2Q,
302 2.321937586453963310008279956042545173930E-1Q,
303 9.640277413988055668692438709376437553804E-1Q,
304 1.911021064710270904508663334033003246028E0Q,
305 1.600811610164341450262992138893970224971E0Q,
306 4.266299218652587901171386591543457861138E-1Q,
307 1.316470424456061252962568223251247207325E-2Q,
309 #define NP2r7_3r2D 8
310 static const __float128 P2r7_3r2D[NP2r7_3r2D + 1] = {
311 3.924508608545520758883457108453520099610E-6Q,
312 4.029707889408829273226495756222078039823E-4Q,
313 1.484629715787703260797886463307469600219E-2Q,
314 2.553136379967180865331706538897231588685E-1Q,
315 2.229457223891676394409880026887106228740E0Q,
316 1.005708903856384091956550845198392117318E1Q,
317 2.277082659664386953166629360352385889558E1Q,
318 2.384726835193630788249826630376533988245E1Q,
319 9.700989749041320895890113781610939632410E0Q,
320 /* 1.000000000000000000000000000000000000000E0 */
323 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
324 Peak relative error 1.7e-36
325 0.3125 <= 1/x <= 0.4375 */
326 #define NP2r3_2r7N 9
327 static const __float128 P2r3_2r7N[NP2r3_2r7N + 1] = {
328 3.916766777108274628543759603786857387402E-6Q,
329 3.212176636756546217390661984304645137013E-4Q,
330 9.255768488524816445220126081207248947118E-3Q,
331 1.214853146369078277453080641911700735354E-1Q,
332 7.855163309847214136198449861311404633665E-1Q,
333 2.520058073282978403655488662066019816540E0Q,
334 3.825136484837545257209234285382183711466E0Q,
335 2.432569427554248006229715163865569506873E0Q,
336 4.877934835018231178495030117729800489743E-1Q,
337 1.109902737860249670981355149101343427885E-2Q,
339 #define NP2r3_2r7D 8
340 static const __float128 P2r3_2r7D[NP2r3_2r7D + 1] = {
341 3.342307880794065640312646341190547184461E-5Q,
342 2.782182891138893201544978009012096558265E-3Q,
343 8.221304931614200702142049236141249929207E-2Q,
344 1.123728246291165812392918571987858010949E0Q,
345 7.740482453652715577233858317133423434590E0Q,
346 2.737624677567945952953322566311201919139E1Q,
347 4.837181477096062403118304137851260715475E1Q,
348 3.941098643468580791437772701093795299274E1Q,
349 1.245821247166544627558323920382547533630E1Q,
350 /* 1.000000000000000000000000000000000000000E0 */
353 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
354 Peak relative error 1.7e-35
355 0.4375 <= 1/x <= 0.5 */
356 #define NP2_2r3N 8
357 static const __float128 P2_2r3N[NP2_2r3N + 1] = {
358 3.397930802851248553545191160608731940751E-4Q,
359 2.104020902735482418784312825637833698217E-2Q,
360 4.442291771608095963935342749477836181939E-1Q,
361 4.131797328716583282869183304291833754967E0Q,
362 1.819920169779026500146134832455189917589E1Q,
363 3.781779616522937565300309684282401791291E1Q,
364 3.459605449728864218972931220783543410347E1Q,
365 1.173594248397603882049066603238568316561E1Q,
366 9.455702270242780642835086549285560316461E-1Q,
368 #define NP2_2r3D 8
369 static const __float128 P2_2r3D[NP2_2r3D + 1] = {
370 2.899568897241432883079888249845707400614E-3Q,
371 1.831107138190848460767699919531132426356E-1Q,
372 3.999350044057883839080258832758908825165E0Q,
373 3.929041535867957938340569419874195303712E1Q,
374 1.884245613422523323068802689915538908291E2Q,
375 4.461469948819229734353852978424629815929E2Q,
376 5.004998753999796821224085972610636347903E2Q,
377 2.386342520092608513170837883757163414100E2Q,
378 3.791322528149347975999851588922424189957E1Q,
379 /* 1.000000000000000000000000000000000000000E0 */
382 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
383 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
384 Peak relative error 8.0e-36
385 0 <= 1/x <= .0625 */
386 #define NQ16_IN 10
387 static const __float128 Q16_IN[NQ16_IN + 1] = {
388 -3.917420835712508001321875734030357393421E-18Q,
389 -4.440311387483014485304387406538069930457E-15Q,
390 -1.951635424076926487780929645954007139616E-12Q,
391 -4.318256438421012555040546775651612810513E-10Q,
392 -5.231244131926180765270446557146989238020E-8Q,
393 -3.540072702902043752460711989234732357653E-6Q,
394 -1.311017536555269966928228052917534882984E-4Q,
395 -2.495184669674631806622008769674827575088E-3Q,
396 -2.141868222987209028118086708697998506716E-2Q,
397 -6.184031415202148901863605871197272650090E-2Q,
398 -1.922298704033332356899546792898156493887E-2Q,
400 #define NQ16_ID 9
401 static const __float128 Q16_ID[NQ16_ID + 1] = {
402 3.820418034066293517479619763498400162314E-17Q,
403 4.340702810799239909648911373329149354911E-14Q,
404 1.914985356383416140706179933075303538524E-11Q,
405 4.262333682610888819476498617261895474330E-9Q,
406 5.213481314722233980346462747902942182792E-7Q,
407 3.585741697694069399299005316809954590558E-5Q,
408 1.366513429642842006385029778105539457546E-3Q,
409 2.745282599850704662726337474371355160594E-2Q,
410 2.637644521611867647651200098449903330074E-1Q,
411 1.006953426110765984590782655598680488746E0Q,
412 /* 1.000000000000000000000000000000000000000E0 */
415 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
416 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
417 Peak relative error 1.9e-36
418 0.0625 <= 1/x <= 0.125 */
419 #define NQ8_16N 11
420 static const __float128 Q8_16N[NQ8_16N + 1] = {
421 -2.028630366670228670781362543615221542291E-17Q,
422 -1.519634620380959966438130374006858864624E-14Q,
423 -4.540596528116104986388796594639405114524E-12Q,
424 -7.085151756671466559280490913558388648274E-10Q,
425 -6.351062671323970823761883833531546885452E-8Q,
426 -3.390817171111032905297982523519503522491E-6Q,
427 -1.082340897018886970282138836861233213972E-4Q,
428 -2.020120801187226444822977006648252379508E-3Q,
429 -2.093169910981725694937457070649605557555E-2Q,
430 -1.092176538874275712359269481414448063393E-1Q,
431 -2.374790947854765809203590474789108718733E-1Q,
432 -1.365364204556573800719985118029601401323E-1Q,
434 #define NQ8_16D 11
435 static const __float128 Q8_16D[NQ8_16D + 1] = {
436 1.978397614733632533581207058069628242280E-16Q,
437 1.487361156806202736877009608336766720560E-13Q,
438 4.468041406888412086042576067133365913456E-11Q,
439 7.027822074821007443672290507210594648877E-9Q,
440 6.375740580686101224127290062867976007374E-7Q,
441 3.466887658320002225888644977076410421940E-5Q,
442 1.138625640905289601186353909213719596986E-3Q,
443 2.224470799470414663443449818235008486439E-2Q,
444 2.487052928527244907490589787691478482358E-1Q,
445 1.483927406564349124649083853892380899217E0Q,
446 4.182773513276056975777258788903489507705E0Q,
447 4.419665392573449746043880892524360870944E0Q,
448 /* 1.000000000000000000000000000000000000000E0 */
451 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
452 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
453 Peak relative error 1.5e-35
454 0.125 <= 1/x <= 0.1875 */
455 #define NQ5_8N 10
456 static const __float128 Q5_8N[NQ5_8N + 1] = {
457 -3.656082407740970534915918390488336879763E-13Q,
458 -1.344660308497244804752334556734121771023E-10Q,
459 -1.909765035234071738548629788698150760791E-8Q,
460 -1.366668038160120210269389551283666716453E-6Q,
461 -5.392327355984269366895210704976314135683E-5Q,
462 -1.206268245713024564674432357634540343884E-3Q,
463 -1.515456784370354374066417703736088291287E-2Q,
464 -1.022454301137286306933217746545237098518E-1Q,
465 -3.373438906472495080504907858424251082240E-1Q,
466 -4.510782522110845697262323973549178453405E-1Q,
467 -1.549000892545288676809660828213589804884E-1Q,
469 #define NQ5_8D 10
470 static const __float128 Q5_8D[NQ5_8D + 1] = {
471 3.565550843359501079050699598913828460036E-12Q,
472 1.321016015556560621591847454285330528045E-9Q,
473 1.897542728662346479999969679234270605975E-7Q,
474 1.381720283068706710298734234287456219474E-5Q,
475 5.599248147286524662305325795203422873725E-4Q,
476 1.305442352653121436697064782499122164843E-2Q,
477 1.750234079626943298160445750078631894985E-1Q,
478 1.311420542073436520965439883806946678491E0Q,
479 5.162757689856842406744504211089724926650E0Q,
480 9.527760296384704425618556332087850581308E0Q,
481 6.604648207463236667912921642545100248584E0Q,
482 /* 1.000000000000000000000000000000000000000E0 */
485 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
486 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
487 Peak relative error 1.3e-35
488 0.1875 <= 1/x <= 0.25 */
489 #define NQ4_5N 10
490 static const __float128 Q4_5N[NQ4_5N + 1] = {
491 -4.079513568708891749424783046520200903755E-11Q,
492 -9.326548104106791766891812583019664893311E-9Q,
493 -8.016795121318423066292906123815687003356E-7Q,
494 -3.372350544043594415609295225664186750995E-5Q,
495 -7.566238665947967882207277686375417983917E-4Q,
496 -9.248861580055565402130441618521591282617E-3Q,
497 -6.033106131055851432267702948850231270338E-2Q,
498 -1.966908754799996793730369265431584303447E-1Q,
499 -2.791062741179964150755788226623462207560E-1Q,
500 -1.255478605849190549914610121863534191666E-1Q,
501 -4.320429862021265463213168186061696944062E-3Q,
503 #define NQ4_5D 9
504 static const __float128 Q4_5D[NQ4_5D + 1] = {
505 3.978497042580921479003851216297330701056E-10Q,
506 9.203304163828145809278568906420772246666E-8Q,
507 8.059685467088175644915010485174545743798E-6Q,
508 3.490187375993956409171098277561669167446E-4Q,
509 8.189109654456872150100501732073810028829E-3Q,
510 1.072572867311023640958725265762483033769E-1Q,
511 7.790606862409960053675717185714576937994E-1Q,
512 3.016049768232011196434185423512777656328E0Q,
513 5.722963851442769787733717162314477949360E0Q,
514 4.510527838428473279647251350931380867663E0Q,
515 /* 1.000000000000000000000000000000000000000E0 */
518 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
519 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
520 Peak relative error 2.1e-35
521 0.25 <= 1/x <= 0.3125 */
522 #define NQ3r2_4N 9
523 static const __float128 Q3r2_4N[NQ3r2_4N + 1] = {
524 -1.087480809271383885936921889040388133627E-8Q,
525 -1.690067828697463740906962973479310170932E-6Q,
526 -9.608064416995105532790745641974762550982E-5Q,
527 -2.594198839156517191858208513873961837410E-3Q,
528 -3.610954144421543968160459863048062977822E-2Q,
529 -2.629866798251843212210482269563961685666E-1Q,
530 -9.709186825881775885917984975685752956660E-1Q,
531 -1.667521829918185121727268867619982417317E0Q,
532 -1.109255082925540057138766105229900943501E0Q,
533 -1.812932453006641348145049323713469043328E-1Q,
535 #define NQ3r2_4D 9
536 static const __float128 Q3r2_4D[NQ3r2_4D + 1] = {
537 1.060552717496912381388763753841473407026E-7Q,
538 1.676928002024920520786883649102388708024E-5Q,
539 9.803481712245420839301400601140812255737E-4Q,
540 2.765559874262309494758505158089249012930E-2Q,
541 4.117921827792571791298862613287549140706E-1Q,
542 3.323769515244751267093378361930279161413E0Q,
543 1.436602494405814164724810151689705353670E1Q,
544 3.163087869617098638064881410646782408297E1Q,
545 3.198181264977021649489103980298349589419E1Q,
546 1.203649258862068431199471076202897823272E1Q,
547 /* 1.000000000000000000000000000000000000000E0 */
550 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
551 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
552 Peak relative error 1.6e-36
553 0.3125 <= 1/x <= 0.375 */
554 #define NQ2r7_3r2N 9
555 static const __float128 Q2r7_3r2N[NQ2r7_3r2N + 1] = {
556 -1.723405393982209853244278760171643219530E-7Q,
557 -2.090508758514655456365709712333460087442E-5Q,
558 -9.140104013370974823232873472192719263019E-4Q,
559 -1.871349499990714843332742160292474780128E-2Q,
560 -1.948930738119938669637865956162512983416E-1Q,
561 -1.048764684978978127908439526343174139788E0Q,
562 -2.827714929925679500237476105843643064698E0Q,
563 -3.508761569156476114276988181329773987314E0Q,
564 -1.669332202790211090973255098624488308989E0Q,
565 -1.930796319299022954013840684651016077770E-1Q,
567 #define NQ2r7_3r2D 9
568 static const __float128 Q2r7_3r2D[NQ2r7_3r2D + 1] = {
569 1.680730662300831976234547482334347983474E-6Q,
570 2.084241442440551016475972218719621841120E-4Q,
571 9.445316642108367479043541702688736295579E-3Q,
572 2.044637889456631896650179477133252184672E-1Q,
573 2.316091982244297350829522534435350078205E0Q,
574 1.412031891783015085196708811890448488865E1Q,
575 4.583830154673223384837091077279595496149E1Q,
576 7.549520609270909439885998474045974122261E1Q,
577 5.697605832808113367197494052388203310638E1Q,
578 1.601496240876192444526383314589371686234E1Q,
579 /* 1.000000000000000000000000000000000000000E0 */
582 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
583 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
584 Peak relative error 9.5e-36
585 0.375 <= 1/x <= 0.4375 */
586 #define NQ2r3_2r7N 9
587 static const __float128 Q2r3_2r7N[NQ2r3_2r7N + 1] = {
588 -8.603042076329122085722385914954878953775E-7Q,
589 -7.701746260451647874214968882605186675720E-5Q,
590 -2.407932004380727587382493696877569654271E-3Q,
591 -3.403434217607634279028110636919987224188E-2Q,
592 -2.348707332185238159192422084985713102877E-1Q,
593 -7.957498841538254916147095255700637463207E-1Q,
594 -1.258469078442635106431098063707934348577E0Q,
595 -8.162415474676345812459353639449971369890E-1Q,
596 -1.581783890269379690141513949609572806898E-1Q,
597 -1.890595651683552228232308756569450822905E-3Q,
599 #define NQ2r3_2r7D 8
600 static const __float128 Q2r3_2r7D[NQ2r3_2r7D + 1] = {
601 8.390017524798316921170710533381568175665E-6Q,
602 7.738148683730826286477254659973968763659E-4Q,
603 2.541480810958665794368759558791634341779E-2Q,
604 3.878879789711276799058486068562386244873E-1Q,
605 3.003783779325811292142957336802456109333E0Q,
606 1.206480374773322029883039064575464497400E1Q,
607 2.458414064785315978408974662900438351782E1Q,
608 2.367237826273668567199042088835448715228E1Q,
609 9.231451197519171090875569102116321676763E0Q,
610 /* 1.000000000000000000000000000000000000000E0 */
613 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
614 Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
615 Peak relative error 1.4e-36
616 0.4375 <= 1/x <= 0.5 */
617 #define NQ2_2r3N 9
618 static const __float128 Q2_2r3N[NQ2_2r3N + 1] = {
619 -5.552507516089087822166822364590806076174E-6Q,
620 -4.135067659799500521040944087433752970297E-4Q,
621 -1.059928728869218962607068840646564457980E-2Q,
622 -1.212070036005832342565792241385459023801E-1Q,
623 -6.688350110633603958684302153362735625156E-1Q,
624 -1.793587878197360221340277951304429821582E0Q,
625 -2.225407682237197485644647380483725045326E0Q,
626 -1.123402135458940189438898496348239744403E0Q,
627 -1.679187241566347077204805190763597299805E-1Q,
628 -1.458550613639093752909985189067233504148E-3Q,
630 #define NQ2_2r3D 8
631 static const __float128 Q2_2r3D[NQ2_2r3D + 1] = {
632 5.415024336507980465169023996403597916115E-5Q,
633 4.179246497380453022046357404266022870788E-3Q,
634 1.136306384261959483095442402929502368598E-1Q,
635 1.422640343719842213484515445393284072830E0Q,
636 8.968786703393158374728850922289204805764E0Q,
637 2.914542473339246127533384118781216495934E1Q,
638 4.781605421020380669870197378210457054685E1Q,
639 3.693865837171883152382820584714795072937E1Q,
640 1.153220502744204904763115556224395893076E1Q,
641 /* 1.000000000000000000000000000000000000000E0 */
645 /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
647 static __float128
648 neval (__float128 x, const __float128 *p, int n)
650 __float128 y;
652 p += n;
653 y = *p--;
656 y = y * x + *p--;
658 while (--n > 0);
659 return y;
663 /* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */
665 static __float128
666 deval (__float128 x, const __float128 *p, int n)
668 __float128 y;
670 p += n;
671 y = x + *p--;
674 y = y * x + *p--;
676 while (--n > 0);
677 return y;
681 /* Bessel function of the first kind, order one. */
683 __float128
684 j1q (__float128 x)
686 __float128 xx, xinv, z, p, q, c, s, cc, ss;
688 if (! finiteq (x))
690 if (x != x)
691 return x + x;
692 else
693 return 0.0Q;
695 if (x == 0.0Q)
696 return x;
697 xx = fabsq (x);
698 if (xx <= 0x1p-58Q)
700 __float128 ret = x * 0.5Q;
701 math_check_force_underflow (ret);
702 if (ret == 0)
703 errno = ERANGE;
704 return ret;
706 if (xx <= 2.0Q)
708 /* 0 <= x <= 2 */
709 z = xx * xx;
710 p = xx * z * neval (z, J0_2N, NJ0_2N) / deval (z, J0_2D, NJ0_2D);
711 p += 0.5Q * xx;
712 if (x < 0)
713 p = -p;
714 return p;
717 /* X = x - 3 pi/4
718 cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
719 = 1/sqrt(2) * (-cos(x) + sin(x))
720 sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
721 = -1/sqrt(2) * (sin(x) + cos(x))
722 cf. Fdlibm. */
723 sincosq (xx, &s, &c);
724 ss = -s - c;
725 cc = s - c;
726 if (xx <= FLT128_MAX / 2.0Q)
728 z = cosq (xx + xx);
729 if ((s * c) > 0)
730 cc = z / ss;
731 else
732 ss = z / cc;
735 if (xx > 0x1p256Q)
737 z = ONEOSQPI * cc / sqrtq (xx);
738 if (x < 0)
739 z = -z;
740 return z;
743 xinv = 1.0Q / xx;
744 z = xinv * xinv;
745 if (xinv <= 0.25)
747 if (xinv <= 0.125)
749 if (xinv <= 0.0625)
751 p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID);
752 q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID);
754 else
756 p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D);
757 q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D);
760 else if (xinv <= 0.1875)
762 p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D);
763 q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D);
765 else
767 p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D);
768 q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D);
770 } /* .25 */
771 else /* if (xinv <= 0.5) */
773 if (xinv <= 0.375)
775 if (xinv <= 0.3125)
777 p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D);
778 q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D);
780 else
782 p = neval (z, P2r7_3r2N, NP2r7_3r2N)
783 / deval (z, P2r7_3r2D, NP2r7_3r2D);
784 q = neval (z, Q2r7_3r2N, NQ2r7_3r2N)
785 / deval (z, Q2r7_3r2D, NQ2r7_3r2D);
788 else if (xinv <= 0.4375)
790 p = neval (z, P2r3_2r7N, NP2r3_2r7N)
791 / deval (z, P2r3_2r7D, NP2r3_2r7D);
792 q = neval (z, Q2r3_2r7N, NQ2r3_2r7N)
793 / deval (z, Q2r3_2r7D, NQ2r3_2r7D);
795 else
797 p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
798 q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
801 p = 1.0Q + z * p;
802 q = z * q;
803 q = q * xinv + 0.375Q * xinv;
804 z = ONEOSQPI * (p * cc - q * ss) / sqrtq (xx);
805 if (x < 0)
806 z = -z;
807 return z;
811 /* Y1(x) = 2/pi * (log(x) * J1(x) - 1/x) + x R(x^2)
812 Peak relative error 6.2e-38
813 0 <= x <= 2 */
814 #define NY0_2N 7
815 static __float128 Y0_2N[NY0_2N + 1] = {
816 -6.804415404830253804408698161694720833249E19Q,
817 1.805450517967019908027153056150465849237E19Q,
818 -8.065747497063694098810419456383006737312E17Q,
819 1.401336667383028259295830955439028236299E16Q,
820 -1.171654432898137585000399489686629680230E14Q,
821 5.061267920943853732895341125243428129150E11Q,
822 -1.096677850566094204586208610960870217970E9Q,
823 9.541172044989995856117187515882879304461E5Q,
825 #define NY0_2D 7
826 static __float128 Y0_2D[NY0_2D + 1] = {
827 3.470629591820267059538637461549677594549E20Q,
828 4.120796439009916326855848107545425217219E18Q,
829 2.477653371652018249749350657387030814542E16Q,
830 9.954678543353888958177169349272167762797E13Q,
831 2.957927997613630118216218290262851197754E11Q,
832 6.748421382188864486018861197614025972118E8Q,
833 1.173453425218010888004562071020305709319E6Q,
834 1.450335662961034949894009554536003377187E3Q,
835 /* 1.000000000000000000000000000000000000000E0 */
839 /* Bessel function of the second kind, order one. */
841 __float128
842 y1q (__float128 x)
844 __float128 xx, xinv, z, p, q, c, s, cc, ss;
846 if (! finiteq (x))
847 return 1 / (x + x * x);
848 if (x <= 0.0Q)
850 if (x < 0.0Q)
851 return (zero / (zero * x));
852 return -1 / zero; /* -inf and divide by zero exception. */
854 xx = fabsq (x);
855 if (xx <= 0x1p-114)
857 z = -TWOOPI / x;
858 if (isinfq (z))
859 errno = ERANGE;
860 return z;
862 if (xx <= 2.0Q)
864 /* 0 <= x <= 2 */
865 /* FIXME: SET_RESTORE_ROUNDL (FE_TONEAREST); */
866 z = xx * xx;
867 p = xx * neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D);
868 p = -TWOOPI / xx + p;
869 p = TWOOPI * logq (x) * j1q (x) + p;
870 return p;
873 /* X = x - 3 pi/4
874 cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
875 = 1/sqrt(2) * (-cos(x) + sin(x))
876 sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
877 = -1/sqrt(2) * (sin(x) + cos(x))
878 cf. Fdlibm. */
879 sincosq (xx, &s, &c);
880 ss = -s - c;
881 cc = s - c;
882 if (xx <= FLT128_MAX / 2.0Q)
884 z = cosq (xx + xx);
885 if ((s * c) > 0)
886 cc = z / ss;
887 else
888 ss = z / cc;
891 if (xx > 0x1p256Q)
892 return ONEOSQPI * ss / sqrtq (xx);
894 xinv = 1.0Q / xx;
895 z = xinv * xinv;
896 if (xinv <= 0.25)
898 if (xinv <= 0.125)
900 if (xinv <= 0.0625)
902 p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID);
903 q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID);
905 else
907 p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D);
908 q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D);
911 else if (xinv <= 0.1875)
913 p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D);
914 q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D);
916 else
918 p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D);
919 q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D);
921 } /* .25 */
922 else /* if (xinv <= 0.5) */
924 if (xinv <= 0.375)
926 if (xinv <= 0.3125)
928 p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D);
929 q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D);
931 else
933 p = neval (z, P2r7_3r2N, NP2r7_3r2N)
934 / deval (z, P2r7_3r2D, NP2r7_3r2D);
935 q = neval (z, Q2r7_3r2N, NQ2r7_3r2N)
936 / deval (z, Q2r7_3r2D, NQ2r7_3r2D);
939 else if (xinv <= 0.4375)
941 p = neval (z, P2r3_2r7N, NP2r3_2r7N)
942 / deval (z, P2r3_2r7D, NP2r3_2r7D);
943 q = neval (z, Q2r3_2r7N, NQ2r3_2r7N)
944 / deval (z, Q2r3_2r7D, NQ2r3_2r7D);
946 else
948 p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
949 q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
952 p = 1.0Q + z * p;
953 q = z * q;
954 q = q * xinv + 0.375Q * xinv;
955 z = ONEOSQPI * (p * ss + q * cc) / sqrtq (xx);
956 return z;