// my_function.h
#include <iostream>
#include <vector>

const double pi = 3.1415926;
const double sqrt_2_pi = sqrt(2 * pi);

const double a0 = 2.50662823884;
const double a1 = -18.61500062529;
const double a2 = 41.39119773534;
const double a3 = -25.44106049637;
const double b1 = -8.47351093090;
const double b2 = 23.08336743743;
const double b3 = -21.06224101826;
const double b4 = 3.13082909833;
const double c0 = 0.3374754822726147;
const double c1 = 0.9761690190917186;
const double c2 = 0.1607979714918209;
const double c3 = 0.0276438810333863;
const double c4 = 0.0038405729373609;
const double c5 = 0.0003951896511919;
const double c6 = 0.0000321767881768;
const double c7 = 0.0000002888167364;
const double c8 = 0.0000003960315187;

int recursiveAdd(int n);

void recursivePrint(int a, int b);

int fibonacciCal(int n);

double normcdf(double x);

double my_normcdf(double x);

double hornerFunction(double x, double a0, double a1);

double hornerFunction(double x, double a0, double a1, double a2);

double hornerFunction(double x, double a0, double a1, double a2, double a3);

double hornerFunction(double x, double a0, double a1, double a2, double a3, double a4);

double hornerFunction(double x, double a0, double a1, double a2, double a3, double a4,
    double a5);

double hornerFunction(double x, double a0, double a1, double a2, double a3, double a4,
    double a5, double a6);

double hornerFunction(double x, double a0, double a1, double a2, double a3, double a4,
    double a5, double a6, double a7);

double hornerFunction(double x, double a0, double a1, double a2, double a3, double a4,
    double a5, double a6, double a7, double a8);

double my_hornerFunction(const double& x, std::vector<double>& a);

double norminv(double x);

double my_norminv(double x);

double my_blackScholesCallPrice(double strike, double timeToMaturity,
    double spot, double volatility, double riskFreeRate);

double blackScholesCallPrice(double strike, double timeToMaturity,
    double spot, double volatility, double riskFreeRate);

my_function.cpp

#include "my_function.h"
#include <cmath>



int recursiveAdd(int n) {
	/*Write a recursive function to compute the sum of the numbers between 1 and n*/
	if (n == 1) {
		return 1;
	}
	else {
		return n + recursiveAdd(n - 1);
	}
}

void recursivePrint(int a, int b) {
	if (a == b) {
		std::cout << a << std::endl;
	}
	else if (a < b){
		std::cout << a << " ";
		recursivePrint(++a, b);
		
	}
	else {
		std::cout << " the input maybe wrong";
	}
}

int fibonacciCal(int n) {
	if (n == 0) {
		return 1;
	}
	else if (n == 1) {
		return 1;
	}
	else if (n < 0) {
		std::cout << "warnings,input n is wrong, please input n>=0";
		return NULL;
	}
	else {
		return fibonacciCal(n - 1) + fibonacciCal(n - 2);
	}
}

double normcdf(double x) {
	return 0.5 * erfc(-x * sqrt(0.5));
}

double my_normcdf(double x) {
	/* 计算标准正太分布的累计概率分布
	*/
	if (x >= 0) {
		x = abs(x);
		double pi = 3.1415926;
		double k = 1 / (1 + 0.2316419 * x);
		double v5 = k * (0.319381530 + k * (-0.356563782 + k * (1.781477937 + k * (-1.821255978 + 1.330274429 * k))));
		double temp = 1.0 / sqrt_2_pi * exp(-0.5 * x * x) * v5;
		return 1 - temp;
	}
	else {
		x = abs(x);
		double pi = 3.1415926;
		double k = 1 / (1 + 0.2316419 * x);
		double v5 = k * (0.319381530 + k * (-0.356563782 + k * (1.781477937 + k * (-1.821255978 + 1.330274429 * k))));
		double temp = 1.0 / sqrt_2_pi * exp(-0.5 * x * x) * v5;
		return temp;
	}
	

}

double my_hornerFunction(const double & x, std::vector<double> &a) {
	/*std::cout << "a value is " << std::endl;
	for (int size_t = 0; size_t < a.size(); size_t++) {
		std::cout << "  " << a[size_t];
	}*/
	std::vector<double> a_copy(a);
	/*std::cout << "a_copy value is " << std::endl;
	for (int size_t = 0; size_t < a_copy.size(); size_t++) {
		std::cout << "  " << a_copy[size_t];
	}*/
	//std::cout << std::endl;
	if (a_copy.size() == 2) {
		return a_copy[0] + x * a_copy[1];
	}
	else if (a_copy.size() <= 1) {
		std::cout << "warnings,input vector or array is wrong" << std::endl;
		return NULL;
	}
	else {
		double a_value = a_copy[0];
		a_copy.erase(a_copy.begin(), a_copy.begin() + 1);
		return a_value + x * my_hornerFunction(x, a_copy);
	}
}


