#pragma once
// Description:
// This is a code for a generic rootfinder for real valued functions
// given an interval in which a root lies, the code finds the root 
// to a desired accuracy.
//
// Algorithm:
// This is a C++ implementation based upon algortihm M from J.C.P. Bus 
// & T.J. Dekker "Two Efficient Algorithms with Guaranteed Convergenece
// for finding a Zero of a Fuction" in ACM Trans. on Math. Software (TOMS)
// Vol 1., No. 4 Dec 1975 pages 330-345.
//
// Copyright:
// This implementation in C++ is Copyright (C) 2005 Alan Bain
// and may be freely used in other applications provided
// that this copyright notice is retained unaltered.
///////////////////////////////////////////////////////////////////////////
//
// C++ notes:
// The following discourse applies equally to e.g. numerical 
// integrators.
//
// Various standard approaches exist in C/C++ to supply the function
// whose root is to be found to a root finder.  None of these
// are used in this code.  Their various disadvantages are explained
// below.
//
// in ISO C, we can pass a function pointer
// double (*f)(double) to the routine.  Many useful routines
// whose zeros we wish to find also need parameters; which must
// generically be added as a void * so the function signature becomes
// double (*f)(double, void*) and internally the function must 
// cast the void* to some useful structure internally (the COMMON
// BLOCK appoach).  This is bug prone like its FORTRAN antecedent
// and ugly; despite this, since it is the only approach possible
// within the C language, it is used heavily e.g. NAg C library.
//
// C++: The most basic C++ "textbook" implementation is 
// using polymorphism i.e. virtual function & abstract base
//
// An abstract base class can be declared e.g. Solvable which
// contains a pure virtual function declaration e.g.
// virtual double penalty(double) =0; thus a derived class must
// provide an implementation which can access the class member data
// which solves the problem of function parameters.  
//
// This approach is bad on two counts; one is that it fixes the
// name of the function (not a good idea as a given object 
// might have several associated functions whose zero we wish to
// find).
// The second problem is that a virtual base in the class forces
// the derived class to have a vtable and thus the function dispatch
// must be done via (in effect) a pointer, which makes inlining
// impossible and carries a performance penalty in its own right.
//
//
// Approach B: member function pointer
// template <class T> 
// double BDSolver (double b, double c, T& obj, double (T::* f)(double) );
//
// The caller must supply must supply an object and the name of
// a class member function for that object.  This means that for every
// call to f the code must look up the pointer and make a function call;
// again the compiler cannot inline the function.
//
// Approach C: templates
// template <double f(double)> 
// double BDSolver(double b, double c)
// 
// This uses the C++ generic function "template" mechanism, thus one calls
// BDSolver<function>(b,c).  This is resolved at compile time so inlining
// is possible, but this approach only works for functions (not member
// functions) so we return to the problem of how to pass arguments to the
// function in a generic fashion.
//
// The approach taken is to use a template for a functor object; that is 
// one which provides double operator() (double). This allows inlining; but
// it appears to have restricted the name of the function (to operator()).
// 
// However the functor classes can be declared as nested classes of the main
// class, and since the 2003 ISO C++ standard allows such classes priviledged
// access to the containing class member elements it is simply necessary
// to pass the ctors of these nested classes the "this" pointer for the main class.  
// 
// Example:
// 
//
//class quadraticWithDeriv { 
//	public:
//		quadraticWithDeriv(double a,double b, double c) : a(a), b(b), c(c) {};		
//		class evalf { 
//		public: 
//			public: 
//				evalf(quadraticWithDeriv* obj) : ptr(obj) {}
//			double operator()(double x) const {
//				return (ptr->a*x*x+ptr->b*x+ptr->c);
//			}
//			private:
//			quadraticWithDeriv *ptr;
//		};
//
//		class evalfdash {
//			public: 
//				evalfdash(quadraticWithDeriv* obj) : ptr(obj) {}
//			double operator()(double x) const {
//				return (2*ptr->a*x+ptr->b);
//			}
//		private:
//			quadraticWithDeriv *ptr;
//		};
//
//		template <class T> T make () {
//			return T(this);
//		}	
//	private:
//	double a,b,c;
//};
// we can access the evalf function as a functor given an object obj of type
// quadraticWithDerivative via obj.make<quadraticWithDerivative::evalf>()
// which can be passes as an argument to BDSolve and this will (in most compilers
// result in code inlining.
//
// END Explanation


// some compilers ignore #pragma once....
#ifndef BUS_DEKKER_SOLVER_H
#define BUS_DEKKER_SOLVER_H



#include <cmath>

namespace NumericalLib {

enum StepType { 
	BISECT=0,
	LINEAR1,
	LINEAR2,
	EXTRAPOLATE_RATIONAL,
	NUMBER_STEP_TYPES};
// Swap two values over
template <class T> void swap(T& a, T&b) {
	T c = a;
	a=b;
	b=c;
}

template <class T> int sgn(T x) {
	return (x>0? 1 : (x<0 ? -1 : 0));
}

static const double DEFAULT_TOL=1e-10;

// Functor class providing a tolerance function which is absolute for
// small values of x, switching to relative for larger values of x
class relativeTol {
public:
		
