home *** CD-ROM | disk | FTP | other *** search
/ InfoMagic Source Code 1993 July / THE_SOURCE_CODE_CD_ROM.iso / bsd_srcs / lib / libm / common_source / jn.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-12-16  |  9.1 KB  |  313 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[] = "@(#)jn.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. /*
  69.  * jn(int n, double x), yn(int n, double x)
  70.  * floating point Bessel's function of the 1st and 2nd kind
  71.  * of order n
  72.  *          
  73.  * Special cases:
  74.  *    y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
  75.  *    y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
  76.  * Note 2. About jn(n,x), yn(n,x)
  77.  *    For n=0, j0(x) is called,
  78.  *    for n=1, j1(x) is called,
  79.  *    for n<x, forward recursion us used starting
  80.  *    from values of j0(x) and j1(x).
  81.  *    for n>x, a continued fraction approximation to
  82.  *    j(n,x)/j(n-1,x) is evaluated and then backward
  83.  *    recursion is used starting from a supposed value
  84.  *    for j(n,x). The resulting value of j(0,x) is
  85.  *    compared with the actual value to correct the
  86.  *    supposed value of j(n,x).
  87.  *
  88.  *    yn(n,x) is similar in all respects, except
  89.  *    that forward recursion is used for all
  90.  *    values of n>1.
  91.  *    
  92.  */
  93.  
  94. #include <math.h>
  95. #include <float.h>
  96. #include <errno.h>
  97.  
  98. #if defined(vax) || defined(tahoe)
  99. #define _IEEE    0
  100. #else
  101. #define _IEEE    1
  102. #define infnan(x) (0.0)
  103. #endif
  104.  
  105. extern double j0(),j1(),log(),fabs(),sqrt(),cos(),sin(),y0(),y1();
  106. static double
  107. invsqrtpi= 5.641895835477562869480794515607725858441e-0001,
  108. two  = 2.0,
  109. zero = 0.0,
  110. one  = 1.0;
  111.  
  112. double jn(n,x)
  113.     int n; double x;
  114. {
  115.     int i, sgn;
  116.     double a, b, temp;
  117.     double z, w;
  118.  
  119.     /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
  120.      * Thus, J(-n,x) = J(n,-x)
  121.      */
  122.     /* if J(n,NaN) is NaN */
  123.     if (_IEEE && isnan(x)) return x+x;
  124.     if (n<0){        
  125.         n = -n;
  126.         x = -x;
  127.     }
  128.     if (n==0) return(j0(x));
  129.     if (n==1) return(j1(x));
  130.     sgn = (n&1)&(x < zero);        /* even n -- 0, odd n -- sign(x) */
  131.     x = fabs(x);
  132.     if (x == 0 || !finite (x))     /* if x is 0 or inf */
  133.         b = zero;
  134.     else if ((double) n <= x) {
  135.             /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
  136.         if (_IEEE && x >= 8.148143905337944345e+090) {
  137.                     /* x >= 2**302 */
  138.     /* (x >> n**2) 
  139.      *        Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
  140.      *        Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
  141.      *        Let s=sin(x), c=cos(x), 
  142.      *        xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
  143.      *
  144.      *           n    sin(xn)*sqt2    cos(xn)*sqt2
  145.      *        ----------------------------------
  146.      *           0     s-c         c+s
  147.      *           1    -s-c         -c+s
  148.      *           2    -s+c        -c-s
  149.      *           3     s+c         c-s
  150.      */
  151.         switch(n&3) {
  152.             case 0: temp =  cos(x)+sin(x); break;
  153.             case 1: temp = -cos(x)+sin(x); break;
  154.             case 2: temp = -cos(x)-sin(x); break;
  155.             case 3: temp =  cos(x)-sin(x); break;
  156.         }
  157.         b = invsqrtpi*temp/sqrt(x);
  158.         } else {    
  159.             a = j0(x);
  160.             b = j1(x);
  161.             for(i=1;i<n;i++){
  162.             temp = b;
  163.             b = b*((double)(i+i)/x) - a; /* avoid underflow */
  164.             a = temp;
  165.             }
  166.         }
  167.     } else {
  168.         if (x < 1.86264514923095703125e-009) { /* x < 2**-29 */
  169.     /* x is tiny, return the first Taylor expansion of J(n,x) 
  170.      * J(n,x) = 1/n!*(x/2)^n  - ...
  171.      */
  172.         if (n > 33)    /* underflow */
  173.             b = zero;
  174.         else {
  175.             temp = x*0.5; b = temp;
  176.             for (a=one,i=2;i<=n;i++) {
  177.             a *= (double)i;        /* a = n! */
  178.             b *= temp;        /* b = (x/2)^n */
  179.             }
  180.             b = b/a;
  181.         }
  182.         } else {
  183.         /* use backward recurrence */
  184.         /*             x      x^2      x^2       
  185.          *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
  186.          *            2n  - 2(n+1) - 2(n+2)
  187.          *
  188.          *             1      1        1       
  189.          *  (for large x)   =  ----  ------   ------   .....
  190.          *            2n   2(n+1)   2(n+2)
  191.          *            -- - ------ - ------ - 
  192.          *             x     x         x
  193.          *
  194.          * Let w = 2n/x and h=2/x, then the above quotient
  195.          * is equal to the continued fraction:
  196.          *            1
  197.          *    = -----------------------
  198.          *               1
  199.          *       w - -----------------
  200.          *              1
  201.          *             w+h - ---------
  202.          *               w+2h - ...
  203.          *
  204.          * To determine how many terms needed, let
  205.          * Q(0) = w, Q(1) = w(w+h) - 1,
  206.          * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
  207.          * When Q(k) > 1e4    good for single 
  208.          * When Q(k) > 1e9    good for double 
  209.          * When Q(k) > 1e17    good for quadruple 
  210.          */
  211.         /* determine k */
  212.         double t,v;
  213.         double q0,q1,h,tmp; int k,m;
  214.         w  = (n+n)/(double)x; h = 2.0/(double)x;
  215.         q0 = w;  z = w+h; q1 = w*z - 1.0; k=1;
  216.         while (q1<1.0e9) {
  217.             k += 1; z += h;
  218.             tmp = z*q1 - q0;
  219.             q0 = q1;
  220.             q1 = tmp;
  221.         }
  222.         m = n+n;
  223.         for(t=zero, i = 2*(n+k); i>=m; i -= 2) t = one/(i/x-t);
  224.         a = t;
  225.         b = one;
  226.         /*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
  227.          *  Hence, if n*(log(2n/x)) > ...
  228.          *  single 8.8722839355e+01
  229.          *  double 7.09782712893383973096e+02
  230.          *  long double 1.1356523406294143949491931077970765006170e+04
  231.          *  then recurrent value may overflow and the result will
  232.          *  likely underflow to zero
  233.          */
  234.         tmp = n;
  235.         v = two/x;
  236.         tmp = tmp*log(fabs(v*tmp));
  237.             for (i=n-1;i>0;i--){
  238.                 temp = b;
  239.                 b = ((i+i)/x)*b - a;
  240.                 a = temp;
  241.             /* scale b to avoid spurious overflow */
  242. #            if defined(vax) || defined(tahoe)
  243. #                define BMAX 1e13
  244. #            else
  245. #                define BMAX 1e100
  246. #            endif /* defined(vax) || defined(tahoe) */
  247.             if (b > BMAX) {
  248.                 a /= b;
  249.                 t /= b;
  250.                 b = one;
  251.             }
  252.         }
  253.             b = (t*j0(x)/b);
  254.         }
  255.     }
  256.     return ((sgn == 1) ? -b : b);
  257. }
  258. double yn(n,x) 
  259.     int n; double x;
  260. {
  261.     int i, sign;
  262.     double a, b, temp;
  263.  
  264.     /* Y(n,NaN), Y(n, x < 0) is NaN */
  265.     if (x <= 0 || (_IEEE && x != x))
  266.         if (_IEEE && x < 0) return zero/zero;
  267.         else if (x < 0)     return (infnan(EDOM));
  268.         else if (_IEEE)     return -one/zero;
  269.         else            return(infnan(-ERANGE));
  270.     else if (!finite(x)) return(0);
  271.     sign = 1;
  272.     if (n<0){
  273.         n = -n;
  274.         sign = 1 - ((n&1)<<2);
  275.     }
  276.     if (n == 0) return(y0(x));
  277.     if (n == 1) return(sign*y1(x));
  278.     if(_IEEE && x >= 8.148143905337944345e+090) { /* x > 2**302 */
  279.     /* (x >> n**2) 
  280.      *        Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
  281.      *        Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
  282.      *        Let s=sin(x), c=cos(x), 
  283.      *        xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
  284.      *
  285.      *           n    sin(xn)*sqt2    cos(xn)*sqt2
  286.      *        ----------------------------------
  287.      *           0     s-c         c+s
  288.      *           1    -s-c         -c+s
  289.      *           2    -s+c        -c-s
  290.      *           3     s+c         c-s
  291.      */
  292.         switch (n&3) {
  293.             case 0: temp =  sin(x)-cos(x); break;
  294.             case 1: temp = -sin(x)-cos(x); break;
  295.             case 2: temp = -sin(x)+cos(x); break;
  296.             case 3: temp =  sin(x)+cos(x); break;
  297.         }
  298.         b = invsqrtpi*temp/sqrt(x);
  299.     } else {
  300.         a = y0(x);
  301.         b = y1(x);
  302.     /* quit if b is -inf */
  303.         for (i = 1; i < n && !finite(b); i++){
  304.         temp = b;
  305.         b = ((double)(i+i)/x)*b - a;
  306.         a = temp;
  307.         }
  308.     }
  309.     if (!_IEEE && !finite(b))
  310.         return (infnan(-sign * ERANGE));
  311.     return ((sign > 0) ? b : -b);
  312. }
  313.