Back to index

glibc  2.9
w_jn.c
Go to the documentation of this 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 "math_private.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 #ifdef NO_LONG_DOUBLE
00066 strong_alias (jn, jnl)
00067 #endif
00068 
00069 
00070 #ifdef __STDC__
00071        double yn(int n, double x)  /* wrapper yn */
00072 #else
00073        double yn(n,x)                     /* wrapper yn */
00074        double x; int n;
00075 #endif
00076 {
00077 #ifdef _IEEE_LIBM
00078        return __ieee754_yn(n,x);
00079 #else
00080        double z;
00081        z = __ieee754_yn(n,x);
00082        if(_LIB_VERSION == _IEEE_ || __isnan(x) ) return z;
00083         if(x <= 0.0){
00084                 if(x==0.0)
00085                     /* d= -one/(x-x); */
00086                     return __kernel_standard((double)n,x,12);
00087                 else
00088                     /* d = zero/(x-x); */
00089                     return __kernel_standard((double)n,x,13);
00090         }
00091        if(x>X_TLOSS) {
00092            return __kernel_standard((double)n,x,39); /* yn(x>X_TLOSS,n) */
00093        } else
00094            return z;
00095 #endif
00096 }
00097 #ifdef NO_LONG_DOUBLE
00098 strong_alias (yn, ynl)
00099 #endif