/*****************************************************************************
 * Project: RooFit                                                           *
 * Package: RooFitCore                                                       *
 * @(#)root/roofitcore:$Id$
 * Authors:                                                                  *
 *   WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu       *
 *   DK, David Kirkby,    UC Irvine,         dkirkby@uci.edu                 *
 *                                                                           *
 * Copyright (c) 2000-2005, Regents of the University of California          *
 *                          and Stanford University. All rights reserved.    *
 *                                                                           *
 * Redistribution and use in source and binary forms,                        *
 * with or without modification, are permitted according to the terms        *
 * listed in LICENSE (http://roofit.sourceforge.net/license.txt)             *
 *****************************************************************************/

//////////////////////////////////////////////////////////////////////////////
// 
// BEGIN_HTML
// Implement the abstract 1-dimensional root finding interface using
// the Brent-Decker method. This implementation is based on the one
// in the GNU scientific library (v0.99).
// END_HTML
//

#include "RooFit.h"

#include "RooBrentRootFinder.h"
#include "RooBrentRootFinder.h"
#include "RooAbsFunc.h"
#include <math.h>
#include "Riostream.h"
#include "RooMsgService.h"

using namespace std;

ClassImp(RooBrentRootFinder)
;


//_____________________________________________________________________________
RooBrentRootFinder::RooBrentRootFinder(const RooAbsFunc& function) :
  RooAbsRootFinder(function),
  _tol(2.2204460492503131e-16)
{
  // Constructor taking function binding as input
}



