ceph/src/boost/libs/math/example/autodiff_multiprecision.cpp

   1 //           Copyright Matthew Pulver 2018 - 2019.
   2 // Distributed under the Boost Software License, Version 1.0.
   3 //      (See accompanying file LICENSE_1_0.txt or copy at
   4 //           https://www.boost.org/LICENSE_1_0.txt)
   5
   6 #include <boost/math/differentiation/autodiff.hpp>
   7 #include <boost/multiprecision/cpp_bin_float.hpp>
   8 #include <iostream>
   9
  10 using namespace boost::math::differentiation;
  11
  12 template <typename W, typename X, typename Y, typename Z>
  13 promote<W, X, Y, Z> f(const W& w, const X& x, const Y& y, const Z& z) {
  14   using namespace std;
  15   return exp(w * sin(x * log(y) / z) + sqrt(w * z / (x * y))) + w * w / tan(z);
  16 }
  17
  18 int main() {
  19   using float50 = boost::multiprecision::cpp_bin_float_50;
  20
  21   constexpr unsigned Nw = 3;  // Max order of derivative to calculate for w
  22   constexpr unsigned Nx = 2;  // Max order of derivative to calculate for x
  23   constexpr unsigned Ny = 4;  // Max order of derivative to calculate for y
  24   constexpr unsigned Nz = 3;  // Max order of derivative to calculate for z
  25   // Declare 4 independent variables together into a std::tuple.
  26   auto const variables = make_ftuple<float50, Nw, Nx, Ny, Nz>(11, 12, 13, 14);
  27   auto const& w = std::get<0>(variables);  // Up to Nw derivatives at w=11
  28   auto const& x = std::get<1>(variables);  // Up to Nx derivatives at x=12
  29   auto const& y = std::get<2>(variables);  // Up to Ny derivatives at y=13
  30   auto const& z = std::get<3>(variables);  // Up to Nz derivatives at z=14
  31   auto const v = f(w, x, y, z);
  32   // Calculated from Mathematica symbolic differentiation.
  33   float50 const answer("1976.319600747797717779881875290418720908121189218755");
  34   std::cout << std::setprecision(std::numeric_limits<float50>::digits10)
  35             << "mathematica   : " << answer << '\n'
  36             << "autodiff      : " << v.derivative(Nw, Nx, Ny, Nz) << '\n'
  37             << std::setprecision(3)
  38             << "relative error: " << (v.derivative(Nw, Nx, Ny, Nz) / answer - 1) << '\n';
  39   return 0;
  40 }
  41 /*
  42 Output:
  43 mathematica   : 1976.3196007477977177798818752904187209081211892188
  44 autodiff      : 1976.3196007477977177798818752904187209081211892188
  45 relative error: 2.67e-50
  46 **/