sysdeps/libm-ieee754/w_jn.c

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