//_____________________________________________________________________________
Bool_t RooBrentRootFinder::findRoot(Double_t &result, Double_t xlo, Double_t xhi, Double_t value) const
{
  // Do the root finding using the Brent-Decker method. Returns a boolean status and
  // loads 'result' with our best guess at the root if true.
  // Prints a warning if the initial interval does not bracket a single
  // root or if the root is not found after a fixed number of iterations.

  _function->saveXVec() ;

  Double_t a(xlo),b(xhi);
  Double_t fa= (*_function)(&a) - value;
  Double_t fb= (*_function)(&b) - value;
  if(fb*fa > 0) {
    oocxcoutD((TObject*)0,Eval) << "RooBrentRootFinder::findRoot(" << _function->getName() << "): initial interval does not bracket a root: ("
				<< a << "," << b << "), value = " << value << " f[xlo] = " << fa << " f[xhi] = " << fb << endl;
    return kFALSE;
  }

  Bool_t ac_equal(kFALSE);
  Double_t fc= fb;
  Double_t c(0),d(0),e(0);
  for(Int_t iter= 0; iter <= MaxIterations; iter++) {

    if ((fb < 0 && fc < 0) || (fb > 0 && fc > 0)) {
      // Rename a,b,c and adjust bounding interval d
      ac_equal = kTRUE;
      c = a;
      fc = fa;
      d = b - a;
      e = b - a;
    }
  
    if (fabs (fc) < fabs (fb)) {
      ac_equal = kTRUE;
      a = b;
      b = c;
      c = a;
      fa = fb;
      fb = fc;
      fc = fa;
    }

    Double_t tol = 0.5 * _tol * fabs(b);
    Double_t m = 0.5 * (c - b);


    if (fb == 0 || fabs(m) <= tol) {
      //cout << "RooBrentRootFinder: iter = " << iter << " m = " << m << " tol = " << tol << endl ;
      result= b;
      _function->restoreXVec() ;      
      return kTRUE;
    }
  
    if (fabs (e) < tol || fabs (fa) <= fabs (fb)) {
      // Bounds decreasing too slowly: use bisection
      d = m;
      e = m;
    }
    else {
      // Attempt inverse cubic interpolation
      Double_t p, q, r;
      Double_t s = fb / fa;
      
      if (ac_equal) {
	p = 2 * m * s;
	q = 1 - s;
      }
      else {
	q = fa / fc;
	r = fb / fc;
	p = s * (2 * m * q * (q - r) - (b - a) * (r - 1));
	q = (q - 1) * (r - 1) * (s - 1);
      }
      // Check whether we are in bounds
      if (p > 0) {
	q = -q;
      }
      else {
	p = -p;
      }
      
      Double_t min1= 3 * m * q - fabs (tol * q);
      Double_t min2= fabs (e * q);
      if (2 * p < (min1 < min2 ? min1 : min2)) {
	// Accept the interpolation
	e = d;
	d = p / q;
      }
      else {
	// Interpolation failed: use bisection.
	d = m;
	e = m;
      }
    }
    // Move last best guess to a
    a = b;
    fa = fb;
    // Evaluate new trial root
    if (fabs (d) > tol) {
      b += d;
    }
    else {
      b += (m > 0 ? +tol : -tol);
    }
    fb= (*_function)(&b) - value;

  }
  // Return our best guess if we run out of iterations
  oocoutE((TObject*)0,Eval) << "RooBrentRootFinder::findRoot(" << _function->getName() << "): maximum iterations exceeded." << endl;
  result= b;

  _function->restoreXVec() ;

  return kFALSE;
}
 RooBrentRootFinder.cxx:1
 RooBrentRootFinder.cxx:2
 RooBrentRootFinder.cxx:3
 RooBrentRootFinder.cxx:4
 RooBrentRootFinder.cxx:5
 RooBrentRootFinder.cxx:6
 RooBrentRootFinder.cxx:7
 RooBrentRootFinder.cxx:8
 RooBrentRootFinder.cxx:9
 RooBrentRootFinder.cxx:10
 RooBrentRootFinder.cxx:11
 RooBrentRootFinder.cxx:12
 RooBrentRootFinder.cxx:13
 RooBrentRootFinder.cxx:14
 RooBrentRootFinder.cxx:15
 RooBrentRootFinder.cxx:16
 RooBrentRootFinder.cxx:17
 RooBrentRootFinder.cxx:18
 RooBrentRootFinder.cxx:19
 RooBrentRootFinder.cxx:20
 RooBrentRootFinder.cxx:21
 RooBrentRootFinder.cxx:22
 RooBrentRootFinder.cxx:23
 RooBrentRootFinder.cxx:24
 RooBrentRootFinder.cxx:25
 RooBrentRootFinder.cxx:26
 RooBrentRootFinder.cxx:27
 RooBrentRootFinder.cxx:28
 RooBrentRootFinder.cxx:29
 RooBrentRootFinder.cxx:30
 RooBrentRootFinder.cxx:31
 RooBrentRootFinder.cxx:32
 RooBrentRootFinder.cxx:33
 RooBrentRootFinder.cxx:34
 RooBrentRootFinder.cxx:35
 RooBrentRootFinder.cxx:36
 RooBrentRootFinder.cxx:37
 RooBrentRootFinder.cxx:38
 RooBrentRootFinder.cxx:39
 RooBrentRootFinder.cxx:40
 RooBrentRootFinder.cxx:41
 RooBrentRootFinder.cxx:42
 RooBrentRootFinder.cxx:43
 RooBrentRootFinder.cxx:44
 RooBrentRootFinder.cxx:45
 RooBrentRootFinder.cxx:46
 RooBrentRootFinder.cxx:47
 RooBrentRootFinder.cxx:48
 RooBrentRootFinder.cxx:49
 RooBrentRootFinder.cxx:50
 RooBrentRootFinder.cxx:51
 RooBrentRootFinder.cxx:52
 RooBrentRootFinder.cxx:53
 RooBrentRootFinder.cxx:54
 RooBrentRootFinder.cxx:55
 RooBrentRootFinder.cxx:56
 RooBrentRootFinder.cxx:57
 RooBrentRootFinder.cxx:58
 RooBrentRootFinder.cxx:59
 RooBrentRootFinder.cxx:60
 RooBrentRootFinder.cxx:61
 RooBrentRootFinder.cxx:62
 RooBrentRootFinder.cxx:63
 RooBrentRootFinder.cxx:64
 RooBrentRootFinder.cxx:65
 RooBrentRootFinder.cxx:66
 RooBrentRootFinder.cxx:67
 RooBrentRootFinder.cxx:68
 RooBrentRootFinder.cxx:69
 RooBrentRootFinder.cxx:70
 RooBrentRootFinder.cxx:71
 RooBrentRootFinder.cxx:72
 RooBrentRootFinder.cxx:73
 RooBrentRootFinder.cxx:74
 RooBrentRootFinder.cxx:75
 RooBrentRootFinder.cxx:76
 RooBrentRootFinder.cxx:77
 RooBrentRootFinder.cxx:78
 RooBrentRootFinder.cxx:79
 RooBrentRootFinder.cxx:80
 RooBrentRootFinder.cxx:81
 RooBrentRootFinder.cxx:82
 RooBrentRootFinder.cxx:83
 RooBrentRootFinder.cxx:84
 RooBrentRootFinder.cxx:85
 RooBrentRootFinder.cxx:86
 RooBrentRootFinder.cxx:87
 RooBrentRootFinder.cxx:88
 RooBrentRootFinder.cxx:89
 RooBrentRootFinder.cxx:90
 RooBrentRootFinder.cxx:91
 RooBrentRootFinder.cxx:92
 RooBrentRootFinder.cxx:93
 RooBrentRootFinder.cxx:94
 RooBrentRootFinder.cxx:95
 RooBrentRootFinder.cxx:96
 RooBrentRootFinder.cxx:97
 RooBrentRootFinder.cxx:98
 RooBrentRootFinder.cxx:99
 RooBrentRootFinder.cxx:100
 RooBrentRootFinder.cxx:101
 RooBrentRootFinder.cxx:102
 RooBrentRootFinder.cxx:103
 RooBrentRootFinder.cxx:104
 RooBrentRootFinder.cxx:105
 RooBrentRootFinder.cxx:106
 RooBrentRootFinder.cxx:107
 RooBrentRootFinder.cxx:108
 RooBrentRootFinder.cxx:109
 RooBrentRootFinder.cxx:110
 RooBrentRootFinder.cxx:111
 RooBrentRootFinder.cxx:112
 RooBrentRootFinder.cxx:113
 RooBrentRootFinder.cxx:114
 RooBrentRootFinder.cxx:115
 RooBrentRootFinder.cxx:116
 RooBrentRootFinder.cxx:117
 RooBrentRootFinder.cxx:118
 RooBrentRootFinder.cxx:119
 RooBrentRootFinder.cxx:120
 RooBrentRootFinder.cxx:121
 RooBrentRootFinder.cxx:122
 RooBrentRootFinder.cxx:123
 RooBrentRootFinder.cxx:124
 RooBrentRootFinder.cxx:125
 RooBrentRootFinder.cxx:126
 RooBrentRootFinder.cxx:127
 RooBrentRootFinder.cxx:128
 RooBrentRootFinder.cxx:129
 RooBrentRootFinder.cxx:130
 RooBrentRootFinder.cxx:131
 RooBrentRootFinder.cxx:132
 RooBrentRootFinder.cxx:133
 RooBrentRootFinder.cxx:134
 RooBrentRootFinder.cxx:135
 RooBrentRootFinder.cxx:136
 RooBrentRootFinder.cxx:137
 RooBrentRootFinder.cxx:138
 RooBrentRootFinder.cxx:139
 RooBrentRootFinder.cxx:140
 RooBrentRootFinder.cxx:141
 RooBrentRootFinder.cxx:142
 RooBrentRootFinder.cxx:143
 RooBrentRootFinder.cxx:144
 RooBrentRootFinder.cxx:145
 RooBrentRootFinder.cxx:146
 RooBrentRootFinder.cxx:147
 RooBrentRootFinder.cxx:148
 RooBrentRootFinder.cxx:149
 RooBrentRootFinder.cxx:150
 RooBrentRootFinder.cxx:151
 RooBrentRootFinder.cxx:152
 RooBrentRootFinder.cxx:153
 RooBrentRootFinder.cxx:154
 RooBrentRootFinder.cxx:155
 RooBrentRootFinder.cxx:156
 RooBrentRootFinder.cxx:157
 RooBrentRootFinder.cxx:158
 RooBrentRootFinder.cxx:159
 RooBrentRootFinder.cxx:160
 RooBrentRootFinder.cxx:161
 RooBrentRootFinder.cxx:162
 RooBrentRootFinder.cxx:163
 RooBrentRootFinder.cxx:164
 RooBrentRootFinder.cxx:165
 RooBrentRootFinder.cxx:166