/*
								+----------------------------------+
								|                                  |
								| ***  Numerical math methods  *** |
								|                                  |
								|   Copyright © -tHE SWINe- 2008   |
								|                                  |
								|           Numerical.h            |
								|                                  |
								+----------------------------------+
*/

#pragma once
#ifndef __NUMERICAL_MATH_INCLUDED
#define __NUMERICAL_MATH_INCLUDED

/**
 *	@file Numerical.h
 *	@date 2008
 *	@author -tHE SWINe-
 *	@brief numerical math methods
 */

/**
 *	@brief simple one-dimensional integrator with max error parameter
 *
 *	@tparam CFunctor is a function object type
 *	@tparam CScalarType is scalar type
 */
template <class CFunctor, class CScalarType>
class CIntegrator {
protected:
	CFunctor m_function; /**< @brief function to be solved */
	CScalarType m_f_min; /**< @brief lower bound of the integration interval */
	CScalarType m_f_scale; /**< @brief size of the integration interval */
	CScalarType m_f_max_error; /**< @brief maximum permitted error */

public:
	/**
	 *	@brief default constructor
	 *
	 *	@param[in] function is the function being integrated
	 *	@param[in] f_min is lower bound of the integration interval
	 *	@param[in] f_max is lower bound of the integration interval
	 *	@param[in] f_max_error is maximal relative error of the solution
	 */
	inline CIntegrator(CFunctor function, CScalarType f_min,
		CScalarType f_max, CScalarType f_max_error)
		:m_function(function), m_f_min(f_min), m_f_scale(f_max - f_min),
		m_f_max_error(f_max_error)
	{}

	/**
	 *	@brief calculates the integral
	 *	@return Returns the value of the integral.
	 *	@note This may take some time if high precision is required
	 *		and / or the integration range is wide.
	 */
	inline operator CScalarType () const
	{
		CScalarType f_scale = m_f_scale;
		CScalarType f_first_integral((m_function(m_f_min) +
			m_function(m_f_min + f_scale)) * f_scale * CScalarType(.5));
		// first integral approximation (trapezoid rule)

		if(m_f_scale == 0)
			return 0;

		for(int n_pass = 1; n_pass < 20; ++ n_pass) {
			CScalarType f_integral(f_first_integral * CScalarType(.5));
			// prepare integration variable

			f_scale *= CScalarType(.5);
			// scale is going down with every refinement

			CScalarType f_refinement(0), f_arg(f_scale + m_f_min);
			for(uint32_t i = 1, n = (uint32_t(1) << n_pass) + 1; i < n; i += 2, f_arg += 2 * f_scale)
				f_refinement += m_function(f_arg); // bloody hell optimized
			f_refinement *= f_scale;
			f_integral += f_refinement;
			// refine integral sollution

			if(fabs((f_integral - f_first_integral) / f_first_integral) <= m_f_max_error)
				return f_integral;
			// is error in required range now?

			f_first_integral = f_integral;
		}
		// adaptive subdivision to fit the error

		return f_first_integral;
	}
};

/**
 *	@brief integrates one-dimensional function
 *
 *	@tparam CFunctor is a function object type
 *	@tparam CScalarType is scalar type
 *
 *	@param[in] function is the function being integrated
 *	@param[in] f_min is lower bound of the integration interval
 *	@param[in] f_max is lower bound of the integration interval
 *	@param[in] f_max_error is maximal relative error of the solution
 *
 *	@return Returns the value of the integral.
 */
template <class CFunctor, class CScalarType>
inline double f_Integrate(CFunctor function, CScalarType f_min, CScalarType f_max, CScalarType f_max_error)
{
	return CIntegrator<CFunctor, CScalarType>(function, f_min, f_max, f_max_error);
}

#if 0 // deprecated, utterly useless

/**
 *	@brief simple one-dimensional numerical sovler with max error parameter
 *
 *	@tparam CFunctor is a function object type
 *	@tparam CScalarType is scalar type
 *
 *	@note Using the golden cut method.
 */
template <class CFunctor, class CScalarType>
class CGoldenCutSover {
protected:
	CFunctor m_function; /**< @brief function to be solved */
	CScalarType m_f_min; /**< @brief lower bound of argument values */
	CScalarType m_f_max; /**< @brief upper bound of argument values */
	CScalarType m_f_max_error; /**< @brief maximal error of the solution */
	CScalarType m_f_root; /**< @brief target value of the function */

public:
	/*
	 *	inline CGoldenCutSover::CGoldenCutSover(CFunctor function, CScalarType f_root,
	 *		CScalarType f_min, CScalarType f_max, CScalarType f_max_error)
	 *		- default constructor; function is root lies in range f_min to f_max (estimate),
	 *		  calculates root with maximal relative error f_max_error
	 */
	inline CGoldenCutSover(CFunctor function, CScalarType f_root,
		CScalarType f_min, CScalarType f_max, CScalarType f_max_error)
		:m_function(function), m_f_min(f_min), m_f_max(f_max),
		m_f_root(f_root), m_f_max_error(f_max_error * (f_max - f_min))
	{
		_ASSERTE(m_f_min < m_f_max);
	}

