mirror of
https://github.com/boostorg/math.git
synced 2026-01-19 04:22:09 +00:00
Add make_ftuple(), digamma(), lgamma(), tgamma(), doc/test updates. (#218)
Improve tests and coverage. C++11/14 support. (@kedarbhat)
This commit is contained in:
pulver
GitHub
@@ -139,6 +139,7 @@ test-suite examples :
|
||||
[ run cubic_b_spline_example.cpp : : : [ requires cxx11_smart_ptr cxx11_hdr_random cxx11_defaulted_functions ] ]
|
||||
[ compile naive_monte_carlo_example.cpp : [ requires cxx11_auto_declarations cxx11_lambdas cxx11_unified_initialization_syntax cxx11_hdr_thread cxx11_hdr_atomic cxx11_decltype cxx11_hdr_future cxx11_hdr_chrono cxx11_hdr_random cxx11_allocator ] ] # requires user input, can't run it, take a long time too!
|
||||
[ run catmull_rom_example.cpp : : : [ requires cxx17_if_constexpr cxx11_auto_declarations cxx17_std_apply ] ] # Actually the C++17 features used is std::size, not if constexpr; looks like there isn't yet a test for it.
|
||||
[ run autodiff_black_scholes_brief.cpp : : : [ requires cxx11_inline_namespaces ] ]
|
||||
[ run autodiff_black_scholes.cpp : : : [ requires cxx11_inline_namespaces ] ]
|
||||
[ run autodiff_fourth_power.cpp : : : [ requires cxx11_inline_namespaces ] ]
|
||||
[ run autodiff_mixed_partials.cpp : : : [ requires cxx11_inline_namespaces ] ]
|
||||
|
||||
@@ -12,118 +12,126 @@ using namespace boost::math::differentiation;
|
||||
// https://en.wikipedia.org/wiki/Greeks_(finance)#Formulas_for_European_option_Greeks
|
||||
|
||||
// Standard normal probability density function
|
||||
template<typename X>
|
||||
X phi(const X& x)
|
||||
{
|
||||
return boost::math::constants::one_div_root_two_pi<X>()*exp(-0.5*x*x);
|
||||
template <typename X>
|
||||
X phi(X const& x) {
|
||||
return boost::math::constants::one_div_root_two_pi<X>() * exp(-0.5 * x * x);
|
||||
}
|
||||
|
||||
// Standard normal cumulative distribution function
|
||||
template<typename X>
|
||||
X Phi(const X& x)
|
||||
{
|
||||
return 0.5*erfc(-boost::math::constants::one_div_root_two<X>()*x);
|
||||
template <typename X>
|
||||
X Phi(X const& x) {
|
||||
return 0.5 * erfc(-boost::math::constants::one_div_root_two<X>() * x);
|
||||
}
|
||||
|
||||
enum CP { call, put };
|
||||
enum class CP { call, put };
|
||||
|
||||
// Assume zero annual dividend yield (q=0).
|
||||
template<typename Price,typename Sigma,typename Tau,typename Rate>
|
||||
promote<Price,Sigma,Tau,Rate>
|
||||
black_scholes_option_price(CP cp, double K, const Price& S, const Sigma& sigma, const Tau& tau, const Rate& r)
|
||||
{
|
||||
template <typename Price, typename Sigma, typename Tau, typename Rate>
|
||||
promote<Price, Sigma, Tau, Rate> black_scholes_option_price(CP cp,
|
||||
double K,
|
||||
Price const& S,
|
||||
Sigma const& sigma,
|
||||
Tau const& tau,
|
||||
Rate const& r) {
|
||||
using namespace std;
|
||||
const auto d1 = (log(S/K) + (r+sigma*sigma/2)*tau) / (sigma*sqrt(tau));
|
||||
const auto d2 = (log(S/K) + (r-sigma*sigma/2)*tau) / (sigma*sqrt(tau));
|
||||
if (cp == call)
|
||||
return S*Phi(d1) - exp(-r*tau)*K*Phi(d2);
|
||||
else
|
||||
return exp(-r*tau)*K*Phi(-d2) - S*Phi(-d1);
|
||||
auto const d1 = (log(S / K) + (r + sigma * sigma / 2) * tau) / (sigma * sqrt(tau));
|
||||
auto const d2 = (log(S / K) + (r - sigma * sigma / 2) * tau) / (sigma * sqrt(tau));
|
||||
switch (cp) {
|
||||
case CP::call:
|
||||
return S * Phi(d1) - exp(-r * tau) * K * Phi(d2);
|
||||
case CP::put:
|
||||
return exp(-r * tau) * K * Phi(-d2) - S * Phi(-d1);
|
||||
}
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
const double K = 100.0; // Strike price.
|
||||
const auto S = make_fvar<double,3>(105); // Stock price.
|
||||
const auto sigma = make_fvar<double,0,3>(5); // Volatility.
|
||||
const auto tau = make_fvar<double,0,0,1>(30.0/365); // Time to expiration in years. (30 days).
|
||||
const auto r = make_fvar<double,0,0,0,1>(1.25/100); // Interest rate.
|
||||
const auto call_price = black_scholes_option_price(call, K, S, sigma, tau, r);
|
||||
const auto put_price = black_scholes_option_price(put, K, S, sigma, tau, r);
|
||||
int main() {
|
||||
double const K = 100.0; // Strike price.
|
||||
auto const variables = make_ftuple<double, 3, 3, 1, 1>(105, 5, 30.0 / 365, 1.25 / 100);
|
||||
auto const& S = std::get<0>(variables); // Stock price.
|
||||
auto const& sigma = std::get<1>(variables); // Volatility.
|
||||
auto const& tau = std::get<2>(variables); // Time to expiration in years. (30 days).
|
||||
auto const& r = std::get<3>(variables); // Interest rate.
|
||||
auto const call_price = black_scholes_option_price(CP::call, K, S, sigma, tau, r);
|
||||
auto const put_price = black_scholes_option_price(CP::put, K, S, sigma, tau, r);
|
||||
|
||||
// Compare automatically calculated greeks by autodiff with formulas for greeks.
|
||||
// https://en.wikipedia.org/wiki/Greeks_(finance)#Formulas_for_European_option_Greeks
|
||||
const double d1 = static_cast<double>((log(S/K) + (r+sigma*sigma/2)*tau) / (sigma*sqrt(tau)));
|
||||
const double d2 = static_cast<double>((log(S/K) + (r-sigma*sigma/2)*tau) / (sigma*sqrt(tau)));
|
||||
const double formula_call_delta = +Phi(+d1);
|
||||
const double formula_put_delta = -Phi(-d1);
|
||||
const double formula_vega = static_cast<double>(S*phi(d1)*sqrt(tau));
|
||||
const double formula_call_theta = static_cast<double>(-S*phi(d1)*sigma/(2*sqrt(tau))-r*K*exp(-r*tau)*Phi(+d2));
|
||||
const double formula_put_theta = static_cast<double>(-S*phi(d1)*sigma/(2*sqrt(tau))+r*K*exp(-r*tau)*Phi(-d2));
|
||||
const double formula_call_rho = static_cast<double>(+K*tau*exp(-r*tau)*Phi(+d2));
|
||||
const double formula_put_rho = static_cast<double>(-K*tau*exp(-r*tau)*Phi(-d2));
|
||||
const double formula_gamma = static_cast<double>(phi(d1)/(S*sigma*sqrt(tau)));
|
||||
const double formula_vanna = static_cast<double>(-phi(d1)*d2/sigma);
|
||||
const double formula_charm = static_cast<double>(phi(d1)*(d2*sigma*sqrt(tau)-2*r*tau)/(2*tau*sigma*sqrt(tau)));
|
||||
const double formula_vomma = static_cast<double>(S*phi(d1)*sqrt(tau)*d1*d2/sigma);
|
||||
const double formula_veta = static_cast<double>(-S*phi(d1)*sqrt(tau)*(r*d1/(sigma*sqrt(tau))-(1+d1*d2)/(2*tau)));
|
||||
const double formula_speed = static_cast<double>(-phi(d1)*(d1/(sigma*sqrt(tau))+1)/(S*S*sigma*sqrt(tau)));
|
||||
const double formula_zomma = static_cast<double>(phi(d1)*(d1*d2-1)/(S*sigma*sigma*sqrt(tau)));
|
||||
const double formula_color =
|
||||
static_cast<double>(-phi(d1)/(2*S*tau*sigma*sqrt(tau))*(1+(2*r*tau-d2*sigma*sqrt(tau))*d1/(sigma*sqrt(tau))));
|
||||
const double formula_ultima = -formula_vega*static_cast<double>((d1*d2*(1-d1*d2)+d1*d1+d2*d2)/(sigma*sigma));
|
||||
double const d1 = static_cast<double>((log(S / K) + (r + sigma * sigma / 2) * tau) / (sigma * sqrt(tau)));
|
||||
double const d2 = static_cast<double>((log(S / K) + (r - sigma * sigma / 2) * tau) / (sigma * sqrt(tau)));
|
||||
double const formula_call_delta = +Phi(+d1);
|
||||
double const formula_put_delta = -Phi(-d1);
|
||||
double const formula_vega = static_cast<double>(S * phi(d1) * sqrt(tau));
|
||||
double const formula_call_theta =
|
||||
static_cast<double>(-S * phi(d1) * sigma / (2 * sqrt(tau)) - r * K * exp(-r * tau) * Phi(+d2));
|
||||
double const formula_put_theta =
|
||||
static_cast<double>(-S * phi(d1) * sigma / (2 * sqrt(tau)) + r * K * exp(-r * tau) * Phi(-d2));
|
||||
double const formula_call_rho = static_cast<double>(+K * tau * exp(-r * tau) * Phi(+d2));
|
||||
double const formula_put_rho = static_cast<double>(-K * tau * exp(-r * tau) * Phi(-d2));
|
||||
double const formula_gamma = static_cast<double>(phi(d1) / (S * sigma * sqrt(tau)));
|
||||
double const formula_vanna = static_cast<double>(-phi(d1) * d2 / sigma);
|
||||
double const formula_charm =
|
||||
static_cast<double>(phi(d1) * (d2 * sigma * sqrt(tau) - 2 * r * tau) / (2 * tau * sigma * sqrt(tau)));
|
||||
double const formula_vomma = static_cast<double>(S * phi(d1) * sqrt(tau) * d1 * d2 / sigma);
|
||||
double const formula_veta = static_cast<double>(-S * phi(d1) * sqrt(tau) *
|
||||
(r * d1 / (sigma * sqrt(tau)) - (1 + d1 * d2) / (2 * tau)));
|
||||
double const formula_speed =
|
||||
static_cast<double>(-phi(d1) * (d1 / (sigma * sqrt(tau)) + 1) / (S * S * sigma * sqrt(tau)));
|
||||
double const formula_zomma = static_cast<double>(phi(d1) * (d1 * d2 - 1) / (S * sigma * sigma * sqrt(tau)));
|
||||
double const formula_color =
|
||||
static_cast<double>(-phi(d1) / (2 * S * tau * sigma * sqrt(tau)) *
|
||||
(1 + (2 * r * tau - d2 * sigma * sqrt(tau)) * d1 / (sigma * sqrt(tau))));
|
||||
double const formula_ultima =
|
||||
-formula_vega * static_cast<double>((d1 * d2 * (1 - d1 * d2) + d1 * d1 + d2 * d2) / (sigma * sigma));
|
||||
|
||||
std::cout << std::setprecision(std::numeric_limits<double>::digits10)
|
||||
<< "autodiff black-scholes call price = " << call_price.derivative(0,0,0,0) << '\n'
|
||||
<< "autodiff black-scholes put price = " << put_price.derivative(0,0,0,0) << '\n'
|
||||
<< "\n## First-order Greeks\n"
|
||||
<< "autodiff call delta = " << call_price.derivative(1,0,0,0) << '\n'
|
||||
<< " formula call delta = " << formula_call_delta << '\n'
|
||||
<< "autodiff call vega = " << call_price.derivative(0,1,0,0) << '\n'
|
||||
<< " formula call vega = " << formula_vega << '\n'
|
||||
<< "autodiff call theta = " << -call_price.derivative(0,0,1,0) << '\n' // minus sign due to tau = T-time
|
||||
<< " formula call theta = " << formula_call_theta << '\n'
|
||||
<< "autodiff call rho = " << call_price.derivative(0,0,0,1) << '\n'
|
||||
<< " formula call rho = " << formula_call_rho << '\n'
|
||||
<< '\n'
|
||||
<< "autodiff put delta = " << put_price.derivative(1,0,0,0) << '\n'
|
||||
<< " formula put delta = " << formula_put_delta << '\n'
|
||||
<< "autodiff put vega = " << put_price.derivative(0,1,0,0) << '\n'
|
||||
<< " formula put vega = " << formula_vega << '\n'
|
||||
<< "autodiff put theta = " << -put_price.derivative(0,0,1,0) << '\n'
|
||||
<< " formula put theta = " << formula_put_theta << '\n'
|
||||
<< "autodiff put rho = " << put_price.derivative(0,0,0,1) << '\n'
|
||||
<< " formula put rho = " << formula_put_rho << '\n'
|
||||
<< "\n## Second-order Greeks\n"
|
||||
<< "autodiff call gamma = " << call_price.derivative(2,0,0,0) << '\n'
|
||||
<< "autodiff put gamma = " << put_price.derivative(2,0,0,0) << '\n'
|
||||
<< " formula gamma = " << formula_gamma << '\n'
|
||||
<< "autodiff call vanna = " << call_price.derivative(1,1,0,0) << '\n'
|
||||
<< "autodiff put vanna = " << put_price.derivative(1,1,0,0) << '\n'
|
||||
<< " formula vanna = " << formula_vanna << '\n'
|
||||
<< "autodiff call charm = " << -call_price.derivative(1,0,1,0) << '\n'
|
||||
<< "autodiff put charm = " << -put_price.derivative(1,0,1,0) << '\n'
|
||||
<< " formula charm = " << formula_charm << '\n'
|
||||
<< "autodiff call vomma = " << call_price.derivative(0,2,0,0) << '\n'
|
||||
<< "autodiff put vomma = " << put_price.derivative(0,2,0,0) << '\n'
|
||||
<< " formula vomma = " << formula_vomma << '\n'
|
||||
<< "autodiff call veta = " << call_price.derivative(0,1,1,0) << '\n'
|
||||
<< "autodiff put veta = " << put_price.derivative(0,1,1,0) << '\n'
|
||||
<< " formula veta = " << formula_veta << '\n'
|
||||
<< "\n## Third-order Greeks\n"
|
||||
<< "autodiff call speed = " << call_price.derivative(3,0,0,0) << '\n'
|
||||
<< "autodiff put speed = " << put_price.derivative(3,0,0,0) << '\n'
|
||||
<< " formula speed = " << formula_speed << '\n'
|
||||
<< "autodiff call zomma = " << call_price.derivative(2,1,0,0) << '\n'
|
||||
<< "autodiff put zomma = " << put_price.derivative(2,1,0,0) << '\n'
|
||||
<< " formula zomma = " << formula_zomma << '\n'
|
||||
<< "autodiff call color = " << call_price.derivative(2,0,1,0) << '\n'
|
||||
<< "autodiff put color = " << put_price.derivative(2,0,1,0) << '\n'
|
||||
<< " formula color = " << formula_color << '\n'
|
||||
<< "autodiff call ultima = " << call_price.derivative(0,3,0,0) << '\n'
|
||||
<< "autodiff put ultima = " << put_price.derivative(0,3,0,0) << '\n'
|
||||
<< " formula ultima = " << formula_ultima << '\n'
|
||||
;
|
||||
<< "autodiff black-scholes call price = " << call_price.derivative(0, 0, 0, 0) << '\n'
|
||||
<< "autodiff black-scholes put price = " << put_price.derivative(0, 0, 0, 0) << '\n'
|
||||
<< "\n## First-order Greeks\n"
|
||||
<< "autodiff call delta = " << call_price.derivative(1, 0, 0, 0) << '\n'
|
||||
<< " formula call delta = " << formula_call_delta << '\n'
|
||||
<< "autodiff call vega = " << call_price.derivative(0, 1, 0, 0) << '\n'
|
||||
<< " formula call vega = " << formula_vega << '\n'
|
||||
<< "autodiff call theta = " << -call_price.derivative(0, 0, 1, 0)
|
||||
<< '\n' // minus sign due to tau = T-time
|
||||
<< " formula call theta = " << formula_call_theta << '\n'
|
||||
<< "autodiff call rho = " << call_price.derivative(0, 0, 0, 1) << '\n'
|
||||
<< " formula call rho = " << formula_call_rho << '\n'
|
||||
<< '\n'
|
||||
<< "autodiff put delta = " << put_price.derivative(1, 0, 0, 0) << '\n'
|
||||
<< " formula put delta = " << formula_put_delta << '\n'
|
||||
<< "autodiff put vega = " << put_price.derivative(0, 1, 0, 0) << '\n'
|
||||
<< " formula put vega = " << formula_vega << '\n'
|
||||
<< "autodiff put theta = " << -put_price.derivative(0, 0, 1, 0) << '\n'
|
||||
<< " formula put theta = " << formula_put_theta << '\n'
|
||||
<< "autodiff put rho = " << put_price.derivative(0, 0, 0, 1) << '\n'
|
||||
<< " formula put rho = " << formula_put_rho << '\n'
|
||||
<< "\n## Second-order Greeks\n"
|
||||
<< "autodiff call gamma = " << call_price.derivative(2, 0, 0, 0) << '\n'
|
||||
<< "autodiff put gamma = " << put_price.derivative(2, 0, 0, 0) << '\n'
|
||||
<< " formula gamma = " << formula_gamma << '\n'
|
||||
<< "autodiff call vanna = " << call_price.derivative(1, 1, 0, 0) << '\n'
|
||||
<< "autodiff put vanna = " << put_price.derivative(1, 1, 0, 0) << '\n'
|
||||
<< " formula vanna = " << formula_vanna << '\n'
|
||||
<< "autodiff call charm = " << -call_price.derivative(1, 0, 1, 0) << '\n'
|
||||
<< "autodiff put charm = " << -put_price.derivative(1, 0, 1, 0) << '\n'
|
||||
<< " formula charm = " << formula_charm << '\n'
|
||||
<< "autodiff call vomma = " << call_price.derivative(0, 2, 0, 0) << '\n'
|
||||
<< "autodiff put vomma = " << put_price.derivative(0, 2, 0, 0) << '\n'
|
||||
<< " formula vomma = " << formula_vomma << '\n'
|
||||
<< "autodiff call veta = " << call_price.derivative(0, 1, 1, 0) << '\n'
|
||||
<< "autodiff put veta = " << put_price.derivative(0, 1, 1, 0) << '\n'
|
||||
<< " formula veta = " << formula_veta << '\n'
|
||||
<< "\n## Third-order Greeks\n"
|
||||
<< "autodiff call speed = " << call_price.derivative(3, 0, 0, 0) << '\n'
|
||||
<< "autodiff put speed = " << put_price.derivative(3, 0, 0, 0) << '\n'
|
||||
<< " formula speed = " << formula_speed << '\n'
|
||||
<< "autodiff call zomma = " << call_price.derivative(2, 1, 0, 0) << '\n'
|
||||
<< "autodiff put zomma = " << put_price.derivative(2, 1, 0, 0) << '\n'
|
||||
<< " formula zomma = " << formula_zomma << '\n'
|
||||
<< "autodiff call color = " << call_price.derivative(2, 0, 1, 0) << '\n'
|
||||
<< "autodiff put color = " << put_price.derivative(2, 0, 1, 0) << '\n'
|
||||
<< " formula color = " << formula_color << '\n'
|
||||
<< "autodiff call ultima = " << call_price.derivative(0, 3, 0, 0) << '\n'
|
||||
<< "autodiff put ultima = " << put_price.derivative(0, 3, 0, 0) << '\n'
|
||||
<< " formula ultima = " << formula_ultima << '\n';
|
||||
return 0;
|
||||
}
|
||||
/*
|
||||
|
||||
67
example/autodiff_black_scholes_brief.cpp
Normal file
67
example/autodiff_black_scholes_brief.cpp
Normal file
@@ -0,0 +1,67 @@
|
||||
// Copyright Matthew Pulver 2018 - 2019.
|
||||
// Distributed under the Boost Software License, Version 1.0.
|
||||
// (See accompanying file LICENSE_1_0.txt or copy at
|
||||
// https://www.boost.org/LICENSE_1_0.txt)
|
||||
|
||||
#include <boost/math/differentiation/autodiff.hpp>
|
||||
#include <iostream>
|
||||
|
||||
using namespace boost::math::constants;
|
||||
using namespace boost::math::differentiation;
|
||||
|
||||
// Equations and function/variable names are from
|
||||
// https://en.wikipedia.org/wiki/Greeks_(finance)#Formulas_for_European_option_Greeks
|
||||
|
||||
// Standard normal cumulative distribution function
|
||||
template <typename X>
|
||||
X Phi(X const& x) {
|
||||
return 0.5 * erfc(-one_div_root_two<X>() * x);
|
||||
}
|
||||
|
||||
enum class CP { call, put };
|
||||
|
||||
// Assume zero annual dividend yield (q=0).
|
||||
template <typename Price, typename Sigma, typename Tau, typename Rate>
|
||||
promote<Price, Sigma, Tau, Rate> black_scholes_option_price(CP cp,
|
||||
double K,
|
||||
Price const& S,
|
||||
Sigma const& sigma,
|
||||
Tau const& tau,
|
||||
Rate const& r) {
|
||||
using namespace std;
|
||||
auto const d1 = (log(S / K) + (r + sigma * sigma / 2) * tau) / (sigma * sqrt(tau));
|
||||
auto const d2 = (log(S / K) + (r - sigma * sigma / 2) * tau) / (sigma * sqrt(tau));
|
||||
switch (cp) {
|
||||
case CP::call:
|
||||
return S * Phi(d1) - exp(-r * tau) * K * Phi(d2);
|
||||
case CP::put:
|
||||
return exp(-r * tau) * K * Phi(-d2) - S * Phi(-d1);
|
||||
}
|
||||
}
|
||||
|
||||
int main() {
|
||||
double const K = 100.0; // Strike price.
|
||||
auto const S = make_fvar<double, 2>(105); // Stock price.
|
||||
double const sigma = 5; // Volatility.
|
||||
double const tau = 30.0 / 365; // Time to expiration in years. (30 days).
|
||||
double const r = 1.25 / 100; // Interest rate.
|
||||
auto const call_price = black_scholes_option_price(CP::call, K, S, sigma, tau, r);
|
||||
auto const put_price = black_scholes_option_price(CP::put, K, S, sigma, tau, r);
|
||||
|
||||
std::cout << "black-scholes call price = " << call_price.derivative(0) << '\n'
|
||||
<< "black-scholes put price = " << put_price.derivative(0) << '\n'
|
||||
<< "call delta = " << call_price.derivative(1) << '\n'
|
||||
<< "put delta = " << put_price.derivative(1) << '\n'
|
||||
<< "call gamma = " << call_price.derivative(2) << '\n'
|
||||
<< "put gamma = " << put_price.derivative(2) << '\n';
|
||||
return 0;
|
||||
}
|
||||
/*
|
||||
Output:
|
||||
black-scholes call price = 56.5136
|
||||
black-scholes put price = 51.4109
|
||||
call delta = 0.773818
|
||||
put delta = -0.226182
|
||||
call gamma = 0.00199852
|
||||
put gamma = 0.00199852
|
||||
**/
|
||||
@@ -6,23 +6,22 @@
|
||||
#include <boost/math/differentiation/autodiff.hpp>
|
||||
#include <iostream>
|
||||
|
||||
template<typename T>
|
||||
T fourth_power(T x)
|
||||
{
|
||||
x *= x;
|
||||
return x *= x;
|
||||
template <typename T>
|
||||
T fourth_power(T const& x) {
|
||||
T x4 = x * x; // retval in operator*() uses x4's memory via NRVO.
|
||||
x4 *= x4; // No copies of x4 are made within operator*=() even when squaring.
|
||||
return x4; // x4 uses y's memory in main() via NRVO.
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
using namespace boost::math::differentiation;
|
||||
int main() {
|
||||
using namespace boost::math::differentiation;
|
||||
|
||||
constexpr int Order=5; // The highest order derivative to be calculated.
|
||||
const autodiff_fvar<double,Order> x = make_fvar<double,Order>(2.0); // Find derivatives at x=2.
|
||||
const autodiff_fvar<double,Order> y = fourth_power(x);
|
||||
for (int i=0 ; i<=Order ; ++i)
|
||||
std::cout << "y.derivative("<<i<<") = " << y.derivative(i) << std::endl;
|
||||
return 0;
|
||||
constexpr unsigned Order = 5; // Highest order derivative to be calculated.
|
||||
auto const x = make_fvar<double, Order>(2.0); // Find derivatives at x=2.
|
||||
auto const y = fourth_power(x);
|
||||
for (unsigned i = 0; i <= Order; ++i)
|
||||
std::cout << "y.derivative(" << i << ") = " << y.derivative(i) << std::endl;
|
||||
return 0;
|
||||
}
|
||||
/*
|
||||
Output:
|
||||
|
||||
File diff suppressed because one or more lines are too long
@@ -7,35 +7,35 @@
|
||||
#include <boost/multiprecision/cpp_bin_float.hpp>
|
||||
#include <iostream>
|
||||
|
||||
template<typename T>
|
||||
T f(const T& w, const T& x, const T& y, const T& z)
|
||||
{
|
||||
using namespace boost::math::differentiation;
|
||||
|
||||
template <typename W, typename X, typename Y, typename Z>
|
||||
promote<W, X, Y, Z> f(const W& w, const X& x, const Y& y, const Z& z) {
|
||||
using namespace std;
|
||||
return exp(w*sin(x*log(y)/z) + sqrt(w*z/(x*y))) + w*w/tan(z);
|
||||
return exp(w * sin(x * log(y) / z) + sqrt(w * z / (x * y))) + w * w / tan(z);
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
using cpp_bin_float_50 = boost::multiprecision::cpp_bin_float_50;
|
||||
using namespace boost::math::differentiation;
|
||||
int main() {
|
||||
using float50 = boost::multiprecision::cpp_bin_float_50;
|
||||
|
||||
constexpr int Nw=3; // Max order of derivative to calculate for w
|
||||
constexpr int Nx=2; // Max order of derivative to calculate for x
|
||||
constexpr int Ny=4; // Max order of derivative to calculate for y
|
||||
constexpr int Nz=3; // Max order of derivative to calculate for z
|
||||
using var = autodiff_fvar<cpp_bin_float_50,Nw,Nx,Ny,Nz>;
|
||||
const var w = make_fvar<cpp_bin_float_50,Nw>(11);
|
||||
const var x = make_fvar<cpp_bin_float_50,0,Nx>(12);
|
||||
const var y = make_fvar<cpp_bin_float_50,0,0,Ny>(13);
|
||||
const var z = make_fvar<cpp_bin_float_50,0,0,0,Nz>(14);
|
||||
const var v = f(w,x,y,z);
|
||||
constexpr unsigned Nw = 3; // Max order of derivative to calculate for w
|
||||
constexpr unsigned Nx = 2; // Max order of derivative to calculate for x
|
||||
constexpr unsigned Ny = 4; // Max order of derivative to calculate for y
|
||||
constexpr unsigned Nz = 3; // Max order of derivative to calculate for z
|
||||
// Declare 4 independent variables together into a std::tuple.
|
||||
auto const variables = make_ftuple<float50, Nw, Nx, Ny, Nz>(11, 12, 13, 14);
|
||||
auto const& w = std::get<0>(variables); // Up to Nw derivatives at w=11
|
||||
auto const& x = std::get<1>(variables); // Up to Nx derivatives at x=12
|
||||
auto const& y = std::get<2>(variables); // Up to Ny derivatives at y=13
|
||||
auto const& z = std::get<3>(variables); // Up to Nz derivatives at z=14
|
||||
auto const v = f(w, x, y, z);
|
||||
// Calculated from Mathematica symbolic differentiation.
|
||||
const cpp_bin_float_50 answer("1976.319600747797717779881875290418720908121189218755");
|
||||
std::cout << std::setprecision(std::numeric_limits<cpp_bin_float_50>::digits10)
|
||||
<< "mathematica : " << answer << '\n'
|
||||
<< "autodiff : " << v.derivative(Nw,Nx,Ny,Nz) << '\n'
|
||||
<< "relative error: " << std::setprecision(3) << (v.derivative(Nw,Nx,Ny,Nz)/answer-1)
|
||||
<< std::endl;
|
||||
float50 const answer("1976.319600747797717779881875290418720908121189218755");
|
||||
std::cout << std::setprecision(std::numeric_limits<float50>::digits10)
|
||||
<< "mathematica : " << answer << '\n'
|
||||
<< "autodiff : " << v.derivative(Nw, Nx, Ny, Nz) << '\n'
|
||||
<< std::setprecision(3)
|
||||
<< "relative error: " << (v.derivative(Nw, Nx, Ny, Nz) / answer - 1) << '\n';
|
||||
return 0;
|
||||
}
|
||||
/*
|
||||
|
||||
Reference in New Issue
Block a user