Back to index

glibc  2.9
e_j1l.c
Go to the documentation of this file.
00001 /*                                               j1l.c
00002  *
00003  *     Bessel function of order one
00004  *
00005  *
00006  *
00007  * SYNOPSIS:
00008  *
00009  * long double x, y, j1l();
00010  *
00011  * y = j1l( x );
00012  *
00013  *
00014  *
00015  * DESCRIPTION:
00016  *
00017  * Returns Bessel function of first kind, order one of the argument.
00018  *
00019  * The domain is divided into two major intervals [0, 2] and
00020  * (2, infinity). In the first interval the rational approximation is
00021  * J1(x) = .5x + x x^2 R(x^2)
00022  *
00023  * The second interval is further partitioned into eight equal segments
00024  * of 1/x.
00025  * J1(x) = sqrt(2/(pi x)) (P1(x) cos(X) - Q1(x) sin(X)),
00026  * X = x - 3 pi / 4,
00027  *
00028  * and the auxiliary functions are given by
00029  *
00030  * J1(x)cos(X) + Y1(x)sin(X) = sqrt( 2/(pi x)) P1(x),
00031  * P1(x) = 1 + 1/x^2 R(1/x^2)
00032  *
00033  * Y1(x)cos(X) - J1(x)sin(X) = sqrt( 2/(pi x)) Q1(x),
00034  * Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)).
00035  *
00036  *
00037  *
00038  * ACCURACY:
00039  *
00040  *                      Absolute error:
00041  * arithmetic   domain      # trials      peak         rms
00042  *    IEEE      0, 30       100000      2.8e-34      2.7e-35
00043  *
00044  *
00045  */
00046 
00047 /*                                               y1l.c
00048  *
00049  *     Bessel function of the second kind, order one
00050  *
00051  *
00052  *
00053  * SYNOPSIS:
00054  *
00055  * double x, y, y1l();
00056  *
00057  * y = y1l( x );
00058  *
00059  *
00060  *
00061  * DESCRIPTION:
00062  *
00063  * Returns Bessel function of the second kind, of order
00064  * one, of the argument.
00065  *
00066  * The domain is divided into two major intervals [0, 2] and
00067  * (2, infinity). In the first interval the rational approximation is
00068  * Y1(x) = 2/pi * (log(x) * J1(x) - 1/x) + x R(x^2) .
00069  * In the second interval the approximation is the same as for J1(x), and
00070  * Y1(x) = sqrt(2/(pi x)) (P1(x) sin(X) + Q1(x) cos(X)),
00071  * X = x - 3 pi / 4.
00072  *
00073  * ACCURACY:
00074  *
00075  *  Absolute error, when y0(x) < 1; else relative error:
00076  *
00077  * arithmetic   domain     # trials      peak         rms
00078  *    IEEE      0, 30       100000      2.7e-34     2.9e-35
00079  *
00080  */
00081 
00082 /* Copyright 2001 by Stephen L. Moshier (moshier@na-net.onrl.gov).
00083 
00084     This library is free software; you can redistribute it and/or
00085     modify it under the terms of the GNU Lesser General Public
00086     License as published by the Free Software Foundation; either
00087     version 2.1 of the License, or (at your option) any later version.
00088 
00089     This library is distributed in the hope that it will be useful,
00090     but WITHOUT ANY WARRANTY; without even the implied warranty of
00091     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00092     Lesser General Public License for more details.
00093 
00094     You should have received a copy of the GNU Lesser General Public
00095     License along with this library; if not, write to the Free Software
00096     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
00097 
00098 #include "math.h"
00099 #include "math_private.h"
00100 
00101 /* 1 / sqrt(pi) */
00102 static const long double ONEOSQPI = 5.6418958354775628694807945156077258584405E-1L;
00103 /* 2 / pi */
00104 static const long double TWOOPI = 6.3661977236758134307553505349005744813784E-1L;
00105 static const long double zero = 0.0L;
00106 
00107 /* J1(x) = .5x + x x^2 R(x^2)
00108    Peak relative error 1.9e-35
00109    0 <= x <= 2  */
00110 #define NJ0_2N 6
00111 static const long double J0_2N[NJ0_2N + 1] = {
00112  -5.943799577386942855938508697619735179660E16L,
00113   1.812087021305009192259946997014044074711E15L,
00114  -2.761698314264509665075127515729146460895E13L,
00115   2.091089497823600978949389109350658815972E11L,
00116  -8.546413231387036372945453565654130054307E8L,
00117   1.797229225249742247475464052741320612261E6L,
00118  -1.559552840946694171346552770008812083969E3L
00119 };
00120 #define NJ0_2D 6
00121 static const long double J0_2D[NJ0_2D + 1] = {
00122   9.510079323819108569501613916191477479397E17L,
00123   1.063193817503280529676423936545854693915E16L,
00124   5.934143516050192600795972192791775226920E13L,
00125   2.168000911950620999091479265214368352883E11L,
00126   5.673775894803172808323058205986256928794E8L,
00127   1.080329960080981204840966206372671147224E6L,
00128   1.411951256636576283942477881535283304912E3L,
00129  /* 1.000000000000000000000000000000000000000E0L */
00130 };
00131 
00132 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
00133    0 <= 1/x <= .0625
00134    Peak relative error 3.6e-36  */
00135 #define NP16_IN 9
00136 static const long double P16_IN[NP16_IN + 1] = {
00137   5.143674369359646114999545149085139822905E-16L,
00138   4.836645664124562546056389268546233577376E-13L,
00139   1.730945562285804805325011561498453013673E-10L,
00140   3.047976856147077889834905908605310585810E-8L,
00141   2.855227609107969710407464739188141162386E-6L,
00142   1.439362407936705484122143713643023998457E-4L,
00143   3.774489768532936551500999699815873422073E-3L,
00144   4.723962172984642566142399678920790598426E-2L,
00145   2.359289678988743939925017240478818248735E-1L,
00146   3.032580002220628812728954785118117124520E-1L,
00147 };
00148 #define NP16_ID 9
00149 static const long double P16_ID[NP16_ID + 1] = {
00150   4.389268795186898018132945193912677177553E-15L,
00151   4.132671824807454334388868363256830961655E-12L,
00152   1.482133328179508835835963635130894413136E-9L,
00153   2.618941412861122118906353737117067376236E-7L,
00154   2.467854246740858470815714426201888034270E-5L,
00155   1.257192927368839847825938545925340230490E-3L,
00156   3.362739031941574274949719324644120720341E-2L,
00157   4.384458231338934105875343439265370178858E-1L,
00158   2.412830809841095249170909628197264854651E0L,
00159   4.176078204111348059102962617368214856874E0L,
00160  /* 1.000000000000000000000000000000000000000E0 */
00161 };
00162 
00163 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
00164     0.0625 <= 1/x <= 0.125
00165     Peak relative error 1.9e-36  */
00166 #define NP8_16N 11
00167 static const long double P8_16N[NP8_16N + 1] = {
00168   2.984612480763362345647303274082071598135E-16L,
00169   1.923651877544126103941232173085475682334E-13L,
00170   4.881258879388869396043760693256024307743E-11L,
00171   6.368866572475045408480898921866869811889E-9L,
00172   4.684818344104910450523906967821090796737E-7L,
00173   2.005177298271593587095982211091300382796E-5L,
00174   4.979808067163957634120681477207147536182E-4L,
00175   6.946005761642579085284689047091173581127E-3L,
00176   5.074601112955765012750207555985299026204E-2L,
00177   1.698599455896180893191766195194231825379E-1L,
00178   1.957536905259237627737222775573623779638E-1L,
00179   2.991314703282528370270179989044994319374E-2L,
00180 };
00181 #define NP8_16D 10
00182 static const long double P8_16D[NP8_16D + 1] = {
00183   2.546869316918069202079580939942463010937E-15L,
00184   1.644650111942455804019788382157745229955E-12L,
00185   4.185430770291694079925607420808011147173E-10L,
00186   5.485331966975218025368698195861074143153E-8L,
00187   4.062884421686912042335466327098932678905E-6L,
00188   1.758139661060905948870523641319556816772E-4L,
00189   4.445143889306356207566032244985607493096E-3L,
00190   6.391901016293512632765621532571159071158E-2L,
00191   4.933040207519900471177016015718145795434E-1L,
00192   1.839144086168947712971630337250761842976E0L,
00193   2.715120873995490920415616716916149586579E0L,
00194  /* 1.000000000000000000000000000000000000000E0 */
00195 };
00196 
00197 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
00198   0.125 <= 1/x <= 0.1875
00199   Peak relative error 1.3e-36  */
00200 #define NP5_8N 10
00201 static const long double P5_8N[NP5_8N + 1] = {
00202   2.837678373978003452653763806968237227234E-12L,
00203   9.726641165590364928442128579282742354806E-10L,
00204   1.284408003604131382028112171490633956539E-7L,
00205   8.524624695868291291250573339272194285008E-6L,
00206   3.111516908953172249853673787748841282846E-4L,
00207   6.423175156126364104172801983096596409176E-3L,
00208   7.430220589989104581004416356260692450652E-2L,
00209   4.608315409833682489016656279567605536619E-1L,
00210   1.396870223510964882676225042258855977512E0L,
00211   1.718500293904122365894630460672081526236E0L,
00212   5.465927698800862172307352821870223855365E-1L
00213 };
00214 #define NP5_8D 10
00215 static const long double P5_8D[NP5_8D + 1] = {
00216   2.421485545794616609951168511612060482715E-11L,
00217   8.329862750896452929030058039752327232310E-9L,
00218   1.106137992233383429630592081375289010720E-6L,
00219   7.405786153760681090127497796448503306939E-5L,
00220   2.740364785433195322492093333127633465227E-3L,
00221   5.781246470403095224872243564165254652198E-2L,
00222   6.927711353039742469918754111511109983546E-1L,
00223   4.558679283460430281188304515922826156690E0L,
00224   1.534468499844879487013168065728837900009E1L,
00225   2.313927430889218597919624843161569422745E1L,
00226   1.194506341319498844336768473218382828637E1L,
00227  /* 1.000000000000000000000000000000000000000E0 */
00228 };
00229 
00230 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
00231    Peak relative error 1.4e-36
00232    0.1875 <= 1/x <= 0.25  */
00233 #define NP4_5N 10
00234 static const long double P4_5N[NP4_5N + 1] = {
00235   1.846029078268368685834261260420933914621E-10L,
00236   3.916295939611376119377869680335444207768E-8L,
00237   3.122158792018920627984597530935323997312E-6L,
00238   1.218073444893078303994045653603392272450E-4L,
00239   2.536420827983485448140477159977981844883E-3L,
00240   2.883011322006690823959367922241169171315E-2L,
00241   1.755255190734902907438042414495469810830E-1L,
00242   5.379317079922628599870898285488723736599E-1L,
00243   7.284904050194300773890303361501726561938E-1L,
00244   3.270110346613085348094396323925000362813E-1L,
00245   1.804473805689725610052078464951722064757E-2L,
00246 };
00247 #define NP4_5D 9
00248 static const long double P4_5D[NP4_5D + 1] = {
00249   1.575278146806816970152174364308980863569E-9L,
00250   3.361289173657099516191331123405675054321E-7L,
00251   2.704692281550877810424745289838790693708E-5L,
00252   1.070854930483999749316546199273521063543E-3L,
00253   2.282373093495295842598097265627962125411E-2L,
00254   2.692025460665354148328762368240343249830E-1L,
00255   1.739892942593664447220951225734811133759E0L,
00256   5.890727576752230385342377570386657229324E0L,
00257   9.517442287057841500750256954117735128153E0L,
00258   6.100616353935338240775363403030137736013E0L,
00259  /* 1.000000000000000000000000000000000000000E0 */
00260 };
00261 
00262 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
00263    Peak relative error 3.0e-36
00264    0.25 <= 1/x <= 0.3125  */
00265 #define NP3r2_4N 9
00266 static const long double P3r2_4N[NP3r2_4N + 1] = {
00267   8.240803130988044478595580300846665863782E-8L,
00268   1.179418958381961224222969866406483744580E-5L,
00269   6.179787320956386624336959112503824397755E-4L,
00270   1.540270833608687596420595830747166658383E-2L,
00271   1.983904219491512618376375619598837355076E-1L,
00272   1.341465722692038870390470651608301155565E0L,
00273   4.617865326696612898792238245990854646057E0L,
00274   7.435574801812346424460233180412308000587E0L,
00275   4.671327027414635292514599201278557680420E0L,
00276   7.299530852495776936690976966995187714739E-1L,
00277 };
00278 #define NP3r2_4D 9
00279 static const long double P3r2_4D[NP3r2_4D + 1] = {
00280   7.032152009675729604487575753279187576521E-7L,
00281   1.015090352324577615777511269928856742848E-4L,
00282   5.394262184808448484302067955186308730620E-3L,
00283   1.375291438480256110455809354836988584325E-1L,
00284   1.836247144461106304788160919310404376670E0L,
00285   1.314378564254376655001094503090935880349E1L,
00286   4.957184590465712006934452500894672343488E1L,
00287   9.287394244300647738855415178790263465398E1L,
00288   7.652563275535900609085229286020552768399E1L,
00289   2.147042473003074533150718117770093209096E1L,
00290  /* 1.000000000000000000000000000000000000000E0 */
00291 };
00292 
00293 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
00294    Peak relative error 1.0e-35
00295    0.3125 <= 1/x <= 0.375  */
00296 #define NP2r7_3r2N 9
00297 static const long double P2r7_3r2N[NP2r7_3r2N + 1] = {
00298   4.599033469240421554219816935160627085991E-7L,
00299   4.665724440345003914596647144630893997284E-5L,
00300   1.684348845667764271596142716944374892756E-3L,
00301   2.802446446884455707845985913454440176223E-2L,
00302   2.321937586453963310008279956042545173930E-1L,
00303   9.640277413988055668692438709376437553804E-1L,
00304   1.911021064710270904508663334033003246028E0L,
00305   1.600811610164341450262992138893970224971E0L,
00306   4.266299218652587901171386591543457861138E-1L,
00307   1.316470424456061252962568223251247207325E-2L,
00308 };
00309 #define NP2r7_3r2D 8
00310 static const long double P2r7_3r2D[NP2r7_3r2D + 1] = {
00311   3.924508608545520758883457108453520099610E-6L,
00312   4.029707889408829273226495756222078039823E-4L,
00313   1.484629715787703260797886463307469600219E-2L,
00314   2.553136379967180865331706538897231588685E-1L,
00315   2.229457223891676394409880026887106228740E0L,
00316   1.005708903856384091956550845198392117318E1L,
00317   2.277082659664386953166629360352385889558E1L,
00318   2.384726835193630788249826630376533988245E1L,
00319   9.700989749041320895890113781610939632410E0L,
00320  /* 1.000000000000000000000000000000000000000E0 */
00321 };
00322 
00323 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
00324    Peak relative error 1.7e-36
00325    0.3125 <= 1/x <= 0.4375  */
00326 #define NP2r3_2r7N 9
00327 static const long double P2r3_2r7N[NP2r3_2r7N + 1] = {
00328   3.916766777108274628543759603786857387402E-6L,
00329   3.212176636756546217390661984304645137013E-4L,
00330   9.255768488524816445220126081207248947118E-3L,
00331   1.214853146369078277453080641911700735354E-1L,
00332   7.855163309847214136198449861311404633665E-1L,
00333   2.520058073282978403655488662066019816540E0L,
00334   3.825136484837545257209234285382183711466E0L,
00335   2.432569427554248006229715163865569506873E0L,
00336   4.877934835018231178495030117729800489743E-1L,
00337   1.109902737860249670981355149101343427885E-2L,
00338 };
00339 #define NP2r3_2r7D 8
00340 static const long double P2r3_2r7D[NP2r3_2r7D + 1] = {
00341   3.342307880794065640312646341190547184461E-5L,
00342   2.782182891138893201544978009012096558265E-3L,
00343   8.221304931614200702142049236141249929207E-2L,
00344   1.123728246291165812392918571987858010949E0L,
00345   7.740482453652715577233858317133423434590E0L,
00346   2.737624677567945952953322566311201919139E1L,
00347   4.837181477096062403118304137851260715475E1L,
00348   3.941098643468580791437772701093795299274E1L,
00349   1.245821247166544627558323920382547533630E1L,
00350  /* 1.000000000000000000000000000000000000000E0 */
00351 };
00352 
00353 /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
00354    Peak relative error 1.7e-35
00355    0.4375 <= 1/x <= 0.5  */
00356 #define NP2_2r3N 8
00357 static const long double P2_2r3N[NP2_2r3N + 1] = {
00358   3.397930802851248553545191160608731940751E-4L,
00359   2.104020902735482418784312825637833698217E-2L,
00360   4.442291771608095963935342749477836181939E-1L,
00361   4.131797328716583282869183304291833754967E0L,
00362   1.819920169779026500146134832455189917589E1L,
00363   3.781779616522937565300309684282401791291E1L,
00364   3.459605449728864218972931220783543410347E1L,
00365   1.173594248397603882049066603238568316561E1L,
00366   9.455702270242780642835086549285560316461E-1L,
00367 };
00368 #define NP2_2r3D 8
00369 static const long double P2_2r3D[NP2_2r3D + 1] = {
00370   2.899568897241432883079888249845707400614E-3L,
00371   1.831107138190848460767699919531132426356E-1L,
00372   3.999350044057883839080258832758908825165E0L,
00373   3.929041535867957938340569419874195303712E1L,
00374   1.884245613422523323068802689915538908291E2L,
00375   4.461469948819229734353852978424629815929E2L,
00376   5.004998753999796821224085972610636347903E2L,
00377   2.386342520092608513170837883757163414100E2L,
00378   3.791322528149347975999851588922424189957E1L,
00379  /* 1.000000000000000000000000000000000000000E0 */
00380 };
00381 
00382 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
00383    Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
00384    Peak relative error 8.0e-36
00385    0 <= 1/x <= .0625  */
00386 #define NQ16_IN 10
00387 static const long double Q16_IN[NQ16_IN + 1] = {
00388   -3.917420835712508001321875734030357393421E-18L,
00389   -4.440311387483014485304387406538069930457E-15L,
00390   -1.951635424076926487780929645954007139616E-12L,
00391   -4.318256438421012555040546775651612810513E-10L,
00392   -5.231244131926180765270446557146989238020E-8L,
00393   -3.540072702902043752460711989234732357653E-6L,
00394   -1.311017536555269966928228052917534882984E-4L,
00395   -2.495184669674631806622008769674827575088E-3L,
00396   -2.141868222987209028118086708697998506716E-2L,
00397   -6.184031415202148901863605871197272650090E-2L,
00398   -1.922298704033332356899546792898156493887E-2L,
00399 };
00400 #define NQ16_ID 9
00401 static const long double Q16_ID[NQ16_ID + 1] = {
00402   3.820418034066293517479619763498400162314E-17L,
00403   4.340702810799239909648911373329149354911E-14L,
00404   1.914985356383416140706179933075303538524E-11L,
00405   4.262333682610888819476498617261895474330E-9L,
00406   5.213481314722233980346462747902942182792E-7L,
00407   3.585741697694069399299005316809954590558E-5L,
00408   1.366513429642842006385029778105539457546E-3L,
00409   2.745282599850704662726337474371355160594E-2L,
00410   2.637644521611867647651200098449903330074E-1L,
00411   1.006953426110765984590782655598680488746E0L,
00412  /* 1.000000000000000000000000000000000000000E0 */
00413  };
00414 
00415 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
00416    Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
00417    Peak relative error 1.9e-36
00418    0.0625 <= 1/x <= 0.125  */
00419 #define NQ8_16N 11
00420 static const long double Q8_16N[NQ8_16N + 1] = {
00421   -2.028630366670228670781362543615221542291E-17L,
00422   -1.519634620380959966438130374006858864624E-14L,
00423   -4.540596528116104986388796594639405114524E-12L,
00424   -7.085151756671466559280490913558388648274E-10L,
00425   -6.351062671323970823761883833531546885452E-8L,
00426   -3.390817171111032905297982523519503522491E-6L,
00427   -1.082340897018886970282138836861233213972E-4L,
00428   -2.020120801187226444822977006648252379508E-3L,
00429   -2.093169910981725694937457070649605557555E-2L,
00430   -1.092176538874275712359269481414448063393E-1L,
00431   -2.374790947854765809203590474789108718733E-1L,
00432   -1.365364204556573800719985118029601401323E-1L,
00433 };
00434 #define NQ8_16D 11
00435 static const long double Q8_16D[NQ8_16D + 1] = {
00436   1.978397614733632533581207058069628242280E-16L,
00437   1.487361156806202736877009608336766720560E-13L,
00438   4.468041406888412086042576067133365913456E-11L,
00439   7.027822074821007443672290507210594648877E-9L,
00440   6.375740580686101224127290062867976007374E-7L,
00441   3.466887658320002225888644977076410421940E-5L,
00442   1.138625640905289601186353909213719596986E-3L,
00443   2.224470799470414663443449818235008486439E-2L,
00444   2.487052928527244907490589787691478482358E-1L,
00445   1.483927406564349124649083853892380899217E0L,
00446   4.182773513276056975777258788903489507705E0L,
00447   4.419665392573449746043880892524360870944E0L,
00448  /* 1.000000000000000000000000000000000000000E0 */
00449 };
00450 
00451 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
00452    Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
00453    Peak relative error 1.5e-35
00454    0.125 <= 1/x <= 0.1875  */
00455 #define NQ5_8N 10
00456 static const long double Q5_8N[NQ5_8N + 1] = {
00457   -3.656082407740970534915918390488336879763E-13L,
00458   -1.344660308497244804752334556734121771023E-10L,
00459   -1.909765035234071738548629788698150760791E-8L,
00460   -1.366668038160120210269389551283666716453E-6L,
00461   -5.392327355984269366895210704976314135683E-5L,
00462   -1.206268245713024564674432357634540343884E-3L,
00463   -1.515456784370354374066417703736088291287E-2L,
00464   -1.022454301137286306933217746545237098518E-1L,
00465   -3.373438906472495080504907858424251082240E-1L,
00466   -4.510782522110845697262323973549178453405E-1L,
00467   -1.549000892545288676809660828213589804884E-1L,
00468 };
00469 #define NQ5_8D 10
00470 static const long double Q5_8D[NQ5_8D + 1] = {
00471   3.565550843359501079050699598913828460036E-12L,
00472   1.321016015556560621591847454285330528045E-9L,
00473   1.897542728662346479999969679234270605975E-7L,
00474   1.381720283068706710298734234287456219474E-5L,
00475   5.599248147286524662305325795203422873725E-4L,
00476   1.305442352653121436697064782499122164843E-2L,
00477   1.750234079626943298160445750078631894985E-1L,
00478   1.311420542073436520965439883806946678491E0L,
00479   5.162757689856842406744504211089724926650E0L,
00480   9.527760296384704425618556332087850581308E0L,
00481   6.604648207463236667912921642545100248584E0L,
00482  /* 1.000000000000000000000000000000000000000E0 */
00483 };
00484 
00485 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
00486    Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
00487    Peak relative error 1.3e-35
00488    0.1875 <= 1/x <= 0.25  */
00489 #define NQ4_5N 10
00490 static const long double Q4_5N[NQ4_5N + 1] = {
00491   -4.079513568708891749424783046520200903755E-11L,
00492   -9.326548104106791766891812583019664893311E-9L,
00493   -8.016795121318423066292906123815687003356E-7L,
00494   -3.372350544043594415609295225664186750995E-5L,
00495   -7.566238665947967882207277686375417983917E-4L,
00496   -9.248861580055565402130441618521591282617E-3L,
00497   -6.033106131055851432267702948850231270338E-2L,
00498   -1.966908754799996793730369265431584303447E-1L,
00499   -2.791062741179964150755788226623462207560E-1L,
00500   -1.255478605849190549914610121863534191666E-1L,
00501   -4.320429862021265463213168186061696944062E-3L,
00502 };
00503 #define NQ4_5D 9
00504 static const long double Q4_5D[NQ4_5D + 1] = {
00505   3.978497042580921479003851216297330701056E-10L,
00506   9.203304163828145809278568906420772246666E-8L,
00507   8.059685467088175644915010485174545743798E-6L,
00508   3.490187375993956409171098277561669167446E-4L,
00509   8.189109654456872150100501732073810028829E-3L,
00510   1.072572867311023640958725265762483033769E-1L,
00511   7.790606862409960053675717185714576937994E-1L,
00512   3.016049768232011196434185423512777656328E0L,
00513   5.722963851442769787733717162314477949360E0L,
00514   4.510527838428473279647251350931380867663E0L,
00515  /* 1.000000000000000000000000000000000000000E0 */
00516 };
00517 
00518 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
00519    Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
00520    Peak relative error 2.1e-35
00521    0.25 <= 1/x <= 0.3125  */
00522 #define NQ3r2_4N 9
00523 static const long double Q3r2_4N[NQ3r2_4N + 1] = {
00524   -1.087480809271383885936921889040388133627E-8L,
00525   -1.690067828697463740906962973479310170932E-6L,
00526   -9.608064416995105532790745641974762550982E-5L,
00527   -2.594198839156517191858208513873961837410E-3L,
00528   -3.610954144421543968160459863048062977822E-2L,
00529   -2.629866798251843212210482269563961685666E-1L,
00530   -9.709186825881775885917984975685752956660E-1L,
00531   -1.667521829918185121727268867619982417317E0L,
00532   -1.109255082925540057138766105229900943501E0L,
00533   -1.812932453006641348145049323713469043328E-1L,
00534 };
00535 #define NQ3r2_4D 9
00536 static const long double Q3r2_4D[NQ3r2_4D + 1] = {
00537   1.060552717496912381388763753841473407026E-7L,
00538   1.676928002024920520786883649102388708024E-5L,
00539   9.803481712245420839301400601140812255737E-4L,
00540   2.765559874262309494758505158089249012930E-2L,
00541   4.117921827792571791298862613287549140706E-1L,
00542   3.323769515244751267093378361930279161413E0L,
00543   1.436602494405814164724810151689705353670E1L,
00544   3.163087869617098638064881410646782408297E1L,
00545   3.198181264977021649489103980298349589419E1L,
00546   1.203649258862068431199471076202897823272E1L,
00547  /* 1.000000000000000000000000000000000000000E0  */
00548 };
00549 
00550 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
00551    Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
00552    Peak relative error 1.6e-36
00553    0.3125 <= 1/x <= 0.375  */
00554 #define NQ2r7_3r2N 9
00555 static const long double Q2r7_3r2N[NQ2r7_3r2N + 1] = {
00556   -1.723405393982209853244278760171643219530E-7L,
00557   -2.090508758514655456365709712333460087442E-5L,
00558   -9.140104013370974823232873472192719263019E-4L,
00559   -1.871349499990714843332742160292474780128E-2L,
00560   -1.948930738119938669637865956162512983416E-1L,
00561   -1.048764684978978127908439526343174139788E0L,
00562   -2.827714929925679500237476105843643064698E0L,
00563   -3.508761569156476114276988181329773987314E0L,
00564   -1.669332202790211090973255098624488308989E0L,
00565   -1.930796319299022954013840684651016077770E-1L,
00566 };
00567 #define NQ2r7_3r2D 9
00568 static const long double Q2r7_3r2D[NQ2r7_3r2D + 1] = {
00569   1.680730662300831976234547482334347983474E-6L,
00570   2.084241442440551016475972218719621841120E-4L,
00571   9.445316642108367479043541702688736295579E-3L,
00572   2.044637889456631896650179477133252184672E-1L,
00573   2.316091982244297350829522534435350078205E0L,
00574   1.412031891783015085196708811890448488865E1L,
00575   4.583830154673223384837091077279595496149E1L,
00576   7.549520609270909439885998474045974122261E1L,
00577   5.697605832808113367197494052388203310638E1L,
00578   1.601496240876192444526383314589371686234E1L,
00579   /* 1.000000000000000000000000000000000000000E0 */
00580 };
00581 
00582 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
00583    Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
00584    Peak relative error 9.5e-36
00585    0.375 <= 1/x <= 0.4375  */
00586 #define NQ2r3_2r7N 9
00587 static const long double Q2r3_2r7N[NQ2r3_2r7N + 1] = {
00588   -8.603042076329122085722385914954878953775E-7L,
00589   -7.701746260451647874214968882605186675720E-5L,
00590   -2.407932004380727587382493696877569654271E-3L,
00591   -3.403434217607634279028110636919987224188E-2L,
00592   -2.348707332185238159192422084985713102877E-1L,
00593   -7.957498841538254916147095255700637463207E-1L,
00594   -1.258469078442635106431098063707934348577E0L,
00595   -8.162415474676345812459353639449971369890E-1L,
00596   -1.581783890269379690141513949609572806898E-1L,
00597   -1.890595651683552228232308756569450822905E-3L,
00598 };
00599 #define NQ2r3_2r7D 8
00600 static const long double Q2r3_2r7D[NQ2r3_2r7D + 1] = {
00601   8.390017524798316921170710533381568175665E-6L,
00602   7.738148683730826286477254659973968763659E-4L,
00603   2.541480810958665794368759558791634341779E-2L,
00604   3.878879789711276799058486068562386244873E-1L,
00605   3.003783779325811292142957336802456109333E0L,
00606   1.206480374773322029883039064575464497400E1L,
00607   2.458414064785315978408974662900438351782E1L,
00608   2.367237826273668567199042088835448715228E1L,
00609   9.231451197519171090875569102116321676763E0L,
00610  /* 1.000000000000000000000000000000000000000E0 */
00611 };
00612 
00613 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
00614    Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
00615    Peak relative error 1.4e-36
00616    0.4375 <= 1/x <= 0.5  */
00617 #define NQ2_2r3N 9
00618 static const long double Q2_2r3N[NQ2_2r3N + 1] = {
00619   -5.552507516089087822166822364590806076174E-6L,
00620   -4.135067659799500521040944087433752970297E-4L,
00621   -1.059928728869218962607068840646564457980E-2L,
00622   -1.212070036005832342565792241385459023801E-1L,
00623   -6.688350110633603958684302153362735625156E-1L,
00624   -1.793587878197360221340277951304429821582E0L,
00625   -2.225407682237197485644647380483725045326E0L,
00626   -1.123402135458940189438898496348239744403E0L,
00627   -1.679187241566347077204805190763597299805E-1L,
00628   -1.458550613639093752909985189067233504148E-3L,
00629 };
00630 #define NQ2_2r3D 8
00631 static const long double Q2_2r3D[NQ2_2r3D + 1] = {
00632   5.415024336507980465169023996403597916115E-5L,
00633   4.179246497380453022046357404266022870788E-3L,
00634   1.136306384261959483095442402929502368598E-1L,
00635   1.422640343719842213484515445393284072830E0L,
00636   8.968786703393158374728850922289204805764E0L,
00637   2.914542473339246127533384118781216495934E1L,
00638   4.781605421020380669870197378210457054685E1L,
00639   3.693865837171883152382820584714795072937E1L,
00640   1.153220502744204904763115556224395893076E1L,
00641   /* 1.000000000000000000000000000000000000000E0 */
00642 };
00643 
00644 
00645 /* Evaluate P[n] x^n  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
00646 
00647 static long double
00648 neval (long double x, const long double *p, int n)
00649 {
00650   long double y;
00651 
00652   p += n;
00653   y = *p--;
00654   do
00655     {
00656       y = y * x + *p--;
00657     }
00658   while (--n > 0);
00659   return y;
00660 }
00661 
00662 
00663 /* Evaluate x^n+1  +  P[n] x^(n)  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
00664 
00665 static long double
00666 deval (long double x, const long double *p, int n)
00667 {
00668   long double y;
00669 
00670   p += n;
00671   y = x + *p--;
00672   do
00673     {
00674       y = y * x + *p--;
00675     }
00676   while (--n > 0);
00677   return y;
00678 }
00679 
00680 
00681 /* Bessel function of the first kind, order one.  */
00682 
00683 long double
00684 __ieee754_j1l (long double x)
00685 {
00686   long double xx, xinv, z, p, q, c, s, cc, ss;
00687 
00688   if (! __finitel (x))
00689     {
00690       if (x != x)
00691        return x;
00692       else
00693        return 0.0L;
00694     }
00695   if (x == 0.0L)
00696     return x;
00697   xx = fabsl (x);
00698   if (xx <= 2.0L)
00699     {
00700       /* 0 <= x <= 2 */
00701       z = xx * xx;
00702       p = xx * z * neval (z, J0_2N, NJ0_2N) / deval (z, J0_2D, NJ0_2D);
00703       p += 0.5L * xx;
00704       if (x < 0)
00705        p = -p;
00706       return p;
00707     }
00708 
00709   xinv = 1.0L / xx;
00710   z = xinv * xinv;
00711   if (xinv <= 0.25)
00712     {
00713       if (xinv <= 0.125)
00714        {
00715          if (xinv <= 0.0625)
00716            {
00717              p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID);
00718              q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID);
00719            }
00720          else
00721            {
00722              p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D);
00723              q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D);
00724            }
00725        }
00726       else if (xinv <= 0.1875)
00727        {
00728          p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D);
00729          q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D);
00730        }
00731       else
00732        {
00733          p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D);
00734          q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D);
00735        }
00736     }                       /* .25 */
00737   else /* if (xinv <= 0.5) */
00738     {
00739       if (xinv <= 0.375)
00740        {
00741          if (xinv <= 0.3125)
00742            {
00743              p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D);
00744              q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D);
00745            }
00746          else
00747            {
00748              p = neval (z, P2r7_3r2N, NP2r7_3r2N)
00749                 / deval (z, P2r7_3r2D, NP2r7_3r2D);
00750              q = neval (z, Q2r7_3r2N, NQ2r7_3r2N)
00751                 / deval (z, Q2r7_3r2D, NQ2r7_3r2D);
00752            }
00753        }
00754       else if (xinv <= 0.4375)
00755        {
00756          p = neval (z, P2r3_2r7N, NP2r3_2r7N)
00757              / deval (z, P2r3_2r7D, NP2r3_2r7D);
00758          q = neval (z, Q2r3_2r7N, NQ2r3_2r7N)
00759              / deval (z, Q2r3_2r7D, NQ2r3_2r7D);
00760        }
00761       else
00762        {
00763          p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
00764          q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
00765        }
00766     }
00767   p = 1.0L + z * p;
00768   q = z * q;
00769   q = q * xinv + 0.375L * xinv;
00770   /* X = x - 3 pi/4
00771      cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
00772      = 1/sqrt(2) * (-cos(x) + sin(x))
00773      sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
00774      = -1/sqrt(2) * (sin(x) + cos(x))
00775      cf. Fdlibm.  */
00776   __sincosl (xx, &s, &c);
00777   ss = -s - c;
00778   cc = s - c;
00779   z = __cosl (xx + xx);
00780   if ((s * c) > 0)
00781     cc = z / ss;
00782   else
00783     ss = z / cc;
00784   z = ONEOSQPI * (p * cc - q * ss) / __ieee754_sqrtl (xx);
00785   if (x < 0)
00786     z = -z;
00787   return z;
00788 }
00789 
00790 
00791 /* Y1(x) = 2/pi * (log(x) * J1(x) - 1/x) + x R(x^2)
00792    Peak relative error 6.2e-38
00793    0 <= x <= 2   */
00794 #define NY0_2N 7
00795 static long double Y0_2N[NY0_2N + 1] = {
00796   -6.804415404830253804408698161694720833249E19L,
00797   1.805450517967019908027153056150465849237E19L,
00798   -8.065747497063694098810419456383006737312E17L,
00799   1.401336667383028259295830955439028236299E16L,
00800   -1.171654432898137585000399489686629680230E14L,
00801   5.061267920943853732895341125243428129150E11L,
00802   -1.096677850566094204586208610960870217970E9L,
00803   9.541172044989995856117187515882879304461E5L,
00804 };
00805 #define NY0_2D 7
00806 static long double Y0_2D[NY0_2D + 1] = {
00807   3.470629591820267059538637461549677594549E20L,
00808   4.120796439009916326855848107545425217219E18L,
00809   2.477653371652018249749350657387030814542E16L,
00810   9.954678543353888958177169349272167762797E13L,
00811   2.957927997613630118216218290262851197754E11L,
00812   6.748421382188864486018861197614025972118E8L,
00813   1.173453425218010888004562071020305709319E6L,
00814   1.450335662961034949894009554536003377187E3L,
00815   /* 1.000000000000000000000000000000000000000E0 */
00816 };
00817 
00818 
00819 /* Bessel function of the second kind, order one.  */
00820 
00821 long double
00822 __ieee754_y1l (long double x)
00823 {
00824   long double xx, xinv, z, p, q, c, s, cc, ss;
00825 
00826   if (! __finitel (x))
00827     {
00828       if (x != x)
00829        return x;
00830       else
00831        return 0.0L;
00832     }
00833   if (x <= 0.0L)
00834     {
00835       if (x < 0.0L)
00836        return (zero / (zero * x));
00837       return -HUGE_VALL + x;
00838     }
00839   xx = fabsl (x);
00840   if (xx <= 2.0L)
00841     {
00842       /* 0 <= x <= 2 */
00843       z = xx * xx;
00844       p = xx * neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D);
00845       p = -TWOOPI / xx + p;
00846       p = TWOOPI * __ieee754_logl (x) * __ieee754_j1l (x) + p;
00847       return p;
00848     }
00849 
00850   xinv = 1.0L / xx;
00851   z = xinv * xinv;
00852   if (xinv <= 0.25)
00853     {
00854       if (xinv <= 0.125)
00855        {
00856          if (xinv <= 0.0625)
00857            {
00858              p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID);
00859              q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID);
00860            }
00861          else
00862            {
00863              p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D);
00864              q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D);
00865            }
00866        }
00867       else if (xinv <= 0.1875)
00868        {
00869          p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D);
00870          q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D);
00871        }
00872       else
00873        {
00874          p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D);
00875          q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D);
00876        }
00877     }                       /* .25 */
00878   else /* if (xinv <= 0.5) */
00879     {
00880       if (xinv <= 0.375)
00881        {
00882          if (xinv <= 0.3125)
00883            {
00884              p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D);
00885              q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D);
00886            }
00887          else
00888            {
00889              p = neval (z, P2r7_3r2N, NP2r7_3r2N)
00890                 / deval (z, P2r7_3r2D, NP2r7_3r2D);
00891              q = neval (z, Q2r7_3r2N, NQ2r7_3r2N)
00892                 / deval (z, Q2r7_3r2D, NQ2r7_3r2D);
00893            }
00894        }
00895       else if (xinv <= 0.4375)
00896        {
00897          p = neval (z, P2r3_2r7N, NP2r3_2r7N)
00898              / deval (z, P2r3_2r7D, NP2r3_2r7D);
00899          q = neval (z, Q2r3_2r7N, NQ2r3_2r7N)
00900              / deval (z, Q2r3_2r7D, NQ2r3_2r7D);
00901        }
00902       else
00903        {
00904          p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
00905          q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
00906        }
00907     }
00908   p = 1.0L + z * p;
00909   q = z * q;
00910   q = q * xinv + 0.375L * xinv;
00911   /* X = x - 3 pi/4
00912      cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
00913      = 1/sqrt(2) * (-cos(x) + sin(x))
00914      sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
00915      = -1/sqrt(2) * (sin(x) + cos(x))
00916      cf. Fdlibm.  */
00917   __sincosl (xx, &s, &c);
00918   ss = -s - c;
00919   cc = s - c;
00920   z = __cosl (xx + xx);
00921   if ((s * c) > 0)
00922     cc = z / ss;
00923   else
00924     ss = z / cc;
00925   z = ONEOSQPI * (p * ss + q * cc) / __ieee754_sqrtl (xx);
00926   return z;
00927 }