RTAI API: base/math/w_jn.c Source File

00001 /* @(#)w_jn.c 5.1 93/09/24 */ 00002 /* 00003 * ==================================================== 00004 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 00005 * 00006 * Developed at SunPro, a Sun Microsystems, Inc. business. 00007 * Permission to use, copy, modify, and distribute this 00008 * software is freely granted, provided that this notice 00009 * is preserved. 00010 * ==================================================== 00011 */ 00012 00013 #if defined(LIBM_SCCS) && !defined(lint) 00014 static char rcsid[] = "$NetBSD: w_jn.c,v 1.6 1995/05/10 20:49:19 jtc Exp $"; 00015 #endif 00016 00017 /* 00018 * wrapper jn(int n, double x), yn(int n, double x) 00019 * floating point Bessel's function of the 1st and 2nd kind 00020 * of order n 00021 * 00022 * Special cases: 00023 * y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal; 00024 * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal. 00025 * Note 2. About jn(n,x), yn(n,x) 00026 * For n=0, j0(x) is called, 00027 * for n=1, j1(x) is called, 00028 * for n<x, forward recursion us used starting 00029 * from values of j0(x) and j1(x). 00030 * for n>x, a continued fraction approximation to 00031 * j(n,x)/j(n-1,x) is evaluated and then backward 00032 * recursion is used starting from a supposed value 00033 * for j(n,x). The resulting value of j(0,x) is 00034 * compared with the actual value to correct the 00035 * supposed value of j(n,x). 00036 * 00037 * yn(n,x) is similar in all respects, except 00038 * that forward recursion is used for all 00039 * values of n>1. 00040 * 00041 */ 00042 00043 #include "math.h" 00044 #include "mathP.h" 00045 00046 #ifdef __STDC__ 00047 double jn(int n, double x) /* wrapper jn */ 00048 #else 00049 double jn(n,x) /* wrapper jn */ 00050 double x; int n; 00051 #endif 00052 { 00053 #ifdef _IEEE_LIBM 00054 return __ieee754_jn(n,x); 00055 #else 00056 double z; 00057 z = __ieee754_jn(n,x); 00058 if(_LIB_VERSION == _IEEE_ || isnan(x) ) return z; 00059 if(fabs(x)>X_TLOSS) { 00060 return __kernel_standard((double)n,x,38); /* jn(|x|>X_TLOSS,n) */ 00061 } else 00062 return z; 00063 #endif 00064 } 00065 00066 #ifdef __STDC__ 00067 double yn(int n, double x) /* wrapper yn */ 00068 #else 00069 double yn(n,x) /* wrapper yn */ 00070 double x; int n; 00071 #endif 00072 { 00073 #ifdef _IEEE_LIBM 00074 return __ieee754_yn(n,x); 00075 #else 00076 double z; 00077 z = __ieee754_yn(n,x); 00078 if(_LIB_VERSION == _IEEE_ || isnan(x) ) return z; 00079 if(x <= 0.0){ 00080 if(x==0.0) 00081 /* d= -one/(x-x); */ 00082 return __kernel_standard((double)n,x,12); 00083 else 00084 /* d = zero/(x-x); */ 00085 return __kernel_standard((double)n,x,13); 00086 } 00087 if(x>X_TLOSS) { 00088 return __kernel_standard((double)n,x,39); /* yn(x>X_TLOSS,n) */ 00089 } else 00090 return z; 00091 #endif 00092 }