//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
template <class Function>
G4bool G4Solver<Function>::Bisection(Function & theFunction)
{
    // Check the interval before start
    if (a > b || std::abs(a-b) <= tolerance) 
    {
	G4cerr << "G4Solver::Bisection: The interval must be properly set." << G4endl;
	return false;
    }
    G4double fa = theFunction(a);
    G4double fb = theFunction(b);
    if (fa*fb > 0.0) 
    {
	G4cerr << "G4Solver::Bisection: The interval must include a root." << G4endl;
	return false;
    }
    
    G4double eps=tolerance*(b-a);
    
    
    // Finding the root
    for (G4int i = 0; i < MaxIter; i++) 
    {
	G4double c = (a+b)/2.0;
	if ((b-a) < eps) 
	{
	    root = c;
	    return true;
	}
	G4double fc = theFunction(c);
	if (fc == 0.0) 
	{
	    root = c;
	    return true;
	}
	if (fa*fc < 0.0) 
	{
	    a=c;
	    fa=fc;
	} 
	else 
	{
	    b=c;
	    fb=fc;
	}
    }
    G4cerr << "G4Solver::Bisection: Excedded maximum number of iterations whithout convegence." << G4endl;
    return false;
}


template <class Function>
G4bool G4Solver<Function>::RegulaFalsi(Function & theFunction)
{
    // Check the interval before start
    if (a > b || std::abs(a-b) <= tolerance) 
    {
	G4cerr << "G4Solver::RegulaFalsi: The interval must be properly set." << G4endl;
	return false;
    }
    G4double fa = theFunction(a);
    G4double fb = theFunction(b);
    if (fa*fb > 0.0) 
    {
	G4cerr << "G4Solver::RegulaFalsi: The interval must include a root." << G4endl;
	return false;
    }
    
    G4double eps=tolerance*(b-a);
	
	
    // Finding the root
    for (G4int i = 0; i < MaxIter; i++) 
    {
	G4double c = (a*fb-b*fa)/(fb-fa);
	G4double delta = std::min(std::abs(c-a),std::abs(b-c));
	if (delta < eps) 
	{
	    root = c;
	    return true;
	}
	G4double fc = theFunction(c);
	if (fc == 0.0) 
	{
	    root = c;
	    return true;
	}
	if (fa*fc < 0.0) 
	{
	    b=c;
	    fb=fc;
	} 
	else 
	{
	    a=c;
	    fa=fc;
	}
    }
    G4cerr << "G4Solver::Bisection: Excedded maximum number of iterations whithout convegence." << G4endl;
    return false;

}	

template <class Function>
G4bool G4Solver<Function>::Brent(Function & theFunction)
{
    
    const G4double precision = 3.0e-8;
    
    // Check the interval before start
    if (a > b || std::abs(a-b) <= tolerance) 
    {
	G4cerr << "G4Solver::Brent: The interval must be properly set." << G4endl;
	return false;
    }
    G4double fa = theFunction(a);
    G4double fb = theFunction(b);
    if (fa*fb > 0.0) 
    {
	G4cerr << "G4Solver::Brent: The interval must include a root." << G4endl;
	return false;
    }
    
    G4double c = b;
    G4double fc = fb;
    G4double d = 0.0;
    G4double e = 0.0;
    
    for (G4int i=0; i < MaxIter; i++) 
    {
	// Rename a,b,c and adjust bounding interval d
	if (fb*fc > 0.0) 
	{
	    c = a;
	    fc = fa;
	    d = b - a;
	    e = d;
	}
	if (std::abs(fc) < std::abs(fb)) 
	{
	    a = b;
	    b = c;
	    c = a;
	    fa = fb;
	    fb = fc;
	    fc = fa;
	}
	G4double Tol1 = 2.0*precision*std::abs(b) + 0.5*tolerance;
	G4double xm = 0.5*(c-b);
	if (std::abs(xm) <= Tol1 || fb == 0.0) 
	{
	    root = b;
	    return true;
	}
	// Inverse quadratic interpolation
	if (std::abs(e) >= Tol1 && std::abs(fa) > std::abs(fb)) 
	{
	    G4double s = fb/fa;
	    G4double p = 0.0;
	    G4double q = 0.0;
	    if (a == c) 
	    {
		p = 2.0*xm*s;
		q = 1.0 - s;
	    } 
	    else 
	    {
		q = fa/fc;
		G4double r = fb/fc;
		p = s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
		q = (q-1.0)*(r-1.0)*(s-1.0);
	    }
	    // Check bounds
	    if (p > 0.0) q = -q;
	    p = std::abs(p);
	    G4double min1 = 3.0*xm*q-std::abs(Tol1*q);
	    G4double min2 = std::abs(e*q);
	    if (2.0*p < std::min(min1,min2)) 
	    {
		// Interpolation
		e = d;
		d = p/q;
	    } 
	    else 
	    {
		// Bisection
		d = xm;
		e = d;
	    }
	} 
	else 
	{
	    // Bounds decreasing too slowly, use bisection
	    d = xm;
	    e = d;
	}
	// Move last guess to a 
	a = b;
	fa = fb;
	if (std::abs(d) > Tol1) b += d;
	else 
	{
	    if (xm >= 0.0) b += std::abs(Tol1);
	    else b -= std::abs(Tol1);
	}
	fb = theFunction(b);
    }
    G4cerr << "G4Solver::Brent: Number of iterations exceeded." << G4endl;
    return false;
}



