e_j1l.c

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

/*							j1l.c
 *
 *	Bessel function of order one
 *
 *
 *
 * SYNOPSIS:
 *
 * long double x, y, j1l();
 *
 * y = j1l( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns Bessel function of first kind, order one of the argument.
 *
 * The domain is divided into two major intervals [0, 2] and
 * (2, infinity). In the first interval the rational approximation is
 * J1(x) = .5x + x x^2 R(x^2)
 *
 * The second interval is further partitioned into eight equal segments
 * of 1/x.
 * J1(x) = sqrt(2/(pi x)) (P1(x) cos(X) - Q1(x) sin(X)),
 * X = x - 3 pi / 4,
 *
 * and the auxiliary functions are given by
 *
 * J1(x)cos(X) + Y1(x)sin(X) = sqrt( 2/(pi x)) P1(x),
 * P1(x) = 1 + 1/x^2 R(1/x^2)
 *
 * Y1(x)cos(X) - J1(x)sin(X) = sqrt( 2/(pi x)) Q1(x),
 * Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)).
 *
 *
 *
 * ACCURACY:
 *
 *                      Absolute error:
 * arithmetic   domain      # trials      peak         rms
 *    IEEE      0, 30       100000      2.8e-34      2.7e-35
 *
 *
 */

/*							y1l.c
 *
 *	Bessel function of the second kind, order one
 *
 *
 *
 * SYNOPSIS:
 *
 * double x, y, y1l();
 *
 * y = y1l( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns Bessel function of the second kind, of order
 * one, of the argument.
 *
 * The domain is divided into two major intervals [0, 2] and
 * (2, infinity). In the first interval the rational approximation is
 * Y1(x) = 2/pi * (log(x) * J1(x) - 1/x) + x R(x^2) .
 * In the second interval the approximation is the same as for J1(x), and
 * Y1(x) = sqrt(2/(pi x)) (P1(x) sin(X) + Q1(x) cos(X)),
 * X = x - 3 pi / 4.
 *
 * ACCURACY:
 *
 *  Absolute error, when y0(x) < 1; else relative error:
 *
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      0, 30       100000      2.7e-34     2.9e-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;

/* J1(x) = .5x + x x^2 R(x^2)
   Peak relative error 1.9e-35
   0 <= x <= 2  */
#define NJ0_2N 6
static const long double J0_2N[NJ0_2N + 1] = {
 -5.943799577386942855938508697619735179660E16L,
  1.812087021305009192259946997014044074711E15L,
 -2.761698314264509665075127515729146460895E13L,
  2.091089497823600978949389109350658815972E11L,
 -8.546413231387036372945453565654130054307E8L,
  1.797229225249742247475464052741320612261E6L,
 -1.559552840946694171346552770008812083969E3L
};
#define NJ0_2D 6
static const long double J0_2D[NJ0_2D + 1] = {
  9.510079323819108569501613916191477479397E17L,
  1.063193817503280529676423936545854693915E16L,
  5.934143516050192600795972192791775226920E13L,
  2.168000911950620999091479265214368352883E11L,
  5.673775894803172808323058205986256928794E8L,
  1.080329960080981204840966206372671147224E6L,
  1.411951256636576283942477881535283304912E3L,
 /* 1.000000000000000000000000000000000000000E0L */
};

/* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
   0 <= 1/x <= .0625
   Peak relative error 3.6e-36  */
#define NP16_IN 9
static const long double P16_IN[NP16_IN + 1] = {
  5.143674369359646114999545149085139822905E-16L,
  4.836645664124562546056389268546233577376E-13L,
  1.730945562285804805325011561498453013673E-10L,
  3.047976856147077889834905908605310585810E-8L,
  2.855227609107969710407464739188141162386E-6L,
  1.439362407936705484122143713643023998457E-4L,
  3.774489768532936551500999699815873422073E-3L,
  4.723962172984642566142399678920790598426E-2L,
  2.359289678988743939925017240478818248735E-1L,
  3.032580002220628812728954785118117124520E-1L,
};
#define NP16_ID 9
static const long double P16_ID[NP16_ID + 1] = {
  4.389268795186898018132945193912677177553E-15L,
  4.132671824807454334388868363256830961655E-12L,
  1.482133328179508835835963635130894413136E-9L,
  2.618941412861122118906353737117067376236E-7L,
  2.467854246740858470815714426201888034270E-5L,
  1.257192927368839847825938545925340230490E-3L,
  3.362739031941574274949719324644120720341E-2L,
  4.384458231338934105875343439265370178858E-1L,
  2.412830809841095249170909628197264854651E0L,
  4.176078204111348059102962617368214856874E0L,
 /* 1.000000000000000000000000000000000000000E0 */
};

/* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
    0.0625 <= 1/x <= 0.125
    Peak relative error 1.9e-36  */
#define NP8_16N 11
static const long double P8_16N[NP8_16N + 1] = {
  2.984612480763362345647303274082071598135E-16L,
  1.923651877544126103941232173085475682334E-13L,
  4.881258879388869396043760693256024307743E-11L,
  6.368866572475045408480898921866869811889E-9L,
  4.684818344104910450523906967821090796737E-7L,
  2.005177298271593587095982211091300382796E-5L,
  4.979808067163957634120681477207147536182E-4L,
  6.946005761642579085284689047091173581127E-3L,
  5.074601112955765012750207555985299026204E-2L,
  1.698599455896180893191766195194231825379E-1L,
  1.957536905259237627737222775573623779638E-1L,
  2.991314703282528370270179989044994319374E-2L,
};
#define NP8_16D 10
static const long double P8_16D[NP8_16D + 1] = {
  2.546869316918069202079580939942463010937E-15L,
  1.644650111942455804019788382157745229955E-12L,
  4.185430770291694079925607420808011147173E-10L,
  5.485331966975218025368698195861074143153E-8L,
  4.062884421686912042335466327098932678905E-6L,
  1.758139661060905948870523641319556816772E-4L,
  4.445143889306356207566032244985607493096E-3L,
  6.391901016293512632765621532571159071158E-2L,
  4.933040207519900471177016015718145795434E-1L,
  1.839144086168947712971630337250761842976E0L,
  2.715120873995490920415616716916149586579E0L,
 /* 1.000000000000000000000000000000000000000E0 */
};

/* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
  0.125 <= 1/x <= 0.1875
  Peak relative error 1.3e-36  */
#define NP5_8N 10
static const long double P5_8N[NP5_8N + 1] = {
  2.837678373978003452653763806968237227234E-12L,
  9.726641165590364928442128579282742354806E-10L,
  1.284408003604131382028112171490633956539E-7L,
  8.524624695868291291250573339272194285008E-6L,
  3.111516908953172249853673787748841282846E-4L,
  6.423175156126364104172801983096596409176E-3L,
  7.430220589989104581004416356260692450652E-2L,
  4.608315409833682489016656279567605536619E-1L,
  1.396870223510964882676225042258855977512E0L,
  1.718500293904122365894630460672081526236E0L,
  5.465927698800862172307352821870223855365E-1L
};
#define NP5_8D 10
static const long double P5_8D[NP5_8D + 1] = {
  2.421485545794616609951168511612060482715E-11L,
  8.329862750896452929030058039752327232310E-9L,
  1.106137992233383429630592081375289010720E-6L,
  7.405786153760681090127497796448503306939E-5L,
  2.740364785433195322492093333127633465227E-3L,
  5.781246470403095224872243564165254652198E-2L,
  6.927711353039742469918754111511109983546E-1L,
  4.558679283460430281188304515922826156690E0L,
  1.534468499844879487013168065728837900009E1L,
  2.313927430889218597919624843161569422745E1L,
  1.194506341319498844336768473218382828637E1L,
 /* 1.000000000000000000000000000000000000000E0 */
};

/* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
   Peak relative error 1.4e-36
   0.1875 <= 1/x <= 0.25  */
#define NP4_5N 10
static const long double P4_5N[NP4_5N + 1] = {
  1.846029078268368685834261260420933914621E-10L,
  3.916295939611376119377869680335444207768E-8L,
  3.122158792018920627984597530935323997312E-6L,
  1.218073444893078303994045653603392272450E-4L,
  2.536420827983485448140477159977981844883E-3L,
  2.883011322006690823959367922241169171315E-2L,
  1.755255190734902907438042414495469810830E-1L,
  5.379317079922628599870898285488723736599E-1L,
  7.284904050194300773890303361501726561938E-1L,
  3.270110346613085348094396323925000362813E-1L,
  1.804473805689725610052078464951722064757E-2L,
};
#define NP4_5D 9
static const long double P4_5D[NP4_5D + 1] = {
  1.575278146806816970152174364308980863569E-9L,
  3.361289173657099516191331123405675054321E-7L,
  2.704692281550877810424745289838790693708E-5L,
  1.070854930483999749316546199273521063543E-3L,
  2.282373093495295842598097265627962125411E-2L,
  2.692025460665354148328762368240343249830E-1L,
  1.739892942593664447220951225734811133759E0L,
  5.890727576752230385342377570386657229324E0L,
  9.517442287057841500750256954117735128153E0L,
  6.100616353935338240775363403030137736013E0L,
 /* 1.000000000000000000000000000000000000000E0 */
};

/* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
   Peak relative error 3.0e-36
   0.25 <= 1/x <= 0.3125  */
#define NP3r2_4N 9
static const long double P3r2_4N[NP3r2_4N + 1] = {
  8.240803130988044478595580300846665863782E-8L,
  1.179418958381961224222969866406483744580E-5L,
  6.179787320956386624336959112503824397755E-4L,
  1.540270833608687596420595830747166658383E-2L,
  1.983904219491512618376375619598837355076E-1L,
  1.341465722692038870390470651608301155565E0L,
  4.617865326696612898792238245990854646057E0L,
  7.435574801812346424460233180412308000587E0L,
  4.671327027414635292514599201278557680420E0L,
  7.299530852495776936690976966995187714739E-1L,
};
#define NP3r2_4D 9
static const long double P3r2_4D[NP3r2_4D + 1] = {
  7.032152009675729604487575753279187576521E-7L,
  1.015090352324577615777511269928856742848E-4L,
  5.394262184808448484302067955186308730620E-3L,
  1.375291438480256110455809354836988584325E-1L,
  1.836247144461106304788160919310404376670E0L,
  1.314378564254376655001094503090935880349E1L,
  4.957184590465712006934452500894672343488E1L,
  9.287394244300647738855415178790263465398E1L,
  7.652563275535900609085229286020552768399E1L,
  2.147042473003074533150718117770093209096E1L,
 /* 1.000000000000000000000000000000000000000E0 */
};

/* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(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] = {
  4.599033469240421554219816935160627085991E-7L,
  4.665724440345003914596647144630893997284E-5L,
  1.684348845667764271596142716944374892756E-3L,
  2.802446446884455707845985913454440176223E-2L,
  2.321937586453963310008279956042545173930E-1L,
  9.640277413988055668692438709376437553804E-1L,
  1.911021064710270904508663334033003246028E0L,
  1.600811610164341450262992138893970224971E0L,
  4.266299218652587901171386591543457861138E-1L,
  1.316470424456061252962568223251247207325E-2L,
};
#define NP2r7_3r2D 8
static const long double P2r7_3r2D[NP2r7_3r2D + 1] = {