lib/libm/chbevl.c

   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 \f/*                                                     chbevl.c        */
  53
  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
  60 double chbevl( x, array, n )
  61 double x;
  62 double array[];
  63 int n;
  64 {
  65 double b0, b1, b2, *p;
  66 int i;
  67
  68 p = array;
  69 b0 = *p++;
  70 b1 = 0.0;
  71 i = n - 1;
  72
  73 do
  74         {
  75         b2 = b1;
  76         b1 = b0;
  77         b0 = x * b1  -  b2  + *p++;
  78         }
  79 while( --i );
  80
  81 return( 0.5*(b0-b2) );
  82 }