#include "option.h"
#define C(i,j)  C[(i)*(2*Nj+1)+(j)]
#define OPTION_DEBUG
#undef OPTION_DEBUG

double exfdm(double K, double T, double S, double r, int N,  
             double sigma, double dx, double dividend, 
             double *rstock, double *roption, int opt, int Nj)
{

 double dt,n,p,discount,invexpdx;
 double dxu,dxd,pu,pm,pd,v;
 double *St,*C;
 int i,j,k;
 double delta,gamma;
 double a1,a2,b2,c1,c2;
 double dpd,dpu,dpm,oneoveru;

 St = (double *) malloc((2*Nj+1)*sizeof(double));
#ifdef SANITY_CHECK
 check_if_null((void *)St,"St","option");
#endif
 C  = (double *) malloc(2*(2*Nj+1)*sizeof(double));
#ifdef SANITY_CHECK
 check_if_null((void *)St,"C","option");
#endif


 dt = (double) T / (double) N;
 n = (double) N;
 v = r - dividend - 0.5*sigma*sigma;
 invexpdx=1.0/exp(dx);
 pu  = 0.5 * dt *((sigma/dx)*(sigma/dx)+v/dx);
 pm  = 1.0 - dt* (sigma/dx)*(sigma/dx) - r*dt;
 pd  = 0.5 * dt *((sigma/dx)*(sigma/dx)-v/dx);
 discount = exp(-r*dt);

#ifdef OPTION_DEBUG
 fprintf(stdout, "dt=%f n=%f v=%f expdx=%f iexpdx=%f (p)=(%f,%f,%f) discount=%f\n",dt,n,v,1.0/invexpdx,invexpdx,pu,pm,pd,discount);
 fflush(stdout);
#endif


 St[0] = S * exp(Nj*dx);
 for(j=1;j<=2*Nj;j++)
   St[j] = St[j-1]*invexpdx;
#ifdef OPTION_DEBUG
 PRINTMAT2("St",16,St,(2*N+1),0);
#endif

 if ((opt == AMERICANCALL) || (opt == EUROPEANCALL)) {
    for(j=0;j<=2*Nj;j++)
      C(0,j) = MAX((0.0), (St[j]-K));
 }    
 else {  /* It is an AMERICANPUT or a EUROPEANPUT */
    for(j=0;j<=2*Nj;j++)
      C(0,j) = MAX((0.0), (K-St[j]));
 }
#ifdef OPTION_DEBUG
/*
 PRINTMATRIX("C ","Init",C,(2*Nj+1)); NEWLINE;
*/
 PRINTMAT2("C ",16,C,2*Nj+1,0); 
#endif

 oneoveru= invexpdx;
 dpd=pd;
 dpu=pu;
 dpm=pm;
if ((opt == EUROPEANCALL) ) {
 for(i=N-1;i>=0;i--) {
   for(j=0;j<2*Nj-1;j++) {
      C(1,j) = dpu*C((0),j) + dpm*C((0),(j+1))+dpd*C((0),(j+2));
      St[j]= St[j]*oneoveru;
   }
   C(1,(2*Nj-1))=C(1, (2*Nj-2));
   C(1,(2*Nj))=C(1, (2*Nj-1));
   for(j=0;j<=2*Nj;j++) {
     C(0,j) = C(1,j); /* Change accordingly for American options*/
   }
#ifdef OPTION_DEBUG
/*
  PRINTMAT("St",i,St,2*i+1); 
*/
  PRINTMAT2("C ",i,C,2*Nj+1,0); 
#endif

 }
}
if ((opt == EUROPEANPUT) ) {
 for(i=N-1;i>=0;i--) {
   for(j=0;j<2*Nj-1;j++) {
      C(1,j) = dpu*C((0),j) + dpm*C((0),(j+1))+dpd*C((0),(j+2));
      St[j]= St[j]*oneoveru;
   }
   C(1,(2*Nj-1))=C(1,(2*Nj-2)) + (St[2*Nj-2]-St[2*Nj-1]);
   C(1,(2*Nj))=C(1,(2*Nj-1)) + (St[2*Nj-1]-St[2*Nj]);
   for(j=0;j<=2*Nj;j++) {
     C(0,j) = C(1,j); /* Change accordingly for American options*/
   }
#ifdef OPTION_DEBUG
/*
  PRINTMAT("St",i,St,2*i+1); 
*/
  PRINTMAT2("C ",i,C,2*Nj+1,0); 
#endif

 }
}
if ((opt == AMERICANCALL) ) {
 for(i=N-1;i>=0;i--) {
   for(j=0;j<2*Nj-1;j++) {
      C(1,j) = dpu*C((0),j) + dpm*C((0),(j+1))+dpd*C((0),(j+2));
      St[j]= St[j]*oneoveru;
   }
   C(1,(2*Nj-1))=C(1, (2*Nj-2));
   C(1,(2*Nj))=C(1, (2*Nj-1));
   for(j=0;j<=2*Nj;j++) {
     C(0,j) = MAX(C(1,j), (St[j]-K)); 
   }
#ifdef OPTION_DEBUG
/*
  PRINTMAT("St",i,St,2*i+1); 
*/
  PRINTMAT2("C ",i,C,2*Nj+1,0); 
#endif

 }
}
if ((opt == AMERICANPUT) ) {
 for(i=N-1;i>=0;i--) {
   for(j=0;j<2*Nj-1;j++) {
      C(1,j) = dpu*C((0),j) + dpm*C((0),(j+1))+dpd*C((0),(j+2));
      St[j]= St[j]*oneoveru;
   }
   C(1,(2*Nj-1))=C(1,(2*Nj-2)) + (St[2*Nj-2]-St[2*Nj-1]);
   C(1,(2*Nj))=C(1,(2*Nj-1)) + (St[2*Nj-1]-St[2*Nj]);
   for(j=0;j<=2*Nj;j++) {
     C(0,j) = MAX(C(1,j), (K-St[j])); 
   }
#ifdef OPTION_DEBUG
/*
  PRINTMAT("St",i,St,2*i+1); 
  PRINTMAT2("C*",i,C,2*Nj+1,(2*Nj+1)); 
*/
  PRINTMAT2("C ",i,C,2*Nj+1,0); 
#endif

 }
}


*rstock= St[0]; *roption= C(0,0);
return(*roption);
}



