chbevl.c 1.6 KB

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