e_j0l.c

Glibc 2.3.2源代码(解压后有100多M)

/*							j0l.c
 *
 *	Bessel function of order zero
 *
 *
 *
 * SYNOPSIS:
 *
 * long double x, y, j0l();
 *
 * y = j0l( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns Bessel function of first kind, order zero of the argument.
 *
 * The domain is divided into two major intervals [0, 2] and
 * (2, infinity). In the first interval the rational approximation
 * is J0(x) = 1 - x^2 / 4 + x^4 R(x^2)
 * The second interval is further partitioned into eight equal segments
 * of 1/x.
 *
 * J0(x) = sqrt(2/(pi x)) (P0(x) cos(X) - Q0(x) sin(X)),
 * X = x - pi/4,
 *
 * and the auxiliary functions are given by
 *
 * J0(x)cos(X) + Y0(x)sin(X) = sqrt( 2/(pi x)) P0(x),
 * P0(x) = 1 + 1/x^2 R(1/x^2)
 *
 * Y0(x)cos(X) - J0(x)sin(X) = sqrt( 2/(pi x)) Q0(x),
 * Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2))
 *
 *
 *
 * ACCURACY:
 *
 *                      Absolute error:
 * arithmetic   domain      # trials      peak         rms
 *    IEEE      0, 30       100000      1.7e-34      2.4e-35
 *
 *
 */

/*							y0l.c
 *
 *	Bessel function of the second kind, order zero
 *
 *
 *
 * SYNOPSIS:
 *
 * double x, y, y0l();
 *
 * y = y0l( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns Bessel function of the second kind, of order
 * zero, of the argument.
 *
 * The approximation is the same as for J0(x), and
 * Y0(x) = sqrt(2/(pi x)) (P0(x) sin(X) + Q0(x) cos(X)).
 *
 * ACCURACY:
 *
 *  Absolute error, when y0(x) < 1; else relative error:
 *
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      0, 30       100000      3.0e-34     2.7e-35
 *
 */

/* Copyright 2001 by Stephen L. Moshier ([e-mail]).

    This library is free software; you can redistribute it and/or
    modify it under the terms of the GNU Lesser General Public
    License as published by the Free Software Foundation; either
    version 2.1 of the License, or (at your option) any later version.

    This library is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public
    License along with this library; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */

#include "math.h"
#include "math_private.h"

/* 1 / sqrt(pi) */
static const long double ONEOSQPI = 5.6418958354775628694807945156077258584405E-1L;
/* 2 / pi */
static const long double TWOOPI = 6.3661977236758134307553505349005744813784E-1L;
static const long double zero = 0.0L;

/* J0(x) = 1 - x^2/4 + x^2 x^2 R(x^2)
   Peak relative error 3.4e-37
   0 <= x <= 2  */
#define NJ0_2N 6
static const long double J0_2N[NJ0_2N + 1] = {
  3.133239376997663645548490085151484674892E16L,
 -5.479944965767990821079467311839107722107E14L,
  6.290828903904724265980249871997551894090E12L,
 -3.633750176832769659849028554429106299915E10L,
  1.207743757532429576399485415069244807022E8L,
 -2.107485999925074577174305650549367415465E5L,
  1.562826808020631846245296572935547005859E2L,
};
#define NJ0_2D 6
static const long double J0_2D[NJ0_2D + 1] = {
  2.005273201278504733151033654496928968261E18L,
  2.063038558793221244373123294054149790864E16L,
  1.053350447931127971406896594022010524994E14L,
  3.496556557558702583143527876385508882310E11L,
  8.249114511878616075860654484367133976306E8L,
  1.402965782449571800199759247964242790589E6L,
  1.619910762853439600957801751815074787351E3L,
 /* 1.000000000000000000000000000000000000000E0 */
};

/* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2),
   0 <= 1/x <= .0625
   Peak relative error 3.3e-36  */
