QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
brent.hpp
Go to the documentation of this file.
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_);
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
68 if ((froot > 0.0 && fxMax_ > 0.0) ||
69 (froot < 0.0 && fxMax_ < 0.0)) {
70
71 // Rename xMin_, root_, xMax_ and adjust bounds
74 e=d=root_-xMin_;
75 }
76 if (std::fabs(fxMax_) < std::fabs(froot)) {
80 fxMin_=froot;
81 froot=fxMax_;
83 }
84 // Convergence check
85 xAcc1=2.0*QL_EPSILON*std::fabs(root_)+0.5*xAccuracy;
86 xMid=(xMax_-root_)/2.0;
87 if (std::fabs(xMid) <= xAcc1 || (close(froot, 0.0))) {
88 f(root_);
90 return root_;
91 }
92 if (std::fabs(e) >= xAcc1 &&
93 std::fabs(fxMin_) > std::fabs(froot)) {
94
95 // Attempt inverse quadratic interpolation
96 s=froot/fxMin_;
97 if (close(xMin_,xMax_)) {
98 p=2.0*xMid*s;
99 q=1.0-s;
100 } else {
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(p);
108 min1=3.0*xMid*q-std::fabs(xAcc1*q);
109 min2=std::fabs(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(d) > xAcc1)
125 root_ += d;
126 else
127 root_ += sign(xAcc1,xMid);
128 froot=f(root_);
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(a)) : Real(-std::fabs(a));
137 }
138 };
139
140}
141
142#endif
Brent 1-D solver
Definition: brent.hpp:37
Real sign(Real a, Real b) const
Definition: brent.hpp:135
Real solveImpl(const F &f, Real xAccuracy) const
Definition: brent.hpp:40
Base class for 1-D solvers.
Definition: solver1d.hpp:67
#define QL_FAIL(message)
throw an error (possibly with file and line information)
Definition: errors.hpp:92
Date d
ext::function< Real(Real)> b
#define QL_EPSILON
Definition: qldefines.hpp:178
QL_REAL Real
real number
Definition: types.hpp:50
Definition: any.hpp:35
bool close(const Quantity &m1, const Quantity &m2, Size n)
Definition: quantity.cpp:163
ext::shared_ptr< YieldTermStructure > q
ext::shared_ptr< YieldTermStructure > r
Real F
Definition: sabr.cpp:200
Abstract 1-D solver class.