Back to index

glibc  2.9
e_j0l.c
Go to the documentation of this file.
00001 /*                                               j0l.c
00002  *
00003  *     Bessel function of order zero
00004  *
00005  *
00006  *
00007  * SYNOPSIS:
00008  *
00009  * long double x, y, j0l();
00010  *
00011  * y = j0l( x );
00012  *
00013  *
00014  *
00015  * DESCRIPTION:
00016  *
00017  * Returns Bessel function of first kind, order zero 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
00021  * is J0(x) = 1 - x^2 / 4 + x^4 R(x^2)
00022  * The second interval is further partitioned into eight equal segments
00023  * of 1/x.
00024  *
00025  * J0(x) = sqrt(2/(pi x)) (P0(x) cos(X) - Q0(x) sin(X)),
00026  * X = x - pi/4,
00027  *
00028  * and the auxiliary functions are given by
00029  *
00030  * J0(x)cos(X) + Y0(x)sin(X) = sqrt( 2/(pi x)) P0(x),
00031  * P0(x) = 1 + 1/x^2 R(1/x^2)
00032  *
00033  * Y0(x)cos(X) - J0(x)sin(X) = sqrt( 2/(pi x)) Q0(x),
00034  * Q0(x) = 1/x (-.125 + 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      1.7e-34      2.4e-35
00043  *
00044  *
00045  */
00046 
00047 /*                                               y0l.c
00048  *
00049  *     Bessel function of the second kind, order zero
00050  *
00051  *
00052  *
00053  * SYNOPSIS:
00054  *
00055  * double x, y, y0l();
00056  *
00057  * y = y0l( x );
00058  *
00059  *
00060  *
00061  * DESCRIPTION:
00062  *
00063  * Returns Bessel function of the second kind, of order
00064  * zero, of the argument.
00065  *
00066  * The approximation is the same as for J0(x), and
00067  * Y0(x) = sqrt(2/(pi x)) (P0(x) sin(X) + Q0(x) cos(X)).
00068  *
00069  * ACCURACY:
00070  *
00071  *  Absolute error, when y0(x) < 1; else relative error:
00072  *
00073  * arithmetic   domain     # trials      peak         rms
00074  *    IEEE      0, 30       100000      3.0e-34     2.7e-35
00075  *
00076  */
00077 
00078 /* Copyright 2001 by Stephen L. Moshier (moshier@na-net.ornl.gov).
00079 
00080     This library is free software; you can redistribute it and/or
00081     modify it under the terms of the GNU Lesser General Public
00082     License as published by the Free Software Foundation; either
00083     version 2.1 of the License, or (at your option) any later version.
00084 
00085     This library is distributed in the hope that it will be useful,
00086     but WITHOUT ANY WARRANTY; without even the implied warranty of
00087     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00088     Lesser General Public License for more details.
00089 
00090     You should have received a copy of the GNU Lesser General Public
00091     License along with this library; if not, write to the Free Software
00092     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
00093 
00094 #include "math.h"
00095 #include "math_private.h"
00096 
00097 /* 1 / sqrt(pi) */
00098 static const long double ONEOSQPI = 5.6418958354775628694807945156077258584405E-1L;
00099 /* 2 / pi */
00100 static const long double TWOOPI = 6.3661977236758134307553505349005744813784E-1L;
00101 static const long double zero = 0.0L;
00102 
00103 /* J0(x) = 1 - x^2/4 + x^2 x^2 R(x^2)
00104    Peak relative error 3.4e-37
00105    0 <= x <= 2  */
00106 #define NJ0_2N 6
00107 static const long double J0_2N[NJ0_2N + 1] = {
00108   3.133239376997663645548490085151484674892E16L,
00109  -5.479944965767990821079467311839107722107E14L,
00110   6.290828903904724265980249871997551894090E12L,
00111  -3.633750176832769659849028554429106299915E10L,
00112   1.207743757532429576399485415069244807022E8L,
00113  -2.107485999925074577174305650549367415465E5L,
00114   1.562826808020631846245296572935547005859E2L,
00115 };
00116 #define NJ0_2D 6
00117 static const long double J0_2D[NJ0_2D + 1] = {
00118   2.005273201278504733151033654496928968261E18L,
00119   2.063038558793221244373123294054149790864E16L,
00120   1.053350447931127971406896594022010524994E14L,
00121   3.496556557558702583143527876385508882310E11L,
00122   8.249114511878616075860654484367133976306E8L,
00123   1.402965782449571800199759247964242790589E6L,
00124   1.619910762853439600957801751815074787351E3L,
00125  /* 1.000000000000000000000000000000000000000E0 */
00126 };
00127 
00128 /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2),
00129    0 <= 1/x <= .0625
00130    Peak relative error 3.3e-36  */
00131 #define NP16_IN 9
00132 static const long double P16_IN[NP16_IN + 1] = {
00133   -1.901689868258117463979611259731176301065E-16L,
00134   -1.798743043824071514483008340803573980931E-13L,
00135   -6.481746687115262291873324132944647438959E-11L,
00136   -1.150651553745409037257197798528294248012E-8L,
00137   -1.088408467297401082271185599507222695995E-6L,
00138   -5.551996725183495852661022587879817546508E-5L,
00139   -1.477286941214245433866838787454880214736E-3L,
00140   -1.882877976157714592017345347609200402472E-2L,
00141   -9.620983176855405325086530374317855880515E-2L,
00142   -1.271468546258855781530458854476627766233E-1L,
00143 };
00144 #define NP16_ID 9
00145 static const long double P16_ID[NP16_ID + 1] = {
00146   2.704625590411544837659891569420764475007E-15L,
00147   2.562526347676857624104306349421985403573E-12L,
00148   9.259137589952741054108665570122085036246E-10L,
00149   1.651044705794378365237454962653430805272E-7L,
00150   1.573561544138733044977714063100859136660E-5L,
00151   8.134482112334882274688298469629884804056E-4L,
00152   2.219259239404080863919375103673593571689E-2L,
00153   2.976990606226596289580242451096393862792E-1L,
00154   1.713895630454693931742734911930937246254E0L,
00155   3.231552290717904041465898249160757368855E0L,
00156   /* 1.000000000000000000000000000000000000000E0 */
00157 };
00158 
00159 /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2)
00160     0.0625 <= 1/x <= 0.125
00161     Peak relative error 2.4e-35  */
00162 #define NP8_16N 10
00163 static const long double P8_16N[NP8_16N + 1] = {
00164   -2.335166846111159458466553806683579003632E-15L,
00165   -1.382763674252402720401020004169367089975E-12L,
00166   -3.192160804534716696058987967592784857907E-10L,
00167   -3.744199606283752333686144670572632116899E-8L,
00168   -2.439161236879511162078619292571922772224E-6L,
00169   -9.068436986859420951664151060267045346549E-5L,
00170   -1.905407090637058116299757292660002697359E-3L,
00171   -2.164456143936718388053842376884252978872E-2L,
00172   -1.212178415116411222341491717748696499966E-1L,
00173   -2.782433626588541494473277445959593334494E-1L,
00174   -1.670703190068873186016102289227646035035E-1L,
00175 };
00176 #define NP8_16D 10
00177 static const long double P8_16D[NP8_16D + 1] = {
00178   3.321126181135871232648331450082662856743E-14L,
00179   1.971894594837650840586859228510007703641E-11L,
00180   4.571144364787008285981633719513897281690E-9L,
00181   5.396419143536287457142904742849052402103E-7L,
00182   3.551548222385845912370226756036899901549E-5L,
00183   1.342353874566932014705609788054598013516E-3L,
00184   2.899133293006771317589357444614157734385E-2L,
00185   3.455374978185770197704507681491574261545E-1L,
00186   2.116616964297512311314454834712634820514E0L,
00187   5.850768316827915470087758636881584174432E0L,
00188   5.655273858938766830855753983631132928968E0L,
00189   /* 1.000000000000000000000000000000000000000E0 */
00190 };
00191 
00192 /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2)
00193   0.125 <= 1/x <= 0.1875
00194   Peak relative error 2.7e-35  */
00195 #define NP5_8N 10
00196 static const long double P5_8N[NP5_8N + 1] = {
00197   -1.270478335089770355749591358934012019596E-12L,
00198   -4.007588712145412921057254992155810347245E-10L,
00199   -4.815187822989597568124520080486652009281E-8L,
00200   -2.867070063972764880024598300408284868021E-6L,
00201   -9.218742195161302204046454768106063638006E-5L,
00202   -1.635746821447052827526320629828043529997E-3L,
00203   -1.570376886640308408247709616497261011707E-2L,
00204   -7.656484795303305596941813361786219477807E-2L,
00205   -1.659371030767513274944805479908858628053E-1L,
00206   -1.185340550030955660015841796219919804915E-1L,
00207   -8.920026499909994671248893388013790366712E-3L,
00208 };
00209 #define NP5_8D 9
00210 static const long double P5_8D[NP5_8D + 1] = {
00211   1.806902521016705225778045904631543990314E-11L,
00212   5.728502760243502431663549179135868966031E-9L,
00213   6.938168504826004255287618819550667978450E-7L,
00214   4.183769964807453250763325026573037785902E-5L,
00215   1.372660678476925468014882230851637878587E-3L,
00216   2.516452105242920335873286419212708961771E-2L,
00217   2.550502712902647803796267951846557316182E-1L,
00218   1.365861559418983216913629123778747617072E0L,
00219   3.523825618308783966723472468855042541407E0L,
00220   3.656365803506136165615111349150536282434E0L,
00221   /* 1.000000000000000000000000000000000000000E0 */
00222 };
00223 
00224 /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2)
00225    Peak relative error 3.5e-35
00226    0.1875 <= 1/x <= 0.25  */
00227 #define NP4_5N 9
00228 static const long double P4_5N[NP4_5N + 1] = {
00229   -9.791405771694098960254468859195175708252E-10L,
00230   -1.917193059944531970421626610188102836352E-7L,
00231   -1.393597539508855262243816152893982002084E-5L,
00232   -4.881863490846771259880606911667479860077E-4L,
00233   -8.946571245022470127331892085881699269853E-3L,
00234   -8.707474232568097513415336886103899434251E-2L,
00235   -4.362042697474650737898551272505525973766E-1L,
00236   -1.032712171267523975431451359962375617386E0L,
00237   -9.630502683169895107062182070514713702346E-1L,
00238   -2.251804386252969656586810309252357233320E-1L,
00239 };
00240 #define NP4_5D 9
00241 static const long double P4_5D[NP4_5D + 1] = {
00242   1.392555487577717669739688337895791213139E-8L,
00243   2.748886559120659027172816051276451376854E-6L,
00244   2.024717710644378047477189849678576659290E-4L,
00245   7.244868609350416002930624752604670292469E-3L,
00246   1.373631762292244371102989739300382152416E-1L,
00247   1.412298581400224267910294815260613240668E0L,
00248   7.742495637843445079276397723849017617210E0L,
00249   2.138429269198406512028307045259503811861E1L,
00250   2.651547684548423476506826951831712762610E1L,
00251   1.167499382465291931571685222882909166935E1L,
00252   /* 1.000000000000000000000000000000000000000E0 */
00253 };
00254 
00255 /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2)
00256    Peak relative error 2.3e-36
00257    0.25 <= 1/x <= 0.3125  */
00258 #define NP3r2_4N 9
00259 static const long double P3r2_4N[NP3r2_4N + 1] = {
00260   -2.589155123706348361249809342508270121788E-8L,
00261   -3.746254369796115441118148490849195516593E-6L,
00262   -1.985595497390808544622893738135529701062E-4L,
00263   -5.008253705202932091290132760394976551426E-3L,
00264   -6.529469780539591572179155511840853077232E-2L,
00265   -4.468736064761814602927408833818990271514E-1L,
00266   -1.556391252586395038089729428444444823380E0L,
00267   -2.533135309840530224072920725976994981638E0L,
00268   -1.605509621731068453869408718565392869560E0L,
00269   -2.518966692256192789269859830255724429375E-1L,
00270 };
00271 #define NP3r2_4D 9
00272 static const long double P3r2_4D[NP3r2_4D + 1] = {
00273   3.682353957237979993646169732962573930237E-7L,
00274   5.386741661883067824698973455566332102029E-5L,
00275   2.906881154171822780345134853794241037053E-3L,
00276   7.545832595801289519475806339863492074126E-2L,
00277   1.029405357245594877344360389469584526654E0L,
00278   7.565706120589873131187989560509757626725E0L,
00279   2.951172890699569545357692207898667665796E1L,
00280   5.785723537170311456298467310529815457536E1L,
00281   5.095621464598267889126015412522773474467E1L,
00282   1.602958484169953109437547474953308401442E1L,
00283   /* 1.000000000000000000000000000000000000000E0 */
00284 };
00285 
00286 /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2)
00287    Peak relative error 1.0e-35
00288    0.3125 <= 1/x <= 0.375  */
00289 #define NP2r7_3r2N 9
00290 static const long double P2r7_3r2N[NP2r7_3r2N + 1] = {
00291   -1.917322340814391131073820537027234322550E-7L,
00292   -1.966595744473227183846019639723259011906E-5L,
00293   -7.177081163619679403212623526632690465290E-4L,
00294   -1.206467373860974695661544653741899755695E-2L,
00295   -1.008656452188539812154551482286328107316E-1L,
00296   -4.216016116408810856620947307438823892707E-1L,
00297   -8.378631013025721741744285026537009814161E-1L,
00298   -6.973895635309960850033762745957946272579E-1L,
00299   -1.797864718878320770670740413285763554812E-1L,
00300   -4.098025357743657347681137871388402849581E-3L,
00301 };
00302 #define NP2r7_3r2D 8
00303 static const long double P2r7_3r2D[NP2r7_3r2D + 1] = {
00304   2.726858489303036441686496086962545034018E-6L,
00305   2.840430827557109238386808968234848081424E-4L,
00306   1.063826772041781947891481054529454088832E-2L,
00307   1.864775537138364773178044431045514405468E-1L,
00308   1.665660052857205170440952607701728254211E0L,
00309   7.723745889544331153080842168958348568395E0L,
00310   1.810726427571829798856428548102077799835E1L,
00311   1.986460672157794440666187503833545388527E1L,
00312   8.645503204552282306364296517220055815488E0L,
00313   /* 1.000000000000000000000000000000000000000E0 */
00314 };
00315 
00316 /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2)
00317    Peak relative error 1.3e-36
00318    0.3125 <= 1/x <= 0.4375  */
00319 #define NP2r3_2r7N 9
00320 static const long double P2r3_2r7N[NP2r3_2r7N + 1] = {
00321   -1.594642785584856746358609622003310312622E-6L,
00322   -1.323238196302221554194031733595194539794E-4L,
00323   -3.856087818696874802689922536987100372345E-3L,
00324   -5.113241710697777193011470733601522047399E-2L,
00325   -3.334229537209911914449990372942022350558E-1L,
00326   -1.075703518198127096179198549659283422832E0L,
00327   -1.634174803414062725476343124267110981807E0L,
00328   -1.030133247434119595616826842367268304880E0L,
00329   -1.989811539080358501229347481000707289391E-1L,
00330   -3.246859189246653459359775001466924610236E-3L,
00331 };
00332 #define NP2r3_2r7D 8
00333 static const long double P2r3_2r7D[NP2r3_2r7D + 1] = {
00334   2.267936634217251403663034189684284173018E-5L,
00335   1.918112982168673386858072491437971732237E-3L,
00336   5.771704085468423159125856786653868219522E-2L,
00337   8.056124451167969333717642810661498890507E-1L,
00338   5.687897967531010276788680634413789328776E0L,
00339   2.072596760717695491085444438270778394421E1L,
00340   3.801722099819929988585197088613160496684E1L,
00341   3.254620235902912339534998592085115836829E1L,
00342   1.104847772130720331801884344645060675036E1L,
00343   /* 1.000000000000000000000000000000000000000E0 */
00344 };
00345 
00346 /* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2)
00347    Peak relative error 1.2e-35
00348    0.4375 <= 1/x <= 0.5  */
00349 #define NP2_2r3N 8
00350 static const long double P2_2r3N[NP2_2r3N + 1] = {
00351   -1.001042324337684297465071506097365389123E-4L,
00352   -6.289034524673365824853547252689991418981E-3L,
00353   -1.346527918018624234373664526930736205806E-1L,
00354   -1.268808313614288355444506172560463315102E0L,
00355   -5.654126123607146048354132115649177406163E0L,
00356   -1.186649511267312652171775803270911971693E1L,
00357   -1.094032424931998612551588246779200724257E1L,
00358   -3.728792136814520055025256353193674625267E0L,
00359   -3.000348318524471807839934764596331810608E-1L,
00360 };
00361 #define NP2_2r3D 8
00362 static const long double P2_2r3D[NP2_2r3D + 1] = {
00363   1.423705538269770974803901422532055612980E-3L,
00364   9.171476630091439978533535167485230575894E-2L,
00365   2.049776318166637248868444600215942828537E0L,
00366   2.068970329743769804547326701946144899583E1L,
00367   1.025103500560831035592731539565060347709E2L,
00368   2.528088049697570728252145557167066708284E2L,
00369   2.992160327587558573740271294804830114205E2L,
00370   1.540193761146551025832707739468679973036E2L,
00371   2.779516701986912132637672140709452502650E1L,
00372   /* 1.000000000000000000000000000000000000000E0 */
00373 };
00374 
00375 /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x),
00376    Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2))
00377    Peak relative error 2.2e-35
00378    0 <= 1/x <= .0625  */
00379 #define NQ16_IN 10
00380 static const long double Q16_IN[NQ16_IN + 1] = {
00381   2.343640834407975740545326632205999437469E-18L,
00382   2.667978112927811452221176781536278257448E-15L,
00383   1.178415018484555397390098879501969116536E-12L,
00384   2.622049767502719728905924701288614016597E-10L,
00385   3.196908059607618864801313380896308968673E-8L,
00386   2.179466154171673958770030655199434798494E-6L,
00387   8.139959091628545225221976413795645177291E-5L,
00388   1.563900725721039825236927137885747138654E-3L,
00389   1.355172364265825167113562519307194840307E-2L,
00390   3.928058355906967977269780046844768588532E-2L,
00391   1.107891967702173292405380993183694932208E-2L,
00392 };
00393 #define NQ16_ID 9
00394 static const long double Q16_ID[NQ16_ID + 1] = {
00395   3.199850952578356211091219295199301766718E-17L,
00396   3.652601488020654842194486058637953363918E-14L,
00397   1.620179741394865258354608590461839031281E-11L,
00398   3.629359209474609630056463248923684371426E-9L,
00399   4.473680923894354600193264347733477363305E-7L,
00400   3.106368086644715743265603656011050476736E-5L,
00401   1.198239259946770604954664925153424252622E-3L,
00402   2.446041004004283102372887804475767568272E-2L,
00403   2.403235525011860603014707768815113698768E-1L,
00404   9.491006790682158612266270665136910927149E-1L,
00405  /* 1.000000000000000000000000000000000000000E0 */
00406  };
00407 
00408 /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x),
00409    Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2))
00410    Peak relative error 5.1e-36
00411    0.0625 <= 1/x <= 0.125  */
00412 #define NQ8_16N 11
00413 static const long double Q8_16N[NQ8_16N + 1] = {
00414   1.001954266485599464105669390693597125904E-17L,
00415   7.545499865295034556206475956620160007849E-15L,
00416   2.267838684785673931024792538193202559922E-12L,
00417   3.561909705814420373609574999542459912419E-10L,
00418   3.216201422768092505214730633842924944671E-8L,
00419   1.731194793857907454569364622452058554314E-6L,
00420   5.576944613034537050396518509871004586039E-5L,
00421   1.051787760316848982655967052985391418146E-3L,
00422   1.102852974036687441600678598019883746959E-2L,
00423   5.834647019292460494254225988766702933571E-2L,
00424   1.290281921604364618912425380717127576529E-1L,
00425   7.598886310387075708640370806458926458301E-2L,
00426 };
00427 #define NQ8_16D 11
00428 static const long double Q8_16D[NQ8_16D + 1] = {
00429   1.368001558508338469503329967729951830843E-16L,
00430   1.034454121857542147020549303317348297289E-13L,
00431   3.128109209247090744354764050629381674436E-11L,
00432   4.957795214328501986562102573522064468671E-9L,
00433   4.537872468606711261992676606899273588899E-7L,
00434   2.493639207101727713192687060517509774182E-5L,
00435   8.294957278145328349785532236663051405805E-4L,
00436   1.646471258966713577374948205279380115839E-2L,
00437   1.878910092770966718491814497982191447073E-1L,
00438   1.152641605706170353727903052525652504075E0L,
00439   3.383550240669773485412333679367792932235E0L,
00440   3.823875252882035706910024716609908473970E0L,
00441  /* 1.000000000000000000000000000000000000000E0 */
00442 };
00443 
00444 /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x),
00445    Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2))
00446    Peak relative error 3.9e-35
00447    0.125 <= 1/x <= 0.1875  */
00448 #define NQ5_8N 10
00449 static const long double Q5_8N[NQ5_8N + 1] = {
00450   1.750399094021293722243426623211733898747E-13L,
00451   6.483426211748008735242909236490115050294E-11L,
00452   9.279430665656575457141747875716899958373E-9L,
00453   6.696634968526907231258534757736576340266E-7L,
00454   2.666560823798895649685231292142838188061E-5L,
00455   6.025087697259436271271562769707550594540E-4L,
00456   7.652807734168613251901945778921336353485E-3L,
00457   5.226269002589406461622551452343519078905E-2L,
00458   1.748390159751117658969324896330142895079E-1L,
00459   2.378188719097006494782174902213083589660E-1L,
00460   8.383984859679804095463699702165659216831E-2L,
00461 };
00462 #define NQ5_8D 10
00463 static const long double Q5_8D[NQ5_8D + 1] = {
00464   2.389878229704327939008104855942987615715E-12L,
00465   8.926142817142546018703814194987786425099E-10L,
00466   1.294065862406745901206588525833274399038E-7L,
00467   9.524139899457666250828752185212769682191E-6L,
00468   3.908332488377770886091936221573123353489E-4L,
00469   9.250427033957236609624199884089916836748E-3L,
00470   1.263420066165922645975830877751588421451E-1L,
00471   9.692527053860420229711317379861733180654E-1L,
00472   3.937813834630430172221329298841520707954E0L,
00473   7.603126427436356534498908111445191312181E0L,
00474   5.670677653334105479259958485084550934305E0L,
00475  /* 1.000000000000000000000000000000000000000E0 */
00476 };
00477 
00478 /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x),
00479    Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2))
00480    Peak relative error 3.2e-35
00481    0.1875 <= 1/x <= 0.25  */
00482 #define NQ4_5N 10
00483 static const long double Q4_5N[NQ4_5N + 1] = {
00484   2.233870042925895644234072357400122854086E-11L,
00485   5.146223225761993222808463878999151699792E-9L,
00486   4.459114531468296461688753521109797474523E-7L,
00487   1.891397692931537975547242165291668056276E-5L,
00488   4.279519145911541776938964806470674565504E-4L,
00489   5.275239415656560634702073291768904783989E-3L,
00490   3.468698403240744801278238473898432608887E-2L,
00491   1.138773146337708415188856882915457888274E-1L,
00492   1.622717518946443013587108598334636458955E-1L,
00493   7.249040006390586123760992346453034628227E-2L,
00494   1.941595365256460232175236758506411486667E-3L,
00495 };
00496 #define NQ4_5D 9
00497 static const long double Q4_5D[NQ4_5D + 1] = {
00498   3.049977232266999249626430127217988047453E-10L,
00499   7.120883230531035857746096928889676144099E-8L,
00500   6.301786064753734446784637919554359588859E-6L,
00501   2.762010530095069598480766869426308077192E-4L,
00502   6.572163250572867859316828886203406361251E-3L,
00503   8.752566114841221958200215255461843397776E-2L,
00504   6.487654992874805093499285311075289932664E-1L,
00505   2.576550017826654579451615283022812801435E0L,
00506   5.056392229924022835364779562707348096036E0L,
00507   4.179770081068251464907531367859072157773E0L,
00508  /* 1.000000000000000000000000000000000000000E0 */
00509 };
00510 
00511 /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x),
00512    Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2))
00513    Peak relative error 1.4e-36
00514    0.25 <= 1/x <= 0.3125  */
00515 #define NQ3r2_4N 10
00516 static const long double Q3r2_4N[NQ3r2_4N + 1] = {
00517   6.126167301024815034423262653066023684411E-10L,
00518   1.043969327113173261820028225053598975128E-7L,
00519   6.592927270288697027757438170153763220190E-6L,
00520   2.009103660938497963095652951912071336730E-4L,
00521   3.220543385492643525985862356352195896964E-3L,
00522   2.774405975730545157543417650436941650990E-2L,
00523   1.258114008023826384487378016636555041129E-1L,
00524   2.811724258266902502344701449984698323860E-1L,
00525   2.691837665193548059322831687432415014067E-1L,
00526   7.949087384900985370683770525312735605034E-2L,
00527   1.229509543620976530030153018986910810747E-3L,
00528 };
00529 #define NQ3r2_4D 9
00530 static const long double Q3r2_4D[NQ3r2_4D + 1] = {
00531   8.364260446128475461539941389210166156568E-9L,
00532   1.451301850638956578622154585560759862764E-6L,
00533   9.431830010924603664244578867057141839463E-5L,
00534   3.004105101667433434196388593004526182741E-3L,
00535   5.148157397848271739710011717102773780221E-2L,
00536   4.901089301726939576055285374953887874895E-1L,
00537   2.581760991981709901216967665934142240346E0L,
00538   7.257105880775059281391729708630912791847E0L,
00539   1.006014717326362868007913423810737369312E1L,
00540   5.879416600465399514404064187445293212470E0L,
00541  /* 1.000000000000000000000000000000000000000E0*/
00542 };
00543 
00544 /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x),
00545    Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2))
00546    Peak relative error 3.8e-36
00547    0.3125 <= 1/x <= 0.375  */
00548 #define NQ2r7_3r2N 9
00549 static const long double Q2r7_3r2N[NQ2r7_3r2N + 1] = {
00550   7.584861620402450302063691901886141875454E-8L,
00551   9.300939338814216296064659459966041794591E-6L,
00552   4.112108906197521696032158235392604947895E-4L,
00553   8.515168851578898791897038357239630654431E-3L,
00554   8.971286321017307400142720556749573229058E-2L,
00555   4.885856732902956303343015636331874194498E-1L,
00556   1.334506268733103291656253500506406045846E0L,
00557   1.681207956863028164179042145803851824654E0L,
00558   8.165042692571721959157677701625853772271E-1L,
00559   9.805848115375053300608712721986235900715E-2L,
00560 };
00561 #define NQ2r7_3r2D 9
00562 static const long double Q2r7_3r2D[NQ2r7_3r2D + 1] = {
00563   1.035586492113036586458163971239438078160E-6L,
00564   1.301999337731768381683593636500979713689E-4L,
00565   5.993695702564527062553071126719088859654E-3L,
00566   1.321184892887881883489141186815457808785E-1L,
00567   1.528766555485015021144963194165165083312E0L,
00568   9.561463309176490874525827051566494939295E0L,
00569   3.203719484883967351729513662089163356911E1L,
00570   5.497294687660930446641539152123568668447E1L,
00571   4.391158169390578768508675452986948391118E1L,
00572   1.347836630730048077907818943625789418378E1L,
00573  /* 1.000000000000000000000000000000000000000E0 */
00574 };
00575 
00576 /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x),
00577    Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2))
00578    Peak relative error 2.2e-35
00579    0.375 <= 1/x <= 0.4375  */
00580 #define NQ2r3_2r7N 9
00581 static const long double Q2r3_2r7N[NQ2r3_2r7N + 1] = {
00582   4.455027774980750211349941766420190722088E-7L,
00583   4.031998274578520170631601850866780366466E-5L,
00584   1.273987274325947007856695677491340636339E-3L,
00585   1.818754543377448509897226554179659122873E-2L,
00586   1.266748858326568264126353051352269875352E-1L,
00587   4.327578594728723821137731555139472880414E-1L,
00588   6.892532471436503074928194969154192615359E-1L,
00589   4.490775818438716873422163588640262036506E-1L,
00590   8.649615949297322440032000346117031581572E-2L,
00591   7.261345286655345047417257611469066147561E-4L,
00592 };
00593 #define NQ2r3_2r7D 8
00594 static const long double Q2r3_2r7D[NQ2r3_2r7D + 1] = {
00595   6.082600739680555266312417978064954793142E-6L,
00596   5.693622538165494742945717226571441747567E-4L,
00597   1.901625907009092204458328768129666975975E-2L,
00598   2.958689532697857335456896889409923371570E-1L,
00599   2.343124711045660081603809437993368799568E0L,
00600   9.665894032187458293568704885528192804376E0L,
00601   2.035273104990617136065743426322454881353E1L,
00602   2.044102010478792896815088858740075165531E1L,
00603   8.445937177863155827844146643468706599304E0L,
00604  /* 1.000000000000000000000000000000000000000E0 */
00605 };
00606 
00607 /* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x),
00608    Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2))
00609    Peak relative error 3.1e-36
00610    0.4375 <= 1/x <= 0.5  */
00611 #define NQ2_2r3N 9
00612 static const long double Q2_2r3N[NQ2_2r3N + 1] = {
00613   2.817566786579768804844367382809101929314E-6L,
00614   2.122772176396691634147024348373539744935E-4L,
00615   5.501378031780457828919593905395747517585E-3L,
00616   6.355374424341762686099147452020466524659E-2L,
00617   3.539652320122661637429658698954748337223E-1L,
00618   9.571721066119617436343740541777014319695E-1L,
00619   1.196258777828426399432550698612171955305E0L,
00620   6.069388659458926158392384709893753793967E-1L,
00621   9.026746127269713176512359976978248763621E-2L,
00622   5.317668723070450235320878117210807236375E-4L,
00623 };
00624 #define NQ2_2r3D 8
00625 static const long double Q2_2r3D[NQ2_2r3D + 1] = {
00626   3.846924354014260866793741072933159380158E-5L,
00627   3.017562820057704325510067178327449946763E-3L,
00628   8.356305620686867949798885808540444210935E-2L,
00629   1.068314930499906838814019619594424586273E0L,
00630   6.900279623894821067017966573640732685233E0L,
00631   2.307667390886377924509090271780839563141E1L,
00632   3.921043465412723970791036825401273528513E1L,
00633   3.167569478939719383241775717095729233436E1L,
00634   1.051023841699200920276198346301543665909E1L,
00635  /* 1.000000000000000000000000000000000000000E0*/
00636 };
00637 
00638 
00639 /* Evaluate P[n] x^n  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
00640 
00641 static long double
00642 neval (long double x, const long double *p, int n)
00643 {
00644   long double y;
00645 
00646   p += n;
00647   y = *p--;
00648   do
00649     {
00650       y = y * x + *p--;
00651     }
00652   while (--n > 0);
00653   return y;
00654 }
00655 
00656 
00657 /* Evaluate x^n+1  +  P[n] x^(n)  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
00658 
00659 static long double
00660 deval (long double x, const long double *p, int n)
00661 {
00662   long double y;
00663 
00664   p += n;
00665   y = x + *p--;
00666   do
00667     {
00668       y = y * x + *p--;
00669     }
00670   while (--n > 0);
00671   return y;
00672 }
00673 
00674 
00675 /* Bessel function of the first kind, order zero.  */
00676 
00677 long double
00678 __ieee754_j0l (long double x)
00679 {
00680   long double xx, xinv, z, p, q, c, s, cc, ss;
00681 
00682   if (! __finitel (x))
00683     {
00684       if (x != x)
00685        return x;
00686       else
00687        return 0.0L;
00688     }
00689   if (x == 0.0L)
00690     return 1.0L;
00691 
00692   xx = fabsl (x);
00693   if (xx <= 2.0L)
00694     {
00695       /* 0 <= x <= 2 */
00696       z = xx * xx;
00697       p = z * z * neval (z, J0_2N, NJ0_2N) / deval (z, J0_2D, NJ0_2D);
00698       p -= 0.25L * z;
00699       p += 1.0L;
00700       return p;
00701     }
00702 
00703   xinv = 1.0L / xx;
00704   z = xinv * xinv;
00705   if (xinv <= 0.25)
00706     {
00707       if (xinv <= 0.125)
00708        {
00709          if (xinv <= 0.0625)
00710            {
00711              p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID);
00712              q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID);
00713            }
00714          else
00715            {
00716              p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D);
00717              q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D);
00718            }
00719        }
00720       else if (xinv <= 0.1875)
00721        {
00722          p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D);
00723          q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D);
00724        }
00725       else
00726        {
00727          p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D);
00728          q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D);
00729        }
00730     }                       /* .25 */
00731   else /* if (xinv <= 0.5) */
00732     {
00733       if (xinv <= 0.375)
00734        {
00735          if (xinv <= 0.3125)
00736            {
00737              p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D);
00738              q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D);
00739            }
00740          else
00741            {
00742              p = neval (z, P2r7_3r2N, NP2r7_3r2N)
00743                 / deval (z, P2r7_3r2D, NP2r7_3r2D);
00744              q = neval (z, Q2r7_3r2N, NQ2r7_3r2N)
00745                 / deval (z, Q2r7_3r2D, NQ2r7_3r2D);
00746            }
00747        }
00748       else if (xinv <= 0.4375)
00749        {
00750          p = neval (z, P2r3_2r7N, NP2r3_2r7N)
00751              / deval (z, P2r3_2r7D, NP2r3_2r7D);
00752          q = neval (z, Q2r3_2r7N, NQ2r3_2r7N)
00753              / deval (z, Q2r3_2r7D, NQ2r3_2r7D);
00754        }
00755       else
00756        {
00757          p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
00758          q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
00759        }
00760     }
00761   p = 1.0L + z * p;
00762   q = z * xinv * q;
00763   q = q - 0.125L * xinv;
00764   /* X = x - pi/4
00765      cos(X) = cos(x) cos(pi/4) + sin(x) sin(pi/4)
00766      = 1/sqrt(2) * (cos(x) + sin(x))
00767      sin(X) = sin(x) cos(pi/4) - cos(x) sin(pi/4)
00768      = 1/sqrt(2) * (sin(x) - cos(x))
00769      sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
00770      cf. Fdlibm.  */
00771   __sincosl (xx, &s, &c);
00772   ss = s - c;
00773   cc = s + c;
00774   z = -__cosl (xx + xx);
00775   if ((s * c) < 0)
00776     cc = z / ss;
00777   else
00778     ss = z / cc;
00779   z = ONEOSQPI * (p * cc - q * ss) / __ieee754_sqrtl (xx);
00780   return z;
00781 }
00782 
00783 
00784 /* Y0(x) = 2/pi * log(x) * J0(x) + R(x^2)
00785    Peak absolute error 1.7e-36 (relative where Y0 > 1)
00786    0 <= x <= 2   */
00787 #define NY0_2N 7
00788 static long double Y0_2N[NY0_2N + 1] = {
00789  -1.062023609591350692692296993537002558155E19L,
00790   2.542000883190248639104127452714966858866E19L,
00791  -1.984190771278515324281415820316054696545E18L,
00792   4.982586044371592942465373274440222033891E16L,
00793  -5.529326354780295177243773419090123407550E14L,
00794   3.013431465522152289279088265336861140391E12L,
00795  -7.959436160727126750732203098982718347785E9L,
00796   8.230845651379566339707130644134372793322E6L,
00797 };
00798 #define NY0_2D 7
00799 static long double Y0_2D[NY0_2D + 1] = {
00800   1.438972634353286978700329883122253752192E20L,
00801   1.856409101981569254247700169486907405500E18L,
00802   1.219693352678218589553725579802986255614E16L,
00803   5.389428943282838648918475915779958097958E13L,
00804   1.774125762108874864433872173544743051653E11L,
00805   4.522104832545149534808218252434693007036E8L,
00806   8.872187401232943927082914504125234454930E5L,
00807   1.251945613186787532055610876304669413955E3L,
00808  /* 1.000000000000000000000000000000000000000E0 */
00809 };
00810 
00811 
00812 /* Bessel function of the second kind, order zero.  */
00813 
00814 long double
00815  __ieee754_y0l(long double x)
00816 {
00817   long double xx, xinv, z, p, q, c, s, cc, ss;
00818 
00819   if (! __finitel (x))
00820     {
00821       if (x != x)
00822        return x;
00823       else
00824        return 0.0L;
00825     }
00826   if (x <= 0.0L)
00827     {
00828       if (x < 0.0L)
00829        return (zero / (zero * x));
00830       return -HUGE_VALL + x;
00831     }
00832   xx = fabsl (x);
00833   if (xx <= 2.0L)
00834     {
00835       /* 0 <= x <= 2 */
00836       z = xx * xx;
00837       p = neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D);
00838       p = TWOOPI * __ieee754_logl (x) * __ieee754_j0l (x) + p;
00839       return p;
00840     }
00841 
00842   xinv = 1.0L / xx;
00843   z = xinv * xinv;
00844   if (xinv <= 0.25)
00845     {
00846       if (xinv <= 0.125)
00847        {
00848          if (xinv <= 0.0625)
00849            {
00850              p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID);
00851              q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID);
00852            }
00853          else
00854            {
00855              p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D);
00856              q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D);
00857            }
00858        }
00859       else if (xinv <= 0.1875)
00860        {
00861          p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D);
00862          q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D);
00863        }
00864       else
00865        {
00866          p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D);
00867          q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D);
00868        }
00869     }                       /* .25 */
00870   else /* if (xinv <= 0.5) */
00871     {
00872       if (xinv <= 0.375)
00873        {
00874          if (xinv <= 0.3125)
00875            {
00876              p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D);
00877              q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D);
00878            }
00879          else
00880            {
00881              p = neval (z, P2r7_3r2N, NP2r7_3r2N)
00882                 / deval (z, P2r7_3r2D, NP2r7_3r2D);
00883              q = neval (z, Q2r7_3r2N, NQ2r7_3r2N)
00884                 / deval (z, Q2r7_3r2D, NQ2r7_3r2D);
00885            }
00886        }
00887       else if (xinv <= 0.4375)
00888        {
00889          p = neval (z, P2r3_2r7N, NP2r3_2r7N)
00890              / deval (z, P2r3_2r7D, NP2r3_2r7D);
00891          q = neval (z, Q2r3_2r7N, NQ2r3_2r7N)
00892              / deval (z, Q2r3_2r7D, NQ2r3_2r7D);
00893        }
00894       else
00895        {
00896          p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
00897          q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
00898        }
00899     }
00900   p = 1.0L + z * p;
00901   q = z * xinv * q;
00902   q = q - 0.125L * xinv;
00903   /* X = x - pi/4
00904      cos(X) = cos(x) cos(pi/4) + sin(x) sin(pi/4)
00905      = 1/sqrt(2) * (cos(x) + sin(x))
00906      sin(X) = sin(x) cos(pi/4) - cos(x) sin(pi/4)
00907      = 1/sqrt(2) * (sin(x) - cos(x))
00908      sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
00909      cf. Fdlibm.  */
00910   __sincosl (x, &s, &c);
00911   ss = s - c;
00912   cc = s + c;
00913   z = -__cosl (x + x);
00914   if ((s * c) < 0)
00915     cc = z / ss;
00916   else
00917     ss = z / cc;
00918   z = ONEOSQPI * (p * ss + q * cc) / __ieee754_sqrtl (x);
00919   return z;
00920 }