yn.c 1.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106
  1. /* yn.c
  2. *
  3. * Bessel function of second kind of integer order
  4. *
  5. *
  6. *
  7. * SYNOPSIS:
  8. *
  9. * double x, y, yn();
  10. * int n;
  11. *
  12. * y = yn( n, x );
  13. *
  14. *
  15. *
  16. * DESCRIPTION:
  17. *
  18. * Returns Bessel function of order n, where n is a
  19. * (possibly negative) integer.
  20. *
  21. * The function is evaluated by forward recurrence on
  22. * n, starting with values computed by the routines
  23. * y0() and y1().
  24. *
  25. * If n = 0 or 1 the routine for y0 or y1 is called
  26. * directly.
  27. *
  28. *
  29. *
  30. * ACCURACY:
  31. *
  32. *
  33. * Absolute error, except relative
  34. * when y > 1:
  35. * arithmetic domain # trials peak rms
  36. * DEC 0, 30 2200 2.9e-16 5.3e-17
  37. * IEEE 0, 30 30000 3.4e-15 4.3e-16
  38. *
  39. *
  40. * ERROR MESSAGES:
  41. *
  42. * message condition value returned
  43. * yn singularity x = 0 MAXNUM
  44. * yn overflow MAXNUM
  45. *
  46. * Spot checked against tables for x, n between 0 and 100.
  47. *
  48. */
  49. /*
  50. Cephes Math Library Release 2.8: June, 2000
  51. Copyright 1984, 1987, 2000 by Stephen L. Moshier
  52. */
  53. #include "mconf.h"
  54. #ifdef ANSIPROT
  55. extern double y0(double);
  56. extern double y1(double);
  57. extern double log(double);
  58. #else
  59. double y0(), y1(), log();
  60. #endif
  61. extern double MAXNUM, MAXLOG;
  62. double yn(n, x) int n;
  63. double x;
  64. {
  65. double an, anm1, anm2, r;
  66. int k, sign;
  67. if (n < 0) {
  68. n = -n;
  69. if ((n & 1) == 0) /* -1**n */
  70. sign = 1;
  71. else
  72. sign = -1;
  73. } else
  74. sign = 1;
  75. if (n == 0)
  76. return (sign * y0(x));
  77. if (n == 1)
  78. return (sign * y1(x));
  79. /* test for overflow */
  80. if (x <= 0.0) {
  81. mtherr("yn", SING);
  82. return (-MAXNUM);
  83. }
  84. /* forward recurrence on n */
  85. anm2 = y0(x);
  86. anm1 = y1(x);
  87. k = 1;
  88. r = 2 * k;
  89. do {
  90. an = r * anm1 / x - anm2;
  91. anm2 = anm1;
  92. anm1 = an;
  93. r += 2.0;
  94. ++k;
  95. } while (k < n);
  96. return (sign * an);
  97. }