chbevl.c 1.6 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980
  1. #include "mconf.h"
  2. /* chbevl.c
  3. *
  4. * Evaluate Chebyshev series
  5. *
  6. *
  7. *
  8. * SYNOPSIS:
  9. *
  10. * int N;
  11. * double x, y, coef[N], chebevl();
  12. *
  13. * y = chbevl( x, coef, N );
  14. *
  15. *
  16. *
  17. * DESCRIPTION:
  18. *
  19. * Evaluates the series
  20. *
  21. * N-1
  22. * - '
  23. * y = > coef[i] T (x/2)
  24. * - i
  25. * i=0
  26. *
  27. * of Chebyshev polynomials Ti at argument x/2.
  28. *
  29. * Coefficients are stored in reverse order, i.e. the zero
  30. * order term is last in the array. Note N is the number of
  31. * coefficients, not the order.
  32. *
  33. * If coefficients are for the interval a to b, x must
  34. * have been transformed to x -> 2(2x - b - a)/(b-a) before
  35. * entering the routine. This maps x from (a, b) to (-1, 1),
  36. * over which the Chebyshev polynomials are defined.
  37. *
  38. * If the coefficients are for the inverted interval, in
  39. * which (a, b) is mapped to (1/b, 1/a), the transformation
  40. * required is x -> 2(2ab/x - b - a)/(b-a). If b is infinity,
  41. * this becomes x -> 4a/x - 1.
  42. *
  43. *
  44. *
  45. * SPEED:
  46. *
  47. * Taking advantage of the recurrence properties of the
  48. * Chebyshev polynomials, the routine requires one more
  49. * addition per loop than evaluating a nested polynomial of
  50. * the same degree.
  51. *
  52. */
  53. /* chbevl.c */
  54. /*
  55. Cephes Math Library Release 2.0: April, 1987
  56. Copyright 1985, 1987 by Stephen L. Moshier
  57. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  58. */
  59. double chbevl(x, array, n) double x;
  60. double array[];
  61. int n;
  62. {
  63. double b0, b1, b2, *p;
  64. int i;
  65. p = array;
  66. b0 = *p++;
  67. b1 = 0.0;
  68. i = n - 1;
  69. do {
  70. b2 = b1;
  71. b1 = b0;
  72. b0 = x * b1 - b2 + *p++;
  73. } while (--i);
  74. return (0.5 * (b0 - b2));
  75. }