39944165411

Committed 23 Oct 2026 11:09AM UTC coverage: 79.009% (+0.5%) from 78.485%

Build # 39944165411

Build Type

Pull #2908

github

Committed by

web-flow

Commit Message

Merge 2415ffbff into 14b91b151

Pull Request Pull Request #2908: Avoid underflow in the non-central chi-square distribution

Coverage Stats

70 of 70 new or added lines in 2 files covered. (100.0%)

235 existing lines in 16 files now uncovered.

64818 of 82039 relevant lines covered (79.01%)

8871140.49 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

40.3

/ql/math/abcdmathfunction.cpp

/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*
 Copyright (C) 2006, 2007, 2015 Ferdinando Ametrano
 Copyright (C) 2006 Cristina Duminuco
 Copyright (C) 2005, 2006 Klaus Spanderen
 Copyright (C) 2007 Giorgio Facchinetti
 Copyright (C) 2015 Paolo Mazzocchi

 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/

 QuantLib is free software: you can redistribute it and/or modify it
 under the terms of the QuantLib license.  You should have received a
 copy of the license along with this program; if not, please email
 <quantlib-dev@lists.sf.net>. The license is also available online at
 <https://www.quantlib.org/license.shtml>.

 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the license for more details.
*/

#include <ql/math/abcdmathfunction.hpp>
#include <utility>

namespace QuantLib {

    void AbcdMathFunction::validate(Real a,
                                    Real b,
                                    Real c,
                                    Real d) {
        QL_REQUIRE(c>0, "c (" << c << ") must be positive");
        QL_REQUIRE(d>=0, "d (" << d << ") must be non negative");
        QL_REQUIRE(a+d>=0,
                   "a+d (" << a << "+" << d << ") must be non negative");

        if (b>=0.0)
            return;

        // the one and only stationary point...
        Time zeroFirstDerivative = 1.0/c-a/b;
        if (zeroFirstDerivative>=0.0) {
            // ... is a minimum
            // must be abcd(zeroFirstDerivative)>=0
            QL_REQUIRE(b>=-(d*c)/std::exp(c*a/b-1.0),
                       "b (" << b << ") less than " <<
                       -(d*c)/std::exp(c*a/b-1.0) << ": negative function"
                       " value at stationary point " << zeroFirstDerivative);
        }

    }

    void AbcdMathFunction::initialize_() {
        validate(a_, b_, c_, d_);
        da_ = b_ - c_*a_;
        db_ = -c_*b_;
        dabcd_[0]=da_;
        dabcd_[1]=db_;
        dabcd_[2]=c_;
        dabcd_[3]=0.0;

        pa_ = -(a_ + b_/c_)/c_;
        pb_ = -b_/c_;
        K_ = 0.0;

        dibc_ = b_/c_;
        diacplusbcc_ = a_/c_ + dibc_/c_;
    }

    AbcdMathFunction::AbcdMathFunction(Real aa, Real bb, Real cc, Real dd)
    : a_(aa), b_(bb), c_(cc), d_(dd), abcd_(4), dabcd_(4) {
        abcd_[0]=a_;
        abcd_[1]=b_;
        abcd_[2]=c_;
        abcd_[3]=d_;
        initialize_();
    }

    AbcdMathFunction::AbcdMathFunction(std::vector<Real> abcd) : abcd_(std::move(abcd)), dabcd_(4) {
        a_=abcd_[0];
        b_=abcd_[1];
        c_=abcd_[2];
        d_=abcd_[3];
        initialize_();
    }

    Time AbcdMathFunction::maximumLocation() const {
        if (b_ == 0.0) {
            // f(t) = a*exp(-c*t) + d: monotonically decreasing (a>0) or increasing (a<0)
            if (a_ >= 0.0)
                return 0.0;
            else
                return QL_MAX_REAL;
        }

        // Stationary point: f'(t*) = 0  =>  t* = 1/c - a/b
        Time t_star = 1.0/c_ - a_/b_;

        if (b_ > 0.0) {
            // f''(t*) = -c*b*exp(-c*t*) < 0  =>  t* is a maximum.
            // If t* <= 0 the function is decreasing on [0, inf), maximum is at t=0.
            return (t_star > 0.0 ? t_star : 0.0);
        } else {
            // b_ < 0: f''(t*) > 0  =>  t* is a minimum, not a maximum.
            if (t_star <= 0.0) {
                // No interior minimum on [0, inf): function is monotonically
                // increasing toward d_; maximum is at t -> +inf.
                return QL_MAX_REAL;
            } else {
                // Interior minimum at t* > 0: function decreases then increases.
                // Maximum is at one of the two boundaries: t=0 (value a_+d_) or
                // t -> +inf (value d_).  validate() guarantees a_+d_ >= 0 and d_ >= 0.
                return (a_ >= 0.0 ? 0.0 : QL_MAX_REAL);
            }
        }
    }

