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