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