#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<vector>
#include<iostream>
using namespace std;


double MAX(double a,double b)
{
	return a>b?a:b;
}

int MIN(int a, int b)
{
    return a>b?b:a;
}

void interpolate(double* target, double max, double min, int k){
	int i;
	for(i = 0; i <= k; i++){
		target[i] =  min*i/k + max*(k-i)/k;
	}
}

int main()
{
	double S,sigma,X,year,delta,m,r;
	int n,k;
	double u,d,R,p;
	int i,j,l,q,tmp;
	double max_aver,min_aver,Au,Ad,Cu,Cd,x,C1,C2;
	
	cout << "enter stock price" << endl;
	cin >> S;
	cout << "enter strike price" << endl;
	cin >> X;
	cout << "enter year" << endl;
	cin >> year;
	cout << "the number of period" << endl;
	cin >> n;
	cout << "enter volatility " << endl;
	cin >> sigma;
	cout << "enter interest rate " << endl;
	cin >> r;
	cout << "enter the number of running average on each node" << endl;
	cin >> k;
	
	sigma/=100;
	r/=100;
	
	u = exp(sigma * sqrt(year/n));
	d = exp(-sigma * sqrt(year/n));
	R = exp(r * year / n);
	p = (R - d) / (u - d);
	
	double** C = new double*[n+1];
	for(int i=0;i<=n;i++)
	{
		C[i] = new double[k+1];
	}
	double** A = new double*[n+1];
	for(int i=0;i<=n;i++)
	{
		A[i] = new double[k+1];
	}
	double* put_bucket = new double[k+1];
	
	for(i = 0; i <= n; i++){
		max_aver = S*(1-pow(u,n-i+1))/(1-u) + S*pow(u,n-i)*d*(1-pow(d,i)  )/(1-d);
		max_aver /= (n+1);
		min_aver = S*(1-pow(d,i+1))/(1-d) + S*pow(d,i)*u*(1-pow(u,n-i))/(1-u);
        min_aver /= (n+1);
        interpolate(A[i],max_aver,min_aver,k);
		for(j = 0; j <=k; j++)
			C[i][j] = MAX((double)0,X-A[i][j]);
	}
	
	for(i = n-1; i >=0; i--){
		for(j = 0; j <= i; j++){
			max_aver = S*(1-pow(u,i-j+1))/(1-u) + S*pow(u,i-j)*d*(1-pow(d,j)  )/(1-d);
			max_aver /= (i+1);
            min_aver = S*(1-pow(d,j+1))  /(1-d) + S*pow(d,j  )*u*(1-pow(u,i-j))/(1-u);
			min_aver /= (i+1);
            double* stock_bucket = new double[k+1];
			interpolate(stock_bucket,max_aver,min_aver,k);
			for(l = 0; l <= k; l++){
                //next step is up
				Au = (stock_bucket[l]*(i+1) + S*pow(u,i+1-j)*pow(d,j)) / (i+2);
                tmp = 0;
                for(q = 0; q <= k; q++){
					if(Au < A[j][q])
                        tmp++; 
                }
                if(tmp == 0 || tmp == k+1 || A[j][tmp] == A[j][tmp+1]){
                    Cu = C[j][MIN(tmp,k)];
                }
                else{
                    x = (Au - A[j][tmp]) / (A[j][tmp-1] - A[j][tmp]);
                    Cu = x*C[j][tmp-1] + (1-x)*C[j][tmp];
                }
                
				//next step is down
				Ad = (stock_bucket[l]*(i+1) + S*pow(u,i-j)*pow(d,j+1)) / (i+2);
                tmp = 0;
                for(q = 0; q <= k; q++){
					if(Ad < A[j+1][q])
                        tmp++; 
                }
				if(tmp == 0 || tmp ==k+1 || Ad == A[j+1][tmp]){
                    Cd = C[j+1][MIN(tmp,k)];
                }
                else{
                    x = (Ad - A[j+1][tmp]) / (A[j+1][tmp-1] - A[j+1][tmp]);
				    Cd = x*C[j+1][tmp-1] + (1-x)*C[j+1][tmp];
                }
				put_bucket[l] = MAX(X-stock_bucket[l] , (p*Cu+(1-p)*Cd)/R);
            }
			
			//copy call_bucket[k] to C[j][k] 
			for(l = 0; l <= k; l++){
				C[j][l] = put_bucket[l];
            }
            //copy stock_bucket[k] to A[j][k]
            for(l = 0; l <= k; l++){
                A[j][l] = stock_bucket[l];
            }
		}
        if( i == 1){
		    C1=0;C2=0;
            for(j=0;j<=k;j++){
                C1+=C[0][j];
                C2+=C[1][j];
            }
            C1/=(k+1);
            C2/=(k+1);
        }
	}
	cout << "Put price is " << C[0][0] << endl;
	delta = (C1-C2)/(S*u - S*d);
	cout << "delta is " << delta << endl;
    system("pause");
}