// finance++

struct IP_Binom2
{
  double S0;
  unsigned int N;
  vector<double> dt;
  double sigma;
  double r;
};

typedef struct IP_Binom2 IP_Binom2;

/* Binom *********************/
class Binom2
{
 protected:
  static const unsigned int Q = 2;
  static const int I = 1;
  static const int J = 1;
  static const int K = 0;
  bool recomb;
  
  vector< vector<double> > M;
  vector< vector<double> > P;
  vector< double > disc;

  double S0;
  unsigned int N;
  vector<double> dt;


 public:
  Binom2(IP_Binom2 inIP)
    {
      // dt can either have just one element or
      // it can have N-1 elements, nothing else is allowed

      if(! (inIP.dt.size() == 1 || inIP.dt.size() == (inIP.N - 1)) )
	{
	  fprintf(stderr, "Binom2, invalid dt size %d, must be 1 or N-1.\n", inIP.dt.size());
	  return;
	}

      if(inIP.sigma < 0.0)
	{
	  fprintf(stderr, "Binom2, invalid sigma %f, cannot be less than 0.\n", inIP.sigma);
	  return;
	}

      if(inIP.dt.size() == 1)
	recomb = true;
      else
	recomb = false;

      S0 = inIP.S0;
      N = inIP.N;

      dt = vector<double>(inIP.dt.size());
      
      for(unsigned int i = 0; i < inIP.dt.size(); i++)
	{
	  dt[i] = inIP.dt[i];

	  vector<double> vd = vector<double>(2);

	  // fill in the M vector

	  vd[0] = pow(M_E, (-inIP.sigma) * sqrt(dt[i]));      // Di
	  vd[1] = pow(M_E, inIP.sigma * sqrt(dt[i]));   // Ui

	  /**
	  vd[0] = 1.0 - inIP.sigma * sqrt(dt[i]) + 1.0/2.0 * pow(inIP.sigma, 2.0) * dt[i];
	  vd[1] = 1.0 + inIP.sigma * sqrt(dt[i]) + 1.0/2.0 * pow(inIP.sigma, 2.0) * dt[i];
	  **/
	  
	  M.push_back(vd);
	  
	  // fill in the P vector
	  
	  vd = vector<double>(2);


	  /**
	  vd[1] = (inIP.r * dt[i] - M[i][0]) / (M[i][1] - M[i][0]); 
	  vd[0] = 1 - vd[1];
	  **/

	  vd[0] = (M[i][1] - pow(M_E, inIP.r * dt[i])) / (M[i][1] - M[i][0]);  // PDi
	  vd[1] = (pow(M_E, inIP.r * dt[i]) - M[i][0]) / (M[i][1] - M[i][0]);  // PUi

	  /**
	  vd[1] = 1.0/2.0 * (1.0 + (inIP.r - (pow(inIP.sigma, 2.0) / 2.0) / (inIP.sigma * dt[i])) +
			     (pow(inIP.sigma, 2.0) * dt[i]) / 6.0);
	  vd[0] = 1 - vd[1];
	  **/
	  
	  P.push_back(vd);
	  
	  // printf("P[%d][0] is %f\n", i, P[i][0]);
	  // printf("P[%d][1] is %f\n", i, P[i][1]);
	  
	  disc.push_back(pow(M_E, (-inIP.r) * dt[i]));
	  
	  // printf("disc[0] is %f\n", disc[0]);
	  
	}


    }
};



