jn.c 2.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123
  1. /*
  2. * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
  3. * See the copyright notice in the ACK home directory, in the file "Copyright".
  4. *
  5. * Author: Ceriel J.H. Jacobs
  6. */
  7. /* $Id$ */
  8. #include <math.h>
  9. #include <errno.h>
  10. extern int errno;
  11. double
  12. yn(n, x)
  13. double x;
  14. {
  15. /* Use y0, y1, and the recurrence relation
  16. y(n+1,x) = 2*n*y(n,x)/x - y(n-1, x).
  17. According to Hart & Cheney, this is stable for all
  18. x, n.
  19. Also use: y(-n,x) = (-1)^n * y(n, x)
  20. */
  21. int negative = 0;
  22. extern double y0(), y1();
  23. double yn1, yn2;
  24. register int i;
  25. if (x <= 0) {
  26. errno = EDOM;
  27. return -HUGE;
  28. }
  29. if (n < 0) {
  30. n = -n;
  31. negative = (n % 2);
  32. }
  33. if (n == 0) return y0(x);
  34. if (n == 1) return y1(x);
  35. yn2 = y0(x);
  36. yn1 = y1(x);
  37. for (i = 1; i < n; i++) {
  38. double tmp = yn1;
  39. yn1 = (i*2)*yn1/x - yn2;
  40. yn2 = tmp;
  41. }
  42. if (negative) return -yn1;
  43. return yn1;
  44. }
  45. double
  46. jn(n, x)
  47. double x;
  48. {
  49. /* Unfortunately, according to Hart & Cheney, the recurrence
  50. j(n+1,x) = 2*n*j(n,x)/x - j(n-1,x) is unstable for
  51. increasing n, except when x > n.
  52. However, j(n,x)/j(n-1,x) = 2/(2*n-x*x/(2*(n+1)-x*x/( ....
  53. (a continued fraction).
  54. We can use this to determine KJn and KJn-1, where K is a
  55. normalization constant not yet known. This enables us
  56. to determine KJn-2, ...., KJ1, KJ0. Now we can use the
  57. J0 or J1 approximation to determine K.
  58. Use: j(-n, x) = (-1)^n * j(n, x)
  59. j(n, -x) = (-1)^n * j(n, x)
  60. */
  61. extern double j0(), j1();
  62. if (n < 0) {
  63. n = -n;
  64. x = -x;
  65. }
  66. if (n == 0) return j0(x);
  67. if (n == 1) return j1(x);
  68. if (x > n) {
  69. /* in this case, the recurrence relation is stable for
  70. increasing n, so we use that.
  71. */
  72. double jn2 = j0(x), jn1 = j1(x);
  73. register int i;
  74. for (i = 1; i < n; i++) {
  75. double tmp = jn1;
  76. jn1 = (2*i)*jn1/x - jn2;
  77. jn2 = tmp;
  78. }
  79. return jn1;
  80. }
  81. {
  82. /* we first compute j(n,x)/j(n-1,x) */
  83. register int i;
  84. double quotient = 0.0;
  85. double xsqr = x*x;
  86. double jn1, jn2;
  87. for (i = 20; /* ??? how many do we need ??? */
  88. i > 0; i--) {
  89. quotient = xsqr/(2*(i+n) - quotient);
  90. }
  91. quotient = x / (2*n - quotient);
  92. jn1 = quotient;
  93. jn2 = 1.0;
  94. for (i = n-1; i > 0; i--) {
  95. /* recurrence relation is stable for decreasing n
  96. */
  97. double tmp = jn2;
  98. jn2 = (2*i)*jn2/x - jn1;
  99. jn1 = tmp;
  100. }
  101. /* So, now we have K*Jn = quotient and K*J0 = jn2.
  102. Now it is easy; compute real j0, this gives K = jn2/j0,
  103. and this then gives Jn = quotient/K = j0 * quotient / jn2.
  104. */
  105. return j0(x)*quotient/jn2;
  106. }
  107. }