double hornerFunction(double x, double a0, double a1) {
	//std::cout << "01 normal value is " << a0 + x * a1 << std::endl;
	return a0 + x * a1;
}

double hornerFunction(double x, double a0, double a1, double a2) {
	//std::cout << "02 a0 = " << a0 << " x = " << x << " result = " << result << std::endl;
	return a0 + x * hornerFunction(x, a1, a2);
}

double hornerFunction(double x, double a0, double a1, double a2, double a3) {
	return a0 + x * hornerFunction(x, a1, a2, a3);
}

double hornerFunction(double x, double a0, double a1, double a2, double a3, double a4) {
	return a0 + x * hornerFunction(x, a1, a2, a3, a4);
}

double hornerFunction(double x, double a0, double a1, double a2, double a3, double a4,
	double a5) {
	return a0 + x * hornerFunction(x, a1, a2, a3, a4, a5);
}

double hornerFunction(double x, double a0, double a1, double a2, double a3, double a4,
	double a5, double a6) {
	return a0 + x * hornerFunction(x, a1, a2, a3, a4, a5, a6);
}

double hornerFunction(double x, double a0, double a1, double a2, double a3, double a4,
	double a5, double a6, double a7) {
	return a0 + x * hornerFunction(x, a1, a2, a3, a4, a5, a6, a7);
}

double hornerFunction(double x, double a0, double a1, double a2, double a3, double a4,
	double a5, double a6, double a7, double a8) {
	return a0 + x * hornerFunction(x, a1, a2, a3, a4, a5, a6, a7, a8);
}

double norminv(double x) {
	// We use Moro's algorithm
	double y = x - 0.5;
	if (y<0.42 && y>-0.42) {
		double r = y * y;
		return y * hornerFunction(r, a0, a1, a2, a3) / hornerFunction(r, 1.0, b1, b2, b3, b4);
	}
	else {
		double r;
		if (y < 0.0) {
			r = x;
		}
		else {
			r = 1.0 - x;
		}
		double s = log(-log(r));
		double t = hornerFunction(s, c0, c1, c2, c3, c4, c5, c6, c7, c8);
		if (x > 0.5) {
			return t;
		}
		else {
			return -t;
		}
	}
}

double my_norminv(double x) {
	// We use Moro's algorithm
	double y = x - 0.5;
	if (y<0.42 && y>-0.42) {
		double r = y * y;
		std::vector<double> a = { a0, a1, a2, a3 };
		std::vector<double> b = { 1.0, b1, b2, b3, b4 };
		return y * my_hornerFunction(r, a) / my_hornerFunction(r, b);
	}
	else {
		double r;
		if (y < 0.0) {
			r = x;
		}
		else {
			r = 1.0 - x;
		}
		double s = log(-log(r));
		std::vector<double> c = { c0, c1, c2, c3, c4, c5, c6, c7, c8 };
		double t = my_hornerFunction(s, c);
		if (x > 0.5) {
			return t;
		}
		else {
			return -t;
		}
	}
}

double my_blackScholesCallPrice(double strike, double timeToMaturity,
	double spot, double volatility, double riskFreeRate) {
	double numetator = log(spot / strike) + (riskFreeRate + volatility * volatility / 2) * timeToMaturity;
	double e_t = volatility * sqrt(timeToMaturity);
	double d1 = numetator / e_t;
	double d2 = d1 - e_t;
	double nd1 = my_normcdf(d1);
	//std::cout << "d1 = " << d1 << " my_normcdf(d1) = " << my_normcdf(d1) << std::endl;
	double nd2 = my_normcdf(d2);
	double c0 = spot * nd1 - strike * nd2 * exp(-riskFreeRate * timeToMaturity);
	return c0;
}

double blackScholesCallPrice(double strike, double timeToMaturity,
	double spot, double volatility, double riskFreeRate) {
	double numerator = log(spot / strike)
		+ (riskFreeRate + volatility * volatility * 0.5) * timeToMaturity;
	double denominator = volatility * sqrt(timeToMaturity);
	double d1 = numerator / denominator;
	double d2 = d1 - denominator;
	//std::cout << "d1 = " << d1 << " normcdf(d1) = " << normcdf(d1) << std::endl;
	double t1 = normcdf(d1) * spot;
	double t2 = normcdf(d2) * strike * exp(-riskFreeRate * timeToMaturity);
	return t1 - t2;
}