home *** CD-ROM | disk | FTP | other *** search
/ InfoMagic Source Code 1993 July / THE_SOURCE_CODE_CD_ROM.iso / bsd_srcs / lib / libm / common_source / j1.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-12-16  |  15.7 KB  |  447 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[] = "@(#)j1.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 j1(double x), y1(double x)
  69.  * Bessel function of the first and second kinds of order zero.
  70.  * Method -- j1(x):
  71.  *    1. For tiny x, we use j1(x) = x/2 - x^3/16 + x^5/384 - ...
  72.  *    2. Reduce x to |x| since j1(x)=-j1(-x),  and
  73.  *       for x in (0,2)
  74.  *        j1(x) = x/2 + x*z*R0/S0,  where z = x*x;
  75.  *       (precision:  |j1/x - 1/2 - R0/S0 |<2**-61.51 )
  76.  *       for x in (2,inf)
  77.  *         j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1))
  78.  *         y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
  79.  *        where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
  80.  *       as follows:
  81.  *        cos(x1) =  cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
  82.  *            =  1/sqrt(2) * (sin(x) - cos(x))
  83.  *        sin(x1) =  sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
  84.  *            = -1/sqrt(2) * (sin(x) + cos(x))
  85.  *        (To avoid cancellation, use
  86.  *        sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
  87.  *         to compute the worse one.)
  88.  *       
  89.  *    3 Special cases
  90.  *        j1(nan)= nan
  91.  *        j1(0) = 0
  92.  *        j1(inf) = 0
  93.  *        
  94.  * Method -- y1(x):
  95.  *    1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN 
  96.  *    2. For x<2.
  97.  *       Since 
  98.  *        y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...)
  99.  *       therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function.
  100.  *       We use the following function to approximate y1,
  101.  *        y1(x) = x*U(z)/V(z) + (2/pi)*(j1(x)*ln(x)-1/x), z= x^2
  102.  *       where for x in [0,2] (abs err less than 2**-65.89)
  103.  *        U(z) = u0 + u1*z + ... + u4*z^4
  104.  *        V(z) = 1  + v1*z + ... + v5*z^5
  105.  *       Note: For tiny x, 1/x dominate y1 and hence
  106.  *        y1(tiny) = -2/pi/tiny, (choose tiny<2**-54)
  107.  *    3. For x>=2.
  108.  *         y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
  109.  *        where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
  110.  *       by method mentioned above.
  111.  */
  112.  
  113. #include <math.h>
  114. #include <float.h>
  115.  
  116. #if defined(vax) || defined(tahoe)
  117. #define _IEEE    0
  118. #else
  119. #define _IEEE    1
  120. #define infnan(x) (0.0)
  121. #endif
  122.  
  123. static double pone(), qone();
  124.  
  125. static double 
  126. huge    = 1e300,
  127. zero    = 0.0,
  128. one    = 1.0,
  129. invsqrtpi= 5.641895835477562869480794515607725858441e-0001,
  130. tpi    = 0.636619772367581343075535053490057448,
  131.  
  132.     /* R0/S0 on [0,2] */
  133. r00 =  -6.250000000000000020842322918309200910191e-0002,
  134. r01 =   1.407056669551897148204830386691427791200e-0003,
  135. r02 =  -1.599556310840356073980727783817809847071e-0005,
  136. r03 =   4.967279996095844750387702652791615403527e-0008,
  137. s01 =   1.915375995383634614394860200531091839635e-0002,
  138. s02 =   1.859467855886309024045655476348872850396e-0004,
  139. s03 =   1.177184640426236767593432585906758230822e-0006,
  140. s04 =   5.046362570762170559046714468225101016915e-0009,
  141. s05 =   1.235422744261379203512624973117299248281e-0011;
  142.  
  143. #define two_129    6.80564733841876926e+038    /* 2^129 */
  144. #define two_m54    5.55111512312578270e-017    /* 2^-54 */
  145. double j1(x) 
  146.     double x;
  147. {
  148.     double z, s,c,ss,cc,r,u,v,y;
  149.     y = fabs(x);
  150.     if (!finite(x))            /* Inf or NaN */
  151.         if (_IEEE && x != x)
  152.             return(x);
  153.         else
  154.             return (copysign(x, zero));
  155.     y = fabs(x);
  156.     if (y >= 2)            /* |x| >= 2.0 */
  157.     {
  158.         s = sin(y);
  159.         c = cos(y);
  160.         ss = -s-c;
  161.         cc = s-c;
  162.         if (y < .5*DBL_MAX) {      /* make sure y+y not overflow */
  163.             z = cos(y+y);
  164.             if ((s*c)<zero) cc = z/ss;
  165.             else         ss = z/cc;
  166.         }
  167.     /*
  168.      * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
  169.      * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)
  170.      */
  171. #if !defined(vax) && !defined(tahoe)
  172.         if (y > two_129)     /* x > 2^129 */
  173.             z = (invsqrtpi*cc)/sqrt(y);
  174.         else
  175. #endif /* defined(vax) || defined(tahoe) */
  176.         {
  177.             u = pone(y); v = qone(y);
  178.             z = invsqrtpi*(u*cc-v*ss)/sqrt(y);
  179.         }
  180.         if (x < 0) return -z;
  181.         else       return  z;
  182.     }
  183.     if (y < 7.450580596923828125e-009) {    /* |x|<2**-27 */
  184.         if(huge+x>one) return 0.5*x;/* inexact if x!=0 necessary */
  185.     }
  186.     z = x*x;
  187.     r =  z*(r00+z*(r01+z*(r02+z*r03)));
  188.     s =  one+z*(s01+z*(s02+z*(s03+z*(s04+z*s05))));
  189.     r *= x;
  190.     return (x*0.5+r/s);
  191. }
  192.  
  193. static double u0[5] = {
  194.   -1.960570906462389484206891092512047539632e-0001,
  195.    5.044387166398112572026169863174882070274e-0002,
  196.   -1.912568958757635383926261729464141209569e-0003,
  197.    2.352526005616105109577368905595045204577e-0005,
  198.    -9.190991580398788465315411784276789663849e-0008,
  199. };
  200. static double v0[5] = {
  201.    1.991673182366499064031901734535479833387e-0002,
  202.    2.025525810251351806268483867032781294682e-0004,
  203.    1.356088010975162198085369545564475416398e-0006,
  204.    6.227414523646214811803898435084697863445e-0009,
  205.    1.665592462079920695971450872592458916421e-0011,
  206. };
  207.  
  208. double y1(x) 
  209.     double x;
  210. {
  211.     double z, s,c,ss,cc,u,v,j1();
  212.     /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
  213.     if (!finite(x))
  214.         if (!_IEEE) return (infnan(EDOM));
  215.         else if (x < 0)
  216.             return(zero/zero);
  217.         else if (x > 0)
  218.             return (0);
  219.         else
  220.             return(x);
  221.     if (x <= 0) {
  222.         if (_IEEE && x == 0) return -one/zero;
  223.         else if(x == 0) return(infnan(-ERANGE));
  224.         else if(_IEEE) return (zero/zero);
  225.         else return(infnan(EDOM));
  226.     }
  227.         if (x >= 2)             /* |x| >= 2.0 */
  228.     {
  229.                 s = sin(x);
  230.                 c = cos(x);
  231.                 ss = -s-c;
  232.                 cc = s-c;
  233.         if (x < .5 * DBL_MAX)    /* make sure x+x not overflow */
  234.         {
  235.                     z = cos(x+x);
  236.                     if ((s*c)>zero) cc = z/ss;
  237.                     else            ss = z/cc;
  238.                 }
  239.         /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0))
  240.          * where x0 = x-3pi/4
  241.          *      Better formula:
  242.          *              cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
  243.          *                      =  1/sqrt(2) * (sin(x) - cos(x))
  244.          *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
  245.          *                      = -1/sqrt(2) * (cos(x) + sin(x))
  246.          * To avoid cancellation, use
  247.          *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
  248.          * to compute the worse one.
  249.          */
  250.                 if (_IEEE && x>two_129)
  251.             z = (invsqrtpi*ss)/sqrt(x);
  252.                 else {
  253.                     u = pone(x); v = qone(x);
  254.                     z = invsqrtpi*(u*ss+v*cc)/sqrt(x);
  255.                 }
  256.                 return z;
  257.         } 
  258.         if (x <= two_m54) {    /* x < 2**-54 */
  259.             return (-tpi/x);
  260.         } 
  261.         z = x*x;
  262.         u = u0[0]+z*(u0[1]+z*(u0[2]+z*(u0[3]+z*u0[4])));
  263.         v = one+z*(v0[0]+z*(v0[1]+z*(v0[2]+z*(v0[3]+z*v0[4]))));
  264.         return (x*(u/v) + tpi*(j1(x)*log(x)-one/x));
  265. }
  266.  
  267. /* For x >= 8, the asymptotic expansions of pone is
  268.  *    1 + 15/128 s^2 - 4725/2^15 s^4 - ...,    where s = 1/x.
  269.  * We approximate pone by
  270.  *     pone(x) = 1 + (R/S)
  271.  * where  R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10
  272.  *       S = 1 + ps0*s^2 + ... + ps4*s^10
  273.  * and
  274.  *    | pone(x)-1-R/S | <= 2  ** ( -60.06)
  275.  */
  276.  
  277. static double pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
  278.    0.0,
  279.    1.171874999999886486643746274751925399540e-0001,
  280.    1.323948065930735690925827997575471527252e+0001,
  281.    4.120518543073785433325860184116512799375e+0002,
  282.    3.874745389139605254931106878336700275601e+0003,
  283.    7.914479540318917214253998253147871806507e+0003,
  284. };
  285. static double ps8[5] = {
  286.    1.142073703756784104235066368252692471887e+0002,
  287.    3.650930834208534511135396060708677099382e+0003,
  288.    3.695620602690334708579444954937638371808e+0004,
  289.    9.760279359349508334916300080109196824151e+0004,
  290.    3.080427206278887984185421142572315054499e+0004,
  291. };
  292.  
  293. static double pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
  294.    1.319905195562435287967533851581013807103e-0011,
  295.    1.171874931906140985709584817065144884218e-0001,
  296.    6.802751278684328781830052995333841452280e+0000,
  297.    1.083081829901891089952869437126160568246e+0002,
  298.    5.176361395331997166796512844100442096318e+0002,
  299.    5.287152013633375676874794230748055786553e+0002,
  300. };
  301. static double ps5[5] = {
  302.    5.928059872211313557747989128353699746120e+0001,
  303.    9.914014187336144114070148769222018425781e+0002,
  304.    5.353266952914879348427003712029704477451e+0003,
  305.    7.844690317495512717451367787640014588422e+0003,
  306.    1.504046888103610723953792002716816255382e+0003,
  307. };
  308.  
  309. static double pr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
  310.    3.025039161373736032825049903408701962756e-0009,
  311.    1.171868655672535980750284752227495879921e-0001,
  312.    3.932977500333156527232725812363183251138e+0000,
  313.    3.511940355916369600741054592597098912682e+0001,
  314.    9.105501107507812029367749771053045219094e+0001,
  315.    4.855906851973649494139275085628195457113e+0001,
  316. };
  317. static double ps3[5] = {
  318.    3.479130950012515114598605916318694946754e+0001,
  319.    3.367624587478257581844639171605788622549e+0002,
  320.    1.046871399757751279180649307467612538415e+0003,
  321.    8.908113463982564638443204408234739237639e+0002,
  322.    1.037879324396392739952487012284401031859e+0002,
  323. };
  324.  
  325. static double pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
  326.    1.077108301068737449490056513753865482831e-0007,
  327.    1.171762194626833490512746348050035171545e-0001,
  328.    2.368514966676087902251125130227221462134e+0000,
  329.    1.224261091482612280835153832574115951447e+0001,
  330.    1.769397112716877301904532320376586509782e+0001,
  331.    5.073523125888185399030700509321145995160e+0000,
  332. };
  333. static double ps2[5] = {
  334.    2.143648593638214170243114358933327983793e+0001,
  335.    1.252902271684027493309211410842525120355e+0002,
  336.    2.322764690571628159027850677565128301361e+0002,
  337.    1.176793732871470939654351793502076106651e+0002,
  338.    8.364638933716182492500902115164881195742e+0000,
  339. };
  340.  
  341. static double pone(x)
  342.     double x;
  343. {
  344.     double *p,*q,z,r,s;
  345.     if (x >= 8.0)                {p = pr8; q= ps8;}
  346.     else if (x >= 4.54545211791992188) {p = pr5; q= ps5;}
  347.     else if (x >= 2.85714149475097656) {p = pr3; q= ps3;}
  348.     else /* if (x >= 2.0) */       {p = pr2; q= ps2;}
  349.     z = one/(x*x);
  350.     r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
  351.     s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
  352.     return (one + r/s);
  353. }
  354.         
  355.  
  356. /* For x >= 8, the asymptotic expansions of qone is
  357.  *    3/8 s - 105/1024 s^3 - ..., where s = 1/x.
  358.  * We approximate pone by
  359.  *     qone(x) = s*(0.375 + (R/S))
  360.  * where  R = qr1*s^2 + qr2*s^4 + ... + qr5*s^10
  361.  *       S = 1 + qs1*s^2 + ... + qs6*s^12
  362.  * and
  363.  *    | qone(x)/s -0.375-R/S | <= 2  ** ( -61.13)
  364.  */
  365.  
  366. static double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
  367.    0.0,
  368.   -1.025390624999927207385863635575804210817e-0001,
  369.   -1.627175345445899724355852152103771510209e+0001,
  370.   -7.596017225139501519843072766973047217159e+0002,
  371.   -1.184980667024295901645301570813228628541e+0004,
  372.   -4.843851242857503225866761992518949647041e+0004,
  373. };
  374. static double qs8[6] = {
  375.    1.613953697007229231029079421446916397904e+0002,
  376.    7.825385999233484705298782500926834217525e+0003,
  377.    1.338753362872495800748094112937868089032e+0005,
  378.    7.196577236832409151461363171617204036929e+0005,
  379.    6.666012326177764020898162762642290294625e+0005,
  380.   -2.944902643038346618211973470809456636830e+0005,
  381. };
  382.  
  383. static double qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
  384.   -2.089799311417640889742251585097264715678e-0011,
  385.   -1.025390502413754195402736294609692303708e-0001,
  386.   -8.056448281239359746193011295417408828404e+0000,
  387.   -1.836696074748883785606784430098756513222e+0002,
  388.   -1.373193760655081612991329358017247355921e+0003,
  389.   -2.612444404532156676659706427295870995743e+0003,
  390. };
  391. static double qs5[6] = {
  392.    8.127655013843357670881559763225310973118e+0001,
  393.    1.991798734604859732508048816860471197220e+0003,
  394.    1.746848519249089131627491835267411777366e+0004,
  395.    4.985142709103522808438758919150738000353e+0004,
  396.    2.794807516389181249227113445299675335543e+0004,
  397.   -4.719183547951285076111596613593553911065e+0003,
  398. };
  399.  
  400. static double qr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
  401.   -5.078312264617665927595954813341838734288e-0009,
  402.   -1.025378298208370901410560259001035577681e-0001,
  403.   -4.610115811394734131557983832055607679242e+0000,
  404.   -5.784722165627836421815348508816936196402e+0001,
  405.   -2.282445407376317023842545937526967035712e+0002,
  406.   -2.192101284789093123936441805496580237676e+0002,
  407. };
  408. static double qs3[6] = {
  409.    4.766515503237295155392317984171640809318e+0001,
  410.    6.738651126766996691330687210949984203167e+0002,
  411.    3.380152866795263466426219644231687474174e+0003,
  412.    5.547729097207227642358288160210745890345e+0003,
  413.    1.903119193388108072238947732674639066045e+0003,
  414.   -1.352011914443073322978097159157678748982e+0002,
  415. };
  416.  
  417. static double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
  418.   -1.783817275109588656126772316921194887979e-0007,
  419.   -1.025170426079855506812435356168903694433e-0001,
  420.   -2.752205682781874520495702498875020485552e+0000,
  421.   -1.966361626437037351076756351268110418862e+0001,
  422.   -4.232531333728305108194363846333841480336e+0001,
  423.   -2.137192117037040574661406572497288723430e+0001,
  424. };
  425. static double qs2[6] = {
  426.    2.953336290605238495019307530224241335502e+0001,
  427.    2.529815499821905343698811319455305266409e+0002,
  428.    7.575028348686454070022561120722815892346e+0002,
  429.    7.393932053204672479746835719678434981599e+0002,
  430.    1.559490033366661142496448853793707126179e+0002,
  431.   -4.959498988226281813825263003231704397158e+0000,
  432. };
  433.  
  434. static double qone(x)
  435.     double x;
  436. {
  437.     double *p,*q, s,r,z;
  438.     if (x >= 8.0)               {p = qr8; q= qs8;}
  439.     else if (x >= 4.54545211791992188) {p = qr5; q= qs5;}
  440.     else if (x >= 2.85714149475097656) {p = qr3; q= qs3;}
  441.     else /* if (x >= 2.0) */       {p = qr2; q= qs2;}
  442.     z = one/(x*x);
  443.     r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
  444.     s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
  445.     return (.375 + r/s)/x;
  446. }
  447.