/****************************************************
*  Program to demonstrate Chebyshev economization   *
* ------------------------------------------------- *
* Reference: BASIC Scientific Subroutines, Vol. II  *
* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1]. *
*                                                   *
*              C++ Version by J.-P. Moreau, Paris.  *
*                      (www.jpmoreau.fr)            *
* ------------------------------------------------- *
* SAMPLE RUN:                                       *
*                                                   *
* What is the degree of the input polynomial: ? 15  *
*                                                   *
* What is the degree of the desired economized      *
* polynomial: ? 9                                   *
*                                                   *
* What is the range of the input polynomial: ? 1.57 *
*                                                   *
* Input the coefficients:                           *
* C( 0) = ? 0                                       *
* C( 1) = ? 1                                       *
* C( 2) = ? 0                                       *
* C( 3) = ? -.166666666                             *
* C( 4) = ? 0                                       *
* C( 5) = ? .00833333333                            *
* C( 6) = ? 0                                       *
* C( 7) = ? -.0001984127                            *
* C( 8) = ? 0                                       *
* C( 9) = ? .000002755732                           *
* C( 10) = ? 0                                      *
* C( 11) = ? -.000000025052109                      *
* C( 12) = ? 0                                      *
* C( 13) = ? .00000000016059045                     *
* C( 14) = ? 0                                      *
* C( 15) = ? -.00000000000076471635                 *
*                                                   *
* The Chebyshev series coefficients are:            *
*                                                   *
* A( 0) =  0.0000000000                             *
* A( 1) =  1.1334708982                             *
* A( 2) =  0.0000000000                             *
* A( 3) = -0.1378841454                             *
* A( 4) =  0.0000000000                             *
* A( 5) =  0.0044798168                             *
* A( 6) =  0.0000000000                             *
* A( 7) = -0.0000674667                             *
* A( 8) =  0.0000000000                             *
* A( 9) =  0.0000005865                             *
* A(10) =  0.0000000000                             *
* A(11) = -0.0000000033                             *
* A(12) =  0.0000000000                             *
* A(13) =  0.0000000000                             *
* A(14) =  0.0000000000                             *
* A(15) =  0.0000000000                             *
*                                                   *
* The economized polynomial coefficients are:       *
*                                                   *
* C( 0) =  0.0000000000                             *
* C( 1) =  0.9999999767                             *
* C( 2) =  0.0000000000                             *
* C( 3) = -1.6666647620                             *
* C( 4) =  0.0000000000                             *
* C( 5) =  0.0083329009                             *
* C( 6) =  0.0000000000                             *
* C( 7) = -0.0001980098                             *
* C( 8) =  0.0000000000                             *
* C( 9) =  0.0000025907                             *
*                                                   *
*****************************************************/
#include <stdio.h>
#include <math.h>

#define SIZE  40

double B[SIZE+1][SIZE+1];
double A[SIZE+1], C[SIZE+1];

int    i,m,m1;
double x0;


/*******************************************************
* Chebyshev series coefficients evaluation subroutine  *
* The order of the polynomial is n. The coefficients   *
* are returned in the array B[i,j), i is the degree of *
* the polynomial, j is the coefficient order.          *
*******************************************************/
void Cheby_Ser(int n)  {
  int i, j;
  // Establish t0 and t1 coefficients
  B[0][0] = 1; B[1][0] = 0; B[1][1] = 1;
  // Return if order is less than two
  if (n < 2) return;
  for (i = 2; i < n+1; i++)  {
    for (j = 1; j < i+1; j++)
      // Basic recursion relation
      B[i][j] = 2 * B[i - 1][j - 1] - B[i - 2][j];
    B[i][0] = -B[i - 2][0];
  }
}


/***********************************************************
* Chebyshev economization subroutine. The program takes    *
* the input polynomial coefficients, C(i), and returns the *
* Chebyshev series coefficients, A(i). The degree of the   *
* series passed to the routine is m. The degree of the     *
* series returned is m1. The maximum range of x is x0 used *
* for scaling. Note that the input series coefficients are *
* nulled during the process, and then set equal to the     *
* economized series coefficients.                          *
***********************************************************/
void Cheby_Econ()  {
  double bb;  int i,j,l,n;
  // Start by scaling the input coefficients according to C(i)
  bb = x0;
  for (i = 1; i < m+1; i++)  {
    C[i] = C[i] * bb; bb = bb * x0;
  }
  // Call Chebyshev series coefficients subroutine
  for (n = m; n > -1; n--)  {
    Cheby_Ser(n);
    A[n] = C[n] / B[n][n];
    for (l = 0; l < n+1; l++)
      //Chebyshev series of order l is substracted out of the polynomial
      C[l] = C[l] - A[n] * B[n][l];
  }
  //Perform truncation
  for (i = 0; i < m1+1; i++)
    for (j = 0; j < i+1; j++)
      C[j] = C[j] + A[i] * B[i][j];
  // Convert back to the interval X0
  bb = 1.0 / x0;
  for (i = 1; i < m1+1; i++)  {
    C[i] = C[i] * bb;
    bb = bb / x0;
  }
}


void main()  {
  printf("\n What is the degree of the input polynomial: "); scanf("%d",&m);
  printf("\n What is the degree of the desired economized polynomial: "); scanf("%d",&m1);
  printf("\n What is the range of input polynomial: "); scanf("%lf",&x0);
  printf("\n\n Input the coefficients:\n\n");
  
  for (i = 0; i < m+1; i++) {
    printf(" C(%d) = ",i); scanf("%lf",&C[i]);
  }
  printf("\n");

  Cheby_Econ();

  printf("\n The Chebyshev series coefficients are:\n\n");
  for (i = 0; i < m+1; i++) 
    printf(" A(%d) = %12.9f\n",i,A[i]); 
  printf("\n");
  
  if (m+m1 > 12) getchar();

  printf(" The economized polynomial coefficients are:\n\n");
  for (i = 0; i < m1+1; i++)
    printf(" C(%d) = %12.9f\n",i,C[i]); 
  printf("\n");
}

// End of file chebecon.cpp