w_jnl.c 2.7 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495
  1. /* w_jnl.c -- long double version of w_jn.c.
  2. * Conversion to long double by Ulrich Drepper,
  3. * Cygnus Support, drepper@cygnus.com.
  4. */
  5. /*
  6. * ====================================================
  7. * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  8. *
  9. * Developed at SunPro, a Sun Microsystems, Inc. business.
  10. * Permission to use, copy, modify, and distribute this
  11. * software is freely granted, provided that this notice
  12. * is preserved.
  13. * ====================================================
  14. */
  15. #if defined(LIBM_SCCS) && !defined(lint)
  16. static char rcsid[] = "$NetBSD: $";
  17. #endif
  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. #include <math.h>
  44. #include "math_private.h"
  45. #if !defined __NO_LONG_DOUBLE_MATH
  46. # ifndef __DO_XSI_MATH__
  47. long double
  48. jnl(int n, long double x) /* wrapper jnl */
  49. {
  50. # if __UCLIBC_HAS_FENV__
  51. long double z;
  52. z = (long double) __ieee754_jn(n, (double) x);
  53. if (_LIB_VERSION == _IEEE_
  54. || _LIB_VERSION == _POSIX_
  55. || isnan(x))
  56. return z;
  57. if(fabsl(x)>X_TLOSS) {
  58. return __kernel_standard_l((double)n,x,238); /* jn(|x|>X_TLOSS,n) */
  59. } else
  60. return z;
  61. # else
  62. return (long double) __ieee754_jn(n, (double) x);
  63. # endif /* __UCLIBC_HAS_FENV__ */
  64. }
  65. long double
  66. ynl(int n, long double x) /* wrapper ynl */
  67. {
  68. # if __UCLIBC_HAS_FENV__
  69. long double z;
  70. z = (long double) __ieee754_yn(n,(double) x);
  71. if(_LIB_VERSION == _IEEE_ || isnan(x) ) return z;
  72. if(x <= 0.0){
  73. if(x==0.0)
  74. /* d= -one/(x-x); */
  75. return __kernel_standard_l((double)n,x,212);
  76. else
  77. /* d = zero/(x-x); */
  78. return __kernel_standard_l((double)n,x,213);
  79. }
  80. if(x>X_TLOSS && _LIB_VERSION != _POSIX_) {
  81. return __kernel_standard_l((double)n,x,239); /* yn(x>X_TLOSS,n) */
  82. } else
  83. return z;
  84. # else
  85. return (long double) __ieee754_yn(n,(double) x);
  86. # endif /* __UCLIBC_HAS_FENV__ */
  87. }
  88. # endif /* __DO_XSI_MATH__ */
  89. #endif /* __NO_LONG_DOUBLE_MATH */