home *** CD-ROM | disk | FTP | other *** search
/ InfoMagic Source Code 1993 July / THE_SOURCE_CODE_CD_ROM.iso / bsd_srcs / lib / libm / common_source / j0.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-12-16  |  15.7 KB  |  440 lines
  1. /*-
  2.  * Copyright (c) 1992 The Regents of the University of California.
  3.  * All rights reserved.
  4.  *
  5.  * Redistribution and use in source and binary forms, with or without
  6.  * modification, are permitted provided that the following conditions
  7.  * are met:
  8.  * 1. Redistributions of source code must retain the above copyright
  9.  *    notice, this list of conditions and the following disclaimer.
  10.  * 2. Redistributions in binary form must reproduce the above copyright
  11.  *    notice, this list of conditions and the following disclaimer in the
  12.  *    documentation and/or other materials provided with the distribution.
  13.  * 3. All advertising materials mentioning features or use of this software
  14.  *    must display the following acknowledgement:
  15.  *    This product includes software developed by the University of
  16.  *    California, Berkeley and its contributors.
  17.  * 4. Neither the name of the University nor the names of its contributors
  18.  *    may be used to endorse or promote products derived from this software
  19.  *    without specific prior written permission.
  20.  *
  21.  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  22.  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  23.  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  24.  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  25.  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  26.  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  27.  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  28.  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  29.  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  30.  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  31.  * SUCH DAMAGE.
  32.  */
  33.  
  34. #ifndef lint
  35. static char sccsid[] = "@(#)j0.c    5.5 (Berkeley) 12/16/92";
  36. #endif /* not lint */
  37.  
  38. /*
  39.  * 16 December 1992
  40.  * Minor modifications by Peter McIlroy to adapt non-IEEE architecture.
  41.  */
  42.  
  43. /*
  44.  * ====================================================
  45.  * Copyright (C) 1992 by Sun Microsystems, Inc.
  46.  *
  47.  * Developed at SunPro, a Sun Microsystems, Inc. business.
  48.  * Permission to use, copy, modify, and distribute this
  49.  * software is freely granted, provided that this notice 
  50.  * is preserved.
  51.  * ====================================================
  52.  *
  53.  * ******************* WARNING ********************
  54.  * This is an alpha version of SunPro's FDLIBM (Freely
  55.  * Distributable Math Library) for IEEE double precision 
  56.  * arithmetic. FDLIBM is a basic math library written
  57.  * in C that runs on machines that conform to IEEE 
  58.  * Standard 754/854. This alpha version is distributed 
  59.  * for testing purpose. Those who use this software 
  60.  * should report any bugs to 
  61.  *
  62.  *        fdlibm-comments@sunpro.eng.sun.com
  63.  *
  64.  * -- K.C. Ng, Oct 12, 1992
  65.  * ************************************************
  66.  */
  67.  
  68. /* double j0(double x), y0(double x)
  69.  * Bessel function of the first and second kinds of order zero.
  70.  * Method -- j0(x):
  71.  *    1. For tiny x, we use j0(x) = 1 - x^2/4 + x^4/64 - ...
  72.  *    2. Reduce x to |x| since j0(x)=j0(-x),  and
  73.  *       for x in (0,2)
  74.  *        j0(x) = 1-z/4+ z^2*R0/S0,  where z = x*x;
  75.  *       (precision:  |j0-1+z/4-z^2R0/S0 |<2**-63.67 )
  76.  *       for x in (2,inf)
  77.  *         j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
  78.  *        where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
  79.  *       as follow:
  80.  *        cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
  81.  *            = 1/sqrt(2) * (cos(x) + sin(x))
  82.  *        sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
  83.  *            = 1/sqrt(2) * (sin(x) - cos(x))
  84.  *        (To avoid cancellation, use
  85.  *        sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
  86.  *         to compute the worse one.)
  87.  *       
  88.  *    3 Special cases
  89.  *        j0(nan)= nan
  90.  *        j0(0) = 1
  91.  *        j0(inf) = 0
  92.  *        
  93.  * Method -- y0(x):
  94.  *    1. For x<2.
  95.  *       Since 
  96.  *        y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x^2/4 - ...)
  97.  *       therefore y0(x)-2/pi*j0(x)*ln(x) is an even function.
  98.  *       We use the following function to approximate y0,
  99.  *        y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x^2
  100.  *       where 
  101.  *        U(z) = u0 + u1*z + ... + u6*z^6
  102.  *        V(z) = 1  + v1*z + ... + v4*z^4
  103.  *       with absolute approximation error bounded by 2**-72.
  104.  *       Note: For tiny x, U/V = u0 and j0(x)~1, hence
  105.  *        y0(tiny) = u0 + (2/pi)*ln(tiny), (choose tiny<2**-27)
  106.  *    2. For x>=2.
  107.  *         y0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)+q0(x)*sin(x0))
  108.  *        where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
  109.  *       by the method mentioned above.
  110.  *    3. Special cases: y0(0)=-inf, y0(x<0)=NaN, y0(inf)=0.
  111.  */
  112.  
  113. #include <math.h>
  114. #include <float.h>
  115. #if defined(vax) || defined(tahoe)
  116. #define _IEEE    0
  117. #else
  118. #define _IEEE    1
  119. #define infnan(x) (0.0)
  120. #endif
  121.  
  122. static double pzero __P((double)), qzero __P((double));
  123.  
  124. static double 
  125. huge     = 1e300,
  126. zero    = 0.0,
  127. one    = 1.0,
  128. invsqrtpi= 5.641895835477562869480794515607725858441e-0001,
  129. tpi    = 0.636619772367581343075535053490057448,
  130.          /* R0/S0 on [0, 2.00] */
  131. r02 =   1.562499999999999408594634421055018003102e-0002,
  132. r03 =  -1.899792942388547334476601771991800712355e-0004,
  133. r04 =   1.829540495327006565964161150603950916854e-0006,
  134. r05 =  -4.618326885321032060803075217804816988758e-0009,
  135. s01 =   1.561910294648900170180789369288114642057e-0002,
  136. s02 =   1.169267846633374484918570613449245536323e-0004,
  137. s03 =   5.135465502073181376284426245689510134134e-0007,
  138. s04 =   1.166140033337900097836930825478674320464e-0009;
  139.  
  140. double
  141. j0(x) 
  142.     double x;
  143. {
  144.     double z, s,c,ss,cc,r,u,v;
  145.  
  146.     if (!finite(x))
  147.         if (_IEEE) return one/(x*x);
  148.         else return (0);
  149.     x = fabs(x);
  150.     if (x >= 2.0) {    /* |x| >= 2.0 */
  151.         s = sin(x);
  152.         c = cos(x);
  153.         ss = s-c;
  154.         cc = s+c;
  155.         if (x < .5 * DBL_MAX) {  /* make sure x+x not overflow */
  156.             z = -cos(x+x);
  157.             if ((s*c)<zero) cc = z/ss;
  158.             else         ss = z/cc;
  159.         }
  160.     /*
  161.      * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
  162.      * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
  163.      */
  164.         if (_IEEE && x> 6.80564733841876927e+38) /* 2^129 */
  165.             z = (invsqrtpi*cc)/sqrt(x);
  166.         else {
  167.             u = pzero(x); v = qzero(x);
  168.             z = invsqrtpi*(u*cc-v*ss)/sqrt(x);
  169.         }
  170.         return z;
  171.     }
  172.     if (x < 1.220703125e-004) {           /* |x| < 2**-13 */
  173.         if (huge+x > one) {               /* raise inexact if x != 0 */
  174.             if (x < 7.450580596923828125e-009) /* |x|<2**-27 */
  175.             return one;
  176.             else return (one - 0.25*x*x);
  177.         }
  178.     }
  179.     z = x*x;
  180.     r =  z*(r02+z*(r03+z*(r04+z*r05)));
  181.     s =  one+z*(s01+z*(s02+z*(s03+z*s04)));
  182.     if (x < one) {            /* |x| < 1.00 */
  183.         return (one + z*(-0.25+(r/s)));
  184.     } else {
  185.         u = 0.5*x;
  186.         return ((one+u)*(one-u)+z*(r/s));
  187.     }
  188. }
  189.  
  190. static double
  191. u00 =  -7.380429510868722527422411862872999615628e-0002,
  192. u01 =   1.766664525091811069896442906220827182707e-0001,
  193. u02 =  -1.381856719455968955440002438182885835344e-0002,
  194. u03 =   3.474534320936836562092566861515617053954e-0004,
  195. u04 =  -3.814070537243641752631729276103284491172e-0006,
  196. u05 =   1.955901370350229170025509706510038090009e-0008,
  197. u06 =  -3.982051941321034108350630097330144576337e-0011,
  198. v01 =   1.273048348341237002944554656529224780561e-0002,
  199. v02 =   7.600686273503532807462101309675806839635e-0005,
  200. v03 =   2.591508518404578033173189144579208685163e-0007,
  201. v04 =   4.411103113326754838596529339004302243157e-0010;
  202.  
  203. double
  204. y0(x) 
  205.     double x;
  206. {
  207.     double z, s,c,ss,cc,u,v,j0();
  208.     /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0  */
  209.     if (!finite(x))
  210.         if (_IEEE)
  211.             return (one/(x+x*x));
  212.         else
  213.             return (0);
  214.         if (x == 0)
  215.         if (_IEEE)    return (-one/zero);
  216.         else        return(infnan(-ERANGE));
  217.         if (x<0)
  218.         if (_IEEE)    return (zero/zero);
  219.         else        return (infnan(EDOM));
  220.         if (x >= 2.00) {    /* |x| >= 2.0 */
  221.         /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
  222.          * where x0 = x-pi/4
  223.          *      Better formula:
  224.          *              cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
  225.          *                      =  1/sqrt(2) * (sin(x) + cos(x))
  226.          *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
  227.          *                      =  1/sqrt(2) * (sin(x) - cos(x))
  228.          * To avoid cancellation, use
  229.          *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
  230.          * to compute the worse one.
  231.          */
  232.                 s = sin(x);
  233.                 c = cos(x);
  234.                 ss = s-c;
  235.                 cc = s+c;
  236.     /*
  237.      * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
  238.      * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
  239.      */
  240.                 if (x < .5 * DBL_MAX) {  /* make sure x+x not overflow */
  241.                     z = -cos(x+x);
  242.                     if ((s*c)<zero) cc = z/ss;
  243.                     else            ss = z/cc;
  244.                 }
  245.                 if (_IEEE && x > 6.80564733841876927e+38) /* > 2^129 */
  246.             z = (invsqrtpi*ss)/sqrt(x);
  247.                 else {
  248.                     u = pzero(x); v = qzero(x);
  249.                     z = invsqrtpi*(u*ss+v*cc)/sqrt(x);
  250.                 }
  251.                 return z;
  252.     }
  253.     if (x <= 7.450580596923828125e-009) {        /* x < 2**-27 */
  254.         return (u00 + tpi*log(x));
  255.     }
  256.     z = x*x;
  257.     u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06)))));
  258.     v = one+z*(v01+z*(v02+z*(v03+z*v04)));
  259.     return (u/v + tpi*(j0(x)*log(x)));
  260. }
  261.  
  262. /* The asymptotic expansions of pzero is
  263.  *    1 - 9/128 s^2 + 11025/98304 s^4 - ...,    where s = 1/x.
  264.  * For x >= 2, We approximate pzero by
  265.  *     pzero(x) = 1 + (R/S)
  266.  * where  R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10
  267.  *       S = 1 + ps0*s^2 + ... + ps4*s^10
  268.  * and
  269.  *    | pzero(x)-1-R/S | <= 2  ** ( -60.26)
  270.  */
  271. static double pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
  272.    0.0,
  273.   -7.031249999999003994151563066182798210142e-0002,
  274.   -8.081670412753498508883963849859423939871e+0000,
  275.   -2.570631056797048755890526455854482662510e+0002,
  276.   -2.485216410094288379417154382189125598962e+0003,
  277.   -5.253043804907295692946647153614119665649e+0003,
  278. };
  279. static double ps8[5] = {
  280.    1.165343646196681758075176077627332052048e+0002,
  281.    3.833744753641218451213253490882686307027e+0003,
  282.    4.059785726484725470626341023967186966531e+0004,
  283.    1.167529725643759169416844015694440325519e+0005,
  284.    4.762772841467309430100106254805711722972e+0004,
  285. };
  286.  
  287. static double pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
  288.   -1.141254646918944974922813501362824060117e-0011,
  289.   -7.031249408735992804117367183001996028304e-0002,
  290.   -4.159610644705877925119684455252125760478e+0000,
  291.   -6.767476522651671942610538094335912346253e+0001,
  292.   -3.312312996491729755731871867397057689078e+0002,
  293.   -3.464333883656048910814187305901796723256e+0002,
  294. };
  295. static double ps5[5] = {
  296.    6.075393826923003305967637195319271932944e+0001,
  297.    1.051252305957045869801410979087427910437e+0003,
  298.    5.978970943338558182743915287887408780344e+0003,
  299.    9.625445143577745335793221135208591603029e+0003,
  300.    2.406058159229391070820491174867406875471e+0003,
  301. };
  302.  
  303. static double pr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
  304.   -2.547046017719519317420607587742992297519e-0009,
  305.   -7.031196163814817199050629727406231152464e-0002,
  306.   -2.409032215495295917537157371488126555072e+0000,
  307.   -2.196597747348830936268718293366935843223e+0001,
  308.   -5.807917047017375458527187341817239891940e+0001,
  309.   -3.144794705948885090518775074177485744176e+0001,
  310. };
  311. static double ps3[5] = {
  312.    3.585603380552097167919946472266854507059e+0001,
  313.    3.615139830503038919981567245265266294189e+0002,
  314.    1.193607837921115243628631691509851364715e+0003,
  315.    1.127996798569074250675414186814529958010e+0003,
  316.    1.735809308133357510239737333055228118910e+0002,
  317. };
  318.  
  319. static double pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
  320.   -8.875343330325263874525704514800809730145e-0008,
  321.   -7.030309954836247756556445443331044338352e-0002,
  322.   -1.450738467809529910662233622603401167409e+0000,
  323.   -7.635696138235277739186371273434739292491e+0000,
  324.   -1.119316688603567398846655082201614524650e+0001,
  325.   -3.233645793513353260006821113608134669030e+0000,
  326. };
  327. static double ps2[5] = {
  328.    2.222029975320888079364901247548798910952e+0001,
  329.    1.362067942182152109590340823043813120940e+0002,
  330.    2.704702786580835044524562897256790293238e+0002,
  331.    1.538753942083203315263554770476850028583e+0002,
  332.    1.465761769482561965099880599279699314477e+0001,
  333. };
  334.  
  335. static double pzero(x)
  336.     double x;
  337. {
  338.     double *p,*q,z,r,s;
  339.     if (x >= 8.00)               {p = pr8; q= ps8;}
  340.     else if (x >= 4.54545211791992188) {p = pr5; q= ps5;}
  341.     else if (x >= 2.85714149475097656) {p = pr3; q= ps3;}
  342.     else if (x >= 2.00)           {p = pr2; q= ps2;}
  343.     z = one/(x*x);
  344.     r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
  345.     s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
  346.     return one+ r/s;
  347. }
  348.         
  349.  
  350. /* For x >= 8, the asymptotic expansions of qzero is
  351.  *    -1/8 s + 75/1024 s^3 - ..., where s = 1/x.
  352.  * We approximate pzero by
  353.  *     qzero(x) = s*(-1.25 + (R/S))
  354.  * where  R = qr0 + qr1*s^2 + qr2*s^4 + ... + qr5*s^10
  355.  *       S = 1 + qs0*s^2 + ... + qs5*s^12
  356.  * and
  357.  *    | qzero(x)/s +1.25-R/S | <= 2  ** ( -61.22)
  358.  */
  359. static double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
  360.    0.0,
  361.    7.324218749999350414479738504551775297096e-0002,
  362.    1.176820646822526933903301695932765232456e+0001,
  363.    5.576733802564018422407734683549251364365e+0002,
  364.    8.859197207564685717547076568608235802317e+0003,
  365.    3.701462677768878501173055581933725704809e+0004,
  366. };
  367. static double qs8[6] = {
  368.    1.637760268956898345680262381842235272369e+0002,
  369.    8.098344946564498460163123708054674227492e+0003,
  370.    1.425382914191204905277585267143216379136e+0005,
  371.    8.033092571195144136565231198526081387047e+0005,
  372.    8.405015798190605130722042369969184811488e+0005,
  373.   -3.438992935378666373204500729736454421006e+0005,
  374. };
  375.  
  376. static double qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
  377.    1.840859635945155400568380711372759921179e-0011,
  378.    7.324217666126847411304688081129741939255e-0002,
  379.    5.835635089620569401157245917610984757296e+0000,
  380.    1.351115772864498375785526599119895942361e+0002,
  381.    1.027243765961641042977177679021711341529e+0003,
  382.    1.989977858646053872589042328678602481924e+0003,
  383. };
  384. static double qs5[6] = {
  385.    8.277661022365377058749454444343415524509e+0001,
  386.    2.077814164213929827140178285401017305309e+0003,
  387.    1.884728877857180787101956800212453218179e+0004,
  388.    5.675111228949473657576693406600265778689e+0004,
  389.    3.597675384251145011342454247417399490174e+0004,
  390.   -5.354342756019447546671440667961399442388e+0003,
  391. };
  392.  
  393. static double qr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
  394.    4.377410140897386263955149197672576223054e-0009,
  395.    7.324111800429115152536250525131924283018e-0002,
  396.    3.344231375161707158666412987337679317358e+0000,
  397.    4.262184407454126175974453269277100206290e+0001,
  398.    1.708080913405656078640701512007621675724e+0002,
  399.    1.667339486966511691019925923456050558293e+0002,
  400. };
  401. static double qs3[6] = {
  402.    4.875887297245871932865584382810260676713e+0001,
  403.    7.096892210566060535416958362640184894280e+0002,
  404.    3.704148226201113687434290319905207398682e+0003,
  405.    6.460425167525689088321109036469797462086e+0003,
  406.    2.516333689203689683999196167394889715078e+0003,
  407.   -1.492474518361563818275130131510339371048e+0002,
  408. };
  409.  
  410. static double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
  411.    1.504444448869832780257436041633206366087e-0007,
  412.    7.322342659630792930894554535717104926902e-0002,
  413.    1.998191740938159956838594407540292600331e+0000,
  414.    1.449560293478857407645853071687125850962e+0001,
  415.    3.166623175047815297062638132537957315395e+0001,
  416.    1.625270757109292688799540258329430963726e+0001,
  417. };
  418. static double qs2[6] = {
  419.    3.036558483552191922522729838478169383969e+0001,
  420.    2.693481186080498724211751445725708524507e+0002,
  421.    8.447837575953201460013136756723746023736e+0002,
  422.    8.829358451124885811233995083187666981299e+0002,
  423.    2.126663885117988324180482985363624996652e+0002,
  424.   -5.310954938826669402431816125780738924463e+0000,
  425. };
  426.  
  427. static double qzero(x)
  428.     double x;
  429. {
  430.     double *p,*q, s,r,z;
  431.     if (x >= 8.00)               {p = qr8; q= qs8;}
  432.     else if (x >= 4.54545211791992188) {p = qr5; q= qs5;}
  433.     else if (x >= 2.85714149475097656) {p = qr3; q= qs3;}
  434.     else if (x >= 2.00)           {p = qr2; q= qs2;}
  435.     z = one/(x*x);
  436.     r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
  437.     s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
  438.     return (-.125 + r/s)/x;
  439. }
  440.