boost/math/differentiation/autodiff.hpp

   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 #ifndef BOOST_MATH_DIFFERENTIATION_AUTODIFF_HPP
   7 #define BOOST_MATH_DIFFERENTIATION_AUTODIFF_HPP
   8
   9 #include <boost/cstdfloat.hpp>
  10 #include <boost/math/constants/constants.hpp>
  11 #include <boost/math/special_functions.hpp>
  12 #include <boost/math/tools/config.hpp>
  13 #include <boost/math/tools/promotion.hpp>
  14
  15 #include <algorithm>
  16 #include <array>
  17 #include <cmath>
  18 #include <functional>
  19 #include <limits>
  20 #include <numeric>
  21 #include <ostream>
  22 #include <tuple>
  23 #include <type_traits>
  24
  25 namespace boost {
  26 namespace math {
  27 namespace differentiation {
  28 // Automatic Differentiation v1
  29 inline namespace autodiff_v1 {
  30 namespace detail {
  31
  32 template <typename RealType, typename... RealTypes>
  33 struct promote_args_n {
  34   using type = typename tools::promote_args_2<RealType, typename promote_args_n<RealTypes...>::type>::type;
  35 };
  36
  37 template <typename RealType>
  38 struct promote_args_n<RealType> {
  39   using type = typename tools::promote_arg<RealType>::type;
  40 };
  41
  42 }  // namespace detail
  43
  44 template <typename RealType, typename... RealTypes>
  45 using promote = typename detail::promote_args_n<RealType, RealTypes...>::type;
  46
  47 namespace detail {
  48
  49 template <typename RealType, size_t Order>
  50 class fvar;
  51
  52 template <typename T>
  53 struct is_fvar_impl : std::false_type {};
  54
  55 template <typename RealType, size_t Order>
  56 struct is_fvar_impl<fvar<RealType, Order>> : std::true_type {};
  57
  58 template <typename T>
  59 using is_fvar = is_fvar_impl<decay_t<T>>;
  60
  61 template <typename RealType, size_t Order, size_t... Orders>
  62 struct nest_fvar {
  63   using type = fvar<typename nest_fvar<RealType, Orders...>::type, Order>;
  64 };
  65
  66 template <typename RealType, size_t Order>
  67 struct nest_fvar<RealType, Order> {
  68   using type = fvar<RealType, Order>;
  69 };
  70
  71 template <typename>
  72 struct get_depth_impl : std::integral_constant<size_t, 0> {};
  73
  74 template <typename RealType, size_t Order>
  75 struct get_depth_impl<fvar<RealType, Order>>
  76     : std::integral_constant<size_t, get_depth_impl<RealType>::value + 1> {};
  77
  78 template <typename T>
  79 using get_depth = get_depth_impl<decay_t<T>>;
  80
  81 template <typename>
  82 struct get_order_sum_t : std::integral_constant<size_t, 0> {};
  83
  84 template <typename RealType, size_t Order>
  85 struct get_order_sum_t<fvar<RealType, Order>>
  86     : std::integral_constant<size_t, get_order_sum_t<RealType>::value + Order> {};
  87
  88 template <typename T>
  89 using get_order_sum = get_order_sum_t<decay_t<T>>;
  90
  91 template <typename RealType>
  92 struct get_root_type {
  93   using type = RealType;
  94 };
  95
  96 template <typename RealType, size_t Order>
  97 struct get_root_type<fvar<RealType, Order>> {
  98   using type = typename get_root_type<RealType>::type;
  99 };
 100
 101 template <typename RealType, size_t Depth>
 102 struct type_at {
 103   using type = RealType;
 104 };
 105
 106 template <typename RealType, size_t Order, size_t Depth>
 107 struct type_at<fvar<RealType, Order>, Depth> {
 108   using type = typename conditional<Depth == 0,
 109                                     fvar<RealType, Order>,
 110                                     typename type_at<RealType, Depth - 1>::type>::type;
 111 };
 112
 113 template <typename RealType, size_t Depth>
 114 using get_type_at = typename type_at<RealType, Depth>::type;
 115
 116 // Satisfies Boost's Conceptual Requirements for Real Number Types.
 117 // https://www.boost.org/libs/math/doc/html/math_toolkit/real_concepts.html
 118 template <typename RealType, size_t Order>
 119 class fvar {
 120   std::array<RealType, Order + 1> v;
 121
 122  public:
 123   using root_type = typename get_root_type<RealType>::type;  // RealType in the root fvar<RealType,Order>.
 124
 125   fvar() = default;
 126
 127   // Initialize a variable or constant.
 128   fvar(root_type const&, bool const is_variable);
 129
 130   // RealType(cr) | RealType | RealType is copy constructible.
 131   fvar(fvar const&) = default;
 132
 133   // Be aware of implicit casting from one fvar<> type to another by this copy constructor.
 134   template <typename RealType2, size_t Order2>
 135   fvar(fvar<RealType2, Order2> const&);
 136
 137   // RealType(ca) | RealType | RealType is copy constructible from the arithmetic types.
 138   explicit fvar(root_type const&);  // Initialize a constant. (No epsilon terms.)
 139
 140   template <typename RealType2>
 141   fvar(RealType2 const& ca);  // Supports any RealType2 for which static_cast<root_type>(ca) compiles.
 142
 143   // r = cr | RealType& | Assignment operator.
 144   fvar& operator=(fvar const&) = default;
 145
 146   // r = ca | RealType& | Assignment operator from the arithmetic types.
 147   // Handled by constructor that takes a single parameter of generic type.
 148   // fvar& operator=(root_type const&); // Set a constant.
 149
 150   // r += cr | RealType& | Adds cr to r.
 151   template <typename RealType2, size_t Order2>
 152   fvar& operator+=(fvar<RealType2, Order2> const&);
 153
 154   // r += ca | RealType& | Adds ar to r.
 155   fvar& operator+=(root_type const&);
 156
 157   // r -= cr | RealType& | Subtracts cr from r.
 158   template <typename RealType2, size_t Order2>
 159   fvar& operator-=(fvar<RealType2, Order2> const&);
 160
 161   // r -= ca | RealType& | Subtracts ca from r.
 162   fvar& operator-=(root_type const&);
 163
 164   // r *= cr | RealType& | Multiplies r by cr.
 165   template <typename RealType2, size_t Order2>
 166   fvar& operator*=(fvar<RealType2, Order2> const&);
 167
 168   // r *= ca | RealType& | Multiplies r by ca.
 169   fvar& operator*=(root_type const&);
 170
 171   // r /= cr | RealType& | Divides r by cr.
 172   template <typename RealType2, size_t Order2>
 173   fvar& operator/=(fvar<RealType2, Order2> const&);
 174
 175   // r /= ca | RealType& | Divides r by ca.
 176   fvar& operator/=(root_type const&);
 177
 178   // -r | RealType | Unary Negation.
 179   fvar operator-() const;
 180
 181   // +r | RealType& | Identity Operation.
 182   fvar const& operator+() const;
 183
 184   // cr + cr2 | RealType | Binary Addition
 185   template <typename RealType2, size_t Order2>
 186   promote<fvar, fvar<RealType2, Order2>> operator+(fvar<RealType2, Order2> const&) const;
 187
 188   // cr + ca | RealType | Binary Addition
 189   fvar operator+(root_type const&) const;
 190
 191   // ca + cr | RealType | Binary Addition
 192   template <typename RealType2, size_t Order2>
 193   friend fvar<RealType2, Order2> operator+(typename fvar<RealType2, Order2>::root_type const&,
 194                                            fvar<RealType2, Order2> const&);
 195
 196   // cr - cr2 | RealType | Binary Subtraction
 197   template <typename RealType2, size_t Order2>
 198   promote<fvar, fvar<RealType2, Order2>> operator-(fvar<RealType2, Order2> const&) const;
 199
 200   // cr - ca | RealType | Binary Subtraction
 201   fvar operator-(root_type const&) const;
 202
 203   // ca - cr | RealType | Binary Subtraction
 204   template <typename RealType2, size_t Order2>
 205   friend fvar<RealType2, Order2> operator-(typename fvar<RealType2, Order2>::root_type const&,
 206                                            fvar<RealType2, Order2> const&);
 207
 208   // cr * cr2 | RealType | Binary Multiplication
 209   template <typename RealType2, size_t Order2>
 210   promote<fvar, fvar<RealType2, Order2>> operator*(fvar<RealType2, Order2> const&)const;
 211
 212   // cr * ca | RealType | Binary Multiplication
 213   fvar operator*(root_type const&)const;
 214
 215   // ca * cr | RealType | Binary Multiplication
 216   template <typename RealType2, size_t Order2>
 217   friend fvar<RealType2, Order2> operator*(typename fvar<RealType2, Order2>::root_type const&,
 218                                            fvar<RealType2, Order2> const&);
 219
 220   // cr / cr2 | RealType | Binary Subtraction
 221   template <typename RealType2, size_t Order2>
 222   promote<fvar, fvar<RealType2, Order2>> operator/(fvar<RealType2, Order2> const&) const;
 223
 224   // cr / ca | RealType | Binary Subtraction
 225   fvar operator/(root_type const&) const;
 226
 227   // ca / cr | RealType | Binary Subtraction
 228   template <typename RealType2, size_t Order2>
 229   friend fvar<RealType2, Order2> operator/(typename fvar<RealType2, Order2>::root_type const&,
 230                                            fvar<RealType2, Order2> const&);
 231
 232   // For all comparison overloads, only the root term is compared.
 233
 234   // cr == cr2 | bool | Equality Comparison
 235   template <typename RealType2, size_t Order2>
 236   bool operator==(fvar<RealType2, Order2> const&) const;
 237
 238   // cr == ca | bool | Equality Comparison
 239   bool operator==(root_type const&) const;
 240
 241   // ca == cr | bool | Equality Comparison
 242   template <typename RealType2, size_t Order2>
 243   friend bool operator==(typename fvar<RealType2, Order2>::root_type const&, fvar<RealType2, Order2> const&);
 244
 245   // cr != cr2 | bool | Inequality Comparison
 246   template <typename RealType2, size_t Order2>
 247   bool operator!=(fvar<RealType2, Order2> const&) const;
 248
 249   // cr != ca | bool | Inequality Comparison
 250   bool operator!=(root_type const&) const;
 251
 252   // ca != cr | bool | Inequality Comparison
 253   template <typename RealType2, size_t Order2>
 254   friend bool operator!=(typename fvar<RealType2, Order2>::root_type const&, fvar<RealType2, Order2> const&);
 255
 256   // cr <= cr2 | bool | Less than equal to.
 257   template <typename RealType2, size_t Order2>
 258   bool operator<=(fvar<RealType2, Order2> const&) const;
 259
 260   // cr <= ca | bool | Less than equal to.
 261   bool operator<=(root_type const&) const;
 262
 263   // ca <= cr | bool | Less than equal to.
 264   template <typename RealType2, size_t Order2>
 265   friend bool operator<=(typename fvar<RealType2, Order2>::root_type const&, fvar<RealType2, Order2> const&);
 266
 267   // cr >= cr2 | bool | Greater than equal to.
 268   template <typename RealType2, size_t Order2>
 269   bool operator>=(fvar<RealType2, Order2> const&) const;
 270
 271   // cr >= ca | bool | Greater than equal to.
 272   bool operator>=(root_type const&) const;
 273
 274   // ca >= cr | bool | Greater than equal to.
 275   template <typename RealType2, size_t Order2>
 276   friend bool operator>=(typename fvar<RealType2, Order2>::root_type const&, fvar<RealType2, Order2> const&);
 277
 278   // cr < cr2 | bool | Less than comparison.
 279   template <typename RealType2, size_t Order2>
 280   bool operator<(fvar<RealType2, Order2> const&) const;
 281
 282   // cr < ca | bool | Less than comparison.
 283   bool operator<(root_type const&) const;
 284
 285   // ca < cr | bool | Less than comparison.
 286   template <typename RealType2, size_t Order2>
 287   friend bool operator<(typename fvar<RealType2, Order2>::root_type const&, fvar<RealType2, Order2> const&);
 288
 289   // cr > cr2 | bool | Greater than comparison.
 290   template <typename RealType2, size_t Order2>
 291   bool operator>(fvar<RealType2, Order2> const&) const;
 292
 293   // cr > ca | bool | Greater than comparison.
 294   bool operator>(root_type const&) const;
 295
 296   // ca > cr | bool | Greater than comparison.
 297   template <typename RealType2, size_t Order2>
 298   friend bool operator>(typename fvar<RealType2, Order2>::root_type const&, fvar<RealType2, Order2> const&);
 299
 300   // Will throw std::out_of_range if Order < order.
 301   template <typename... Orders>
 302   get_type_at<RealType, sizeof...(Orders)> at(size_t order, Orders... orders) const;
 303
 304   template <typename... Orders>
 305   get_type_at<fvar, sizeof...(Orders)> derivative(Orders... orders) const;
 306
 307   const RealType& operator[](size_t) const;
 308
 309   fvar inverse() const;  // Multiplicative inverse.
 310
 311   fvar& negate();  // Negate and return reference to *this.
 312
 313   static constexpr size_t depth = get_depth<fvar>::value;  // Number of nested std::array<RealType,Order>.
 314
 315   static constexpr size_t order_sum = get_order_sum<fvar>::value;
 316
 317   explicit operator root_type() const;  // Must be explicit, otherwise overloaded operators are ambiguous.
 318
 319   template <typename T, typename = typename boost::enable_if<boost::is_arithmetic<decay_t<T>>>::type>
 320   explicit operator T() const;  // Must be explicit; multiprecision has trouble without the std::enable_if
 321
 322   fvar& set_root(root_type const&);
 323
 324   // Apply coefficients using horner method.
 325   template <typename Func, typename Fvar, typename... Fvars>
 326   promote<fvar<RealType, Order>, Fvar, Fvars...> apply_coefficients(size_t const order,
 327                                                                     Func const& f,
 328                                                                     Fvar const& cr,
 329                                                                     Fvars&&... fvars) const;
 330
 331   template <typename Func>
 332   fvar apply_coefficients(size_t const order, Func const& f) const;
 333
 334   // Use when function returns derivative(i)/factorial(i) and may have some infinite derivatives.
 335   template <typename Func, typename Fvar, typename... Fvars>
 336   promote<fvar<RealType, Order>, Fvar, Fvars...> apply_coefficients_nonhorner(size_t const order,
 337                                                                               Func const& f,
 338                                                                               Fvar const& cr,
 339                                                                               Fvars&&... fvars) const;
 340
 341   template <typename Func>
 342   fvar apply_coefficients_nonhorner(size_t const order, Func const& f) const;
 343
 344   // Apply derivatives using horner method.
 345   template <typename Func, typename Fvar, typename... Fvars>
 346   promote<fvar<RealType, Order>, Fvar, Fvars...> apply_derivatives(size_t const order,
 347                                                                    Func const& f,
 348                                                                    Fvar const& cr,
 349                                                                    Fvars&&... fvars) const;
 350
 351   template <typename Func>
 352   fvar apply_derivatives(size_t const order, Func const& f) const;
 353
 354   // Use when function returns derivative(i) and may have some infinite derivatives.
 355   template <typename Func, typename Fvar, typename... Fvars>
 356   promote<fvar<RealType, Order>, Fvar, Fvars...> apply_derivatives_nonhorner(size_t const order,
 357                                                                              Func const& f,
 358                                                                              Fvar const& cr,
 359                                                                              Fvars&&... fvars) const;
 360
 361   template <typename Func>
 362   fvar apply_derivatives_nonhorner(size_t const order, Func const& f) const;
 363
 364  private:
 365   RealType epsilon_inner_product(size_t z0,
 366                                  size_t isum0,
 367                                  size_t m0,
 368                                  fvar const& cr,
 369                                  size_t z1,
 370                                  size_t isum1,
 371                                  size_t m1,
 372                                  size_t j) const;
 373
 374   fvar epsilon_multiply(size_t z0, size_t isum0, fvar const& cr, size_t z1, size_t isum1) const;
 375
 376   fvar epsilon_multiply(size_t z0, size_t isum0, root_type const& ca) const;
 377
 378   fvar inverse_apply() const;
 379
 380   fvar& multiply_assign_by_root_type(bool is_root, root_type const&);
 381
 382   template <typename RealType2, size_t Orders2>
 383   friend class fvar;
 384
 385   template <typename RealType2, size_t Order2>
 386   friend std::ostream& operator<<(std::ostream&, fvar<RealType2, Order2> const&);
 387
 388   // C++11 Compatibility
 389 #ifdef BOOST_NO_CXX17_IF_CONSTEXPR
 390   template <typename RootType>
 391   void fvar_cpp11(std::true_type, RootType const& ca, bool const is_variable);
 392
 393   template <typename RootType>
 394   void fvar_cpp11(std::false_type, RootType const& ca, bool const is_variable);
 395
 396   template <typename... Orders>
 397   get_type_at<RealType, sizeof...(Orders)> at_cpp11(std::true_type, size_t order, Orders... orders) const;
 398
 399   template <typename... Orders>
 400   get_type_at<RealType, sizeof...(Orders)> at_cpp11(std::false_type, size_t order, Orders... orders) const;
 401
 402   template <typename SizeType>
 403   fvar epsilon_multiply_cpp11(std::true_type,
 404                               SizeType z0,
 405                               size_t isum0,
 406                               fvar const& cr,
 407                               size_t z1,
 408                               size_t isum1) const;
 409
 410   template <typename SizeType>
 411   fvar epsilon_multiply_cpp11(std::false_type,
 412                               SizeType z0,
 413                               size_t isum0,
 414                               fvar const& cr,
 415                               size_t z1,
 416                               size_t isum1) const;
 417
 418   template <typename SizeType>
 419   fvar epsilon_multiply_cpp11(std::true_type, SizeType z0, size_t isum0, root_type const& ca) const;
 420
 421   template <typename SizeType>
 422   fvar epsilon_multiply_cpp11(std::false_type, SizeType z0, size_t isum0, root_type const& ca) const;
 423
 424   template <typename RootType>
 425   fvar& multiply_assign_by_root_type_cpp11(std::true_type, bool is_root, RootType const& ca);
 426
 427   template <typename RootType>
 428   fvar& multiply_assign_by_root_type_cpp11(std::false_type, bool is_root, RootType const& ca);
 429
 430   template <typename RootType>
 431   fvar& negate_cpp11(std::true_type, RootType const&);
 432
 433   template <typename RootType>
 434   fvar& negate_cpp11(std::false_type, RootType const&);
 435
 436   template <typename RootType>
 437   fvar& set_root_cpp11(std::true_type, RootType const& root);
 438
 439   template <typename RootType>
 440   fvar& set_root_cpp11(std::false_type, RootType const& root);
 441 #endif
 442 };
 443
 444 // C++11 compatibility
 445 #ifdef BOOST_NO_CXX17_IF_CONSTEXPR
 446 #define BOOST_AUTODIFF_IF_CONSTEXPR
 447 #else
 448 #define BOOST_AUTODIFF_IF_CONSTEXPR constexpr
 449 #endif
 450
 451 // Standard Library Support Requirements
 452
 453 // fabs(cr1) | RealType
 454 template <typename RealType, size_t Order>
 455 fvar<RealType, Order> fabs(fvar<RealType, Order> const&);
 456
 457 // abs(cr1) | RealType
 458 template <typename RealType, size_t Order>
 459 fvar<RealType, Order> abs(fvar<RealType, Order> const&);
 460
 461 // ceil(cr1) | RealType
 462 template <typename RealType, size_t Order>
 463 fvar<RealType, Order> ceil(fvar<RealType, Order> const&);
 464
 465 // floor(cr1) | RealType
 466 template <typename RealType, size_t Order>
 467 fvar<RealType, Order> floor(fvar<RealType, Order> const&);
 468
 469 // exp(cr1) | RealType
 470 template <typename RealType, size_t Order>
 471 fvar<RealType, Order> exp(fvar<RealType, Order> const&);
 472
 473 // pow(cr, ca) | RealType
 474 template <typename RealType, size_t Order>
 475 fvar<RealType, Order> pow(fvar<RealType, Order> const&, typename fvar<RealType, Order>::root_type const&);
 476
 477 // pow(ca, cr) | RealType
 478 template <typename RealType, size_t Order>
 479 fvar<RealType, Order> pow(typename fvar<RealType, Order>::root_type const&, fvar<RealType, Order> const&);
 480
 481 // pow(cr1, cr2) | RealType
 482 template <typename RealType1, size_t Order1, typename RealType2, size_t Order2>
 483 promote<fvar<RealType1, Order1>, fvar<RealType2, Order2>> pow(fvar<RealType1, Order1> const&,
 484                                                               fvar<RealType2, Order2> const&);
 485
 486 // sqrt(cr1) | RealType
 487 template <typename RealType, size_t Order>
 488 fvar<RealType, Order> sqrt(fvar<RealType, Order> const&);
 489
 490 // log(cr1) | RealType
 491 template <typename RealType, size_t Order>
 492 fvar<RealType, Order> log(fvar<RealType, Order> const&);
 493
 494 // frexp(cr1, &i) | RealType
 495 template <typename RealType, size_t Order>
 496 fvar<RealType, Order> frexp(fvar<RealType, Order> const&, int*);
 497
 498 // ldexp(cr1, i) | RealType
 499 template <typename RealType, size_t Order>
 500 fvar<RealType, Order> ldexp(fvar<RealType, Order> const&, int);
 501
 502 // cos(cr1) | RealType
 503 template <typename RealType, size_t Order>
 504 fvar<RealType, Order> cos(fvar<RealType, Order> const&);
 505
 506 // sin(cr1) | RealType
 507 template <typename RealType, size_t Order>
 508 fvar<RealType, Order> sin(fvar<RealType, Order> const&);
 509
 510 // asin(cr1) | RealType
 511 template <typename RealType, size_t Order>
 512 fvar<RealType, Order> asin(fvar<RealType, Order> const&);
 513
 514 // tan(cr1) | RealType
 515 template <typename RealType, size_t Order>
 516 fvar<RealType, Order> tan(fvar<RealType, Order> const&);
 517
 518 // atan(cr1) | RealType
 519 template <typename RealType, size_t Order>
 520 fvar<RealType, Order> atan(fvar<RealType, Order> const&);
 521
 522 // atan2(cr, ca) | RealType
 523 template <typename RealType, size_t Order>
 524 fvar<RealType, Order> atan2(fvar<RealType, Order> const&, typename fvar<RealType, Order>::root_type const&);
 525
 526 // atan2(ca, cr) | RealType
 527 template <typename RealType, size_t Order>
 528 fvar<RealType, Order> atan2(typename fvar<RealType, Order>::root_type const&, fvar<RealType, Order> const&);
 529
 530 // atan2(cr1, cr2) | RealType
 531 template <typename RealType1, size_t Order1, typename RealType2, size_t Order2>
 532 promote<fvar<RealType1, Order1>, fvar<RealType2, Order2>> atan2(fvar<RealType1, Order1> const&,
 533                                                                 fvar<RealType2, Order2> const&);
 534
 535 // fmod(cr1,cr2) | RealType
 536 template <typename RealType1, size_t Order1, typename RealType2, size_t Order2>
 537 promote<fvar<RealType1, Order1>, fvar<RealType2, Order2>> fmod(fvar<RealType1, Order1> const&,
 538                                                                fvar<RealType2, Order2> const&);
 539
 540 // round(cr1) | RealType
 541 template <typename RealType, size_t Order>
 542 fvar<RealType, Order> round(fvar<RealType, Order> const&);
 543
 544 // iround(cr1) | int
 545 template <typename RealType, size_t Order>
 546 int iround(fvar<RealType, Order> const&);
 547
 548 template <typename RealType, size_t Order>
 549 long lround(fvar<RealType, Order> const&);
 550
 551 template <typename RealType, size_t Order>
 552 long long llround(fvar<RealType, Order> const&);
 553
 554 // trunc(cr1) | RealType
 555 template <typename RealType, size_t Order>
 556 fvar<RealType, Order> trunc(fvar<RealType, Order> const&);
 557
 558 template <typename RealType, size_t Order>
 559 long double truncl(fvar<RealType, Order> const&);
 560
 561 // itrunc(cr1) | int
 562 template <typename RealType, size_t Order>
 563 int itrunc(fvar<RealType, Order> const&);
 564
 565 template <typename RealType, size_t Order>
 566 long long lltrunc(fvar<RealType, Order> const&);
 567
 568 // Additional functions
 569 template <typename RealType, size_t Order>
 570 fvar<RealType, Order> acos(fvar<RealType, Order> const&);
 571
 572 template <typename RealType, size_t Order>
 573 fvar<RealType, Order> acosh(fvar<RealType, Order> const&);
 574
 575 template <typename RealType, size_t Order>
 576 fvar<RealType, Order> asinh(fvar<RealType, Order> const&);
 577
 578 template <typename RealType, size_t Order>
 579 fvar<RealType, Order> atanh(fvar<RealType, Order> const&);
 580
 581 template <typename RealType, size_t Order>
 582 fvar<RealType, Order> cosh(fvar<RealType, Order> const&);
 583
 584 template <typename RealType, size_t Order>
 585 fvar<RealType, Order> digamma(fvar<RealType, Order> const&);
 586
 587 template <typename RealType, size_t Order>
 588 fvar<RealType, Order> erf(fvar<RealType, Order> const&);
 589
 590 template <typename RealType, size_t Order>
 591 fvar<RealType, Order> erfc(fvar<RealType, Order> const&);
 592
 593 template <typename RealType, size_t Order>
 594 fvar<RealType, Order> lambert_w0(fvar<RealType, Order> const&);
 595
 596 template <typename RealType, size_t Order>
 597 fvar<RealType, Order> lgamma(fvar<RealType, Order> const&);
 598
 599 template <typename RealType, size_t Order>
 600 fvar<RealType, Order> sinc(fvar<RealType, Order> const&);
 601
 602 template <typename RealType, size_t Order>
 603 fvar<RealType, Order> sinh(fvar<RealType, Order> const&);
 604
 605 template <typename RealType, size_t Order>
 606 fvar<RealType, Order> tanh(fvar<RealType, Order> const&);
 607
 608 template <typename RealType, size_t Order>
 609 fvar<RealType, Order> tgamma(fvar<RealType, Order> const&);
 610
 611 template <size_t>
 612 struct zero : std::integral_constant<size_t, 0> {};
 613
 614 }  // namespace detail
 615
 616 template <typename RealType, size_t Order, size_t... Orders>
 617 using autodiff_fvar = typename detail::nest_fvar<RealType, Order, Orders...>::type;
 618
 619 template <typename RealType, size_t Order, size_t... Orders>
 620 autodiff_fvar<RealType, Order, Orders...> make_fvar(RealType const& ca) {
 621   return autodiff_fvar<RealType, Order, Orders...>(ca, true);
 622 }
 623
 624 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
 625 namespace detail {
 626
 627 template <typename RealType, size_t Order, size_t... Is>
 628 auto make_fvar_for_tuple(std::index_sequence<Is...>, RealType const& ca) {
 629   return make_fvar<RealType, zero<Is>::value..., Order>(ca);
 630 }
 631
 632 template <typename RealType, size_t... Orders, size_t... Is, typename... RealTypes>
 633 auto make_ftuple_impl(std::index_sequence<Is...>, RealTypes const&... ca) {
 634   return std::make_tuple(make_fvar_for_tuple<RealType, Orders>(std::make_index_sequence<Is>{}, ca)...);
 635 }
 636
 637 }  // namespace detail
 638
 639 template <typename RealType, size_t... Orders, typename... RealTypes>
 640 auto make_ftuple(RealTypes const&... ca) {
 641   static_assert(sizeof...(Orders) == sizeof...(RealTypes),
 642                 "Number of Orders must match number of function parameters.");
 643   return detail::make_ftuple_impl<RealType, Orders...>(std::index_sequence_for<RealTypes...>{}, ca...);
 644 }
 645 #endif
 646
 647 namespace detail {
 648
 649 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
 650 template <typename RealType, size_t Order>
 651 fvar<RealType, Order>::fvar(root_type const& ca, bool const is_variable) {
 652   if constexpr (is_fvar<RealType>::value) {
 653     v.front() = RealType(ca, is_variable);
 654     if constexpr (0 < Order)
 655       std::fill(v.begin() + 1, v.end(), static_cast<RealType>(0));
 656   } else {
 657     v.front() = ca;
 658     if constexpr (0 < Order)
 659       v[1] = static_cast<root_type>(static_cast<int>(is_variable));
 660     if constexpr (1 < Order)
 661       std::fill(v.begin() + 2, v.end(), static_cast<RealType>(0));
 662   }
 663 }
 664 #endif
 665
 666 template <typename RealType, size_t Order>
 667 template <typename RealType2, size_t Order2>
 668 fvar<RealType, Order>::fvar(fvar<RealType2, Order2> const& cr) {
 669   for (size_t i = 0; i <= (std::min)(Order, Order2); ++i)
 670     v[i] = static_cast<RealType>(cr.v[i]);
 671   if BOOST_AUTODIFF_IF_CONSTEXPR (Order2 < Order)
 672     std::fill(v.begin() + (Order2 + 1), v.end(), static_cast<RealType>(0));
 673 }
 674
 675 template <typename RealType, size_t Order>
 676 fvar<RealType, Order>::fvar(root_type const& ca) : v{{static_cast<RealType>(ca)}} {}
 677
 678 // Can cause compiler error if RealType2 cannot be cast to root_type.
 679 template <typename RealType, size_t Order>
 680 template <typename RealType2>
 681 fvar<RealType, Order>::fvar(RealType2 const& ca) : v{{static_cast<RealType>(ca)}} {}
 682
 683 /*
 684 template<typename RealType, size_t Order>
 685 fvar<RealType,Order>& fvar<RealType,Order>::operator=(root_type const& ca)
 686 {
 687     v.front() = static_cast<RealType>(ca);
 688     if constexpr (0 < Order)
 689         std::fill(v.begin()+1, v.end(), static_cast<RealType>(0));
 690     return *this;
 691 }
 692 */
 693
 694 template <typename RealType, size_t Order>
 695 template <typename RealType2, size_t Order2>
 696 fvar<RealType, Order>& fvar<RealType, Order>::operator+=(fvar<RealType2, Order2> const& cr) {
 697   for (size_t i = 0; i <= (std::min)(Order, Order2); ++i)
 698     v[i] += cr.v[i];
 699   return *this;
 700 }
 701
 702 template <typename RealType, size_t Order>
 703 fvar<RealType, Order>& fvar<RealType, Order>::operator+=(root_type const& ca) {
 704   v.front() += ca;
 705   return *this;
 706 }
 707
 708 template <typename RealType, size_t Order>
 709 template <typename RealType2, size_t Order2>
 710 fvar<RealType, Order>& fvar<RealType, Order>::operator-=(fvar<RealType2, Order2> const& cr) {
 711   for (size_t i = 0; i <= Order; ++i)
 712     v[i] -= cr.v[i];
 713   return *this;
 714 }
 715
 716 template <typename RealType, size_t Order>
 717 fvar<RealType, Order>& fvar<RealType, Order>::operator-=(root_type const& ca) {
 718   v.front() -= ca;
 719   return *this;
 720 }
 721
 722 template <typename RealType, size_t Order>
 723 template <typename RealType2, size_t Order2>
 724 fvar<RealType, Order>& fvar<RealType, Order>::operator*=(fvar<RealType2, Order2> const& cr) {
 725   using diff_t = typename std::array<RealType, Order + 1>::difference_type;
 726   promote<RealType, RealType2> const zero(0);
 727   if BOOST_AUTODIFF_IF_CONSTEXPR (Order <= Order2)
 728     for (size_t i = 0, j = Order; i <= Order; ++i, --j)
 729       v[j] = std::inner_product(v.cbegin(), v.cend() - diff_t(i), cr.v.crbegin() + diff_t(i), zero);
 730   else {
 731     for (size_t i = 0, j = Order; i <= Order - Order2; ++i, --j)
 732       v[j] = std::inner_product(cr.v.cbegin(), cr.v.cend(), v.crbegin() + diff_t(i), zero);
 733     for (size_t i = Order - Order2 + 1, j = Order2 - 1; i <= Order; ++i, --j)
 734       v[j] = std::inner_product(cr.v.cbegin(), cr.v.cbegin() + diff_t(j + 1), v.crbegin() + diff_t(i), zero);
 735   }
 736   return *this;
 737 }
 738
 739 template <typename RealType, size_t Order>
 740 fvar<RealType, Order>& fvar<RealType, Order>::operator*=(root_type const& ca) {
 741   return multiply_assign_by_root_type(true, ca);
 742 }
 743
 744 template <typename RealType, size_t Order>
 745 template <typename RealType2, size_t Order2>
 746 fvar<RealType, Order>& fvar<RealType, Order>::operator/=(fvar<RealType2, Order2> const& cr) {
 747   using diff_t = typename std::array<RealType, Order + 1>::difference_type;
 748   RealType const zero(0);
 749   v.front() /= cr.v.front();
 750   if BOOST_AUTODIFF_IF_CONSTEXPR (Order < Order2)
 751     for (size_t i = 1, j = Order2 - 1, k = Order; i <= Order; ++i, --j, --k)
 752       (v[i] -= std::inner_product(
 753            cr.v.cbegin() + 1, cr.v.cend() - diff_t(j), v.crbegin() + diff_t(k), zero)) /= cr.v.front();
 754   else if BOOST_AUTODIFF_IF_CONSTEXPR (0 < Order2)
 755     for (size_t i = 1, j = Order2 - 1, k = Order; i <= Order; ++i, j && --j, --k)
 756       (v[i] -= std::inner_product(
 757            cr.v.cbegin() + 1, cr.v.cend() - diff_t(j), v.crbegin() + diff_t(k), zero)) /= cr.v.front();
 758   else
 759     for (size_t i = 1; i <= Order; ++i)
 760       v[i] /= cr.v.front();
 761   return *this;
 762 }
 763
 764 template <typename RealType, size_t Order>
 765 fvar<RealType, Order>& fvar<RealType, Order>::operator/=(root_type const& ca) {
 766   std::for_each(v.begin(), v.end(), [&ca](RealType& x) { x /= ca; });
 767   return *this;
 768 }
 769
 770 template <typename RealType, size_t Order>
 771 fvar<RealType, Order> fvar<RealType, Order>::operator-() const {
 772   fvar<RealType, Order> retval(*this);
 773   retval.negate();
 774   return retval;
 775 }
 776
 777 template <typename RealType, size_t Order>
 778 fvar<RealType, Order> const& fvar<RealType, Order>::operator+() const {
 779   return *this;
 780 }
 781
 782 template <typename RealType, size_t Order>
 783 template <typename RealType2, size_t Order2>
 784 promote<fvar<RealType, Order>, fvar<RealType2, Order2>> fvar<RealType, Order>::operator+(
 785     fvar<RealType2, Order2> const& cr) const {
 786   promote<fvar<RealType, Order>, fvar<RealType2, Order2>> retval;
 787   for (size_t i = 0; i <= (std::min)(Order, Order2); ++i)
 788     retval.v[i] = v[i] + cr.v[i];
 789   if BOOST_AUTODIFF_IF_CONSTEXPR (Order < Order2)
 790     for (size_t i = Order + 1; i <= Order2; ++i)
 791       retval.v[i] = cr.v[i];
 792   else if BOOST_AUTODIFF_IF_CONSTEXPR (Order2 < Order)
 793     for (size_t i = Order2 + 1; i <= Order; ++i)
 794       retval.v[i] = v[i];
 795   return retval;
 796 }
 797
 798 template <typename RealType, size_t Order>
 799 fvar<RealType, Order> fvar<RealType, Order>::operator+(root_type const& ca) const {
 800   fvar<RealType, Order> retval(*this);
 801   retval.v.front() += ca;
 802   return retval;
 803 }
 804
 805 template <typename RealType, size_t Order>
 806 fvar<RealType, Order> operator+(typename fvar<RealType, Order>::root_type const& ca,
 807                                 fvar<RealType, Order> const& cr) {
 808   return cr + ca;
 809 }
 810
 811 template <typename RealType, size_t Order>
 812 template <typename RealType2, size_t Order2>
 813 promote<fvar<RealType, Order>, fvar<RealType2, Order2>> fvar<RealType, Order>::operator-(
 814     fvar<RealType2, Order2> const& cr) const {
 815   promote<fvar<RealType, Order>, fvar<RealType2, Order2>> retval;
 816   for (size_t i = 0; i <= (std::min)(Order, Order2); ++i)
 817     retval.v[i] = v[i] - cr.v[i];
 818   if BOOST_AUTODIFF_IF_CONSTEXPR (Order < Order2)
 819     for (auto i = Order + 1; i <= Order2; ++i)
 820       retval.v[i] = -cr.v[i];
 821   else if BOOST_AUTODIFF_IF_CONSTEXPR (Order2 < Order)
 822     for (auto i = Order2 + 1; i <= Order; ++i)
 823       retval.v[i] = v[i];
 824   return retval;
 825 }
 826
 827 template <typename RealType, size_t Order>
 828 fvar<RealType, Order> fvar<RealType, Order>::operator-(root_type const& ca) const {
 829   fvar<RealType, Order> retval(*this);
 830   retval.v.front() -= ca;
 831   return retval;
 832 }
 833
 834 template <typename RealType, size_t Order>
 835 fvar<RealType, Order> operator-(typename fvar<RealType, Order>::root_type const& ca,
 836                                 fvar<RealType, Order> const& cr) {
 837   fvar<RealType, Order> mcr = -cr;  // Has same address as retval in operator-() due to NRVO.
 838   mcr += ca;
 839   return mcr;  // <-- This allows for NRVO. The following does not. --> return mcr += ca;
 840 }
 841
 842 template <typename RealType, size_t Order>
 843 template <typename RealType2, size_t Order2>
 844 promote<fvar<RealType, Order>, fvar<RealType2, Order2>> fvar<RealType, Order>::operator*(
 845     fvar<RealType2, Order2> const& cr) const {
 846   using diff_t = typename std::array<RealType, Order + 1>::difference_type;
 847   promote<RealType, RealType2> const zero(0);
 848   promote<fvar<RealType, Order>, fvar<RealType2, Order2>> retval;
 849   if BOOST_AUTODIFF_IF_CONSTEXPR (Order < Order2)
 850     for (size_t i = 0, j = Order, k = Order2; i <= Order2; ++i, j && --j, --k)
 851       retval.v[i] = std::inner_product(v.cbegin(), v.cend() - diff_t(j), cr.v.crbegin() + diff_t(k), zero);
 852   else
 853     for (size_t i = 0, j = Order2, k = Order; i <= Order; ++i, j && --j, --k)
 854       retval.v[i] = std::inner_product(cr.v.cbegin(), cr.v.cend() - diff_t(j), v.crbegin() + diff_t(k), zero);
 855   return retval;
 856 }
 857
 858 template <typename RealType, size_t Order>
 859 fvar<RealType, Order> fvar<RealType, Order>::operator*(root_type const& ca) const {
 860   fvar<RealType, Order> retval(*this);
 861   retval *= ca;
 862   return retval;
 863 }
 864
 865 template <typename RealType, size_t Order>
 866 fvar<RealType, Order> operator*(typename fvar<RealType, Order>::root_type const& ca,
 867                                 fvar<RealType, Order> const& cr) {
 868   return cr * ca;
 869 }
 870
 871 template <typename RealType, size_t Order>
 872 template <typename RealType2, size_t Order2>
 873 promote<fvar<RealType, Order>, fvar<RealType2, Order2>> fvar<RealType, Order>::operator/(
 874     fvar<RealType2, Order2> const& cr) const {
 875   using diff_t = typename std::array<RealType, Order + 1>::difference_type;
 876   promote<RealType, RealType2> const zero(0);
 877   promote<fvar<RealType, Order>, fvar<RealType2, Order2>> retval;
 878   retval.v.front() = v.front() / cr.v.front();
 879   if BOOST_AUTODIFF_IF_CONSTEXPR (Order < Order2) {
 880     for (size_t i = 1, j = Order2 - 1; i <= Order; ++i, --j)
 881       retval.v[i] =
 882           (v[i] - std::inner_product(
 883                       cr.v.cbegin() + 1, cr.v.cend() - diff_t(j), retval.v.crbegin() + diff_t(j + 1), zero)) /
 884           cr.v.front();
 885     for (size_t i = Order + 1, j = Order2 - Order - 1; i <= Order2; ++i, --j)
 886       retval.v[i] =
 887           -std::inner_product(
 888               cr.v.cbegin() + 1, cr.v.cend() - diff_t(j), retval.v.crbegin() + diff_t(j + 1), zero) /
 889           cr.v.front();
 890   } else if BOOST_AUTODIFF_IF_CONSTEXPR (0 < Order2)
 891     for (size_t i = 1, j = Order2 - 1, k = Order; i <= Order; ++i, j && --j, --k)
 892       retval.v[i] =
 893           (v[i] - std::inner_product(
 894                       cr.v.cbegin() + 1, cr.v.cend() - diff_t(j), retval.v.crbegin() + diff_t(k), zero)) /
 895           cr.v.front();
 896   else
 897     for (size_t i = 1; i <= Order; ++i)
 898       retval.v[i] = v[i] / cr.v.front();
 899   return retval;
 900 }
 901
 902 template <typename RealType, size_t Order>
 903 fvar<RealType, Order> fvar<RealType, Order>::operator/(root_type const& ca) const {
 904   fvar<RealType, Order> retval(*this);
 905   retval /= ca;
 906   return retval;
 907 }
 908
 909 template <typename RealType, size_t Order>
 910 fvar<RealType, Order> operator/(typename fvar<RealType, Order>::root_type const& ca,
 911                                 fvar<RealType, Order> const& cr) {
 912   using diff_t = typename std::array<RealType, Order + 1>::difference_type;
 913   fvar<RealType, Order> retval;
 914   retval.v.front() = ca / cr.v.front();
 915   if BOOST_AUTODIFF_IF_CONSTEXPR (0 < Order) {
 916     RealType const zero(0);
 917     for (size_t i = 1, j = Order - 1; i <= Order; ++i, --j)
 918       retval.v[i] =
 919           -std::inner_product(
 920               cr.v.cbegin() + 1, cr.v.cend() - diff_t(j), retval.v.crbegin() + diff_t(j + 1), zero) /
 921           cr.v.front();
 922   }
 923   return retval;
 924 }
 925
 926 template <typename RealType, size_t Order>
 927 template <typename RealType2, size_t Order2>
 928 bool fvar<RealType, Order>::operator==(fvar<RealType2, Order2> const& cr) const {
 929   return v.front() == cr.v.front();
 930 }
 931
 932 template <typename RealType, size_t Order>
 933 bool fvar<RealType, Order>::operator==(root_type const& ca) const {
 934   return v.front() == ca;
 935 }
 936
 937 template <typename RealType, size_t Order>
 938 bool operator==(typename fvar<RealType, Order>::root_type const& ca, fvar<RealType, Order> const& cr) {
 939   return ca == cr.v.front();
 940 }
 941
 942 template <typename RealType, size_t Order>
 943 template <typename RealType2, size_t Order2>
 944 bool fvar<RealType, Order>::operator!=(fvar<RealType2, Order2> const& cr) const {
 945   return v.front() != cr.v.front();
 946 }
 947
 948 template <typename RealType, size_t Order>
 949 bool fvar<RealType, Order>::operator!=(root_type const& ca) const {
 950   return v.front() != ca;
 951 }
 952
 953 template <typename RealType, size_t Order>
 954 bool operator!=(typename fvar<RealType, Order>::root_type const& ca, fvar<RealType, Order> const& cr) {
 955   return ca != cr.v.front();
 956 }
 957
 958 template <typename RealType, size_t Order>
 959 template <typename RealType2, size_t Order2>
 960 bool fvar<RealType, Order>::operator<=(fvar<RealType2, Order2> const& cr) const {
 961   return v.front() <= cr.v.front();
 962 }
 963
 964 template <typename RealType, size_t Order>
 965 bool fvar<RealType, Order>::operator<=(root_type const& ca) const {
 966   return v.front() <= ca;
 967 }
 968
 969 template <typename RealType, size_t Order>
 970 bool operator<=(typename fvar<RealType, Order>::root_type const& ca, fvar<RealType, Order> const& cr) {
 971   return ca <= cr.v.front();
 972 }
 973
 974 template <typename RealType, size_t Order>
 975 template <typename RealType2, size_t Order2>
 976 bool fvar<RealType, Order>::operator>=(fvar<RealType2, Order2> const& cr) const {
 977   return v.front() >= cr.v.front();
 978 }
 979
 980 template <typename RealType, size_t Order>
 981 bool fvar<RealType, Order>::operator>=(root_type const& ca) const {
 982   return v.front() >= ca;
 983 }
 984
 985 template <typename RealType, size_t Order>
 986 bool operator>=(typename fvar<RealType, Order>::root_type const& ca, fvar<RealType, Order> const& cr) {
 987   return ca >= cr.v.front();
 988 }
 989
 990 template <typename RealType, size_t Order>
 991 template <typename RealType2, size_t Order2>
 992 bool fvar<RealType, Order>::operator<(fvar<RealType2, Order2> const& cr) const {
 993   return v.front() < cr.v.front();
 994 }
 995
 996 template <typename RealType, size_t Order>
 997 bool fvar<RealType, Order>::operator<(root_type const& ca) const {
 998   return v.front() < ca;
 999 }
1000
1001 template <typename RealType, size_t Order>
1002 bool operator<(typename fvar<RealType, Order>::root_type const& ca, fvar<RealType, Order> const& cr) {
1003   return ca < cr.v.front();
1004 }
1005
1006 template <typename RealType, size_t Order>
1007 template <typename RealType2, size_t Order2>
1008 bool fvar<RealType, Order>::operator>(fvar<RealType2, Order2> const& cr) const {
1009   return v.front() > cr.v.front();
1010 }
1011
1012 template <typename RealType, size_t Order>
1013 bool fvar<RealType, Order>::operator>(root_type const& ca) const {
1014   return v.front() > ca;
1015 }
1016
1017 template <typename RealType, size_t Order>
1018 bool operator>(typename fvar<RealType, Order>::root_type const& ca, fvar<RealType, Order> const& cr) {
1019   return ca > cr.v.front();
1020 }
1021
1022   /*** Other methods and functions ***/
1023
1024 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1025 // f : order -> derivative(order)/factorial(order)
1026 // Use this when you have the polynomial coefficients, rather than just the derivatives. E.g. See atan2().
1027 template <typename RealType, size_t Order>
1028 template <typename Func, typename Fvar, typename... Fvars>
1029 promote<fvar<RealType, Order>, Fvar, Fvars...> fvar<RealType, Order>::apply_coefficients(
1030     size_t const order,
1031     Func const& f,
1032     Fvar const& cr,
1033     Fvars&&... fvars) const {
1034   fvar<RealType, Order> const epsilon = fvar<RealType, Order>(*this).set_root(0);
1035   size_t i = (std::min)(order, order_sum);
1036   promote<fvar<RealType, Order>, Fvar, Fvars...> accumulator = cr.apply_coefficients(
1037       order - i, [&f, i](auto... indices) { return f(i, indices...); }, std::forward<Fvars>(fvars)...);
1038   while (i--)
1039     (accumulator *= epsilon) += cr.apply_coefficients(
1040         order - i, [&f, i](auto... indices) { return f(i, indices...); }, std::forward<Fvars>(fvars)...);
1041   return accumulator;
1042 }
1043 #endif
1044
1045 // f : order -> derivative(order)/factorial(order)
1046 // Use this when you have the polynomial coefficients, rather than just the derivatives. E.g. See atan().
1047 template <typename RealType, size_t Order>
1048 template <typename Func>
1049 fvar<RealType, Order> fvar<RealType, Order>::apply_coefficients(size_t const order, Func const& f) const {
1050   fvar<RealType, Order> const epsilon = fvar<RealType, Order>(*this).set_root(0);
1051 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1052   size_t i = (std::min)(order, order_sum);
1053 #else  // ODR-use of static constexpr
1054   size_t i = order < order_sum ? order : order_sum;
1055 #endif
1056   fvar<RealType, Order> accumulator = f(i);
1057   while (i--)
1058     (accumulator *= epsilon) += f(i);
1059   return accumulator;
1060 }
1061
1062 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1063 // f : order -> derivative(order)
1064 template <typename RealType, size_t Order>
1065 template <typename Func, typename Fvar, typename... Fvars>
1066 promote<fvar<RealType, Order>, Fvar, Fvars...> fvar<RealType, Order>::apply_coefficients_nonhorner(
1067     size_t const order,
1068     Func const& f,
1069     Fvar const& cr,
1070     Fvars&&... fvars) const {
1071   fvar<RealType, Order> const epsilon = fvar<RealType, Order>(*this).set_root(0);
1072   fvar<RealType, Order> epsilon_i = fvar<RealType, Order>(1);  // epsilon to the power of i
1073   promote<fvar<RealType, Order>, Fvar, Fvars...> accumulator = cr.apply_coefficients_nonhorner(
1074       order,
1075       [&f](auto... indices) { return f(0, static_cast<std::size_t>(indices)...); },
1076       std::forward<Fvars>(fvars)...);
1077   size_t const i_max = (std::min)(order, order_sum);
1078   for (size_t i = 1; i <= i_max; ++i) {
1079     epsilon_i = epsilon_i.epsilon_multiply(i - 1, 0, epsilon, 1, 0);
1080     accumulator += epsilon_i.epsilon_multiply(
1081         i,
1082         0,
1083         cr.apply_coefficients_nonhorner(
1084             order - i,
1085             [&f, i](auto... indices) { return f(i, static_cast<std::size_t>(indices)...); },
1086             std::forward<Fvars>(fvars)...),
1087         0,
1088         0);
1089   }
1090   return accumulator;
1091 }
1092 #endif
1093
1094 // f : order -> coefficient(order)
1095 template <typename RealType, size_t Order>
1096 template <typename Func>
1097 fvar<RealType, Order> fvar<RealType, Order>::apply_coefficients_nonhorner(size_t const order,
1098                                                                           Func const& f) const {
1099   fvar<RealType, Order> const epsilon = fvar<RealType, Order>(*this).set_root(0);
1100   fvar<RealType, Order> epsilon_i = fvar<RealType, Order>(1);  // epsilon to the power of i
1101   fvar<RealType, Order> accumulator = fvar<RealType, Order>(f(0u));
1102 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1103   size_t const i_max = (std::min)(order, order_sum);
1104 #else  // ODR-use of static constexpr
1105   size_t const i_max = order < order_sum ? order : order_sum;
1106 #endif
1107   for (size_t i = 1; i <= i_max; ++i) {
1108     epsilon_i = epsilon_i.epsilon_multiply(i - 1, 0, epsilon, 1, 0);
1109     accumulator += epsilon_i.epsilon_multiply(i, 0, f(i));
1110   }
1111   return accumulator;
1112 }
1113
1114 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1115 // f : order -> derivative(order)
1116 template <typename RealType, size_t Order>
1117 template <typename Func, typename Fvar, typename... Fvars>
1118 promote<fvar<RealType, Order>, Fvar, Fvars...> fvar<RealType, Order>::apply_derivatives(
1119     size_t const order,
1120     Func const& f,
1121     Fvar const& cr,
1122     Fvars&&... fvars) const {
1123   fvar<RealType, Order> const epsilon = fvar<RealType, Order>(*this).set_root(0);
1124   size_t i = (std::min)(order, order_sum);
1125   promote<fvar<RealType, Order>, Fvar, Fvars...> accumulator =
1126       cr.apply_derivatives(
1127           order - i, [&f, i](auto... indices) { return f(i, indices...); }, std::forward<Fvars>(fvars)...) /
1128       factorial<root_type>(static_cast<unsigned>(i));
1129   while (i--)
1130     (accumulator *= epsilon) +=
1131         cr.apply_derivatives(
1132             order - i, [&f, i](auto... indices) { return f(i, indices...); }, std::forward<Fvars>(fvars)...) /
1133         factorial<root_type>(static_cast<unsigned>(i));
1134   return accumulator;
1135 }
1136 #endif
1137
1138 // f : order -> derivative(order)
1139 template <typename RealType, size_t Order>
1140 template <typename Func>
1141 fvar<RealType, Order> fvar<RealType, Order>::apply_derivatives(size_t const order, Func const& f) const {
1142   fvar<RealType, Order> const epsilon = fvar<RealType, Order>(*this).set_root(0);
1143 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1144   size_t i = (std::min)(order, order_sum);
1145 #else  // ODR-use of static constexpr
1146   size_t i = order < order_sum ? order : order_sum;
1147 #endif
1148   fvar<RealType, Order> accumulator = f(i) / factorial<root_type>(static_cast<unsigned>(i));
1149   while (i--)
1150     (accumulator *= epsilon) += f(i) / factorial<root_type>(static_cast<unsigned>(i));
1151   return accumulator;
1152 }
1153
1154 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1155 // f : order -> derivative(order)
1156 template <typename RealType, size_t Order>
1157 template <typename Func, typename Fvar, typename... Fvars>
1158 promote<fvar<RealType, Order>, Fvar, Fvars...> fvar<RealType, Order>::apply_derivatives_nonhorner(
1159     size_t const order,
1160     Func const& f,
1161     Fvar const& cr,
1162     Fvars&&... fvars) const {
1163   fvar<RealType, Order> const epsilon = fvar<RealType, Order>(*this).set_root(0);
1164   fvar<RealType, Order> epsilon_i = fvar<RealType, Order>(1);  // epsilon to the power of i
1165   promote<fvar<RealType, Order>, Fvar, Fvars...> accumulator = cr.apply_derivatives_nonhorner(
1166       order,
1167       [&f](auto... indices) { return f(0, static_cast<std::size_t>(indices)...); },
1168       std::forward<Fvars>(fvars)...);
1169   size_t const i_max = (std::min)(order, order_sum);
1170   for (size_t i = 1; i <= i_max; ++i) {
1171     epsilon_i = epsilon_i.epsilon_multiply(i - 1, 0, epsilon, 1, 0);
1172     accumulator += epsilon_i.epsilon_multiply(
1173         i,
1174         0,
1175         cr.apply_derivatives_nonhorner(
1176             order - i,
1177             [&f, i](auto... indices) { return f(i, static_cast<std::size_t>(indices)...); },
1178             std::forward<Fvars>(fvars)...) /
1179             factorial<root_type>(static_cast<unsigned>(i)),
1180         0,
1181         0);
1182   }
1183   return accumulator;
1184 }
1185 #endif
1186
1187 // f : order -> derivative(order)
1188 template <typename RealType, size_t Order>
1189 template <typename Func>
1190 fvar<RealType, Order> fvar<RealType, Order>::apply_derivatives_nonhorner(size_t const order,
1191                                                                          Func const& f) const {
1192   fvar<RealType, Order> const epsilon = fvar<RealType, Order>(*this).set_root(0);
1193   fvar<RealType, Order> epsilon_i = fvar<RealType, Order>(1);  // epsilon to the power of i
1194   fvar<RealType, Order> accumulator = fvar<RealType, Order>(f(0u));
1195 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1196   size_t const i_max = (std::min)(order, order_sum);
1197 #else  // ODR-use of static constexpr
1198   size_t const i_max = order < order_sum ? order : order_sum;
1199 #endif
1200   for (size_t i = 1; i <= i_max; ++i) {
1201     epsilon_i = epsilon_i.epsilon_multiply(i - 1, 0, epsilon, 1, 0);
1202     accumulator += epsilon_i.epsilon_multiply(i, 0, f(i) / factorial<root_type>(static_cast<unsigned>(i)));
1203   }
1204   return accumulator;
1205 }
1206
1207 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1208 // Can throw "std::out_of_range: array::at: __n (which is 7) >= _Nm (which is 7)"
1209 template <typename RealType, size_t Order>
1210 template <typename... Orders>
1211 get_type_at<RealType, sizeof...(Orders)> fvar<RealType, Order>::at(size_t order, Orders... orders) const {
1212   if constexpr (0 < sizeof...(Orders))
1213     return v.at(order).at(static_cast<std::size_t>(orders)...);
1214   else
1215     return v.at(order);
1216 }
1217 #endif
1218
1219 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1220 // Can throw "std::out_of_range: array::at: __n (which is 7) >= _Nm (which is 7)"
1221 template <typename RealType, size_t Order>
1222 template <typename... Orders>
1223 get_type_at<fvar<RealType, Order>, sizeof...(Orders)> fvar<RealType, Order>::derivative(
1224     Orders... orders) const {
1225   static_assert(sizeof...(Orders) <= depth,
1226                 "Number of parameters to derivative(...) cannot exceed fvar::depth.");
1227   return at(static_cast<std::size_t>(orders)...) *
1228          (... * factorial<root_type>(static_cast<unsigned>(orders)));
1229 }
1230 #endif
1231
1232 template <typename RealType, size_t Order>
1233 const RealType& fvar<RealType, Order>::operator[](size_t i) const {
1234   return v[i];
1235 }
1236
1237 template <typename RealType, size_t Order>
1238 RealType fvar<RealType, Order>::epsilon_inner_product(size_t z0,
1239                                                       size_t const isum0,
1240                                                       size_t const m0,
1241                                                       fvar<RealType, Order> const& cr,
1242                                                       size_t z1,
1243                                                       size_t const isum1,
1244                                                       size_t const m1,
1245                                                       size_t const j) const {
1246   static_assert(is_fvar<RealType>::value, "epsilon_inner_product() must have 1 < depth.");
1247   RealType accumulator = RealType();
1248   auto const i0_max = m1 < j ? j - m1 : 0;
1249   for (auto i0 = m0, i1 = j - m0; i0 <= i0_max; ++i0, --i1)
1250     accumulator += v[i0].epsilon_multiply(z0, isum0 + i0, cr.v[i1], z1, isum1 + i1);
1251   return accumulator;
1252 }
1253
1254 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1255 template <typename RealType, size_t Order>
1256 fvar<RealType, Order> fvar<RealType, Order>::epsilon_multiply(size_t z0,
1257                                                               size_t isum0,
1258                                                               fvar<RealType, Order> const& cr,
1259                                                               size_t z1,
1260                                                               size_t isum1) const {
1261   using diff_t = typename std::array<RealType, Order + 1>::difference_type;
1262   RealType const zero(0);
1263   size_t const m0 = order_sum + isum0 < Order + z0 ? Order + z0 - (order_sum + isum0) : 0;
1264   size_t const m1 = order_sum + isum1 < Order + z1 ? Order + z1 - (order_sum + isum1) : 0;
1265   size_t const i_max = m0 + m1 < Order ? Order - (m0 + m1) : 0;
1266   fvar<RealType, Order> retval = fvar<RealType, Order>();
1267   if constexpr (is_fvar<RealType>::value)
1268     for (size_t i = 0, j = Order; i <= i_max; ++i, --j)
1269       retval.v[j] = epsilon_inner_product(z0, isum0, m0, cr, z1, isum1, m1, j);
1270   else
1271     for (size_t i = 0, j = Order; i <= i_max; ++i, --j)
1272       retval.v[j] = std::inner_product(
1273           v.cbegin() + diff_t(m0), v.cend() - diff_t(i + m1), cr.v.crbegin() + diff_t(i + m0), zero);
1274   return retval;
1275 }
1276 #endif
1277
1278 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1279 // When called from outside this method, z0 should be non-zero. Otherwise if z0=0 then it will give an
1280 // incorrect result of 0 when the root value is 0 and ca=inf, when instead the correct product is nan.
1281 // If z0=0 then use the regular multiply operator*() instead.
1282 template <typename RealType, size_t Order>
1283 fvar<RealType, Order> fvar<RealType, Order>::epsilon_multiply(size_t z0,
1284                                                               size_t isum0,
1285                                                               root_type const& ca) const {
1286   fvar<RealType, Order> retval(*this);
1287   size_t const m0 = order_sum + isum0 < Order + z0 ? Order + z0 - (order_sum + isum0) : 0;
1288   if constexpr (is_fvar<RealType>::value)
1289     for (size_t i = m0; i <= Order; ++i)
1290       retval.v[i] = retval.v[i].epsilon_multiply(z0, isum0 + i, ca);
1291   else
1292     for (size_t i = m0; i <= Order; ++i)
1293       if (retval.v[i] != static_cast<RealType>(0))
1294         retval.v[i] *= ca;
1295   return retval;
1296 }
1297 #endif
1298
1299 template <typename RealType, size_t Order>
1300 fvar<RealType, Order> fvar<RealType, Order>::inverse() const {
1301   return static_cast<root_type>(*this) == 0 ? inverse_apply() : 1 / *this;
1302 }
1303
1304 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1305 template <typename RealType, size_t Order>
1306 fvar<RealType, Order>& fvar<RealType, Order>::negate() {
1307   if constexpr (is_fvar<RealType>::value)
1308     std::for_each(v.begin(), v.end(), [](RealType& r) { r.negate(); });
1309   else
1310     std::for_each(v.begin(), v.end(), [](RealType& a) { a = -a; });
1311   return *this;
1312 }
1313 #endif
1314
1315 // This gives log(0.0) = depth(1)(-inf,inf,-inf,inf,-inf,inf)
1316 // 1 / *this: log(0.0) = depth(1)(-inf,inf,-inf,-nan,-nan,-nan)
1317 template <typename RealType, size_t Order>
1318 fvar<RealType, Order> fvar<RealType, Order>::inverse_apply() const {
1319   root_type derivatives[order_sum + 1];  // LCOV_EXCL_LINE This causes a false negative on lcov coverage test.
1320   root_type const x0 = static_cast<root_type>(*this);
1321   *derivatives = 1 / x0;
1322   for (size_t i = 1; i <= order_sum; ++i)
1323     derivatives[i] = -derivatives[i - 1] * i / x0;
1324   return apply_derivatives_nonhorner(order_sum, [&derivatives](size_t j) { return derivatives[j]; });
1325 }
1326
1327 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1328 template <typename RealType, size_t Order>
1329 fvar<RealType, Order>& fvar<RealType, Order>::multiply_assign_by_root_type(bool is_root,
1330                                                                            root_type const& ca) {
1331   auto itr = v.begin();
1332   if constexpr (is_fvar<RealType>::value) {
1333     itr->multiply_assign_by_root_type(is_root, ca);
1334     for (++itr; itr != v.end(); ++itr)
1335       itr->multiply_assign_by_root_type(false, ca);
1336   } else {
1337     if (is_root || *itr != 0)
1338       *itr *= ca;  // Skip multiplication of 0 by ca=inf to avoid nan, except when is_root.
1339     for (++itr; itr != v.end(); ++itr)
1340       if (*itr != 0)
1341         *itr *= ca;
1342   }
1343   return *this;
1344 }
1345 #endif
1346
1347 template <typename RealType, size_t Order>
1348 fvar<RealType, Order>::operator root_type() const {
1349   return static_cast<root_type>(v.front());
1350 }
1351
1352 template <typename RealType, size_t Order>
1353 template <typename T, typename>
1354 fvar<RealType, Order>::operator T() const {
1355   return static_cast<T>(static_cast<root_type>(v.front()));
1356 }
1357
1358 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1359 template <typename RealType, size_t Order>
1360 fvar<RealType, Order>& fvar<RealType, Order>::set_root(root_type const& root) {
1361   if constexpr (is_fvar<RealType>::value)
1362     v.front().set_root(root);
1363   else
1364     v.front() = root;
1365   return *this;
1366 }
1367 #endif
1368
1369 // Standard Library Support Requirements
1370
1371 template <typename RealType, size_t Order>
1372 fvar<RealType, Order> fabs(fvar<RealType, Order> const& cr) {
1373   typename fvar<RealType, Order>::root_type const zero(0);
1374   return cr < zero ? -cr
1375                    : cr == zero ? fvar<RealType, Order>()  // Canonical fabs'(0) = 0.
1376                                 : cr;                      // Propagate NaN.
1377 }
1378
1379 template <typename RealType, size_t Order>
1380 fvar<RealType, Order> abs(fvar<RealType, Order> const& cr) {
1381   return fabs(cr);
1382 }
1383
1384 template <typename RealType, size_t Order>
1385 fvar<RealType, Order> ceil(fvar<RealType, Order> const& cr) {
1386   using std::ceil;
1387   return fvar<RealType, Order>(ceil(static_cast<typename fvar<RealType, Order>::root_type>(cr)));
1388 }
1389
1390 template <typename RealType, size_t Order>
1391 fvar<RealType, Order> floor(fvar<RealType, Order> const& cr) {
1392   using std::floor;
1393   return fvar<RealType, Order>(floor(static_cast<typename fvar<RealType, Order>::root_type>(cr)));
1394 }
1395
1396 template <typename RealType, size_t Order>
1397 fvar<RealType, Order> exp(fvar<RealType, Order> const& cr) {
1398   using std::exp;
1399   constexpr size_t order = fvar<RealType, Order>::order_sum;
1400   using root_type = typename fvar<RealType, Order>::root_type;
1401   root_type const d0 = exp(static_cast<root_type>(cr));
1402   return cr.apply_derivatives(order, [&d0](size_t) { return d0; });
1403 }
1404
1405 template <typename RealType, size_t Order>
1406 fvar<RealType, Order> pow(fvar<RealType, Order> const& x,
1407                           typename fvar<RealType, Order>::root_type const& y) {
1408   using std::pow;
1409   using root_type = typename fvar<RealType, Order>::root_type;
1410   constexpr size_t order = fvar<RealType, Order>::order_sum;
1411   root_type const x0 = static_cast<root_type>(x);
1412   root_type derivatives[order + 1]{pow(x0, y)};
1413   for (size_t i = 0; i < order && y - i != 0; ++i)
1414     derivatives[i + 1] = (y - i) * derivatives[i] / x0;
1415   return x.apply_derivatives(order, [&derivatives](size_t i) { return derivatives[i]; });
1416 }
1417
1418 template <typename RealType, size_t Order>
1419 fvar<RealType, Order> pow(typename fvar<RealType, Order>::root_type const& x,
1420                           fvar<RealType, Order> const& y) {
1421   BOOST_MATH_STD_USING
1422   using root_type = typename fvar<RealType, Order>::root_type;
1423   constexpr size_t order = fvar<RealType, Order>::order_sum;
1424   root_type const y0 = static_cast<root_type>(y);
1425   root_type derivatives[order + 1];
1426   *derivatives = pow(x, y0);
1427   root_type const logx = log(x);
1428   for (size_t i = 0; i < order; ++i)
1429     derivatives[i + 1] = derivatives[i] * logx;
1430   return y.apply_derivatives(order, [&derivatives](size_t i) { return derivatives[i]; });
1431 }
1432
1433 template <typename RealType1, size_t Order1, typename RealType2, size_t Order2>
1434 promote<fvar<RealType1, Order1>, fvar<RealType2, Order2>> pow(fvar<RealType1, Order1> const& x,
1435                                                               fvar<RealType2, Order2> const& y) {
1436   BOOST_MATH_STD_USING
1437   using return_type = promote<fvar<RealType1, Order1>, fvar<RealType2, Order2>>;
1438   using root_type = typename return_type::root_type;
1439   constexpr size_t order = return_type::order_sum;
1440   root_type const x0 = static_cast<root_type>(x);
1441   root_type const y0 = static_cast<root_type>(y);
1442   root_type dxydx[order + 1]{pow(x0, y0)};
1443   if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0)
1444     return return_type(*dxydx);
1445   else {
1446     for (size_t i = 0; i < order && y0 - i != 0; ++i)
1447       dxydx[i + 1] = (y0 - i) * dxydx[i] / x0;
1448     std::array<fvar<root_type, order>, order + 1> lognx;
1449     lognx.front() = fvar<root_type, order>(1);
1450 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1451     lognx[1] = log(make_fvar<root_type, order>(x0));
1452 #else  // for compilers that compile this branch when order=0.
1453     lognx[(std::min)(size_t(1), order)] = log(make_fvar<root_type, order>(x0));
1454 #endif
1455     for (size_t i = 1; i < order; ++i)
1456       lognx[i + 1] = lognx[i] * lognx[1];
1457     auto const f = [&dxydx, &lognx](size_t i, size_t j) {
1458       size_t binomial = 1;
1459       root_type sum = dxydx[i] * static_cast<root_type>(lognx[j]);
1460       for (size_t k = 1; k <= i; ++k) {
1461         (binomial *= (i - k + 1)) /= k;  // binomial_coefficient(i,k)
1462         sum += binomial * dxydx[i - k] * lognx[j].derivative(k);
1463       }
1464       return sum;
1465     };
1466     if (fabs(x0) < std::numeric_limits<root_type>::epsilon())
1467       return x.apply_derivatives_nonhorner(order, f, y);
1468     return x.apply_derivatives(order, f, y);
1469   }
1470 }
1471
1472 template <typename RealType, size_t Order>
1473 fvar<RealType, Order> sqrt(fvar<RealType, Order> const& cr) {
1474   using std::sqrt;
1475   using root_type = typename fvar<RealType, Order>::root_type;
1476   constexpr size_t order = fvar<RealType, Order>::order_sum;
1477   root_type derivatives[order + 1];
1478   root_type const x = static_cast<root_type>(cr);
1479   *derivatives = sqrt(x);
1480   if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0)
1481     return fvar<RealType, Order>(*derivatives);
1482   else {
1483     root_type numerator = 0.5;
1484     root_type powers = 1;
1485 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1486     derivatives[1] = numerator / *derivatives;
1487 #else  // for compilers that compile this branch when order=0.
1488     derivatives[(std::min)(size_t(1), order)] = numerator / *derivatives;
1489 #endif
1490     using diff_t = typename std::array<RealType, Order + 1>::difference_type;
1491     for (size_t i = 2; i <= order; ++i) {
1492       numerator *= static_cast<root_type>(-0.5) * ((static_cast<diff_t>(i) << 1) - 3);
1493       powers *= x;
1494       derivatives[i] = numerator / (powers * *derivatives);
1495     }
1496     auto const f = [&derivatives](size_t i) { return derivatives[i]; };
1497     if (cr < std::numeric_limits<root_type>::epsilon())
1498       return cr.apply_derivatives_nonhorner(order, f);
1499     return cr.apply_derivatives(order, f);
1500   }
1501 }
1502
1503 // Natural logarithm. If cr==0 then derivative(i) may have nans due to nans from inverse().
1504 template <typename RealType, size_t Order>
1505 fvar<RealType, Order> log(fvar<RealType, Order> const& cr) {
1506   using std::log;
1507   using root_type = typename fvar<RealType, Order>::root_type;
1508   constexpr size_t order = fvar<RealType, Order>::order_sum;
1509   root_type const d0 = log(static_cast<root_type>(cr));
1510   if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0)
1511     return fvar<RealType, Order>(d0);
1512   else {
1513     auto const d1 = make_fvar<root_type, order - 1>(static_cast<root_type>(cr)).inverse();  // log'(x) = 1 / x
1514     return cr.apply_coefficients_nonhorner(order, [&d0, &d1](size_t i) { return i ? d1[i - 1] / i : d0; });
1515   }
1516 }
1517
1518 template <typename RealType, size_t Order>
1519 fvar<RealType, Order> frexp(fvar<RealType, Order> const& cr, int* exp) {
1520   using std::exp2;
1521   using std::frexp;
1522   using root_type = typename fvar<RealType, Order>::root_type;
1523   frexp(static_cast<root_type>(cr), exp);
1524   return cr * static_cast<root_type>(exp2(-*exp));
1525 }
1526
1527 template <typename RealType, size_t Order>
1528 fvar<RealType, Order> ldexp(fvar<RealType, Order> const& cr, int exp) {
1529   // argument to std::exp2 must be casted to root_type, otherwise std::exp2 returns double (always)
1530   using std::exp2;
1531   return cr * exp2(static_cast<typename fvar<RealType, Order>::root_type>(exp));
1532 }
1533
1534 template <typename RealType, size_t Order>
1535 fvar<RealType, Order> cos(fvar<RealType, Order> const& cr) {
1536   BOOST_MATH_STD_USING
1537   using root_type = typename fvar<RealType, Order>::root_type;
1538   constexpr size_t order = fvar<RealType, Order>::order_sum;
1539   root_type const d0 = cos(static_cast<root_type>(cr));
1540   if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0)
1541     return fvar<RealType, Order>(d0);
1542   else {
1543     root_type const d1 = -sin(static_cast<root_type>(cr));
1544     root_type const derivatives[4]{d0, d1, -d0, -d1};
1545     return cr.apply_derivatives(order, [&derivatives](size_t i) { return derivatives[i & 3]; });
1546   }
1547 }
1548
1549 template <typename RealType, size_t Order>
1550 fvar<RealType, Order> sin(fvar<RealType, Order> const& cr) {
1551   BOOST_MATH_STD_USING
1552   using root_type = typename fvar<RealType, Order>::root_type;
1553   constexpr size_t order = fvar<RealType, Order>::order_sum;
1554   root_type const d0 = sin(static_cast<root_type>(cr));
1555   if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0)
1556     return fvar<RealType, Order>(d0);
1557   else {
1558     root_type const d1 = cos(static_cast<root_type>(cr));
1559     root_type const derivatives[4]{d0, d1, -d0, -d1};
1560     return cr.apply_derivatives(order, [&derivatives](size_t i) { return derivatives[i & 3]; });
1561   }
1562 }
1563
1564 template <typename RealType, size_t Order>
1565 fvar<RealType, Order> asin(fvar<RealType, Order> const& cr) {
1566   using std::asin;
1567   using root_type = typename fvar<RealType, Order>::root_type;
1568   constexpr size_t order = fvar<RealType, Order>::order_sum;
1569   root_type const d0 = asin(static_cast<root_type>(cr));
1570   if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0)
1571     return fvar<RealType, Order>(d0);
1572   else {
1573     auto x = make_fvar<root_type, order - 1>(static_cast<root_type>(cr));
1574     auto const d1 = sqrt((x *= x).negate() += 1).inverse();  // asin'(x) = 1 / sqrt(1-x*x).
1575     return cr.apply_coefficients_nonhorner(order, [&d0, &d1](size_t i) { return i ? d1[i - 1] / i : d0; });
1576   }
1577 }
1578
1579 template <typename RealType, size_t Order>
1580 fvar<RealType, Order> tan(fvar<RealType, Order> const& cr) {
1581   using std::tan;
1582   using root_type = typename fvar<RealType, Order>::root_type;
1583   constexpr size_t order = fvar<RealType, Order>::order_sum;
1584   root_type const d0 = tan(static_cast<root_type>(cr));
1585   if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0)
1586     return fvar<RealType, Order>(d0);
1587   else {
1588     auto c = cos(make_fvar<root_type, order - 1>(static_cast<root_type>(cr)));
1589     auto const d1 = (c *= c).inverse();  // tan'(x) = 1 / cos(x)^2
1590     return cr.apply_coefficients_nonhorner(order, [&d0, &d1](size_t i) { return i ? d1[i - 1] / i : d0; });
1591   }
1592 }
1593
1594 template <typename RealType, size_t Order>
1595 fvar<RealType, Order> atan(fvar<RealType, Order> const& cr) {
1596   using std::atan;
1597   using root_type = typename fvar<RealType, Order>::root_type;
1598   constexpr size_t order = fvar<RealType, Order>::order_sum;
1599   root_type const d0 = atan(static_cast<root_type>(cr));
1600   if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0)
1601     return fvar<RealType, Order>(d0);
1602   else {
1603     auto x = make_fvar<root_type, order - 1>(static_cast<root_type>(cr));
1604     auto const d1 = ((x *= x) += 1).inverse();  // atan'(x) = 1 / (x*x+1).
1605     return cr.apply_coefficients(order, [&d0, &d1](size_t i) { return i ? d1[i - 1] / i : d0; });
1606   }
1607 }
1608
1609 template <typename RealType, size_t Order>
1610 fvar<RealType, Order> atan2(fvar<RealType, Order> const& cr,
1611                             typename fvar<RealType, Order>::root_type const& ca) {
1612   using std::atan2;
1613   using root_type = typename fvar<RealType, Order>::root_type;
1614   constexpr size_t order = fvar<RealType, Order>::order_sum;
1615   root_type const d0 = atan2(static_cast<root_type>(cr), ca);
1616   if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0)
1617     return fvar<RealType, Order>(d0);
1618   else {
1619     auto y = make_fvar<root_type, order - 1>(static_cast<root_type>(cr));
1620     auto const d1 = ca / ((y *= y) += (ca * ca));  // (d/dy)atan2(y,x) = x / (y*y+x*x)
1621     return cr.apply_coefficients(order, [&d0, &d1](size_t i) { return i ? d1[i - 1] / i : d0; });
1622   }
1623 }
1624
1625 template <typename RealType, size_t Order>
1626 fvar<RealType, Order> atan2(typename fvar<RealType, Order>::root_type const& ca,
1627                             fvar<RealType, Order> const& cr) {
1628   using std::atan2;
1629   using root_type = typename fvar<RealType, Order>::root_type;
1630   constexpr size_t order = fvar<RealType, Order>::order_sum;
1631   root_type const d0 = atan2(ca, static_cast<root_type>(cr));
1632   if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0)
1633     return fvar<RealType, Order>(d0);
1634   else {
1635     auto x = make_fvar<root_type, order - 1>(static_cast<root_type>(cr));
1636     auto const d1 = -ca / ((x *= x) += (ca * ca));  // (d/dx)atan2(y,x) = -y / (x*x+y*y)
1637     return cr.apply_coefficients(order, [&d0, &d1](size_t i) { return i ? d1[i - 1] / i : d0; });
1638   }
1639 }
1640
1641 template <typename RealType1, size_t Order1, typename RealType2, size_t Order2>
1642 promote<fvar<RealType1, Order1>, fvar<RealType2, Order2>> atan2(fvar<RealType1, Order1> const& cr1,
1643                                                                 fvar<RealType2, Order2> const& cr2) {
1644   using std::atan2;
1645   using return_type = promote<fvar<RealType1, Order1>, fvar<RealType2, Order2>>;
1646   using root_type = typename return_type::root_type;
1647   constexpr size_t order = return_type::order_sum;
1648   root_type const y = static_cast<root_type>(cr1);
1649   root_type const x = static_cast<root_type>(cr2);
1650   root_type const d00 = atan2(y, x);
1651   if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0)
1652     return return_type(d00);
1653   else {
1654     constexpr size_t order1 = fvar<RealType1, Order1>::order_sum;
1655     constexpr size_t order2 = fvar<RealType2, Order2>::order_sum;
1656     auto x01 = make_fvar<typename fvar<RealType2, Order2>::root_type, order2 - 1>(x);
1657     auto const d01 = -y / ((x01 *= x01) += (y * y));
1658     auto y10 = make_fvar<typename fvar<RealType1, Order1>::root_type, order1 - 1>(y);
1659     auto x10 = make_fvar<typename fvar<RealType2, Order2>::root_type, 0, order2>(x);
1660     auto const d10 = x10 / ((x10 * x10) + (y10 *= y10));
1661     auto const f = [&d00, &d01, &d10](size_t i, size_t j) {
1662       return i ? d10[i - 1][j] / i : j ? d01[j - 1] / j : d00;
1663     };
1664     return cr1.apply_coefficients(order, f, cr2);
1665   }
1666 }
1667
1668 template <typename RealType1, size_t Order1, typename RealType2, size_t Order2>
1669 promote<fvar<RealType1, Order1>, fvar<RealType2, Order2>> fmod(fvar<RealType1, Order1> const& cr1,
1670                                                                fvar<RealType2, Order2> const& cr2) {
1671   using boost::math::trunc;
1672   auto const numer = static_cast<typename fvar<RealType1, Order1>::root_type>(cr1);
1673   auto const denom = static_cast<typename fvar<RealType2, Order2>::root_type>(cr2);
1674   return cr1 - cr2 * trunc(numer / denom);
1675 }
1676
1677 template <typename RealType, size_t Order>
1678 fvar<RealType, Order> round(fvar<RealType, Order> const& cr) {
1679   using boost::math::round;
1680   return fvar<RealType, Order>(round(static_cast<typename fvar<RealType, Order>::root_type>(cr)));
1681 }
1682
1683 template <typename RealType, size_t Order>
1684 int iround(fvar<RealType, Order> const& cr) {
1685   using boost::math::iround;
1686   return iround(static_cast<typename fvar<RealType, Order>::root_type>(cr));
1687 }
1688
1689 template <typename RealType, size_t Order>
1690 long lround(fvar<RealType, Order> const& cr) {
1691   using boost::math::lround;
1692   return lround(static_cast<typename fvar<RealType, Order>::root_type>(cr));
1693 }
1694
1695 template <typename RealType, size_t Order>
1696 long long llround(fvar<RealType, Order> const& cr) {
1697   using boost::math::llround;
1698   return llround(static_cast<typename fvar<RealType, Order>::root_type>(cr));
1699 }
1700
1701 template <typename RealType, size_t Order>
1702 fvar<RealType, Order> trunc(fvar<RealType, Order> const& cr) {
1703   using boost::math::trunc;
1704   return fvar<RealType, Order>(trunc(static_cast<typename fvar<RealType, Order>::root_type>(cr)));
1705 }
1706
1707 template <typename RealType, size_t Order>
1708 long double truncl(fvar<RealType, Order> const& cr) {
1709   using std::truncl;
1710   return truncl(static_cast<typename fvar<RealType, Order>::root_type>(cr));
1711 }
1712
1713 template <typename RealType, size_t Order>
1714 int itrunc(fvar<RealType, Order> const& cr) {
1715   using boost::math::itrunc;
1716   return itrunc(static_cast<typename fvar<RealType, Order>::root_type>(cr));
1717 }
1718
1719 template <typename RealType, size_t Order>
1720 long long lltrunc(fvar<RealType, Order> const& cr) {
1721   using boost::math::lltrunc;
1722   return lltrunc(static_cast<typename fvar<RealType, Order>::root_type>(cr));
1723 }
1724
1725 template <typename RealType, size_t Order>
1726 std::ostream& operator<<(std::ostream& out, fvar<RealType, Order> const& cr) {
1727   out << "depth(" << cr.depth << ")(" << cr.v.front();
1728   for (size_t i = 1; i <= Order; ++i)
1729     out << ',' << cr.v[i];
1730   return out << ')';
1731 }
1732
1733 // Additional functions
1734
1735 template <typename RealType, size_t Order>
1736 fvar<RealType, Order> acos(fvar<RealType, Order> const& cr) {
1737   using std::acos;
1738   using root_type = typename fvar<RealType, Order>::root_type;
1739   constexpr size_t order = fvar<RealType, Order>::order_sum;
1740   root_type const d0 = acos(static_cast<root_type>(cr));
1741   if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0)
1742     return fvar<RealType, Order>(d0);
1743   else {
1744     auto x = make_fvar<root_type, order - 1>(static_cast<root_type>(cr));
1745     auto const d1 = sqrt((x *= x).negate() += 1).inverse().negate();  // acos'(x) = -1 / sqrt(1-x*x).
1746     return cr.apply_coefficients(order, [&d0, &d1](size_t i) { return i ? d1[i - 1] / i : d0; });
1747   }
1748 }
1749
1750 template <typename RealType, size_t Order>
1751 fvar<RealType, Order> acosh(fvar<RealType, Order> const& cr) {
1752   using boost::math::acosh;
1753   using root_type = typename fvar<RealType, Order>::root_type;
1754   constexpr size_t order = fvar<RealType, Order>::order_sum;
1755   root_type const d0 = acosh(static_cast<root_type>(cr));
1756   if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0)
1757     return fvar<RealType, Order>(d0);
1758   else {
1759     auto x = make_fvar<root_type, order - 1>(static_cast<root_type>(cr));
1760     auto const d1 = sqrt((x *= x) -= 1).inverse();  // acosh'(x) = 1 / sqrt(x*x-1).
1761     return cr.apply_coefficients(order, [&d0, &d1](size_t i) { return i ? d1[i - 1] / i : d0; });
1762   }
1763 }
1764
1765 template <typename RealType, size_t Order>
1766 fvar<RealType, Order> asinh(fvar<RealType, Order> const& cr) {
1767   using boost::math::asinh;
1768   using root_type = typename fvar<RealType, Order>::root_type;
1769   constexpr size_t order = fvar<RealType, Order>::order_sum;
1770   root_type const d0 = asinh(static_cast<root_type>(cr));
1771   if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0)
1772     return fvar<RealType, Order>(d0);
1773   else {
1774     auto x = make_fvar<root_type, order - 1>(static_cast<root_type>(cr));
1775     auto const d1 = sqrt((x *= x) += 1).inverse();  // asinh'(x) = 1 / sqrt(x*x+1).
1776     return cr.apply_coefficients(order, [&d0, &d1](size_t i) { return i ? d1[i - 1] / i : d0; });
1777   }
1778 }
1779
1780 template <typename RealType, size_t Order>
1781 fvar<RealType, Order> atanh(fvar<RealType, Order> const& cr) {
1782   using boost::math::atanh;
1783   using root_type = typename fvar<RealType, Order>::root_type;
1784   constexpr size_t order = fvar<RealType, Order>::order_sum;
1785   root_type const d0 = atanh(static_cast<root_type>(cr));
1786   if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0)
1787     return fvar<RealType, Order>(d0);
1788   else {
1789     auto x = make_fvar<root_type, order - 1>(static_cast<root_type>(cr));
1790     auto const d1 = ((x *= x).negate() += 1).inverse();  // atanh'(x) = 1 / (1-x*x)
1791     return cr.apply_coefficients(order, [&d0, &d1](size_t i) { return i ? d1[i - 1] / i : d0; });
1792   }
1793 }
1794
1795 template <typename RealType, size_t Order>
1796 fvar<RealType, Order> cosh(fvar<RealType, Order> const& cr) {
1797   BOOST_MATH_STD_USING
1798   using root_type = typename fvar<RealType, Order>::root_type;
1799   constexpr size_t order = fvar<RealType, Order>::order_sum;
1800   root_type const d0 = cosh(static_cast<root_type>(cr));
1801   if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0)
1802     return fvar<RealType, Order>(d0);
1803   else {
1804     root_type const derivatives[2]{d0, sinh(static_cast<root_type>(cr))};
1805     return cr.apply_derivatives(order, [&derivatives](size_t i) { return derivatives[i & 1]; });
1806   }
1807 }
1808
1809 template <typename RealType, size_t Order>
1810 fvar<RealType, Order> digamma(fvar<RealType, Order> const& cr) {
1811   using boost::math::digamma;
1812   using root_type = typename fvar<RealType, Order>::root_type;
1813   constexpr size_t order = fvar<RealType, Order>::order_sum;
1814   root_type const x = static_cast<root_type>(cr);
1815   root_type const d0 = digamma(x);
1816   if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0)
1817     return fvar<RealType, Order>(d0);
1818   else {
1819     static_assert(order <= static_cast<size_t>((std::numeric_limits<int>::max)()),
1820                   "order exceeds maximum derivative for boost::math::polygamma().");
1821     return cr.apply_derivatives(
1822         order, [&x, &d0](size_t i) { return i ? boost::math::polygamma(static_cast<int>(i), x) : d0; });
1823   }
1824 }
1825
1826 template <typename RealType, size_t Order>
1827 fvar<RealType, Order> erf(fvar<RealType, Order> const& cr) {
1828   using boost::math::erf;
1829   using root_type = typename fvar<RealType, Order>::root_type;
1830   constexpr size_t order = fvar<RealType, Order>::order_sum;
1831   root_type const d0 = erf(static_cast<root_type>(cr));
1832   if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0)
1833     return fvar<RealType, Order>(d0);
1834   else {
1835     auto x = make_fvar<root_type, order - 1>(static_cast<root_type>(cr));  // d1 = 2/sqrt(pi)*exp(-x*x)
1836     auto const d1 = 2 * constants::one_div_root_pi<root_type>() * exp((x *= x).negate());
1837     return cr.apply_coefficients(order, [&d0, &d1](size_t i) { return i ? d1[i - 1] / i : d0; });
1838   }
1839 }
1840
1841 template <typename RealType, size_t Order>
1842 fvar<RealType, Order> erfc(fvar<RealType, Order> const& cr) {
1843   using boost::math::erfc;
1844   using root_type = typename fvar<RealType, Order>::root_type;
1845   constexpr size_t order = fvar<RealType, Order>::order_sum;
1846   root_type const d0 = erfc(static_cast<root_type>(cr));
1847   if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0)
1848     return fvar<RealType, Order>(d0);
1849   else {
1850     auto x = make_fvar<root_type, order - 1>(static_cast<root_type>(cr));  // erfc'(x) = -erf'(x)
1851     auto const d1 = -2 * constants::one_div_root_pi<root_type>() * exp((x *= x).negate());
1852     return cr.apply_coefficients(order, [&d0, &d1](size_t i) { return i ? d1[i - 1] / i : d0; });
1853   }
1854 }
1855
1856 template <typename RealType, size_t Order>
1857 fvar<RealType, Order> lambert_w0(fvar<RealType, Order> const& cr) {
1858   using std::exp;
1859   using boost::math::lambert_w0;
1860   using root_type = typename fvar<RealType, Order>::root_type;
1861   constexpr size_t order = fvar<RealType, Order>::order_sum;
1862   root_type derivatives[order + 1];
1863   *derivatives = lambert_w0(static_cast<root_type>(cr));
1864   if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0)
1865     return fvar<RealType, Order>(*derivatives);
1866   else {
1867     root_type const expw = exp(*derivatives);
1868     derivatives[1] = 1 / (static_cast<root_type>(cr) + expw);
1869     if BOOST_AUTODIFF_IF_CONSTEXPR (order == 1)
1870       return cr.apply_derivatives_nonhorner(order, [&derivatives](size_t i) { return derivatives[i]; });
1871     else {
1872       using diff_t = typename std::array<RealType, Order + 1>::difference_type;
1873       root_type d1powers = derivatives[1] * derivatives[1];
1874       root_type const x = derivatives[1] * expw;
1875       derivatives[2] = d1powers * (-1 - x);
1876       std::array<root_type, order> coef{{-1, -1}};  // as in derivatives[2].
1877       for (size_t n = 3; n <= order; ++n) {
1878         coef[n - 1] = coef[n - 2] * -static_cast<root_type>(2 * n - 3);
1879         for (size_t j = n - 2; j != 0; --j)
1880           (coef[j] *= -static_cast<root_type>(n - 1)) -= (n + j - 2) * coef[j - 1];
1881         coef[0] *= -static_cast<root_type>(n - 1);
1882         d1powers *= derivatives[1];
1883         derivatives[n] =
1884             d1powers * std::accumulate(coef.crend() - diff_t(n - 1),
1885                                        coef.crend(),
1886                                        coef[n - 1],
1887                                        [&x](root_type const& a, root_type const& b) { return a * x + b; });
1888       }
1889       return cr.apply_derivatives_nonhorner(order, [&derivatives](size_t i) { return derivatives[i]; });
1890     }
1891   }
1892 }
1893
1894 template <typename RealType, size_t Order>
1895 fvar<RealType, Order> lgamma(fvar<RealType, Order> const& cr) {
1896   using std::lgamma;
1897   using root_type = typename fvar<RealType, Order>::root_type;
1898   constexpr size_t order = fvar<RealType, Order>::order_sum;
1899   root_type const x = static_cast<root_type>(cr);
1900   root_type const d0 = lgamma(x);
1901   if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0)
1902     return fvar<RealType, Order>(d0);
1903   else {
1904     static_assert(order <= static_cast<size_t>((std::numeric_limits<int>::max)()) + 1,
1905                   "order exceeds maximum derivative for boost::math::polygamma().");
1906     return cr.apply_derivatives(
1907         order, [&x, &d0](size_t i) { return i ? boost::math::polygamma(static_cast<int>(i - 1), x) : d0; });
1908   }
1909 }
1910
1911 template <typename RealType, size_t Order>
1912 fvar<RealType, Order> sinc(fvar<RealType, Order> const& cr) {
1913   if (cr != 0)
1914     return sin(cr) / cr;
1915   using root_type = typename fvar<RealType, Order>::root_type;
1916   constexpr size_t order = fvar<RealType, Order>::order_sum;
1917   root_type taylor[order + 1]{1};  // sinc(0) = 1
1918   if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0)
1919     return fvar<RealType, Order>(*taylor);
1920   else {
1921     for (size_t n = 2; n <= order; n += 2)
1922       taylor[n] = (1 - static_cast<int>(n & 2)) / factorial<root_type>(static_cast<unsigned>(n + 1));
1923     return cr.apply_coefficients_nonhorner(order, [&taylor](size_t i) { return taylor[i]; });
1924   }
1925 }
1926
1927 template <typename RealType, size_t Order>
1928 fvar<RealType, Order> sinh(fvar<RealType, Order> const& cr) {
1929   BOOST_MATH_STD_USING
1930   using root_type = typename fvar<RealType, Order>::root_type;
1931   constexpr size_t order = fvar<RealType, Order>::order_sum;
1932   root_type const d0 = sinh(static_cast<root_type>(cr));
1933   if BOOST_AUTODIFF_IF_CONSTEXPR (fvar<RealType, Order>::order_sum == 0)
1934     return fvar<RealType, Order>(d0);
1935   else {
1936     root_type const derivatives[2]{d0, cosh(static_cast<root_type>(cr))};
1937     return cr.apply_derivatives(order, [&derivatives](size_t i) { return derivatives[i & 1]; });
1938   }
1939 }
1940
1941 template <typename RealType, size_t Order>
1942 fvar<RealType, Order> tanh(fvar<RealType, Order> const& cr) {
1943   fvar<RealType, Order> retval = exp(cr * 2);
1944   fvar<RealType, Order> const denom = retval + 1;
1945   (retval -= 1) /= denom;
1946   return retval;
1947 }
1948
1949 template <typename RealType, size_t Order>
1950 fvar<RealType, Order> tgamma(fvar<RealType, Order> const& cr) {
1951   using std::tgamma;
1952   using root_type = typename fvar<RealType, Order>::root_type;
1953   constexpr size_t order = fvar<RealType, Order>::order_sum;
1954   if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0)
1955     return fvar<RealType, Order>(tgamma(static_cast<root_type>(cr)));
1956   else {
1957     if (cr < 0)
1958       return constants::pi<root_type>() / (sin(constants::pi<root_type>() * cr) * tgamma(1 - cr));
1959     return exp(lgamma(cr)).set_root(tgamma(static_cast<root_type>(cr)));
1960   }
1961 }
1962
1963 }  // namespace detail
1964 }  // namespace autodiff_v1
1965 }  // namespace differentiation
1966 }  // namespace math
1967 }  // namespace boost
1968
1969 namespace std {
1970
1971 // boost::math::tools::digits<RealType>() is handled by this std::numeric_limits<> specialization,
1972 // and similarly for max_value, min_value, log_max_value, log_min_value, and epsilon.
1973 template <typename RealType, size_t Order>
1974 class numeric_limits<boost::math::differentiation::detail::fvar<RealType, Order>>
1975     : public numeric_limits<typename boost::math::differentiation::detail::fvar<RealType, Order>::root_type> {
1976 };
1977
1978 }  // namespace std
1979
1980 namespace boost {
1981 namespace math {
1982 namespace tools {
1983 namespace detail {
1984
1985 template <typename RealType, std::size_t Order>
1986 using autodiff_fvar_type = differentiation::detail::fvar<RealType, Order>;
1987
1988 template <typename RealType, std::size_t Order>
1989 using autodiff_root_type = typename autodiff_fvar_type<RealType, Order>::root_type;
1990 }  // namespace detail
1991
1992 // See boost/math/tools/promotion.hpp
1993 template <typename RealType0, size_t Order0, typename RealType1, size_t Order1>
1994 struct promote_args_2<detail::autodiff_fvar_type<RealType0, Order0>,
1995                       detail::autodiff_fvar_type<RealType1, Order1>> {
1996   using type = detail::autodiff_fvar_type<typename promote_args_2<RealType0, RealType1>::type,
1997 #ifndef BOOST_NO_CXX14_CONSTEXPR
1998                                           (std::max)(Order0, Order1)>;
1999 #else
2000         Order0<Order1 ? Order1 : Order0>;
2001 #endif
2002 };
2003
2004 template <typename RealType, size_t Order>
2005 struct promote_args<detail::autodiff_fvar_type<RealType, Order>> {
2006   using type = detail::autodiff_fvar_type<typename promote_args<RealType>::type, Order>;
2007 };
2008
2009 template <typename RealType0, size_t Order0, typename RealType1>
2010 struct promote_args_2<detail::autodiff_fvar_type<RealType0, Order0>, RealType1> {
2011   using type = detail::autodiff_fvar_type<typename promote_args_2<RealType0, RealType1>::type, Order0>;
2012 };
2013
2014 template <typename RealType0, typename RealType1, size_t Order1>
2015 struct promote_args_2<RealType0, detail::autodiff_fvar_type<RealType1, Order1>> {
2016   using type = detail::autodiff_fvar_type<typename promote_args_2<RealType0, RealType1>::type, Order1>;
2017 };
2018
2019 template <typename destination_t, typename RealType, std::size_t Order>
2020 inline BOOST_MATH_CONSTEXPR destination_t real_cast(detail::autodiff_fvar_type<RealType, Order> const& from_v)
2021     BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(destination_t) && BOOST_MATH_IS_FLOAT(RealType)) {
2022   return real_cast<destination_t>(static_cast<detail::autodiff_root_type<RealType, Order>>(from_v));
2023 }
2024
2025 }  // namespace tools
2026
2027 namespace policies {
2028
2029 template <class Policy, std::size_t Order>
2030 using fvar_t = differentiation::detail::fvar<Policy, Order>;
2031 template <class Policy, std::size_t Order>
2032 struct evaluation<fvar_t<float, Order>, Policy> {
2033   using type = fvar_t<typename conditional<Policy::promote_float_type::value, double, float>::type, Order>;
2034 };
2035
2036 template <class Policy, std::size_t Order>
2037 struct evaluation<fvar_t<double, Order>, Policy> {
2038   using type =
2039       fvar_t<typename conditional<Policy::promote_double_type::value, long double, double>::type, Order>;
2040 };
2041
2042 }  // namespace policies
2043 }  // namespace math
2044 }  // namespace boost
2045
2046 #ifdef BOOST_NO_CXX17_IF_CONSTEXPR
2047 #include "autodiff_cpp11.hpp"
2048 #endif
2049
2050 #endif  // BOOST_MATH_DIFFERENTIATION_AUTODIFF_HPP