/src/quantlib/ql/math/solvers1d/brent.hpp
Line | Count | Source |
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 | | <https://www.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 | | |
29 | | namespace 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 | 0 | 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 | 0 | Real min1, min2; |
50 | 0 | 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 | 0 | froot = f(root_); |
56 | 0 | ++evaluationNumber_; |
57 | 0 | if (froot * fxMin_ < 0) { |
58 | 0 | xMax_ = xMin_; |
59 | 0 | fxMax_ = fxMin_; |
60 | 0 | } else { |
61 | 0 | xMin_ = xMax_; |
62 | 0 | fxMin_ = fxMax_; |
63 | 0 | } |
64 | 0 | Real d = root_- xMax_; |
65 | 0 | Real e = d; |
66 | |
|
67 | 0 | while (evaluationNumber_<=maxEvaluations_) { |
68 | 0 | if ((froot > 0.0 && fxMax_ > 0.0) || |
69 | 0 | (froot < 0.0 && fxMax_ < 0.0)) { |
70 | | |
71 | | // Rename xMin_, root_, xMax_ and adjust bounds |
72 | 0 | xMax_=xMin_; |
73 | 0 | fxMax_=fxMin_; |
74 | 0 | e=d=root_-xMin_; |
75 | 0 | } |
76 | 0 | if (std::fabs(fxMax_) < std::fabs(froot)) { |
77 | 0 | xMin_=root_; |
78 | 0 | root_=xMax_; |
79 | 0 | xMax_=xMin_; |
80 | 0 | fxMin_=froot; |
81 | 0 | froot=fxMax_; |
82 | 0 | fxMax_=fxMin_; |
83 | 0 | } |
84 | | // Convergence check |
85 | 0 | xAcc1=2.0*QL_EPSILON*std::fabs(root_)+0.5*xAccuracy; |
86 | 0 | xMid=(xMax_-root_)/2.0; |
87 | 0 | if (std::fabs(xMid) <= xAcc1 || (close(froot, 0.0))) { |
88 | 0 | f(root_); |
89 | 0 | ++evaluationNumber_; |
90 | 0 | return root_; |
91 | 0 | } |
92 | 0 | if (std::fabs(e) >= xAcc1 && |
93 | 0 | std::fabs(fxMin_) > std::fabs(froot)) { |
94 | | |
95 | | // Attempt inverse quadratic interpolation |
96 | 0 | s=froot/fxMin_; |
97 | 0 | if (close(xMin_,xMax_)) { |
98 | 0 | p=2.0*xMid*s; |
99 | 0 | q=1.0-s; |
100 | 0 | } else { |
101 | 0 | q=fxMin_/fxMax_; |
102 | 0 | r=froot/fxMax_; |
103 | 0 | p=s*(2.0*xMid*q*(q-r)-(root_-xMin_)*(r-1.0)); |
104 | 0 | q=(q-1.0)*(r-1.0)*(s-1.0); |
105 | 0 | } |
106 | 0 | if (p > 0.0) q = -q; // Check whether in bounds |
107 | 0 | p=std::fabs(p); |
108 | 0 | min1=3.0*xMid*q-std::fabs(xAcc1*q); |
109 | 0 | min2=std::fabs(e*q); |
110 | 0 | if (2.0*p < (min1 < min2 ? min1 : min2)) { |
111 | 0 | e=d; // Accept interpolation |
112 | 0 | d=p/q; |
113 | 0 | } else { |
114 | 0 | d=xMid; // Interpolation failed, use bisection |
115 | 0 | e=d; |
116 | 0 | } |
117 | 0 | } else { |
118 | | // Bounds decreasing too slowly, use bisection |
119 | 0 | d=xMid; |
120 | 0 | e=d; |
121 | 0 | } |
122 | 0 | xMin_=root_; |
123 | 0 | fxMin_=froot; |
124 | 0 | if (std::fabs(d) > xAcc1) |
125 | 0 | root_ += d; |
126 | 0 | else |
127 | 0 | root_ += sign(xAcc1,xMid); |
128 | 0 | froot=f(root_); |
129 | 0 | ++evaluationNumber_; |
130 | 0 | } |
131 | 0 | QL_FAIL("maximum number of function evaluations (" |
132 | 0 | << maxEvaluations_ << ") exceeded"); |
133 | 0 | } Unexecuted instantiation: cashflows.cpp:double QuantLib::Brent::solveImpl<QuantLib::CashFlows::zSpread(std::__1::vector<boost::shared_ptr<QuantLib::CashFlow>, std::__1::allocator<boost::shared_ptr<QuantLib::CashFlow> > > const&, double, boost::shared_ptr<QuantLib::YieldTermStructure> const&, QuantLib::Compounding, QuantLib::Frequency, std::__1::optional<bool> const&, QuantLib::Date, QuantLib::Date, double, unsigned long, double)::$_0>(QuantLib::CashFlows::zSpread(std::__1::vector<boost::shared_ptr<QuantLib::CashFlow>, std::__1::allocator<boost::shared_ptr<QuantLib::CashFlow> > > const&, double, boost::shared_ptr<QuantLib::YieldTermStructure> const&, QuantLib::Compounding, QuantLib::Frequency, std::__1::optional<bool> const&, QuantLib::Date, QuantLib::Date, double, unsigned long, double)::$_0 const&, double) const Unexecuted instantiation: double QuantLib::Brent::solveImpl<QuantLib::LinearTsrPricer::VegaRatioHelper>(QuantLib::LinearTsrPricer::VegaRatioHelper const&, double) const Unexecuted instantiation: double QuantLib::Brent::solveImpl<QuantLib::LinearTsrPricer::PriceHelper>(QuantLib::LinearTsrPricer::PriceHelper const&, double) const Unexecuted instantiation: double QuantLib::Brent::solveImpl<QuantLib::CallableBond::ImpliedVolHelper>(QuantLib::CallableBond::ImpliedVolHelper const&, double) const Unexecuted instantiation: callablebond.cpp:double QuantLib::Brent::solveImpl<QuantLib::(anonymous namespace)::OASHelper>(QuantLib::(anonymous namespace)::OASHelper const&, double) const Unexecuted instantiation: cdsoption.cpp:double QuantLib::Brent::solveImpl<QuantLib::(anonymous namespace)::ImpliedVolHelper>(QuantLib::(anonymous namespace)::ImpliedVolHelper const&, double) const Unexecuted instantiation: randomdefaultmodel.cpp:double QuantLib::Brent::solveImpl<QuantLib::(anonymous namespace)::Root>(QuantLib::(anonymous namespace)::Root const&, double) const Unexecuted instantiation: syntheticcdo.cpp:double QuantLib::Brent::solveImpl<QuantLib::(anonymous namespace)::ObjectiveFunction>(QuantLib::(anonymous namespace)::ObjectiveFunction const&, double) const Unexecuted instantiation: convolvedstudentt.cpp:double QuantLib::Brent::solveImpl<QuantLib::InverseCumulativeBehrensFisher::operator()(double) const::$_0>(QuantLib::InverseCumulativeBehrensFisher::operator()(double) const::$_0 const&, double) const Unexecuted instantiation: double QuantLib::Brent::solveImpl<QuantLib::GeneralizedHullWhite::Helper>(QuantLib::GeneralizedHullWhite::Helper const&, double) const Unexecuted instantiation: noarbsabr.cpp:double QuantLib::Brent::solveImpl<QuantLib::NoArbSabrModel::NoArbSabrModel(double, double, double, double, double, double)::$_0>(QuantLib::NoArbSabrModel::NoArbSabrModel(double, double, double, double, double, double)::$_0 const&, double) const Unexecuted instantiation: creditdefaultswap.cpp:double QuantLib::Brent::solveImpl<QuantLib::(anonymous namespace)::ObjectiveFunction>(QuantLib::(anonymous namespace)::ObjectiveFunction const&, double) const Unexecuted instantiation: impliedvolatility.cpp:double QuantLib::Brent::solveImpl<QuantLib::(anonymous namespace)::PriceError>(QuantLib::(anonymous namespace)::PriceError const&, double) const Unexecuted instantiation: chisquaredistribution.cpp:double QuantLib::Brent::solveImpl<QuantLib::InverseNonCentralCumulativeChiSquareDistribution::operator()(double) const::$_0>(QuantLib::InverseNonCentralCumulativeChiSquareDistribution::operator()(double) const::$_0 const&, double) const Unexecuted instantiation: richardsonextrapolation.cpp:double QuantLib::Brent::solveImpl<QuantLib::(anonymous namespace)::RichardsonEqn>(QuantLib::(anonymous namespace)::RichardsonEqn const&, double) const Unexecuted instantiation: concentrating1dmesher.cpp:double QuantLib::Brent::solveImpl<QuantLib::Concentrating1dMesher::Concentrating1dMesher(double, double, unsigned long, std::__1::vector<std::__1::tuple<double, double, bool>, std::__1::allocator<std::__1::tuple<double, double, bool> > > const&, double)::$_0>(QuantLib::Concentrating1dMesher::Concentrating1dMesher(double, double, unsigned long, std::__1::vector<std::__1::tuple<double, double, bool>, std::__1::allocator<std::__1::tuple<double, double, bool> > > const&, double)::$_0 const&, double) const Unexecuted instantiation: concentrating1dmesher.cpp:double QuantLib::Brent::solveImpl<QuantLib::Concentrating1dMesher::Concentrating1dMesher(double, double, unsigned long, std::__1::vector<std::__1::tuple<double, double, bool>, std::__1::allocator<std::__1::tuple<double, double, bool> > > const&, double)::$_1>(QuantLib::Concentrating1dMesher::Concentrating1dMesher(double, double, unsigned long, std::__1::vector<std::__1::tuple<double, double, bool>, std::__1::allocator<std::__1::tuple<double, double, bool> > > const&, double)::$_1 const&, double) const Unexecuted instantiation: cevrndcalculator.cpp:double QuantLib::Brent::solveImpl<QuantLib::CEVRNDCalculator::invcdf(double, double) const::$_0>(QuantLib::CEVRNDCalculator::invcdf(double, double) const::$_0 const&, double) const Unexecuted instantiation: gbsmrndcalculator.cpp:double QuantLib::Brent::solveImpl<QuantLib::GBSMRNDCalculator::invcdf(double, double) const::$_0>(QuantLib::GBSMRNDCalculator::invcdf(double, double) const::$_0 const&, double) const Unexecuted instantiation: riskneutraldensitycalculator.cpp:double QuantLib::Brent::solveImpl<QuantLib::RiskNeutralDensityCalculator::InvCDFHelper::inverseCDF(double, double) const::$_0>(QuantLib::RiskNeutralDensityCalculator::InvCDFHelper::inverseCDF(double, double) const::$_0 const&, double) const Unexecuted instantiation: calibrationhelper.cpp:double QuantLib::Brent::solveImpl<QuantLib::BlackCalibrationHelper::impliedVolatility(double, double, unsigned long, double, double) const::$_0>(QuantLib::BlackCalibrationHelper::impliedVolatility(double, double, unsigned long, double, double) const::$_0 const&, double) const Unexecuted instantiation: swaptionpseudojacobian.cpp:double QuantLib::Brent::solveImpl<QuantLib::(anonymous namespace)::QuickCap>(QuantLib::(anonymous namespace)::QuickCap const&, double) const Unexecuted instantiation: double QuantLib::Brent::solveImpl<QuantLib::OneFactorModel::ShortRateTree::Helper>(QuantLib::OneFactorModel::ShortRateTree::Helper const&, double) const Unexecuted instantiation: double QuantLib::Brent::solveImpl<QuantLib::BlackKarasinski::Helper>(QuantLib::BlackKarasinski::Helper const&, double) const Unexecuted instantiation: double QuantLib::Brent::solveImpl<QuantLib::MarkovFunctional::ZeroHelper>(QuantLib::MarkovFunctional::ZeroHelper const&, double) const Unexecuted instantiation: double QuantLib::Brent::solveImpl<QuantLib::G2::SwaptionPricingFunction::SolvingFunction>(QuantLib::G2::SwaptionPricingFunction::SolvingFunction const&, double) const Unexecuted instantiation: double QuantLib::Brent::solveImpl<QuantLib::detail::SumExponentialsRootSolver>(QuantLib::detail::SumExponentialsRootSolver const&, double) const Unexecuted instantiation: blackdeltacalculator.cpp:double QuantLib::Brent::solveImpl<QuantLib::BlackDeltaCalculator::strikeFromDelta(double, QuantLib::DeltaVolQuote::DeltaType) const::$_0>(QuantLib::BlackDeltaCalculator::strikeFromDelta(double, QuantLib::DeltaVolQuote::DeltaType) const::$_0 const&, double) const Unexecuted instantiation: blackdeltacalculator.cpp:double QuantLib::Brent::solveImpl<QuantLib::BlackDeltaCalculator::strikeFromDelta(double, QuantLib::DeltaVolQuote::DeltaType) const::$_1>(QuantLib::BlackDeltaCalculator::strikeFromDelta(double, QuantLib::DeltaVolQuote::DeltaType) const::$_1 const&, double) const Unexecuted instantiation: analyticcompoundoptionengine.cpp:double QuantLib::Brent::solveImpl<QuantLib::(anonymous namespace)::ImpliedSpotHelper>(QuantLib::(anonymous namespace)::ImpliedSpotHelper const&, double) const Unexecuted instantiation: double QuantLib::Brent::solveImpl<QuantLib::Gaussian1dJamshidianSwaptionEngine::rStarFinder>(QuantLib::Gaussian1dJamshidianSwaptionEngine::rStarFinder const&, double) const Unexecuted instantiation: double QuantLib::Brent::solveImpl<QuantLib::JamshidianSwaptionEngine::rStarFinder>(QuantLib::JamshidianSwaptionEngine::rStarFinder const&, double) const Unexecuted instantiation: analytichestonengine.cpp:double QuantLib::Brent::solveImpl<QuantLib::AnalyticHestonEngine::OptimalAlpha::alphaMax(double) const::$_0>(QuantLib::AnalyticHestonEngine::OptimalAlpha::alphaMax(double) const::$_0 const&, double) const Unexecuted instantiation: analytichestonengine.cpp:double QuantLib::Brent::solveImpl<QuantLib::AnalyticHestonEngine::OptimalAlpha::alphaMin(double) const::$_0>(QuantLib::AnalyticHestonEngine::OptimalAlpha::alphaMin(double) const::$_0 const&, double) const Unexecuted instantiation: analytichestonengine.cpp:double QuantLib::Brent::solveImpl<QuantLib::(anonymous namespace)::u_Max>(QuantLib::(anonymous namespace)::u_Max const&, double) const Unexecuted instantiation: analytichestonengine.cpp:double QuantLib::Brent::solveImpl<QuantLib::(anonymous namespace)::uHat_Max>(QuantLib::(anonymous namespace)::uHat_Max const&, double) const Unexecuted instantiation: double QuantLib::Brent::solveImpl<QuantLib::QdPlusBoundaryEvaluator>(QuantLib::QdPlusBoundaryEvaluator const&, double) const Unexecuted instantiation: hestonprocess.cpp:double QuantLib::Brent::solveImpl<QuantLib::HestonProcess::evolve(double, QuantLib::Array const&, double, QuantLib::Array const&) const::$_0>(QuantLib::HestonProcess::evolve(double, QuantLib::Array const&, double, QuantLib::Array const&) const::$_0 const&, double) const Unexecuted instantiation: hestonblackvolsurface.cpp:double QuantLib::Brent::solveImpl<QuantLib::HestonBlackVolSurface::blackVolImpl(double, double) const::$_0>(QuantLib::HestonBlackVolSurface::blackVolImpl(double, double) const::$_0 const&, double) const Unexecuted instantiation: double QuantLib::Brent::solveImpl<QuantLib::KahaleSmileSection::sHelper1>(QuantLib::KahaleSmileSection::sHelper1 const&, double) const Unexecuted instantiation: double QuantLib::Brent::solveImpl<QuantLib::KahaleSmileSection::aHelper>(QuantLib::KahaleSmileSection::aHelper const&, double) const Unexecuted instantiation: double QuantLib::Brent::solveImpl<QuantLib::KahaleSmileSection::sHelper>(QuantLib::KahaleSmileSection::sHelper const&, double) const Unexecuted instantiation: double QuantLib::Brent::solveImpl<QuantLib::OptionletStripper2::ObjectiveFunction>(QuantLib::OptionletStripper2::ObjectiveFunction const&, double) const |
134 | | private: |
135 | 0 | Real sign(Real a, Real b) const { |
136 | 0 | return b >= 0.0 ? Real(std::fabs(a)) : Real(-std::fabs(a)); |
137 | 0 | } |
138 | | }; |
139 | | |
140 | | } |
141 | | |
142 | | #endif |