#define NP16_IN 9
static const long double P16_IN[NP16_IN + 1] = {
  -1.901689868258117463979611259731176301065E-16L,
  -1.798743043824071514483008340803573980931E-13L,
  -6.481746687115262291873324132944647438959E-11L,
  -1.150651553745409037257197798528294248012E-8L,
  -1.088408467297401082271185599507222695995E-6L,
  -5.551996725183495852661022587879817546508E-5L,
  -1.477286941214245433866838787454880214736E-3L,
  -1.882877976157714592017345347609200402472E-2L,
  -9.620983176855405325086530374317855880515E-2L,
  -1.271468546258855781530458854476627766233E-1L,
};
#define NP16_ID 9
static const long double P16_ID[NP16_ID + 1] = {
  2.704625590411544837659891569420764475007E-15L,
  2.562526347676857624104306349421985403573E-12L,
  9.259137589952741054108665570122085036246E-10L,
  1.651044705794378365237454962653430805272E-7L,
  1.573561544138733044977714063100859136660E-5L,
  8.134482112334882274688298469629884804056E-4L,
  2.219259239404080863919375103673593571689E-2L,
  2.976990606226596289580242451096393862792E-1L,
  1.713895630454693931742734911930937246254E0L,
  3.231552290717904041465898249160757368855E0L,
  /* 1.000000000000000000000000000000000000000E0 */
};

/* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2)
    0.0625 <= 1/x <= 0.125
    Peak relative error 2.4e-35  */
#define NP8_16N 10
static const long double P8_16N[NP8_16N + 1] = {
  -2.335166846111159458466553806683579003632E-15L,
  -1.382763674252402720401020004169367089975E-12L,
  -3.192160804534716696058987967592784857907E-10L,
  -3.744199606283752333686144670572632116899E-8L,
  -2.439161236879511162078619292571922772224E-6L,
  -9.068436986859420951664151060267045346549E-5L,
  -1.905407090637058116299757292660002697359E-3L,
  -2.164456143936718388053842376884252978872E-2L,
  -1.212178415116411222341491717748696499966E-1L,
  -2.782433626588541494473277445959593334494E-1L,
  -1.670703190068873186016102289227646035035E-1L,
};
#define NP8_16D 10
static const long double P8_16D[NP8_16D + 1] = {
  3.321126181135871232648331450082662856743E-14L,
  1.971894594837650840586859228510007703641E-11L,
  4.571144364787008285981633719513897281690E-9L,
  5.396419143536287457142904742849052402103E-7L,
  3.551548222385845912370226756036899901549E-5L,
  1.342353874566932014705609788054598013516E-3L,
  2.899133293006771317589357444614157734385E-2L,
  3.455374978185770197704507681491574261545E-1L,
  2.116616964297512311314454834712634820514E0L,
  5.850768316827915470087758636881584174432E0L,
  5.655273858938766830855753983631132928968E0L,
  /* 1.000000000000000000000000000000000000000E0 */
};

/* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2)
  0.125 <= 1/x <= 0.1875
  Peak relative error 2.7e-35  */
#define NP5_8N 10
static const long double P5_8N[NP5_8N + 1] = {
  -1.270478335089770355749591358934012019596E-12L,
  -4.007588712145412921057254992155810347245E-10L,
  -4.815187822989597568124520080486652009281E-8L,
  -2.867070063972764880024598300408284868021E-6L,
  -9.218742195161302204046454768106063638006E-5L,
  -1.635746821447052827526320629828043529997E-3L,
  -1.570376886640308408247709616497261011707E-2L,
  -7.656484795303305596941813361786219477807E-2L,
  -1.659371030767513274944805479908858628053E-1L,
  -1.185340550030955660015841796219919804915E-1L,
  -8.920026499909994671248893388013790366712E-3L,
};
#define NP5_8D 9
static const long double P5_8D[NP5_8D + 1] = {
  1.806902521016705225778045904631543990314E-11L,
  5.728502760243502431663549179135868966031E-9L,
  6.938168504826004255287618819550667978450E-7L,
  4.183769964807453250763325026573037785902E-5L,
  1.372660678476925468014882230851637878587E-3L,
  2.516452105242920335873286419212708961771E-2L,
  2.550502712902647803796267951846557316182E-1L,
  1.365861559418983216913629123778747617072E0L,
  3.523825618308783966723472468855042541407E0L,
  3.656365803506136165615111349150536282434E0L,
  /* 1.000000000000000000000000000000000000000E0 */
};

/* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2)
   Peak relative error 3.5e-35
   0.1875 <= 1/x <= 0.25  */
#define NP4_5N 9
static const long double P4_5N[NP4_5N + 1] = {
  -9.791405771694098960254468859195175708252E-10L,
  -1.917193059944531970421626610188102836352E-7L,
  -1.393597539508855262243816152893982002084E-5L,
  -4.881863490846771259880606911667479860077E-4L,
  -8.946571245022470127331892085881699269853E-3L,
  -8.707474232568097513415336886103899434251E-2L,
  -4.362042697474650737898551272505525973766E-1L,
  -1.032712171267523975431451359962375617386E0L,
  -9.630502683169895107062182070514713702346E-1L,
  -2.251804386252969656586810309252357233320E-1L,
};
#define NP4_5D 9
static const long double P4_5D[NP4_5D + 1] = {
  1.392555487577717669739688337895791213139E-8L,
  2.748886559120659027172816051276451376854E-6L,
  2.024717710644378047477189849678576659290E-4L,
  7.244868609350416002930624752604670292469E-3L,
  1.373631762292244371102989739300382152416E-1L,
  1.412298581400224267910294815260613240668E0L,
  7.742495637843445079276397723849017617210E0L,
  2.138429269198406512028307045259503811861E1L,
  2.651547684548423476506826951831712762610E1L,
  1.167499382465291931571685222882909166935E1L,
  /* 1.000000000000000000000000000000000000000E0 */
};

/* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2)
   Peak relative error 2.3e-36
   0.25 <= 1/x <= 0.3125  */
#define NP3r2_4N 9
static const long double P3r2_4N[NP3r2_4N + 1] = {
  -2.589155123706348361249809342508270121788E-8L,
  -3.746254369796115441118148490849195516593E-6L,
  -1.985595497390808544622893738135529701062E-4L,
  -5.008253705202932091290132760394976551426E-3L,
  -6.529469780539591572179155511840853077232E-2L,
  -4.468736064761814602927408833818990271514E-1L,
  -1.556391252586395038089729428444444823380E0L,
  -2.533135309840530224072920725976994981638E0L,
  -1.605509621731068453869408718565392869560E0L,
  -2.518966692256192789269859830255724429375E-1L,
};
#define NP3r2_4D 9
static const long double P3r2_4D[NP3r2_4D + 1] = {
  3.682353957237979993646169732962573930237E-7L,
  5.386741661883067824698973455566332102029E-5L,
  2.906881154171822780345134853794241037053E-3L,
  7.545832595801289519475806339863492074126E-2L,
  1.029405357245594877344360389469584526654E0L,
  7.565706120589873131187989560509757626725E0L,
  2.951172890699569545357692207898667665796E1L,
  5.785723537170311456298467310529815457536E1L,
  5.095621464598267889126015412522773474467E1L,
  1.602958484169953109437547474953308401442E1L,
  /* 1.000000000000000000000000000000000000000E0 */
};

/* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2)
   Peak relative error 1.0e-35
   0.3125 <= 1/x <= 0.375  */
#define NP2r7_3r2N 9
static const long double P2r7_3r2N[NP2r7_3r2N + 1] = {
  -1.917322340814391131073820537027234322550E-7L,
  -1.966595744473227183846019639723259011906E-5L,
  -7.177081163619679403212623526632690465290E-4L,
  -1.206467373860974695661544653741899755695E-2L,
  -1.008656452188539812154551482286328107316E-1L,
  -4.216016116408810856620947307438823892707E-1L,
  -8.378631013025721741744285026537009814161E-1L,
  -6.973895635309960850033762745957946272579E-1L,
  -1.797864718878320770670740413285763554812E-1L,
  -4.098025357743657347681137871388402849581E-3L,
};
#define NP2r7_3r2D 8
static const long double P2r7_3r2D[NP2r7_3r2D + 1] = {
  2.726858489303036441686496086962545034018E-6L,
  2.840430827557109238386808968234848081424E-4L,
  1.063826772041781947891481054529454088832E-2L,
  1.864775537138364773178044431045514405468E-1L,
  1.665660052857205170440952607701728254211E0L,