	relativeTol(double tol=DEFAULT_TOL) : TOL(tol) { }

// simple rel/abs tolerance function
	double operator()(double x) const {
		return TOL+TOL*fabs(x);
	}

private:
	double TOL;
};

// BDSolve(b,c, function, tolerance)
//
// API:
// [b,c] should be an interval bracketing the root; i.e. f(b) 
// should have the opposite sign from f(c). 
//
// function is a functor object (i.e. one which implement 
// operator() (for speed reasons this is not enforced by
// a virtual base class).
//
// tolerance is a functor object specifying how closely 
// the root interval must be before we conside the routine
//to have converged.
//
// The routine is guaranteed to converge in at most 4 x (number of
// bisection steps) and since bisection will always terminate in 
// a boundeed number of steps, this routine guarantees to terminate
// in a bounded number of steps.  It has an asymptotic order of
// convergence close to the root of 1.618...
//
// Internal description:
//
// variables:
// b aims to be best-approximation to root so far
// c is the most recently computed point such that [b,c]
// is an interval containing a sign change (and hence
// bracketing the root).  bm1 and bm2 are previous
// values of b which are used in extrapolation.
//
// stepType is the next type of step to be taken
// (bisection, linear interpolation/extrapolation
// or rational extrapolation).
//
// At each step we force the move to be at least tol; since small
// moves would cause numerical instability for the extrapolation
// routines.  Two cases may occur:
//
// a) we do better than bisection (i.e shorten interval by 
//    more than 0.5).  In this case we replace the endpoint c
//    and interpolate.
//
// b) we do worse than bisection.  In this case endpoint c 
//    is fixed and we update b, keeping a record of the previous
//    value.  We try up to two steps of linear extrapolation
//    by which point we have 3 values bm2, bm1 and b which we 
//    use in rational extrapolation.  After this we force a 
//    bisection step; so at worst one step in four is a bisect.


template <class T, class S>
double BDsolve(double b, double c, T function, S tolerance) {
	StepType stepType=LINEAR1;
	// [b,c] is a sign change interval which brackets
	// the root.  We keep b as the value which has the
	// lowest absolute residual.
	double bm2,fbm2;

	double fb=function(b);
	double fc=function(c);

	// bm1 is previous value of b; this is used in extrapolation
	// as the second point (along with b); on entry we want to
	// interpolate
	double bm1=c;
	double fbm1=fc;
	
	double tol,midPoint,deltaB;

	// Interpolation is just a special case of extrapolation
	// where the initial endpoints are used as the extrapolands; 
	// hence interpolation will be described as LINEAR1;
	//
	while(midPoint=0.5*(b+c), 
		deltaB=midPoint-b,
		tol=tolerance(b),
		fabs(deltaB)>tol) {

				if (stepType!=BISECT) {
				// Extrapolate routines form a sequence of values of the best
				// endpoint which aim to converge while keeping the opposite sign
				// endpoint fixed.  

				if (fabs(fc)<fabs(fb)) {
					// We are not close enough to the root and
					// have encountered a minimum in the residual;
					// we mustn't extrapolate as this would take us away from the root
					// so we interpolate using the interval endpoints.

					swap(b,c); swap(fb,fc);
					deltaB=midPoint-b;
					
					// This only matters in RATIONAL_EXTRAPOLATE
					if (bm1!=b) { 
						// we need to keep a distinct values of b, bm1 and bm2
						bm2=bm1; fbm2=fbm1;
					}

					// This matters in LINEAR2 or RATIONAL_EXTRAPOLATE
					bm1=c; fbm1=fc;
				}

				

		
				tol=tol*sgn(deltaB); // so b+tol is inside interval [b,c]

				double p=(b-bm1)*fb;
				double q;

				if (stepType <=LINEAR2) {
					// Not enough points to do a 
					// rational extrapolation
					q=fbm1-fb;
					// so b+p/q= b +(b-bm1)*fb/(fb-fbm1)
				} else {
					// rational extrapolation from 
					// pts b, bm1 and bm2
					double fbm2b=(fbm2-fb)/(bm2-b);
					double fbm2bm1=(fbm2-fbm1)/(bm2-bm1);
					p=fbm2bm1*p; 
					q=fbm2b*fbm1-fbm2bm1*fb;
				}

				// Decide if we want to move; only sign of p/
				// matters so make p>0 to make tests simpler
				if (p<0.0) { 
					p=-p;
					q=-q;
				}

				// default here is bisection value
				// change if we prefer extrapolate/interpolated
				// value
				if (p==0.0 || p <= q*tol) { 
					// For stability we must make sure that
					// the point moves by at least tol each time
					deltaB=tol;
				} else if (p< q *deltaB) {
					// must land in half closer to best value
					deltaB=p/q;
				} else {
					// reject and use bisection , but EXTRAPOLATE NEXT time
					stepType=BISECT;
				}
			} // end switch


			// Shift old b values down stack
			bm2=bm1;
			fbm2=fbm1;

			bm1=b;
			fbm1=fb;

			b=b+deltaB;
			fb=function(b); // only place where called

			// Decide what solution method to use for next step

			if ( (fc>=0)? fb>=0: fb<=0)  {
				// better than bisection; we have reduced the
				// interval by more than 0.5 and we replace c
				// by old b endpoint
					c=bm1; fc=fbm1;
					// don't advance state
					if (stepType==BISECT) stepType=LINEAR1;
			}  else { 
					// worse than bisect so advance state
					stepType=static_cast<StepType> ((stepType+1) % NUMBER_STEP_TYPES);
			}
	
		} while (fabs(midPoint-b)>tol);

		return (0.5*(b+c));
	
	}
	




	


}
#endif