template <class Function>
G4bool G4Solver<Function>::Crenshaw(Function & theFunction)
{
    // Check the interval before start
    if (a > b || std::abs(a-b) <= tolerance) 
    {
	G4cerr << "G4Solver::Crenshaw: The interval must be properly set." << G4endl;
	return false;
    }

    G4double fa = theFunction(a);
    if (fa == 0.0) 
    {
	root = a;
	return true;
    }

    G4double Mlast = a;

    G4double fb = theFunction(b);
    if (fb == 0.0)
    {
	root = b;
	return true;
    }

    if (fa*fb > 0.0) 
    {
	G4cerr << "G4Solver::Crenshaw: The interval must include a root." << G4endl;
	return false;
    }

    
    for (G4int i=0; i < MaxIter; i++) 
      {
	G4double c = 0.5 * (b + a);
	G4double fc = theFunction(c);
	if (fc == 0.0 || std::abs(c - a) < tolerance) 
	  {
	    root = c;
	    return true;
	  }

	if (fc * fa  > 0.0)
	  {
	    G4double tmp = a;
	    a = b;
	    b = tmp;
	    tmp = fa;
	    fa = fb;
	    fb = tmp;
	  }
	
	G4double fc0 = fc - fa;
	G4double fb1 = fb - fc;
	G4double fb0 = fb - fa;
	if (fb * fb0 < 2.0 * fc * fc0)
	  {
	    b = c;
	    fb = fc;
	  }
	else 
	  {
	    G4double B = (c - a) / fc0;
	    G4double C = (fc0 - fb1) / (fb1 * fb0);
	    G4double M = a - B * fa * (1.0 - C * fc);
	    G4double fM = theFunction(M);
	    if (fM == 0.0 || std::abs(M - Mlast) < tolerance)
	      {
		root = M;
		return true;
	      }
	    Mlast = M;
	    if (fM * fa < 0.0)
	      {
		b = M;
		fb = fM;
	      }
	    else 
	    {
		a = M;
		fa = fM;
		b = c;
		fb = fc;
	    }
	}
    }
    return false;
}




template <class Function> 	
void G4Solver<Function>::SetIntervalLimits(const G4double Limit1, const G4double Limit2)
{
    if (std::abs(Limit1-Limit2) <= tolerance) 
    {
	G4cerr << "G4Solver::SetIntervalLimits: Interval must be wider than tolerance." << G4endl;
	return;
    }
    if (Limit1 < Limit2) 
    {
	a = Limit1;
	b = Limit2;
    } 
    else 
    {
	a = Limit2;
	b = Limit1;
    }
    return;
}