Subversion Repositories Kolibri OS

Rev

Blame | Last modification | View Log | Download | RSS feed

  1. /* Copyright (C) 1994 DJ Delorie, see COPYING.DJ for details */
  2. /* @(#)w_jn.c 5.1 93/09/24 */
  3. /*
  4.  * ====================================================
  5.  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  6.  *
  7.  * Developed at SunPro, a Sun Microsystems, Inc. business.
  8.  * Permission to use, copy, modify, and distribute this
  9.  * software is freely granted, provided that this notice
  10.  * is preserved.
  11.  * ====================================================
  12.  */
  13.  
  14. #if defined(LIBM_SCCS) && !defined(lint)
  15. static char rcsid[] = "$Id: w_jn.c,v 1.4 1994/08/10 20:34:42 jtc Exp $";
  16. #endif
  17.  
  18. /*
  19.  * wrapper jn(int n, double x), yn(int n, double x)
  20.  * floating point Bessel's function of the 1st and 2nd kind
  21.  * of order n
  22.  *          
  23.  * Special cases:
  24.  *      y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
  25.  *      y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
  26.  * Note 2. About jn(n,x), yn(n,x)
  27.  *      For n=0, j0(x) is called,
  28.  *      for n=1, j1(x) is called,
  29.  *      for n<x, forward recursion us used starting
  30.  *      from values of j0(x) and j1(x).
  31.  *      for n>x, a continued fraction approximation to
  32.  *      j(n,x)/j(n-1,x) is evaluated and then backward
  33.  *      recursion is used starting from a supposed value
  34.  *      for j(n,x). The resulting value of j(0,x) is
  35.  *      compared with the actual value to correct the
  36.  *      supposed value of j(n,x).
  37.  *
  38.  *      yn(n,x) is similar in all respects, except
  39.  *      that forward recursion is used for all
  40.  *      values of n>1.
  41.  *     
  42.  */
  43.  
  44. #include "math.h"
  45. #include "math_private.h"
  46.  
  47. #ifdef __STDC__
  48.         double jn(int n, double x)      /* wrapper jn */
  49. #else
  50.         double jn(n,x)                  /* wrapper jn */
  51.         double x; int n;
  52. #endif
  53. {
  54. #ifdef _IEEE_LIBM
  55.         return __ieee754_jn(n,x);
  56. #else
  57.         double z;
  58.         z = __ieee754_jn(n,x);
  59.         if(_LIB_VERSION == _IEEE_ || isnan(x) ) return z;
  60.         if(fabs(x)>X_TLOSS) {
  61.             return __kernel_standard((double)n,x,38); /* jn(|x|>X_TLOSS,n) */
  62.         } else
  63.             return z;
  64. #endif
  65. }
  66.  
  67. #ifdef __STDC__
  68.         double yn(int n, double x)      /* wrapper yn */
  69. #else
  70.         double yn(n,x)                  /* wrapper yn */
  71.         double x; int n;
  72. #endif
  73. {
  74. #ifdef _IEEE_LIBM
  75.         return __ieee754_yn(n,x);
  76. #else
  77.         double z;
  78.         z = __ieee754_yn(n,x);
  79.         if(_LIB_VERSION == _IEEE_ || isnan(x) ) return z;
  80.         if(x <= 0.0){
  81.                 if(x==0.0)
  82.                     /* d= -one/(x-x); */
  83.                     return __kernel_standard((double)n,x,12);
  84.                 else
  85.                     /* d = zero/(x-x); */
  86.                     return __kernel_standard((double)n,x,13);
  87.         }
  88.         if(x>X_TLOSS) {
  89.             return __kernel_standard((double)n,x,39); /* yn(x>X_TLOSS,n) */
  90.         } else
  91.             return z;
  92. #endif
  93. }
  94.