1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2000, 2001, 2002, 2003 RiskMap srl
5
6 This file is part of QuantLib, a free-software/open-source library
7 for financial quantitative analysts and developers - http://quantlib.org/
8
9 QuantLib is free software: you can redistribute it and/or modify it
10 under the terms of the QuantLib license. You should have received a
11 copy of the license along with this program; if not, please email
12 <quantlib-dev@lists.sf.net>. The license is also available online at
13 <http://quantlib.org/license.shtml>.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the license for more details.
18*/
19
20/*! \file brent.hpp
21 \brief Brent 1-D solver
22*/
23
24#ifndef quantlib_solver1d_brent_h
25#define quantlib_solver1d_brent_h
26
27#include <ql/math/solver1d.hpp>
28
29namespace QuantLib {
30
31 //! %Brent 1-D solver
32 /*! \test the correctness of the returned values is tested by
33 checking them against known good results.
34
35 \ingroup solvers
36 */
37 class Brent : public Solver1D<Brent> {
38 public:
39 template <class F>
40 Real solveImpl(const F& f,
41 Real xAccuracy) const {
42
43 /* The implementation of the algorithm was inspired by
44 Press, Teukolsky, Vetterling, and Flannery,
45 "Numerical Recipes in C", 2nd edition, Cambridge
46 University Press
47 */
48
49 Real min1, min2;
50 Real froot, p, q, r, s, xAcc1, xMid;
51
52 // we want to start with root_ (which equals the guess) on
53 // one side of the bracket and both xMin_ and xMax_ on the
54 // other.
55 froot = f(root_);
56 ++evaluationNumber_;
57 if (froot * fxMin_ < 0) {
58 xMax_ = xMin_;
59 fxMax_ = fxMin_;
60 } else {
61 xMin_ = xMax_;
62 fxMin_ = fxMax_;
63 }
64 Real d = root_- xMax_;
65 Real e = d;
66
67 while (evaluationNumber_<=maxEvaluations_) {
68 if ((froot > 0.0 && fxMax_ > 0.0) ||
69 (froot < 0.0 && fxMax_ < 0.0)) {
70
71 // Rename xMin_, root_, xMax_ and adjust bounds
72 xMax_=xMin_;
73 fxMax_=fxMin_;
74 e=d=root_-xMin_;
75 }
76 if (std::fabs(x: fxMax_) < std::fabs(x: froot)) {
77 xMin_=root_;
78 root_=xMax_;
79 xMax_=xMin_;
80 fxMin_=froot;
81 froot=fxMax_;
82 fxMax_=fxMin_;
83 }
84 // Convergence check
85 xAcc1=2.0*QL_EPSILON*std::fabs(x: root_)+0.5*xAccuracy;
86 xMid=(xMax_-root_)/2.0;
87 if (std::fabs(x: xMid) <= xAcc1 || (close(x: froot, y: 0.0))) {
88 f(root_);
89 ++evaluationNumber_;
90 return root_;
91 }
92 if (std::fabs(x: e) >= xAcc1 &&
93 std::fabs(x: fxMin_) > std::fabs(x: froot)) {
94
95 // Attempt inverse quadratic interpolation
96 s=froot/fxMin_;
97 if (close(x: xMin_,y: xMax_)) {
98 p=2.0*xMid*s;
99 q=1.0-s;
100 } else {
101 q=fxMin_/fxMax_;
102 r=froot/fxMax_;
103 p=s*(2.0*xMid*q*(q-r)-(root_-xMin_)*(r-1.0));
104 q=(q-1.0)*(r-1.0)*(s-1.0);
105 }
106 if (p > 0.0) q = -q; // Check whether in bounds
107 p=std::fabs(x: p);
108 min1=3.0*xMid*q-std::fabs(x: xAcc1*q);
109 min2=std::fabs(x: e*q);
110 if (2.0*p < (min1 < min2 ? min1 : min2)) {
111 e=d; // Accept interpolation
112 d=p/q;
113 } else {
114 d=xMid; // Interpolation failed, use bisection
115 e=d;
116 }
117 } else {
118 // Bounds decreasing too slowly, use bisection
119 d=xMid;
120 e=d;
121 }
122 xMin_=root_;
123 fxMin_=froot;
124 if (std::fabs(x: d) > xAcc1)
125 root_ += d;
126 else
127 root_ += sign(a: xAcc1,b: xMid);
128 froot=f(root_);
129 ++evaluationNumber_;
130 }
131 QL_FAIL("maximum number of function evaluations ("
132 << maxEvaluations_ << ") exceeded");
133 }
134 private:
135 Real sign(Real a, Real b) const {
136 return b >= 0.0 ? Real(std::fabs(x: a)) : Real(-std::fabs(x: a));
137 }
138 };
139
140}
141
142#endif
143

source code of quantlib/ql/math/solvers1d/brent.hpp