    Real AbcdMathFunction::definiteIntegral(Time t1,
                                            Time t2) const {
        return primitive(t2)-primitive(t1);
    }

    std::vector<Real>
    AbcdMathFunction::definiteIntegralCoefficients(Time t,
                                                   Time t2) const {
        Time dt = t2 - t;
        Real expcdt = std::exp(-c_*dt);
        std::vector<Real> result(4);
        result[0] = diacplusbcc_ - (diacplusbcc_ + dibc_*dt)*expcdt;
        result[1] = dibc_ * (1.0 - expcdt);
        result[2] = c_;
        result[3] = d_*dt;
        return result;
    }

    std::vector<Real>
    AbcdMathFunction::definiteDerivativeCoefficients(Time t,
                                                     Time t2) const {
        Time dt = t2 - t;
        Real expcdt = std::exp(-c_*dt);
        std::vector<Real> result(4);
        result[1] = b_*c_/(1.0-expcdt);
        result[0] = a_*c_ - b_ + result[1]*dt*expcdt;
        result[0] /= 1.0-expcdt;
        result[2] = c_;
        result[3] = d_/dt;
        return result;
    }

}

1	/* -- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -- */
2
3	/*
4	Copyright (C) 2006, 2007, 2015 Ferdinando Ametrano
5	Copyright (C) 2006 Cristina Duminuco
6	Copyright (C) 2005, 2006 Klaus Spanderen
7	Copyright (C) 2007 Giorgio Facchinetti
8	Copyright (C) 2015 Paolo Mazzocchi
9
10	This file is part of QuantLib, a free-software/open-source library
11	for financial quantitative analysts and developers - http://quantlib.org/
12
13	QuantLib is free software: you can redistribute it and/or modify it
14	under the terms of the QuantLib license. You should have received a
15	copy of the license along with this program; if not, please email
16	<quantlib-dev@lists.sf.net>. The license is also available online at
17	<https://www.quantlib.org/license.shtml>.
18
19	This program is distributed in the hope that it will be useful, but WITHOUT
20	ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
21	FOR A PARTICULAR PURPOSE. See the license for more details.
22	*/
23
24	#include <ql/math/abcdmathfunction.hpp>
25	#include <utility>
26
27	namespace QuantLib {
28
29	void AbcdMathFunction::validate(Real a,	3,465 ✔
30	Real b,
31	Real c,
32	Real d) {
33	QL_REQUIRE(c>0, "c (" << c << ") must be positive");	3,465 ✔
34	QL_REQUIRE(d>=0, "d (" << d << ") must be non negative");	3,465 ✔
35	QL_REQUIRE(a+d>=0,	3,465 ✔
36	"a+d (" << a << "+" << d << ") must be non negative");
37
38	if (b>=0.0)	3,465 ✔
39	return;
40
41	// the one and only stationary point...
42	Time zeroFirstDerivative = 1.0/c-a/b;	×
43	if (zeroFirstDerivative>=0.0) {	×
44	// ... is a minimum
45	// must be abcd(zeroFirstDerivative)>=0
46	QL_REQUIRE(b>=-(dc)/std::exp(ca/b-1.0),	×
47	"b (" << b << ") less than " <<
48	-(dc)/std::exp(ca/b-1.0) << ": negative function"
49	" value at stationary point " << zeroFirstDerivative);
50	}
51
52	}
53
54	void AbcdMathFunction::initialize_() {	3,463 ✔
55	validate(a_, b_, c_, d_);	3,463 ✔
56	da_ = b_ - c_*a_;	3,463 ✔
57	db_ = -c_*b_;	3,463 ✔
58	dabcd_[0]=da_;	3,463 ✔
59	dabcd_[1]=db_;	3,463 ✔
60	dabcd_[2]=c_;	3,463 ✔
61	dabcd_[3]=0.0;	3,463 ✔
62
63	pa_ = -(a_ + b_/c_)/c_;	3,463 ✔
64	pb_ = -b_/c_;	3,463 ✔
65	K_ = 0.0;	3,463 ✔
66
67	dibc_ = b_/c_;	3,463 ✔
68	diacplusbcc_ = a_/c_ + dibc_/c_;	3,463 ✔
69	}	3,463 ✔
70
71	AbcdMathFunction::AbcdMathFunction(Real aa, Real bb, Real cc, Real dd)	3,463 ✔
72	: a_(aa), b_(bb), c_(cc), d_(dd), abcd_(4), dabcd_(4) {	3,463 ✔
73	abcd_[0]=a_;	3,463 ✔
74	abcd_[1]=b_;	3,463 ✔
75	abcd_[2]=c_;	3,463 ✔
76	abcd_[3]=d_;	3,463 ✔
77	initialize_();	3,463 ✔
78	}	3,463 ✔
79
80	AbcdMathFunction::AbcdMathFunction(std::vector<Real> abcd) : abcd_(std::move(abcd)), dabcd_(4) {	×
81	a_=abcd_[0];	×
82	b_=abcd_[1];	×
83	c_=abcd_[2];	×
84	d_=abcd_[3];	×
85	initialize_();	×
86	}	×
87
88	Time AbcdMathFunction::maximumLocation() const {	×
89	if (b_ == 0.0) {	×
90	// f(t) = aexp(-ct) + d: monotonically decreasing (a>0) or increasing (a<0)
91	if (a_ >= 0.0)	×
UNCOV 92	return 0.0;	×
93	else
94	return QL_MAX_REAL;
95	}
96
97	// Stationary point: f'(t) = 0 => t = 1/c - a/b
UNCOV 98	Time t_star = 1.0/c_ - a_/b_;	×
99
100	if (b_ > 0.0) {	×
101	// f''(t) = -cbexp(-ct) < 0 => t is a maximum.
102	// If t* <= 0 the function is decreasing on [0, inf), maximum is at t=0.
103	return (t_star > 0.0 ? t_star : 0.0);	×
104	} else {
105	// b_ < 0: f''(t) > 0 => t is a minimum, not a maximum.
UNCOV 106	if (t_star <= 0.0) {	×
107	// No interior minimum on [0, inf): function is monotonically
108	// increasing toward d_; maximum is at t -> +inf.
109	return QL_MAX_REAL;
110	} else {
111	// Interior minimum at t* > 0: function decreases then increases.
112	// Maximum is at one of the two boundaries: t=0 (value a_+d_) or
113	// t -> +inf (value d_). validate() guarantees a_+d_ >= 0 and d_ >= 0.
114	return (a_ >= 0.0 ? 0.0 : QL_MAX_REAL);	×
115	}
116	}
117	}
118
UNCOV 119	Real AbcdMathFunction::definiteIntegral(Time t1,	×
120	Time t2) const {
UNCOV 121	return primitive(t2)-primitive(t1);	×
122	}
123
124	std::vector<Real>
125	AbcdMathFunction::definiteIntegralCoefficients(Time t,	×
126	Time t2) const {
127	Time dt = t2 - t;	×
128	Real expcdt = std::exp(-c_*dt);	×
129	std::vector<Real> result(4);	×
130	result[0] = diacplusbcc_ - (diacplusbcc_ + dibc_dt)expcdt;	×
131	result[1] = dibc_ * (1.0 - expcdt);	×
132	result[2] = c_;	×
UNCOV 133	result[3] = d_*dt;	×
UNCOV 134	return result;	×
135	}
136
137	std::vector<Real>
UNCOV 138	AbcdMathFunction::definiteDerivativeCoefficients(Time t,	×
139	Time t2) const {
UNCOV 140	Time dt = t2 - t;	×
UNCOV 141	Real expcdt = std::exp(-c_*dt);	×
UNCOV 142	std::vector<Real> result(4);	×
UNCOV 143	result[1] = b_*c_/(1.0-expcdt);	×
UNCOV 144	result[0] = a_c_ - b_ + result[1]dt*expcdt;	×
UNCOV 145	result[0] /= 1.0-expcdt;	×
UNCOV 146	result[2] = c_;	×
UNCOV 147	result[3] = d_/dt;	×
UNCOV 148	return result;	×
149	}
150
151	}

lballabio / QuantLib / 39944165411

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous