brintos / llvm-project-archived public Read only

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