Coverage Report

Created: 2026-09-28 06:23

next uncovered line (L), next uncovered region (R), next uncovered branch (B)
/src/quantlib/ql/math/solvers1d/newtonsafe.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 newtonsafe.hpp
21
    \brief Safe (bracketed) Newton 1-D solver
22
*/
23
24
#ifndef quantlib_solver1d_newtonsafe_h
25
#define quantlib_solver1d_newtonsafe_h
26
27
#include <ql/math/solver1d.hpp>
28
29
namespace QuantLib {
30
31
    //! safe %Newton 1-D solver
32
    /*! \note This solver requires that the passed function object
33
              implement a method <tt>Real derivative(Real)</tt>.
34
35
        \test the correctness of the returned values is tested by
36
              checking them against known good results.
37
38
        \ingroup solvers
39
    */
40
    class NewtonSafe : public Solver1D<NewtonSafe> {
41
      public:
42
        template <class F>
43
        Real solveImpl(const F& f,
44
539
                       Real xAccuracy) const {
45
46
            /* The implementation of the algorithm was inspired by
47
               Press, Teukolsky, Vetterling, and Flannery,
48
               "Numerical Recipes in C", 2nd edition,
49
               Cambridge University Press
50
            */
51
52
539
            Real froot, dfroot, dx, dxold;
53
539
            Real xh, xl;
54
55
            // Orient the search so that f(xl) < 0
56
539
            if (fxMin_ < 0.0) {
57
143
                xl = xMin_;
58
143
                xh = xMax_;
59
396
            } else {
60
396
                xh = xMin_;
61
396
                xl = xMax_;
62
396
            }
63
64
            // the "stepsize before last"
65
539
            dxold = xMax_-xMin_;
66
            // it was dxold=std::fabs(xMax_-xMin_); in Numerical Recipes
67
            // here (xMax_-xMin_ > 0) is verified in the constructor
68
69
            // and the last step
70
539
            dx = dxold;
71
72
539
            froot = f(root_);
73
539
            dfroot = f.derivative(root_);
74
539
            QL_REQUIRE(dfroot != Null<Real>(),
75
539
                       "NewtonSafe requires function's derivative");
76
539
            ++evaluationNumber_;
77
78
7.77k
            while (evaluationNumber_<=maxEvaluations_) {
79
                // Bisect if (out of range || not decreasing fast enough)
80
7.73k
                if ((((root_-xh)*dfroot-froot)*
81
7.73k
                     ((root_-xl)*dfroot-froot) > 0.0)
82
4.94k
                    || (std::fabs(2.0*froot) > std::fabs(dxold*dfroot))) {
83
84
4.94k
                    dxold = dx;
85
4.94k
                    dx = (xh-xl)/2.0;
86
4.94k
                    root_=xl+dx;
87
4.94k
                } else {
88
2.78k
                    dxold = dx;
89
2.78k
                    dx = froot/dfroot;
90
2.78k
                    root_ -= dx;
91
2.78k
                }
92
                // Convergence criterion
93
7.73k
                if (std::fabs(dx) < xAccuracy) {
94
494
                    f(root_);
95
494
                    ++evaluationNumber_;
96
494
                    return root_;
97
494
                }
98
7.23k
                froot = f(root_);
99
7.23k
                dfroot = f.derivative(root_);
100
7.23k
                ++evaluationNumber_;
101
7.23k
                if (froot < 0.0)
102
3.08k
                    xl=root_;
103
4.15k
                else
104
4.15k
                    xh=root_;
105
7.23k
            }
106
107
539
            QL_FAIL("maximum number of function evaluations ("
108
45
                    << maxEvaluations_ << ") exceeded");
109
45
        }
double QuantLib::NewtonSafe::solveImpl<QuantLib::CashFlows::IrrFinder>(QuantLib::CashFlows::IrrFinder const&, double) const
Line
Count
Source
44
396
                       Real xAccuracy) const {
45
46
            /* The implementation of the algorithm was inspired by
47
               Press, Teukolsky, Vetterling, and Flannery,
48
               "Numerical Recipes in C", 2nd edition,
49
               Cambridge University Press
50
            */
51
52
396
            Real froot, dfroot, dx, dxold;
53
396
            Real xh, xl;
54
55
            // Orient the search so that f(xl) < 0
56
396
            if (fxMin_ < 0.0) {
57
0
                xl = xMin_;
58
0
                xh = xMax_;
59
396
            } else {
60
396
                xh = xMin_;
61
396
                xl = xMax_;
62
396
            }
63
64
            // the "stepsize before last"
65
396
            dxold = xMax_-xMin_;
66
            // it was dxold=std::fabs(xMax_-xMin_); in Numerical Recipes
67
            // here (xMax_-xMin_ > 0) is verified in the constructor
68
69
            // and the last step
70
396
            dx = dxold;
71
72
396
            froot = f(root_);
73
396
            dfroot = f.derivative(root_);
74
396
            QL_REQUIRE(dfroot != Null<Real>(),
75
396
                       "NewtonSafe requires function's derivative");
76
396
            ++evaluationNumber_;
77
78
5.01k
            while (evaluationNumber_<=maxEvaluations_) {
79
                // Bisect if (out of range || not decreasing fast enough)
80
4.97k
                if ((((root_-xh)*dfroot-froot)*
81
4.97k
                     ((root_-xl)*dfroot-froot) > 0.0)
82
2.73k
                    || (std::fabs(2.0*froot) > std::fabs(dxold*dfroot))) {
83
84
2.61k
                    dxold = dx;
85
2.61k
                    dx = (xh-xl)/2.0;
86
2.61k
                    root_=xl+dx;
87
2.61k
                } else {
88
2.36k
                    dxold = dx;
89
2.36k
                    dx = froot/dfroot;
90
2.36k
                    root_ -= dx;
91
2.36k
                }
92
                // Convergence criterion
93
4.97k
                if (std::fabs(dx) < xAccuracy) {
94
351
                    f(root_);
95
351
                    ++evaluationNumber_;
96
351
                    return root_;
97
351
                }
98
4.62k
                froot = f(root_);
99
4.62k
                dfroot = f.derivative(root_);
100
4.62k
                ++evaluationNumber_;
101
4.62k
                if (froot < 0.0)
102
1.87k
                    xl=root_;
103
2.74k
                else
104
2.74k
                    xh=root_;
105
4.62k
            }
106
107
396
            QL_FAIL("maximum number of function evaluations ("
108
45
                    << maxEvaluations_ << ") exceeded");
109
45
        }
Unexecuted instantiation: double QuantLib::NewtonSafe::solveImpl<QuantLib::GFunctionFactory::GFunctionWithShifts::ObjectiveFunction>(QuantLib::GFunctionFactory::GFunctionWithShifts::ObjectiveFunction const&, double) const
Unexecuted instantiation: irregularswaption.cpp:double QuantLib::NewtonSafe::solveImpl<QuantLib::(anonymous namespace)::IrregularImpliedVolHelper>(QuantLib::(anonymous namespace)::IrregularImpliedVolHelper const&, double) const
Unexecuted instantiation: capfloor.cpp:double QuantLib::NewtonSafe::solveImpl<QuantLib::(anonymous namespace)::ImpliedCapVolHelper>(QuantLib::(anonymous namespace)::ImpliedCapVolHelper const&, double) const
Unexecuted instantiation: swaption.cpp:double QuantLib::NewtonSafe::solveImpl<QuantLib::(anonymous namespace)::ImpliedSwaptionVolHelper>(QuantLib::(anonymous namespace)::ImpliedSwaptionVolHelper const&, double) const
Unexecuted instantiation: double QuantLib::NewtonSafe::solveImpl<QuantLib::detail::SumExponentialsRootSolver>(QuantLib::detail::SumExponentialsRootSolver const&, double) const
double QuantLib::NewtonSafe::solveImpl<QuantLib::BlackImpliedStdDevHelper>(QuantLib::BlackImpliedStdDevHelper const&, double) const
Line
Count
Source
44
143
                       Real xAccuracy) const {
45
46
            /* The implementation of the algorithm was inspired by
47
               Press, Teukolsky, Vetterling, and Flannery,
48
               "Numerical Recipes in C", 2nd edition,
49
               Cambridge University Press
50
            */
51
52
143
            Real froot, dfroot, dx, dxold;
53
143
            Real xh, xl;
54
55
            // Orient the search so that f(xl) < 0
56
143
            if (fxMin_ < 0.0) {
57
143
                xl = xMin_;
58
143
                xh = xMax_;
59
143
            } else {
60
0
                xh = xMin_;
61
0
                xl = xMax_;
62
0
            }
63
64
            // the "stepsize before last"
65
143
            dxold = xMax_-xMin_;
66
            // it was dxold=std::fabs(xMax_-xMin_); in Numerical Recipes
67
            // here (xMax_-xMin_ > 0) is verified in the constructor
68
69
            // and the last step
70
143
            dx = dxold;
71
72
143
            froot = f(root_);
73
143
            dfroot = f.derivative(root_);
74
143
            QL_REQUIRE(dfroot != Null<Real>(),
75
143
                       "NewtonSafe requires function's derivative");
76
143
            ++evaluationNumber_;
77
78
2.75k
            while (evaluationNumber_<=maxEvaluations_) {
79
                // Bisect if (out of range || not decreasing fast enough)
80
2.75k
                if ((((root_-xh)*dfroot-froot)*
81
2.75k
                     ((root_-xl)*dfroot-froot) > 0.0)
82
2.33k
                    || (std::fabs(2.0*froot) > std::fabs(dxold*dfroot))) {
83
84
2.33k
                    dxold = dx;
85
2.33k
                    dx = (xh-xl)/2.0;
86
2.33k
                    root_=xl+dx;
87
2.33k
                } else {
88
420
                    dxold = dx;
89
420
                    dx = froot/dfroot;
90
420
                    root_ -= dx;
91
420
                }
92
                // Convergence criterion
93
2.75k
                if (std::fabs(dx) < xAccuracy) {
94
143
                    f(root_);
95
143
                    ++evaluationNumber_;
96
143
                    return root_;
97
143
                }
98
2.61k
                froot = f(root_);
99
2.61k
                dfroot = f.derivative(root_);
100
2.61k
                ++evaluationNumber_;
101
2.61k
                if (froot < 0.0)
102
1.20k
                    xl=root_;
103
1.40k
                else
104
1.40k
                    xh=root_;
105
2.61k
            }
106
107
143
            QL_FAIL("maximum number of function evaluations ("
108
0
                    << maxEvaluations_ << ") exceeded");
109
0
        }
Unexecuted instantiation: double QuantLib::NewtonSafe::solveImpl<QuantLib::QdPlusBoundaryEvaluator>(QuantLib::QdPlusBoundaryEvaluator const&, double) const
Unexecuted instantiation: double QuantLib::NewtonSafe::solveImpl<QuantLib::Gaussian1dSwaptionVolatility::DateHelper>(QuantLib::Gaussian1dSwaptionVolatility::DateHelper const&, double) const
110
    };
111
112
}
113
114
#endif