#include "option.h"

double option(double K, int T, double S, double r, int N,  
             double u, double d, double *rstock, double *roption,
             int opt)
{

 double dt,n,p,discount;
 double *St,*C;
 int i,j,k;
 double oneoveru,dp,odp;

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


 n = (double) N;
 dt = (double) T / (double) N;
 p  = (exp(r *dt) -d)/ (u-d);
 discount = exp(-r*dt);

#ifdef SANITY_CHECK
 fprintf(stdout, "K=%f T=%d S=%f r=%f N=%d, u=%f d=%f\ndt=%f\
                  p=%f discount=%f\n",K,T,S,r,N,u,d,dt,p,discount);
 fflush(stdout);
#endif


 St[0] = S * pow(u,n);
 for(j=1;j<=N;j++)
   St[j] = St[j-1]*(d/u);
#ifdef OPTION_DEBUG
 PRINTMATRIX("St","Init",St,(N+1));
#endif

if ((opt == AMERICANCALL) || (opt == EUROPEANCALL)) { /* CALL option */
 for(j=0;j<=N;j++)
   C[j] = MAX((0.0), (St[j]-K));
}
else {  /* PUT option */
 for(j=0;j<=N;j++)
   C[j] = MAX((0.0), (K-St[j]));
}
#ifdef OPTION_DEBUG
 PRINTMATRIX("C ","Init",C,(N+1)); 
#endif

 oneoveru=1.0/u;
 dp=discount*p;
 odp=discount*(1.0-p);
if ((opt==AMERICANCALL)) {
 for(i=N-1;i>=0;i--) {
   for(j=0;j<=i;j++) {
      C[j] = dp*C[j] + odp*C[j+1];
      St[j]=St[j] *oneoveru;
      C[j] = MAX( C[j], (St[j]-K));
   }
#ifdef AMERICANOPTION_DEBUG
 PRINTMATRIX("C ","   I",C,(i+1));
#endif
 }
}
if ((opt==EUROPEANCALL)) {
 for(i=N-1;i>=0;i--) {
   for(j=0;j<=i;j++) {
      C[j] = dp*C[j] + odp*C[j+1];
      St[j]=St[j] *oneoveru;
   }
#ifdef EUROPEANOPTION_DEBUG
 PRINTMATRIX("C ","   I",C,(i+1));
#endif
 }
}
if ((opt==EUROPEANPUT )) {
 for(i=N-1;i>=0;i--) {
   for(j=0;j<=i;j++) {
      C[j] = dp*C[j] + odp*C[j+1];
      St[j]=St[j] *oneoveru;
   }
#ifdef EUROPEANOPTION_DEBUG
 PRINTMATRIX("C ","   I",C,(i+1));
#endif
 }
}
if ((opt==AMERICANPUT )) {
 for(i=N-1;i>=0;i--) {
   for(j=0;j<=i;j++) {
      C[j] = dp*C[j] + odp*C[j+1];
      St[j]=St[j] *oneoveru;
      C[j] = MAX( C[j], (K-St[j]));
   }
#ifdef AMERICANOPTION_DEBUG
 PRINTMATRIX("C ","   I",C,(i+1));
#endif
 }
}

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

