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