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)
6 #ifndef BOOST_MATH_DIFFERENTIATION_AUTODIFF_HPP
7 #define BOOST_MATH_DIFFERENTIATION_AUTODIFF_HPP
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>
31 #include <type_traits>
35 namespace differentiation {
36 // Automatic Differentiation v1
37 inline namespace autodiff_v1 {
40 template <typename RealType, typename... RealTypes>
41 struct promote_args_n {
42 using type = typename tools::promote_args_2<RealType, typename promote_args_n<RealTypes...>::type>::type;
45 template <typename RealType>
46 struct promote_args_n<RealType> {
47 using type = typename tools::promote_arg<RealType>::type;
52 template <typename RealType, typename... RealTypes>
53 using promote = typename detail::promote_args_n<RealType, RealTypes...>::type;
57 template <typename RealType, size_t Order>
61 struct is_fvar_impl : std::false_type {};
63 template <typename RealType, size_t Order>
64 struct is_fvar_impl<fvar<RealType, Order>> : std::true_type {};
67 using is_fvar = is_fvar_impl<typename std::decay<T>::type>;
69 template <typename RealType, size_t Order, size_t... Orders>
71 using type = fvar<typename nest_fvar<RealType, Orders...>::type, Order>;
74 template <typename RealType, size_t Order>
75 struct nest_fvar<RealType, Order> {
76 using type = fvar<RealType, Order>;
80 struct get_depth_impl : std::integral_constant<size_t, 0> {};
82 template <typename RealType, size_t Order>
83 struct get_depth_impl<fvar<RealType, Order>>
84 : std::integral_constant<size_t, get_depth_impl<RealType>::value + 1> {};
87 using get_depth = get_depth_impl<typename std::decay<T>::type>;
90 struct get_order_sum_t : std::integral_constant<size_t, 0> {};
92 template <typename RealType, size_t Order>
93 struct get_order_sum_t<fvar<RealType, Order>>
94 : std::integral_constant<size_t, get_order_sum_t<RealType>::value + Order> {};
97 using get_order_sum = get_order_sum_t<typename std::decay<T>::type>;
99 template <typename RealType>
100 struct get_root_type {
101 using type = RealType;
104 template <typename RealType, size_t Order>
105 struct get_root_type<fvar<RealType, Order>> {
106 using type = typename get_root_type<RealType>::type;
109 template <typename RealType, size_t Depth>
111 using type = RealType;
114 template <typename RealType, size_t Order, size_t Depth>
115 struct 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;
121 template <typename RealType, size_t Depth>
122 using get_type_at = typename type_at<RealType, Depth>::type;
124 // Satisfies Boost's Conceptual Requirements for Real Number Types.
125 // https://www.boost.org/libs/math/doc/html/math_toolkit/real_concepts.html
126 template <typename RealType, size_t Order>
129 std::array<RealType, Order + 1> v;
132 using root_type = typename get_root_type<RealType>::type; // RealType in the root fvar<RealType,Order>.
136 // Initialize a variable or constant.
137 fvar(root_type const&, bool const is_variable);
139 // RealType(cr) | RealType | RealType is copy constructible.
140 fvar(fvar const&) = default;
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&);
146 // RealType(ca) | RealType | RealType is copy constructible from the arithmetic types.
147 explicit fvar(root_type const&); // Initialize a constant. (No epsilon terms.)
149 template <typename RealType2>
150 fvar(RealType2 const& ca); // Supports any RealType2 for which static_cast<root_type>(ca) compiles.
152 // r = cr | RealType& | Assignment operator.
153 fvar& operator=(fvar const&) = default;
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.
159 // r += cr | RealType& | Adds cr to r.
160 template <typename RealType2, size_t Order2>
161 fvar& operator+=(fvar<RealType2, Order2> const&);
163 // r += ca | RealType& | Adds ar to r.
164 fvar& operator+=(root_type const&);
166 // r -= cr | RealType& | Subtracts cr from r.
167 template <typename RealType2, size_t Order2>
168 fvar& operator-=(fvar<RealType2, Order2> const&);
170 // r -= ca | RealType& | Subtracts ca from r.
171 fvar& operator-=(root_type const&);
173 // r *= cr | RealType& | Multiplies r by cr.
174 template <typename RealType2, size_t Order2>
175 fvar& operator*=(fvar<RealType2, Order2> const&);
177 // r *= ca | RealType& | Multiplies r by ca.
178 fvar& operator*=(root_type const&);
180 // r /= cr | RealType& | Divides r by cr.
181 template <typename RealType2, size_t Order2>
182 fvar& operator/=(fvar<RealType2, Order2> const&);
184 // r /= ca | RealType& | Divides r by ca.
185 fvar& operator/=(root_type const&);
187 // -r | RealType | Unary Negation.
188 fvar operator-() const;
190 // +r | RealType& | Identity Operation.
191 fvar const& operator+() const;
193 // cr + cr2 | RealType | Binary Addition
194 template <typename RealType2, size_t Order2>
195 promote<fvar, fvar<RealType2, Order2>> operator+(fvar<RealType2, Order2> const&) const;
197 // cr + ca | RealType | Binary Addition
198 fvar operator+(root_type const&) const;
200 // ca + cr | RealType | Binary Addition
201 template <typename RealType2, size_t Order2>
202 friend fvar<RealType2, Order2> operator+(typename fvar<RealType2, Order2>::root_type const&,
203 fvar<RealType2, Order2> const&);
205 // cr - cr2 | RealType | Binary Subtraction
206 template <typename RealType2, size_t Order2>
207 promote<fvar, fvar<RealType2, Order2>> operator-(fvar<RealType2, Order2> const&) const;
209 // cr - ca | RealType | Binary Subtraction
210 fvar operator-(root_type const&) const;
212 // ca - cr | RealType | Binary Subtraction
213 template <typename RealType2, size_t Order2>
214 friend fvar<RealType2, Order2> operator-(typename fvar<RealType2, Order2>::root_type const&,
215 fvar<RealType2, Order2> const&);
217 // cr * cr2 | RealType | Binary Multiplication
218 template <typename RealType2, size_t Order2>
219 promote<fvar, fvar<RealType2, Order2>> operator*(fvar<RealType2, Order2> const&)const;
221 // cr * ca | RealType | Binary Multiplication
222 fvar operator*(root_type const&)const;
224 // ca * cr | RealType | Binary Multiplication
225 template <typename RealType2, size_t Order2>
226 friend fvar<RealType2, Order2> operator*(typename fvar<RealType2, Order2>::root_type const&,
227 fvar<RealType2, Order2> const&);
229 // cr / cr2 | RealType | Binary Subtraction
230 template <typename RealType2, size_t Order2>
231 promote<fvar, fvar<RealType2, Order2>> operator/(fvar<RealType2, Order2> const&) const;
233 // cr / ca | RealType | Binary Subtraction
234 fvar operator/(root_type const&) const;
236 // ca / cr | RealType | Binary Subtraction
237 template <typename RealType2, size_t Order2>
238 friend fvar<RealType2, Order2> operator/(typename fvar<RealType2, Order2>::root_type const&,
239 fvar<RealType2, Order2> const&);
241 // For all comparison overloads, only the root term is compared.
243 // cr == cr2 | bool | Equality Comparison
244 template <typename RealType2, size_t Order2>
245 bool operator==(fvar<RealType2, Order2> const&) const;
247 // cr == ca | bool | Equality Comparison
248 bool operator==(root_type const&) const;
250 // ca == cr | bool | Equality Comparison
251 template <typename RealType2, size_t Order2>
252 friend bool operator==(typename fvar<RealType2, Order2>::root_type const&, fvar<RealType2, Order2> const&);
254 // cr != cr2 | bool | Inequality Comparison
255 template <typename RealType2, size_t Order2>
256 bool operator!=(fvar<RealType2, Order2> const&) const;
258 // cr != ca | bool | Inequality Comparison
259 bool operator!=(root_type const&) const;
261 // ca != cr | bool | Inequality Comparison
262 template <typename RealType2, size_t Order2>
263 friend bool operator!=(typename fvar<RealType2, Order2>::root_type const&, fvar<RealType2, Order2> const&);
265 // cr <= cr2 | bool | Less than equal to.
266 template <typename RealType2, size_t Order2>
267 bool operator<=(fvar<RealType2, Order2> const&) const;
269 // cr <= ca | bool | Less than equal to.
270 bool operator<=(root_type const&) const;
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&);
276 // cr >= cr2 | bool | Greater than equal to.
277 template <typename RealType2, size_t Order2>
278 bool operator>=(fvar<RealType2, Order2> const&) const;
280 // cr >= ca | bool | Greater than equal to.
281 bool operator>=(root_type const&) const;
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&);
287 // cr < cr2 | bool | Less than comparison.
288 template <typename RealType2, size_t Order2>
289 bool operator<(fvar<RealType2, Order2> const&) const;
291 // cr < ca | bool | Less than comparison.
292 bool operator<(root_type const&) const;
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&);
298 // cr > cr2 | bool | Greater than comparison.
299 template <typename RealType2, size_t Order2>
300 bool operator>(fvar<RealType2, Order2> const&) const;
302 // cr > ca | bool | Greater than comparison.
303 bool operator>(root_type const&) const;
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&);
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;
313 template <typename... Orders>
314 get_type_at<fvar, sizeof...(Orders)> derivative(Orders... orders) const;
316 const RealType& operator[](size_t) const;
318 fvar inverse() const; // Multiplicative inverse.
320 fvar& negate(); // Negate and return reference to *this.
322 static constexpr size_t depth = get_depth<fvar>::value; // Number of nested std::array<RealType,Order>.
324 static constexpr size_t order_sum = get_order_sum<fvar>::value;
326 explicit operator root_type() const; // Must be explicit, otherwise overloaded operators are ambiguous.
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_if
331 fvar& set_root(root_type const&);
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,
338 Fvars&&... fvars) const;
340 template <typename Func>
341 fvar apply_coefficients(size_t const order, Func const& f) const;
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,
348 Fvars&&... fvars) const;
350 template <typename Func>
351 fvar apply_coefficients_nonhorner(size_t const order, Func const& f) const;
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,
358 Fvars&&... fvars) const;
360 template <typename Func>
361 fvar apply_derivatives(size_t const order, Func const& f) const;
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,
368 Fvars&&... fvars) const;
370 template <typename Func>
371 fvar apply_derivatives_nonhorner(size_t const order, Func const& f) const;
374 RealType epsilon_inner_product(size_t z0,
383 fvar epsilon_multiply(size_t z0, size_t isum0, fvar const& cr, size_t z1, size_t isum1) const;
385 fvar epsilon_multiply(size_t z0, size_t isum0, root_type const& ca) const;
387 fvar inverse_apply() const;
389 fvar& multiply_assign_by_root_type(bool is_root, root_type const&);
391 template <typename RealType2, size_t Orders2>
394 template <typename RealType2, size_t Order2>
395 friend std::ostream& operator<<(std::ostream&, fvar<RealType2, Order2> const&);
397 // C++11 Compatibility
398 #ifdef BOOST_NO_CXX17_IF_CONSTEXPR
399 template <typename RootType>
400 void fvar_cpp11(std::true_type, RootType const& ca, bool const is_variable);
402 template <typename RootType>
403 void fvar_cpp11(std::false_type, RootType const& ca, bool const is_variable);
405 template <typename... Orders>
406 get_type_at<RealType, sizeof...(Orders)> at_cpp11(std::true_type, size_t order, Orders... orders) const;
408 template <typename... Orders>
409 get_type_at<RealType, sizeof...(Orders)> at_cpp11(std::false_type, size_t order, Orders... orders) const;
411 template <typename SizeType>
412 fvar epsilon_multiply_cpp11(std::true_type,
419 template <typename SizeType>
420 fvar epsilon_multiply_cpp11(std::false_type,
427 template <typename SizeType>
428 fvar epsilon_multiply_cpp11(std::true_type, SizeType z0, size_t isum0, root_type const& ca) const;
430 template <typename SizeType>
431 fvar epsilon_multiply_cpp11(std::false_type, SizeType z0, size_t isum0, root_type const& ca) const;
433 template <typename RootType>
434 fvar& multiply_assign_by_root_type_cpp11(std::true_type, bool is_root, RootType const& ca);
436 template <typename RootType>
437 fvar& multiply_assign_by_root_type_cpp11(std::false_type, bool is_root, RootType const& ca);
439 template <typename RootType>
440 fvar& negate_cpp11(std::true_type, RootType const&);
442 template <typename RootType>
443 fvar& negate_cpp11(std::false_type, RootType const&);
445 template <typename RootType>
446 fvar& set_root_cpp11(std::true_type, RootType const& root);
448 template <typename RootType>
449 fvar& set_root_cpp11(std::false_type, RootType const& root);
453 // Standard Library Support Requirements
455 // fabs(cr1) | RealType
456 template <typename RealType, size_t Order>
457 fvar<RealType, Order> fabs(fvar<RealType, Order> const&);
459 // abs(cr1) | RealType
460 template <typename RealType, size_t Order>
461 fvar<RealType, Order> abs(fvar<RealType, Order> const&);
463 // ceil(cr1) | RealType
464 template <typename RealType, size_t Order>
465 fvar<RealType, Order> ceil(fvar<RealType, Order> const&);
467 // floor(cr1) | RealType
468 template <typename RealType, size_t Order>
469 fvar<RealType, Order> floor(fvar<RealType, Order> const&);
471 // exp(cr1) | RealType
472 template <typename RealType, size_t Order>
473 fvar<RealType, Order> exp(fvar<RealType, Order> const&);
475 // pow(cr, ca) | RealType
476 template <typename RealType, size_t Order>
477 fvar<RealType, Order> pow(fvar<RealType, Order> const&, typename fvar<RealType, Order>::root_type const&);
479 // pow(ca, cr) | RealType
480 template <typename RealType, size_t Order>
481 fvar<RealType, Order> pow(typename fvar<RealType, Order>::root_type const&, fvar<RealType, Order> const&);
483 // pow(cr1, cr2) | RealType
484 template <typename RealType1, size_t Order1, typename RealType2, size_t Order2>
485 promote<fvar<RealType1, Order1>, fvar<RealType2, Order2>> pow(fvar<RealType1, Order1> const&,
486 fvar<RealType2, Order2> const&);
488 // sqrt(cr1) | RealType
489 template <typename RealType, size_t Order>
490 fvar<RealType, Order> sqrt(fvar<RealType, Order> const&);
492 // log(cr1) | RealType
493 template <typename RealType, size_t Order>
494 fvar<RealType, Order> log(fvar<RealType, Order> const&);
496 // frexp(cr1, &i) | RealType
497 template <typename RealType, size_t Order>
498 fvar<RealType, Order> frexp(fvar<RealType, Order> const&, int*);
500 // ldexp(cr1, i) | RealType
501 template <typename RealType, size_t Order>
502 fvar<RealType, Order> ldexp(fvar<RealType, Order> const&, int);
504 // cos(cr1) | RealType
505 template <typename RealType, size_t Order>
506 fvar<RealType, Order> cos(fvar<RealType, Order> const&);
508 // sin(cr1) | RealType
509 template <typename RealType, size_t Order>
510 fvar<RealType, Order> sin(fvar<RealType, Order> const&);
512 // asin(cr1) | RealType
513 template <typename RealType, size_t Order>
514 fvar<RealType, Order> asin(fvar<RealType, Order> const&);
516 // tan(cr1) | RealType
517 template <typename RealType, size_t Order>
518 fvar<RealType, Order> tan(fvar<RealType, Order> const&);
520 // atan(cr1) | RealType
521 template <typename RealType, size_t Order>
522 fvar<RealType, Order> atan(fvar<RealType, Order> const&);
524 // atan2(cr, ca) | RealType
525 template <typename RealType, size_t Order>
526 fvar<RealType, Order> atan2(fvar<RealType, Order> const&, typename fvar<RealType, Order>::root_type const&);
528 // atan2(ca, cr) | RealType
529 template <typename RealType, size_t Order>
530 fvar<RealType, Order> atan2(typename fvar<RealType, Order>::root_type const&, fvar<RealType, Order> const&);
532 // atan2(cr1, cr2) | RealType
533 template <typename RealType1, size_t Order1, typename RealType2, size_t Order2>
534 promote<fvar<RealType1, Order1>, fvar<RealType2, Order2>> atan2(fvar<RealType1, Order1> const&,
535 fvar<RealType2, Order2> const&);
537 // fmod(cr1,cr2) | RealType
538 template <typename RealType1, size_t Order1, typename RealType2, size_t Order2>
539 promote<fvar<RealType1, Order1>, fvar<RealType2, Order2>> fmod(fvar<RealType1, Order1> const&,
540 fvar<RealType2, Order2> const&);
542 // round(cr1) | RealType
543 template <typename RealType, size_t Order>
544 fvar<RealType, Order> round(fvar<RealType, Order> const&);
547 template <typename RealType, size_t Order>
548 int iround(fvar<RealType, Order> const&);
550 template <typename RealType, size_t Order>
551 long lround(fvar<RealType, Order> const&);
553 template <typename RealType, size_t Order>
554 long long llround(fvar<RealType, Order> const&);
556 // trunc(cr1) | RealType
557 template <typename RealType, size_t Order>
558 fvar<RealType, Order> trunc(fvar<RealType, Order> const&);
560 template <typename RealType, size_t Order>
561 long double truncl(fvar<RealType, Order> const&);
564 template <typename RealType, size_t Order>
565 int itrunc(fvar<RealType, Order> const&);
567 template <typename RealType, size_t Order>
568 long long lltrunc(fvar<RealType, Order> const&);
570 // Additional functions
571 template <typename RealType, size_t Order>
572 fvar<RealType, Order> acos(fvar<RealType, Order> const&);
574 template <typename RealType, size_t Order>
575 fvar<RealType, Order> acosh(fvar<RealType, Order> const&);
577 template <typename RealType, size_t Order>
578 fvar<RealType, Order> asinh(fvar<RealType, Order> const&);
580 template <typename RealType, size_t Order>
581 fvar<RealType, Order> atanh(fvar<RealType, Order> const&);
583 template <typename RealType, size_t Order>
584 fvar<RealType, Order> cosh(fvar<RealType, Order> const&);
586 template <typename RealType, size_t Order>
587 fvar<RealType, Order> digamma(fvar<RealType, Order> const&);
589 template <typename RealType, size_t Order>
590 fvar<RealType, Order> erf(fvar<RealType, Order> const&);
592 template <typename RealType, size_t Order>
593 fvar<RealType, Order> erfc(fvar<RealType, Order> const&);
595 template <typename RealType, size_t Order>
596 fvar<RealType, Order> lambert_w0(fvar<RealType, Order> const&);
598 template <typename RealType, size_t Order>
599 fvar<RealType, Order> lgamma(fvar<RealType, Order> const&);
601 template <typename RealType, size_t Order>
602 fvar<RealType, Order> sinc(fvar<RealType, Order> const&);
604 template <typename RealType, size_t Order>
605 fvar<RealType, Order> sinh(fvar<RealType, Order> const&);
607 template <typename RealType, size_t Order>
608 fvar<RealType, Order> tanh(fvar<RealType, Order> const&);
610 template <typename RealType, size_t Order>
611 fvar<RealType, Order> tgamma(fvar<RealType, Order> const&);
614 struct zero : std::integral_constant<size_t, 0> {};
616 } // namespace detail
618 template <typename RealType, size_t Order, size_t... Orders>
619 using autodiff_fvar = typename detail::nest_fvar<RealType, Order, Orders...>::type;
621 template <typename RealType, size_t Order, size_t... Orders>
622 autodiff_fvar<RealType, Order, Orders...> make_fvar(RealType const& ca) {
623 return autodiff_fvar<RealType, Order, Orders...>(ca, true);
626 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
629 template <typename RealType, size_t Order, size_t... Is>
630 auto make_fvar_for_tuple(std::index_sequence<Is...>, RealType const& ca) {
631 return make_fvar<RealType, zero<Is>::value..., Order>(ca);
634 template <typename RealType, size_t... Orders, size_t... Is, typename... RealTypes>
635 auto 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)...);
639 } // namespace detail
641 template <typename RealType, size_t... Orders, typename... RealTypes>
642 auto 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...);
651 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
652 template <typename RealType, size_t Order>
653 fvar<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));
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));
668 template <typename RealType, size_t Order>
669 template <typename RealType2, size_t Order2>
670 fvar<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_IF_CONSTEXPR (Order2 < Order)
674 std::fill(v.begin() + (Order2 + 1), v.end(), static_cast<RealType>(0));
677 template <typename RealType, size_t Order>
678 fvar<RealType, Order>::fvar(root_type const& ca) : v{{static_cast<RealType>(ca)}} {}
680 // Can cause compiler error if RealType2 cannot be cast to root_type.
681 template <typename RealType, size_t Order>
682 template <typename RealType2>
683 fvar<RealType, Order>::fvar(RealType2 const& ca) : v{{static_cast<RealType>(ca)}} {}
686 template<typename RealType, size_t Order>
687 fvar<RealType,Order>& fvar<RealType,Order>::operator=(root_type const& ca)
689 v.front() = static_cast<RealType>(ca);
690 if constexpr (0 < Order)
691 std::fill(v.begin()+1, v.end(), static_cast<RealType>(0));
696 template <typename RealType, size_t Order>
697 template <typename RealType2, size_t Order2>
698 fvar<RealType, Order>& fvar<RealType, Order>::operator+=(fvar<RealType2, Order2> const& cr) {
699 for (size_t i = 0; i <= (std::min)(Order, Order2); ++i)
704 template <typename RealType, size_t Order>
705 fvar<RealType, Order>& fvar<RealType, Order>::operator+=(root_type const& ca) {
710 template <typename RealType, size_t Order>
711 template <typename RealType2, size_t Order2>
712 fvar<RealType, Order>& fvar<RealType, Order>::operator-=(fvar<RealType2, Order2> const& cr) {
713 for (size_t i = 0; i <= Order; ++i)
718 template <typename RealType, size_t Order>
719 fvar<RealType, Order>& fvar<RealType, Order>::operator-=(root_type const& ca) {
724 template <typename RealType, size_t Order>
725 template <typename RealType2, size_t Order2>
726 fvar<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_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);
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);
741 template <typename RealType, size_t Order>
742 fvar<RealType, Order>& fvar<RealType, Order>::operator*=(root_type const& ca) {
743 return multiply_assign_by_root_type(true, ca);
746 template <typename RealType, size_t Order>
747 template <typename RealType2, size_t Order2>
748 fvar<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_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_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();
761 for (size_t i = 1; i <= Order; ++i)
762 v[i] /= cr.v.front();
766 template <typename RealType, size_t Order>
767 fvar<RealType, Order>& fvar<RealType, Order>::operator/=(root_type const& ca) {
768 std::for_each(v.begin(), v.end(), [&ca](RealType& x) { x /= ca; });
772 template <typename RealType, size_t Order>
773 fvar<RealType, Order> fvar<RealType, Order>::operator-() const {
774 fvar<RealType, Order> retval(*this);
779 template <typename RealType, size_t Order>
780 fvar<RealType, Order> const& fvar<RealType, Order>::operator+() const {
784 template <typename RealType, size_t Order>
785 template <typename RealType2, size_t Order2>
786 promote<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_IF_CONSTEXPR (Order < Order2)
792 for (size_t i = Order + 1; i <= Order2; ++i)
793 retval.v[i] = cr.v[i];
794 else BOOST_IF_CONSTEXPR (Order2 < Order)
795 for (size_t i = Order2 + 1; i <= Order; ++i)
800 template <typename RealType, size_t Order>
801 fvar<RealType, Order> fvar<RealType, Order>::operator+(root_type const& ca) const {
802 fvar<RealType, Order> retval(*this);
803 retval.v.front() += ca;
807 template <typename RealType, size_t Order>
808 fvar<RealType, Order> operator+(typename fvar<RealType, Order>::root_type const& ca,
809 fvar<RealType, Order> const& cr) {
813 template <typename RealType, size_t Order>
814 template <typename RealType2, size_t Order2>
815 promote<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_IF_CONSTEXPR (Order < Order2)
821 for (auto i = Order + 1; i <= Order2; ++i)
822 retval.v[i] = -cr.v[i];
823 else BOOST_IF_CONSTEXPR (Order2 < Order)
824 for (auto i = Order2 + 1; i <= Order; ++i)
829 template <typename RealType, size_t Order>
830 fvar<RealType, Order> fvar<RealType, Order>::operator-(root_type const& ca) const {
831 fvar<RealType, Order> retval(*this);
832 retval.v.front() -= ca;
836 template <typename RealType, size_t Order>
837 fvar<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.
841 return mcr; // <-- This allows for NRVO. The following does not. --> return mcr += ca;
844 template <typename RealType, size_t Order>
845 template <typename RealType2, size_t Order2>
846 promote<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_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);
855 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);
860 template <typename RealType, size_t Order>
861 fvar<RealType, Order> fvar<RealType, Order>::operator*(root_type const& ca) const {
862 fvar<RealType, Order> retval(*this);
867 template <typename RealType, size_t Order>
868 fvar<RealType, Order> operator*(typename fvar<RealType, Order>::root_type const& ca,
869 fvar<RealType, Order> const& cr) {
873 template <typename RealType, size_t Order>
874 template <typename RealType2, size_t Order2>
875 promote<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_IF_CONSTEXPR (Order < Order2) {
882 for (size_t i = 1, j = Order2 - 1; i <= Order; ++i, --j)
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)) /
887 for (size_t i = Order + 1, j = Order2 - Order - 1; i <= Order2; ++i, --j)
890 cr.v.cbegin() + 1, cr.v.cend() - diff_t(j), retval.v.crbegin() + diff_t(j + 1), zero) /
892 } else BOOST_IF_CONSTEXPR (0 < Order2)
893 for (size_t i = 1, j = Order2 - 1, k = Order; i <= Order; ++i, j && --j, --k)
895 (v[i] - std::inner_product(
896 cr.v.cbegin() + 1, cr.v.cend() - diff_t(j), retval.v.crbegin() + diff_t(k), zero)) /
899 for (size_t i = 1; i <= Order; ++i)
900 retval.v[i] = v[i] / cr.v.front();
904 template <typename RealType, size_t Order>
905 fvar<RealType, Order> fvar<RealType, Order>::operator/(root_type const& ca) const {
906 fvar<RealType, Order> retval(*this);
911 template <typename RealType, size_t Order>
912 fvar<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_IF_CONSTEXPR (0 < Order) {
918 RealType const zero(0);
919 for (size_t i = 1, j = Order - 1; i <= Order; ++i, --j)
922 cr.v.cbegin() + 1, cr.v.cend() - diff_t(j), retval.v.crbegin() + diff_t(j + 1), zero) /
928 template <typename RealType, size_t Order>
929 template <typename RealType2, size_t Order2>
930 bool fvar<RealType, Order>::operator==(fvar<RealType2, Order2> const& cr) const {
931 return v.front() == cr.v.front();
934 template <typename RealType, size_t Order>
935 bool fvar<RealType, Order>::operator==(root_type const& ca) const {
936 return v.front() == ca;
939 template <typename RealType, size_t Order>
940 bool operator==(typename fvar<RealType, Order>::root_type const& ca, fvar<RealType, Order> const& cr) {
941 return ca == cr.v.front();
944 template <typename RealType, size_t Order>
945 template <typename RealType2, size_t Order2>
946 bool fvar<RealType, Order>::operator!=(fvar<RealType2, Order2> const& cr) const {
947 return v.front() != cr.v.front();
950 template <typename RealType, size_t Order>
951 bool fvar<RealType, Order>::operator!=(root_type const& ca) const {
952 return v.front() != ca;
955 template <typename RealType, size_t Order>
956 bool operator!=(typename fvar<RealType, Order>::root_type const& ca, fvar<RealType, Order> const& cr) {
957 return ca != cr.v.front();
960 template <typename RealType, size_t Order>
961 template <typename RealType2, size_t Order2>
962 bool fvar<RealType, Order>::operator<=(fvar<RealType2, Order2> const& cr) const {
963 return v.front() <= cr.v.front();
966 template <typename RealType, size_t Order>
967 bool fvar<RealType, Order>::operator<=(root_type const& ca) const {
968 return v.front() <= ca;
971 template <typename RealType, size_t Order>
972 bool operator<=(typename fvar<RealType, Order>::root_type const& ca, fvar<RealType, Order> const& cr) {
973 return ca <= cr.v.front();
976 template <typename RealType, size_t Order>
977 template <typename RealType2, size_t Order2>
978 bool fvar<RealType, Order>::operator>=(fvar<RealType2, Order2> const& cr) const {
979 return v.front() >= cr.v.front();
982 template <typename RealType, size_t Order>
983 bool fvar<RealType, Order>::operator>=(root_type const& ca) const {
984 return v.front() >= ca;
987 template <typename RealType, size_t Order>
988 bool operator>=(typename fvar<RealType, Order>::root_type const& ca, fvar<RealType, Order> const& cr) {
989 return ca >= cr.v.front();
992 template <typename RealType, size_t Order>
993 template <typename RealType2, size_t Order2>
994 bool fvar<RealType, Order>::operator<(fvar<RealType2, Order2> const& cr) const {
995 return v.front() < cr.v.front();
998 template <typename RealType, size_t Order>
999 bool fvar<RealType, Order>::operator<(root_type const& ca) const {
1000 return v.front() < ca;
1003 template <typename RealType, size_t Order>
1004 bool operator<(typename fvar<RealType, Order>::root_type const& ca, fvar<RealType, Order> const& cr) {
1005 return ca < cr.v.front();
1008 template <typename RealType, size_t Order>
1009 template <typename RealType2, size_t Order2>
1010 bool fvar<RealType, Order>::operator>(fvar<RealType2, Order2> const& cr) const {
1011 return v.front() > cr.v.front();
1014 template <typename RealType, size_t Order>
1015 bool fvar<RealType, Order>::operator>(root_type const& ca) const {
1016 return v.front() > ca;
1019 template <typename RealType, size_t Order>
1020 bool operator>(typename fvar<RealType, Order>::root_type const& ca, fvar<RealType, Order> const& cr) {
1021 return ca > cr.v.front();
1024 /*** Other methods and functions ***/
1026 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1027 // f : order -> derivative(order)/factorial(order)
1028 // Use this when you have the polynomial coefficients, rather than just the derivatives. E.g. See atan2().
1029 template <typename RealType, size_t Order>
1030 template <typename Func, typename Fvar, typename... Fvars>
1031 promote<fvar<RealType, Order>, Fvar, Fvars...> fvar<RealType, Order>::apply_coefficients(
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)...);
1041 (accumulator *= epsilon) += cr.apply_coefficients(
1042 order - i, [&f, i](auto... indices) { return f(i, indices...); }, std::forward<Fvars>(fvars)...);
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().
1049 template <typename RealType, size_t Order>
1050 template <typename Func>
1051 fvar<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_NO_CXX17_IF_CONSTEXPR
1054 size_t i = (std::min)(order, order_sum);
1055 #else // ODR-use of static constexpr
1056 size_t i = order < order_sum ? order : order_sum;
1058 fvar<RealType, Order> accumulator = f(i);
1060 (accumulator *= epsilon) += f(i);
1064 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1065 // f : order -> derivative(order)
1066 template <typename RealType, size_t Order>
1067 template <typename Func, typename Fvar, typename... Fvars>
1068 promote<fvar<RealType, Order>, Fvar, Fvars...> fvar<RealType, Order>::apply_coefficients_nonhorner(
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 i
1075 promote<fvar<RealType, Order>, Fvar, Fvars...> accumulator = cr.apply_coefficients_nonhorner(
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(
1085 cr.apply_coefficients_nonhorner(
1087 [&f, i](auto... indices) { return f(i, static_cast<std::size_t>(indices)...); },
1088 std::forward<Fvars>(fvars)...),
1096 // f : order -> coefficient(order)
1097 template <typename RealType, size_t Order>
1098 template <typename Func>
1099 fvar<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 i
1103 fvar<RealType, Order> accumulator = fvar<RealType, Order>(f(0u));
1104 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1105 size_t const i_max = (std::min)(order, order_sum);
1106 #else // ODR-use of static constexpr
1107 size_t const i_max = order < order_sum ? order : order_sum;
1109 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));
1116 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1117 // f : order -> derivative(order)
1118 template <typename RealType, size_t Order>
1119 template <typename Func, typename Fvar, typename... Fvars>
1120 promote<fvar<RealType, Order>, Fvar, Fvars...> fvar<RealType, Order>::apply_derivatives(
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));
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));
1140 // f : order -> derivative(order)
1141 template <typename RealType, size_t Order>
1142 template <typename Func>
1143 fvar<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_NO_CXX17_IF_CONSTEXPR
1146 size_t i = (std::min)(order, order_sum);
1147 #else // ODR-use of static constexpr
1148 size_t i = order < order_sum ? order : order_sum;
1150 fvar<RealType, Order> accumulator = f(i) / factorial<root_type>(static_cast<unsigned>(i));
1152 (accumulator *= epsilon) += f(i) / factorial<root_type>(static_cast<unsigned>(i));
1156 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1157 // f : order -> derivative(order)
1158 template <typename RealType, size_t Order>
1159 template <typename Func, typename Fvar, typename... Fvars>
1160 promote<fvar<RealType, Order>, Fvar, Fvars...> fvar<RealType, Order>::apply_derivatives_nonhorner(
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 i
1167 promote<fvar<RealType, Order>, Fvar, Fvars...> accumulator = cr.apply_derivatives_nonhorner(
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(
1177 cr.apply_derivatives_nonhorner(
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)),
1189 // f : order -> derivative(order)
1190 template <typename RealType, size_t Order>
1191 template <typename Func>
1192 fvar<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 i
1196 fvar<RealType, Order> accumulator = fvar<RealType, Order>(f(0u));
1197 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1198 size_t const i_max = (std::min)(order, order_sum);
1199 #else // ODR-use of static constexpr
1200 size_t const i_max = order < order_sum ? order : order_sum;
1202 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)));
1209 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1210 // Can throw "std::out_of_range: array::at: __n (which is 7) >= _Nm (which is 7)"
1211 template <typename RealType, size_t Order>
1212 template <typename... Orders>
1213 get_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)...);
1221 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1222 // Can throw "std::out_of_range: array::at: __n (which is 7) >= _Nm (which is 7)"
1223 template <typename RealType, size_t Order>
1224 template <typename... Orders>
1225 get_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)));
1234 template <typename RealType, size_t Order>
1235 const RealType& fvar<RealType, Order>::operator[](size_t i) const {
1239 template <typename RealType, size_t Order>
1240 RealType fvar<RealType, Order>::epsilon_inner_product(size_t z0,
1243 fvar<RealType, Order> const& cr,
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);
1256 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1257 template <typename RealType, size_t Order>
1258 fvar<RealType, Order> fvar<RealType, Order>::epsilon_multiply(size_t z0,
1260 fvar<RealType, Order> const& cr,
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);
1273 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);
1280 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1281 // When called from outside this method, z0 should be non-zero. Otherwise if z0=0 then it will give an
1282 // 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.
1284 template <typename RealType, size_t Order>
1285 fvar<RealType, Order> fvar<RealType, Order>::epsilon_multiply(size_t z0,
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);
1294 for (size_t i = m0; i <= Order; ++i)
1295 if (retval.v[i] != static_cast<RealType>(0))
1301 template <typename RealType, size_t Order>
1302 fvar<RealType, Order> fvar<RealType, Order>::inverse() const {
1303 return static_cast<root_type>(*this) == 0 ? inverse_apply() : 1 / *this;
1306 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1307 template <typename RealType, size_t Order>
1308 fvar<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(); });
1312 std::for_each(v.begin(), v.end(), [](RealType& a) { a = -a; });
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)
1319 template <typename RealType, size_t Order>
1320 fvar<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]; });
1329 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1330 template <typename RealType, size_t Order>
1331 fvar<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);
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)
1349 template <typename RealType, size_t Order>
1350 fvar<RealType, Order>::operator root_type() const {
1351 return static_cast<root_type>(v.front());
1354 template <typename RealType, size_t Order>
1355 template <typename T, typename>
1356 fvar<RealType, Order>::operator T() const {
1357 return static_cast<T>(static_cast<root_type>(v.front()));
1360 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1361 template <typename RealType, size_t Order>
1362 fvar<RealType, Order>& fvar<RealType, Order>::set_root(root_type const& root) {
1363 if constexpr (is_fvar<RealType>::value)
1364 v.front().set_root(root);
1371 // Standard Library Support Requirements
1373 template <typename RealType, size_t Order>
1374 fvar<RealType, Order> fabs(fvar<RealType, Order> const& cr) {
1375 typename fvar<RealType, Order>::root_type const zero(0);
1376 return cr < zero ? -cr
1377 : cr == zero ? fvar<RealType, Order>() // Canonical fabs'(0) = 0.
1378 : cr; // Propagate NaN.
1381 template <typename RealType, size_t Order>
1382 fvar<RealType, Order> abs(fvar<RealType, Order> const& cr) {
1386 template <typename RealType, size_t Order>
1387 fvar<RealType, Order> ceil(fvar<RealType, Order> const& cr) {
1389 return fvar<RealType, Order>(ceil(static_cast<typename fvar<RealType, Order>::root_type>(cr)));
1392 template <typename RealType, size_t Order>
1393 fvar<RealType, Order> floor(fvar<RealType, Order> const& cr) {
1395 return fvar<RealType, Order>(floor(static_cast<typename fvar<RealType, Order>::root_type>(cr)));
1398 template <typename RealType, size_t Order>
1399 fvar<RealType, Order> exp(fvar<RealType, Order> const& cr) {
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; });
1407 template <typename RealType, size_t Order>
1408 fvar<RealType, Order> pow(fvar<RealType, Order> const& x,
1409 typename fvar<RealType, Order>::root_type const& y) {
1410 BOOST_MATH_STD_USING
1411 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()) {
1417 for (size_t i = 0; i < order && y - i != 0; ++i) {
1419 derivatives[i + 1] = coef * pow(x0, y - (i + 1));
1421 return x.apply_derivatives_nonhorner(order, [&derivatives](size_t i) { return derivatives[i]; });
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]; });
1429 template <typename RealType, size_t Order>
1430 fvar<RealType, Order> pow(typename fvar<RealType, Order>::root_type const& x,
1431 fvar<RealType, Order> const& y) {
1432 BOOST_MATH_STD_USING
1433 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]; });
1444 template <typename RealType1, size_t Order1, typename RealType2, size_t Order2>
1445 promote<fvar<RealType1, Order1>, fvar<RealType2, Order2>> pow(fvar<RealType1, Order1> const& x,
1446 fvar<RealType2, Order2> const& y) {
1447 BOOST_MATH_STD_USING
1448 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_IF_CONSTEXPR (order == 0)
1455 return return_type(*dxydx);
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_NO_CXX17_IF_CONSTEXPR
1462 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));
1466 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);
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);
1483 template <typename RealType, size_t Order>
1484 fvar<RealType, Order> sqrt(fvar<RealType, Order> const& cr) {
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_IF_CONSTEXPR (order == 0)
1492 return fvar<RealType, Order>(*derivatives);
1494 root_type numerator = 0.5;
1495 root_type powers = 1;
1496 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
1497 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;
1501 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);
1505 derivatives[i] = numerator / (powers * *derivatives);
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);
1514 // Natural logarithm. If cr==0 then derivative(i) may have nans due to nans from inverse().
1515 template <typename RealType, size_t Order>
1516 fvar<RealType, Order> log(fvar<RealType, Order> const& cr) {
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_IF_CONSTEXPR (order == 0)
1522 return fvar<RealType, Order>(d0);
1524 auto const d1 = make_fvar<root_type, bool(order) ? order - 1 : 0>(static_cast<root_type>(cr)).inverse(); // log'(x) = 1 / x
1525 return cr.apply_coefficients_nonhorner(order, [&d0, &d1](size_t i) { return i ? d1[i - 1] / i : d0; });
1529 template <typename RealType, size_t Order>
1530 fvar<RealType, Order> frexp(fvar<RealType, Order> const& cr, int* exp) {
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));
1538 template <typename RealType, size_t Order>
1539 fvar<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)
1542 return cr * exp2(static_cast<typename fvar<RealType, Order>::root_type>(exp));
1545 template <typename RealType, size_t Order>
1546 fvar<RealType, Order> cos(fvar<RealType, Order> const& cr) {
1547 BOOST_MATH_STD_USING
1548 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_IF_CONSTEXPR (order == 0)
1552 return fvar<RealType, Order>(d0);
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]; });
1560 template <typename RealType, size_t Order>
1561 fvar<RealType, Order> sin(fvar<RealType, Order> const& cr) {
1562 BOOST_MATH_STD_USING
1563 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_IF_CONSTEXPR (order == 0)
1567 return fvar<RealType, Order>(d0);
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]; });
1575 template <typename RealType, size_t Order>
1576 fvar<RealType, Order> asin(fvar<RealType, Order> const& cr) {
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_IF_CONSTEXPR (order == 0)
1582 return fvar<RealType, Order>(d0);
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; });
1590 template <typename RealType, size_t Order>
1591 fvar<RealType, Order> tan(fvar<RealType, Order> const& cr) {
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_IF_CONSTEXPR (order == 0)
1597 return fvar<RealType, Order>(d0);
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)^2
1601 return cr.apply_coefficients_nonhorner(order, [&d0, &d1](size_t i) { return i ? d1[i - 1] / i : d0; });
1605 template <typename RealType, size_t Order>
1606 fvar<RealType, Order> atan(fvar<RealType, Order> const& cr) {
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_IF_CONSTEXPR (order == 0)
1612 return fvar<RealType, Order>(d0);
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; });
1620 template <typename RealType, size_t Order>
1621 fvar<RealType, Order> atan2(fvar<RealType, Order> const& cr,
1622 typename fvar<RealType, Order>::root_type const& ca) {
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_IF_CONSTEXPR (order == 0)
1628 return fvar<RealType, Order>(d0);
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; });
1636 template <typename RealType, size_t Order>
1637 fvar<RealType, Order> atan2(typename fvar<RealType, Order>::root_type const& ca,
1638 fvar<RealType, Order> const& cr) {
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_IF_CONSTEXPR (order == 0)
1644 return fvar<RealType, Order>(d0);
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; });
1652 template <typename RealType1, size_t Order1, typename RealType2, size_t Order2>
1653 promote<fvar<RealType1, Order1>, fvar<RealType2, Order2>> atan2(fvar<RealType1, Order1> const& cr1,
1654 fvar<RealType2, Order2> const& cr2) {
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_IF_CONSTEXPR (order == 0)
1663 return return_type(d00);
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;
1675 return cr1.apply_coefficients(order, f, cr2);
1679 template <typename RealType1, size_t Order1, typename RealType2, size_t Order2>
1680 promote<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);
1688 template <typename RealType, size_t Order>
1689 fvar<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)));
1694 template <typename RealType, size_t Order>
1695 int iround(fvar<RealType, Order> const& cr) {
1696 using boost::math::iround;
1697 return iround(static_cast<typename fvar<RealType, Order>::root_type>(cr));
1700 template <typename RealType, size_t Order>
1701 long lround(fvar<RealType, Order> const& cr) {
1702 using boost::math::lround;
1703 return lround(static_cast<typename fvar<RealType, Order>::root_type>(cr));
1706 template <typename RealType, size_t Order>
1707 long long llround(fvar<RealType, Order> const& cr) {
1708 using boost::math::llround;
1709 return llround(static_cast<typename fvar<RealType, Order>::root_type>(cr));
1712 template <typename RealType, size_t Order>
1713 fvar<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)));
1718 template <typename RealType, size_t Order>
1719 long double truncl(fvar<RealType, Order> const& cr) {
1721 return truncl(static_cast<typename fvar<RealType, Order>::root_type>(cr));
1724 template <typename RealType, size_t Order>
1725 int itrunc(fvar<RealType, Order> const& cr) {
1726 using boost::math::itrunc;
1727 return itrunc(static_cast<typename fvar<RealType, Order>::root_type>(cr));
1730 template <typename RealType, size_t Order>
1731 long long lltrunc(fvar<RealType, Order> const& cr) {
1732 using boost::math::lltrunc;
1733 return lltrunc(static_cast<typename fvar<RealType, Order>::root_type>(cr));
1736 template <typename RealType, size_t Order>
1737 std::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];
1744 // Additional functions
1746 template <typename RealType, size_t Order>
1747 fvar<RealType, Order> acos(fvar<RealType, Order> const& cr) {
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_IF_CONSTEXPR (order == 0)
1753 return fvar<RealType, Order>(d0);
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; });
1761 template <typename RealType, size_t Order>
1762 fvar<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_IF_CONSTEXPR (order == 0)
1768 return fvar<RealType, Order>(d0);
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; });
1776 template <typename RealType, size_t Order>
1777 fvar<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_IF_CONSTEXPR (order == 0)
1783 return fvar<RealType, Order>(d0);
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; });
1791 template <typename RealType, size_t Order>
1792 fvar<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_IF_CONSTEXPR (order == 0)
1798 return fvar<RealType, Order>(d0);
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; });
1806 template <typename RealType, size_t Order>
1807 fvar<RealType, Order> cosh(fvar<RealType, Order> const& cr) {
1808 BOOST_MATH_STD_USING
1809 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_IF_CONSTEXPR (order == 0)
1813 return fvar<RealType, Order>(d0);
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]; });
1820 template <typename RealType, size_t Order>
1821 fvar<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_IF_CONSTEXPR (order == 0)
1828 return fvar<RealType, Order>(d0);
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; });
1837 template <typename RealType, size_t Order>
1838 fvar<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_IF_CONSTEXPR (order == 0)
1844 return fvar<RealType, Order>(d0);
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; });
1852 template <typename RealType, size_t Order>
1853 fvar<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_IF_CONSTEXPR (order == 0)
1859 return fvar<RealType, Order>(d0);
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; });
1867 template <typename RealType, size_t Order>
1868 fvar<RealType, Order> lambert_w0(fvar<RealType, Order> const& cr) {
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_IF_CONSTEXPR (order == 0)
1876 return fvar<RealType, Order>(*derivatives);
1878 root_type const expw = exp(*derivatives);
1879 derivatives[1] = 1 / (static_cast<root_type>(cr) + expw);
1880 BOOST_IF_CONSTEXPR (order == 1)
1881 return cr.apply_derivatives_nonhorner(order, [&derivatives](size_t i) { return derivatives[i]; });
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];
1895 d1powers * std::accumulate(coef.crend() - diff_t(n - 1),
1898 [&x](root_type const& a, root_type const& b) { return a * x + b; });
1900 return cr.apply_derivatives_nonhorner(order, [&derivatives](size_t i) { return derivatives[i]; });
1905 template <typename RealType, size_t Order>
1906 fvar<RealType, Order> lgamma(fvar<RealType, Order> const& cr) {
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_IF_CONSTEXPR (order == 0)
1913 return fvar<RealType, Order>(d0);
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; });
1922 template <typename RealType, size_t Order>
1923 fvar<RealType, Order> sinc(fvar<RealType, Order> const& cr) {
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) = 1
1929 BOOST_IF_CONSTEXPR (order == 0)
1930 return fvar<RealType, Order>(*taylor);
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]; });
1938 template <typename RealType, size_t Order>
1939 fvar<RealType, Order> sinh(fvar<RealType, Order> const& cr) {
1940 BOOST_MATH_STD_USING
1941 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_IF_CONSTEXPR (fvar<RealType, Order>::order_sum == 0)
1945 return fvar<RealType, Order>(d0);
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]; });
1952 template <typename RealType, size_t Order>
1953 fvar<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;
1960 template <typename RealType, size_t Order>
1961 fvar<RealType, Order> tgamma(fvar<RealType, Order> const& cr) {
1963 using root_type = typename fvar<RealType, Order>::root_type;
1964 constexpr size_t order = fvar<RealType, Order>::order_sum;
1965 BOOST_IF_CONSTEXPR (order == 0)
1966 return fvar<RealType, Order>(tgamma(static_cast<root_type>(cr)));
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)));
1974 } // namespace detail
1975 } // namespace autodiff_v1
1976 } // namespace differentiation
1978 } // namespace boost
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.
1984 template <typename RealType, size_t Order>
1985 class numeric_limits<boost::math::differentiation::detail::fvar<RealType, Order>>
1986 : public numeric_limits<typename boost::math::differentiation::detail::fvar<RealType, Order>::root_type> {
1996 template <typename RealType, std::size_t Order>
1997 using autodiff_fvar_type = differentiation::detail::fvar<RealType, Order>;
1999 template <typename RealType, std::size_t Order>
2000 using autodiff_root_type = typename autodiff_fvar_type<RealType, Order>::root_type;
2001 } // namespace detail
2003 // See boost/math/tools/promotion.hpp
2004 template <typename RealType0, size_t Order0, typename RealType1, size_t Order1>
2005 struct promote_args_2<detail::autodiff_fvar_type<RealType0, Order0>,
2006 detail::autodiff_fvar_type<RealType1, Order1>> {
2007 using type = detail::autodiff_fvar_type<typename promote_args_2<RealType0, RealType1>::type,
2008 #ifndef BOOST_NO_CXX14_CONSTEXPR
2009 (std::max)(Order0, Order1)>;
2011 Order0<Order1 ? Order1 : Order0>;
2015 template <typename RealType, size_t Order>
2016 struct promote_args<detail::autodiff_fvar_type<RealType, Order>> {
2017 using type = detail::autodiff_fvar_type<typename promote_args<RealType>::type, Order>;
2020 template <typename RealType0, size_t Order0, typename RealType1>
2021 struct promote_args_2<detail::autodiff_fvar_type<RealType0, Order0>, RealType1> {
2022 using type = detail::autodiff_fvar_type<typename promote_args_2<RealType0, RealType1>::type, Order0>;
2025 template <typename RealType0, typename RealType1, size_t Order1>
2026 struct promote_args_2<RealType0, detail::autodiff_fvar_type<RealType1, Order1>> {
2027 using type = detail::autodiff_fvar_type<typename promote_args_2<RealType0, RealType1>::type, Order1>;
2030 template <typename destination_t, typename RealType, std::size_t Order>
2031 inline 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));
2036 } // namespace tools
2038 namespace policies {
2040 template <class Policy, std::size_t Order>
2041 using fvar_t = differentiation::detail::fvar<Policy, Order>;
2042 template <class Policy, std::size_t Order>
2043 struct evaluation<fvar_t<float, Order>, Policy> {
2044 using type = fvar_t<typename std::conditional<Policy::promote_float_type::value, double, float>::type, Order>;
2047 template <class Policy, std::size_t Order>
2048 struct evaluation<fvar_t<double, Order>, Policy> {
2050 fvar_t<typename std::conditional<Policy::promote_double_type::value, long double, double>::type, Order>;
2053 } // namespace policies
2055 } // namespace boost
2057 #ifdef BOOST_NO_CXX17_IF_CONSTEXPR
2058 #include "autodiff_cpp11.hpp"
2061 #endif // BOOST_MATH_DIFFERENTIATION_AUTODIFF_HPP