3 Nikolai Amelin, Ludmila Malinina, Timur Pocheptsov (C) JINR/Dubna
4 amelin@sunhe.jinr.ru, malinina@sunhe.jinr.ru, pocheptsov@sunhe.jinr.ru
8 //This equation solver class is taken from GEANT4 and modified!!
13 template <class Function>
14 Bool_t EquationSolver<Function>::Brent(Function& theFunction) {
15 const Double_t precision = 3.0e-8;
17 // Check the interval before start
18 if(fA > fB || TMath::Abs(fA-fB) <= fTolerance) {
19 std::cerr << "EquationSolver::Brent: The interval must be properly set." << std::endl;
22 Double_t fa = theFunction(fA);
23 Double_t fb = theFunction(fB);
25 std::cerr << "EquationSolver::Brent: The interval must include a root." << std::endl;
36 for(Int_t i=0; i < fMaxIter; i++) {
37 // Rename a,b,c and adjust bounding interval d
44 if(TMath::Abs(fc) < TMath::Abs(fb)) {
52 Double_t Tol1 = 2.0*precision*TMath::Abs(fB) + 0.5*fTolerance;
53 Double_t xm = 0.5*(c-fB);
54 if(TMath::Abs(xm) <= Tol1 || fb == 0.0) {
58 // Inverse quadratic interpolation
59 if(TMath::Abs(e) >= Tol1 && TMath::Abs(fa) > TMath::Abs(fb)) {
70 p = s*(2.0*xm*q*(q-r)-(fB-fA)*(r-1.0));
71 q = (q-1.0)*(r-1.0)*(s-1.0);
76 Double_t min1 = 3.0*xm*q-TMath::Abs(Tol1*q);
77 Double_t min2 = TMath::Abs(e*q);
78 if (2.0*p < TMath::Min(min1,min2)) {
90 // Bounds decreasing too slowly, use bisection
94 // Move last guess to a
97 if (TMath::Abs(d) > Tol1) fB += d;
99 if (xm >= 0.0) fB += TMath::Abs(Tol1);
100 else fB -= TMath::Abs(Tol1);
102 fb = theFunction(fB);
104 std::cerr << "EquationSolver::Brent: Number of iterations exceeded." << std::endl;
108 template <class Function>
109 void EquationSolver<Function>::SetIntervalLimits(const Double_t Limit1, const Double_t Limit2) {
110 if(TMath::Abs(Limit1-Limit2) <= fTolerance) {
111 std::cerr << "EquationSolver::SetIntervalLimits: Interval must be wider than tolerance." << std::endl;
114 if(Limit1 < Limit2) {