	/**
	 *	@brief finds root position
	 *	@note This might get stuck in an infinite loop if unlucky.
	 */
	inline operator CScalarType() const
	{
		CScalarType f_left_x = m_f_min;
		CScalarType f_right_x = m_f_max;
		// start in min / max

		CScalarType f_left_y = m_function(f_left_x) - m_f_root;
		CScalarType f_right_y = m_function(f_right_x) - m_f_root;
		// calculate function in both points

		const CScalarType f_theta = CScalarType(.6180339889579);
		// golden ratio

		do {
			if(fabs(f_left_y) < fabs(f_right_y)) {
				f_right_x = f_left_x + (f_right_x - f_left_x) * f_theta;
				f_right_y = m_function(f_right_x) - m_f_root;
			} else {
				f_left_x = f_left_x + (f_right_x - f_left_x) * (1 - f_theta);
				f_left_y = m_function(f_left_x) - m_f_root;
			}
			// shift one of points closer to root and recalculate function values
		} while(f_right_x - f_left_x > m_f_max_error);

		return (f_left_x + f_right_x) * .5;
	}
};

/*
 *	template <class CFunctor, class CScalarType>
 *	class CGoldenCutSover2
 *		- simple one-dimensional numerical sovler with max error parameter
 *		- difference between CGoldenCutSover and CGoldenCutSover2
 *		  is CGoldenCutSover2 finds intersection of function with zero
 *		  whereas CGoldenCutSover finds intersection with given constant
 *		- using the golden cut method
 */
template <class CFunctor, class CScalarType>
class CGoldenCutSover2 {
protected:
	CFunctor m_function;
	CScalarType m_f_min, m_f_max;
	CScalarType m_f_max_error;
	CScalarType m_f_root;

public:
	/*
	 *	inline CGoldenCutSover::CGoldenCutSover2(CFunctor function, CScalarType f_min,
	 *		CScalarType f_max, CScalarType f_max_error)
	 *		- default constructor; function is root lies in range f_min to f_max (estimate),
	 *		  calculates root with maximal relative error f_max_error
	 */
	inline CGoldenCutSover2(CFunctor function, CScalarType f_min,
		CScalarType f_max, CScalarType f_max_error)
		:m_function(function), m_f_min(f_min), m_f_max(f_max),
		m_f_max_error(f_max_error * (f_max - f_min))
	{
		_ASSERTE(m_f_min < m_f_max);
	}

	/*
	 *	inline CGoldenCutSover::operator CScalarType () const
	 *		- evaluates root position
	 *		- note this might get stuck in infinite loop if unlucky
	 */
	inline operator CScalarType() const
	{
		CScalarType f_left_x = m_f_min;
		CScalarType f_right_x = m_f_max;
		// start in min / max

		CScalarType f_left_y = m_function(f_left_x);
		CScalarType f_right_y = m_function(f_right_x);
		// calculate function in both points

		do {
			if(fabs(f_left_y) < fabs(f_right_y)) {
				f_right_x = f_left_x + (f_right_x - f_left_x) * f_theta;
				f_right_y = m_function(f_right_x);
			} else {
				f_left_x = f_left_x + (f_right_x - f_left_x) * (1 - f_theta);
				f_left_y = m_function(f_left_x);
			}
			// shift one of points closer to root and recalculate function values
		} while(f_right_x - f_left_x > m_f_max_error);
		return (f_left_x + f_right_x) / 2;
	}
};

/*
 *	template <class CFunctor>
 *	inline float f_Solve(CFunctor function, float f_root, float f_min, float f_max, float f_max_error)
 *		- finds function root (such x where function(x) == f_root) using golden cut method
 *		- f_min to f_max is estimate range where root should be
 *		- returns root with maximal relative error f_max_error
 *		  (or ends in infinite loop if it can't find one)
 */
template <class CFunctor>
inline float f_Solve(CFunctor function, float f_root, float f_min, float f_max, float f_max_error)
{
	return CGoldenCutSover<CFunctor, float>(function, f_root, f_min, f_max, f_max_error);
}

/*
 *	template <class CFunctor>
 *	inline float f_Solve(CFunctor function, float f_min, float f_max, float f_max_error)
 *		- finds function root (such x where function(x) == 0) using golden cut method
 *		- f_min to f_max is estimate range where root should be
 *		- returns root with maximal relative error f_max_error
 *		  (or ends in infinite loop if it can't find one)
 */
template <class CFunctor>
inline float f_Solve(CFunctor function, float f_min, float f_max, float f_max_error)
{
	return CGoldenCutSover2<CFunctor, float>(function, f_min, f_max, f_max_error);
}

#endif // 0

#endif // !__NUMERICAL_MATH_INCLUDED
