/*
    This file is part of Mitsuba, a physically based rendering system.

    Copyright (c) 2007-2011 by Wenzel Jakob and others.

    Mitsuba is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License Version 3
    as published by the Free Software Foundation.

    Mitsuba is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

#if !defined(__BRENT_H)
#define __BRENT_H

#include <mitsuba/mitsuba.h>
#include <boost/function.hpp>

MTS_NAMESPACE_BEGIN

/**
 * \brief Brent's method nonlinear zero finder
 *
 * The implementation is transcribed from the Apache Commons 
 * Java implementation. The supplied function is required to be
 * continuous, but not necessarily smooth.
 *
 * \ingroup libcore
 */
class MTS_EXPORT_CORE BrentSolver {
public:
	/// Return value of \ref BrentSolver::solve()
	struct Result {
		bool success;
		size_t iterations;
		Float x;
		Float y;

		/// Create a new result instance
		inline Result(bool success, size_t iterations, Float x, Float y)
			: success(success), iterations(iterations), x(x), y(y) { }

		/// Return a string representation of the result
		inline std::string toString() const {
			std::ostringstream oss;
			oss << "BrentSolver::Result["
				<< "success=" << success << ", "
				<< "iterations=" << iterations << ", "
				<< "x=" << x << ", "
				<< "y=" << y << "]";
			return oss.str();
		}
	};

	/**
	 * \brief Create a new Brent-style solver with the 
	 * specified accuracy requirements
	 *
	 * \param maxIterations
	 *      Max. number of successive iterations (default: 100)
	 * \param absAccuracy
	 *      Absolute accuracy requirement -- the iterations will stop
	 *      when |f(x)| < absAccuracy.
	 * \param absAccuracyPos
	 *      Absolute accuracy requirement of the position -- the 
	 *      iterations will stop when |minX-maxX| < absAccuracyPos.
	 * \param absAccuracyPos
	 *      Absolute accuracy requirement of the position -- the
	 *      iterations will stop when |minX-maxX|/minX < relAccuracyPos.
	 */
	inline BrentSolver(size_t maxIterations = 100, 
			Float absAccuracy = 1e-6f,
			Float absAccuracyPos = 1e-6f,
			Float relAccuracyPos = 1e-6f) 
		: m_maxIterations(maxIterations),
		  m_absAccuracy(absAccuracy),
		  m_absAccuracyPos(absAccuracyPos),
		  m_relAccuracyPos(relAccuracyPos) { }

	/**
	 * \brief Find a zero in the given interval.
	 *
	 * Requires that the values of the function at the endpoints 
	 * have opposite signs. 
	 *
	 * \param min the lower bound for the interval.
	 * \param max the upper bound for the interval.
	 * \return the value where the function is zero
	 */
	Result solve(const boost::function<Float (Float)> &func,
			Float min, Float max) const;

	/**
	 * \brief Find a zero in the given interval with an initial guess
	 *
	 * Requires that the values of the function at the endpoints 
	 * have opposite signs (note that it is allowed to have endpoints 
	 * with the same sign if the initial point has opposite sign 
	 * function-wise).
	 *
	 * \param min the lower bound for the interval.
	 * \param max the upper bound for the interval.
	 * \param initial the start value to use (must be set to min 
	 *    if no initial point is known)
	 * \return the value where the function is zero
	 */
	Result solve(const boost::function<Float (Float)> &func,
			Float min, Float max, Float initial) const;

	/**
	 * Find a zero starting search according to the three provided points.
	 *
	 * \param x0 old approximation for the root
	 * \param y0 function value at the approximation for the root
	 * \param x1 last calculated approximation for the root
	 * \param y1 function value at the last calculated approximation
	 *    for the root
	 * \param x2 bracket point (must be set to x0 if no bracket point is
	 *   known, this will force starting with linear interpolation)
	 * \param y2 function value at the bracket point.
	 * \return the value where the function is zero
	 */
	Result solve(const boost::function<Float (Float)> &func,
				 Float x0, Float y0,
				 Float x1, Float y1,
				 Float x2, Float y2) const;
protected:
	size_t m_maxIterations;
	Float m_absAccuracy;
	Float m_absAccuracyPos;
	Float m_relAccuracyPos;
}; 

MTS_NAMESPACE_END

#endif /* __BRENT_H */
