#ifndef GELIB_IMPL_ALGORITHM_H_
#define GELIB_IMPL_ALGORITHM_H_

#define _USE_MATH_DEFINES

#include <cmath>
#include <numeric>
#include <algorithm>
#include <vector>
#include <sstream>
#include <functional>
#include "../1_dimensional/Function1d.h"

namespace gelib_impl
{

	namespace Algorithm
	{

		const double inv_sqrt_2pi = 0.3989422804014327;

		const double BRENT_RTOL = 1e-15;

		/*
		 * Maximum allowed number of iterations for brent algorithm
		 */
		const int BRENT_ITMAX = 100;

		/*
		 * Calculate the probability z < x  for z ~ N(0, 1)
		 */
		double normal_CDF(double x);

		/*
		* Evaluate probability density of standard normal distribution at x
		*/
		double normal_PDF(double x);

		/*
		* Evaluate probability density of normal
		* distribution x ~ N(mu, sigma^2) at x
		*/
		double normal_PDF(double x, double mu, double sigma);

		/*
		* Evaluate probability density of normal
		* distribution x ~ N(mu, sigma^2 = variance) at x
		*/
		double normal_PDF(double x, double mu, double sigma, double variance);

		/*
		 * calculate the sum of n elements in an array
		 * i.e: a[0] + a[1] + ... + a[n - 1]
		 */
		double array_sum(double* arr, int n);

		/*
		 * calculate the product of n elements in an array
		 * i.e: a[0] * a[1] * ... * a[n - 1]
		 */
		template<typename T>
		T array_product(T* arr, int n)
		{
			T ret = 1;
			for (int i = 0; i < n; ++i)
				ret *= arr[i];
			return ret;
		}

		/*
		* calculate the product n * (n+1) * ... * (m-1) * m
		* for m > n
		*/
		template<typename T>
		T product(T n, T m, T running = 1)
		{
			return n < m ? product(n + 1, m, n * running) : n * running;
		}

		/*
		 * calculate the binomial coefficient (n C k) n,k > 0 and integers
		 */
		int binomial_coefficient(unsigned int n, unsigned int k);

		/*
		 * calculate x^n for integer n >= 0 and real x
		 */
		double power(double x, size_t n);

		/*
		 * count the number of nan elements in vector v
		 */
		__int64 count_nan(const std::vector<double>& v);

		// TODO move copy_array to Containers
		/*
		 * Copy (double) array a in to array b
		 * where size = number of elements in and b
		 */
		void copy_array(double* arr_a, double* arr_b, size_t size);

		/*
		 * Find iterator pointing at index i of vector v, where v[i - 1] <= x < v[i]
		 */
		std::vector<double>::iterator it_in_range(double x, std::vector<double>& v);

		/*
	 * vector style element-wise operation on arr_a and arr_b, of length n_a and n_b respectively
	 * The result is stored in arr_c, which should be length n_c = max(n_a, n_b)
	 */
		void array_elem_op(double* arr_a, size_t n_a, double* arr_b, size_t n_b, double* arr_c, size_t n_c, std::function<double(double, double)> op);

		/*
		 * vector style element-wise add of arr_a and arr_b, of length n_a and n_b respectively
		 * The result is stored in arr_c, which should be length n_c = max(n_a, n_b)
		 */
		void array_elem_add(double* arr_a, size_t n_a, double* arr_b, size_t n_b, double* arr_c, size_t n_c);

		/*
		 * vector style element-wise subtract of arr_a and arr_b, of length n_a and n_b respectively
		 * The result is stored in arr_c, which should be length n_c = max(n_a, n_b)
		 */
		void array_elem_sub(double* arr_a, size_t n_a, double* arr_b, size_t n_b, double* arr_c, size_t n_c);

		/* Using Brent's method, find the root of a function f bracketed between x1 and x2.
		 * The root will be refined until its accuracy is tol.
		 * NOTE: it is assumed that sign(f(x1)) != sign(f(x2))
		 */
		double brent(Function1d& f, double x1, double x2, double f_x1 = NAN, double f_x2 = NAN, double tol = 1e-15);

		/*
		* Use Simpson's rule to calculate integral to arbitrary accuracy
		* i.e increase strips until prescribed accuracy reached
		*/
		double integrate(Function1d& f, double a, double b, double tol = 1e-9);

	}

}

#endif /* ALGORITHM_H_ */
