chbevlf.c 1.6 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586
  1. /* chbevlf.c
  2. *
  3. * Evaluate Chebyshev series
  4. *
  5. *
  6. *
  7. * SYNOPSIS:
  8. *
  9. * int N;
  10. * float x, y, coef[N], chebevlf();
  11. *
  12. * y = chbevlf( 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. #ifdef ANSIC
  59. float chbevlf( float x, float *array, int n )
  60. #else
  61. float chbevlf( x, array, n )
  62. float x;
  63. float *array;
  64. int n;
  65. #endif
  66. {
  67. float b0, b1, b2, *p;
  68. int i;
  69. p = array;
  70. b0 = *p++;
  71. b1 = 0.0;
  72. i = n - 1;
  73. do
  74. {
  75. b2 = b1;
  76. b1 = b0;
  77. b0 = x * b1 - b2 + *p++;
  78. }
  79. while( --i );
  80. return( 0.5*(b0-b2) );
  81. }