// L-42 MCS 360 Wed 1 Dec 2010 : chebpoly.cpp 

/* We evaluate the n-th Chebyshev polynomial in four ways:
   (1) with the plain recursive definition
   (2) memoization with a static vector
   (3) using a stack of function calls with coefficients
   (4) the iterative version */

#define verbose 1

#include<iostream>
#include<iomanip>
#include<vector>
#include<utility>
#include<stack>
using namespace std;

double rC ( int n, double x );  /* recursive version */
double mC ( int n, double x );  /* with memoization up to n = 50 */
double sC ( int n, double x );  /* version with a stack */
double iC ( int n, double x );  /* iterative version */

int main ( void )
{
   cout << "Give n : "; int n;    cin >> n;
   cout << "Give x : "; double x; cin >> x;
   mC(-1,0.0);

   cout << scientific << setprecision(15);
   cout << " recursive C(" << n << "," << x << ") : " << rC(n,x) << endl;
   cout << "with stack C(" << n << "," << x << ") : " << sC(n,x) << endl;
   cout << "  memoized C(" << n << "," << x << ") : " << mC(n,x) << endl;
   cout << " iterative C(" << n << "," << x << ") : " << iC(n,x) << endl;

   return 0;
}

double rC ( int n, double x )
{
   if(n==0)
      return 1.0;
   else if(n==1)
      return x;
   else
      return 2.0*x*rC(n-1,x) - rC(n-2,x);
}

double mC ( int n, double x )
{
   static vector<double> v;
   if(n==-1)
   {
      v.reserve(51);
      for(int i=0; i<51; i++) v[i] = 2.0;
      return 0.0;
   }
   else
   {
      if(v[n] == 2)
      {
         if(n==0)
            v[n] = 1.0;
         else if(n==1)
            v[n] = x;
         else
            v[n] = 2*x*mC(n-1,x) - mC(n-2,x);
      }
      return v[n];
   }
}

double sC ( int n, double x )
{
   double result = 0.0;
   stack< pair<int,double> > s;
  
   pair<int,double> p;
   p.first = n; p.second = 1.0;
   s.push(p);

   while(!s.empty())
   {
      pair<int,double> q = s.top(); s.pop();

      if(verbose) cout << "<<< popping degree " << q.first 
                       << " and coefficient " << q.second << endl;

      if(q.first == 0)
         result = result + q.second*1.0;
      else if(q.first == 1)
         result = result + q.second*x;
      else // (q.first > 1)
      {
         pair<int,double> q1,q2;
         q1.first = q.first-1; q1.second = q.second*(2*x);
         q2.first = q.first-2; q2.second = q.second*(-1);
         if(verbose) cout << ">>> pushing degree " << q1.first
                          << " and coefficient " << q1.second << endl;
         s.push(q1);
         if(verbose) cout << ">>> pushing degree " << q2.first
                          << " and coefficient " << q2.second << endl;
         s.push(q2);
      }
   }
   return result;
}     

double iC ( int n, double x )
{
   if(n==0)
      return 1;
   else if(n==1)
      return x;
   else
   {
      int i;
      double a,b,c;

      a = 1.0;
      b = x;
      for(i=2; i<=n; i++)
      {  
         c = 2.0*x*b - a;
         a = b;
         b = c;
      }
      return b;
   }
}

