]>
Commit | Line | Data |
---|---|---|
92f5a8d4 TL |
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> | |
f67539c2 TL |
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> | |
92f5a8d4 TL |
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> | |
f67539c2 | 67 | using is_fvar = is_fvar_impl<typename std::decay<T>::type>; |
92f5a8d4 TL |
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> | |
f67539c2 | 87 | using get_depth = get_depth_impl<typename std::decay<T>::type>; |
92f5a8d4 TL |
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> | |
f67539c2 | 97 | using get_order_sum = get_order_sum_t<typename std::decay<T>::type>; |
92f5a8d4 TL |
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 | ||
f67539c2 | 327 | template <typename T, typename = typename std::enable_if<std::is_arithmetic<typename std::decay<T>::type>::value>> |
92f5a8d4 TL |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
779 | 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> | |
786 | 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> | |
792 | 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> | |
807 | 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> | |
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> | |
821 | 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> | |
836 | 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> | |
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> | |
852 | 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> | |
867 | 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> | |
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> | |
881 | 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> | |
911 | 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> | |
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> | |
936 | 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> | |
941 | bool fvar<RealType, Order>::operator==(root_type const& ca) const { | |
942 | return v.front() == ca; | |
943 | } | |
944 | ||
945 | template <typename RealType, size_t Order> | |
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> | |
952 | 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> | |
957 | bool fvar<RealType, Order>::operator!=(root_type const& ca) const { | |
958 | return v.front() != ca; | |
959 | } | |
960 | ||
961 | template <typename RealType, size_t Order> | |
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> | |
968 | 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> | |
973 | bool fvar<RealType, Order>::operator<=(root_type const& ca) const { | |
974 | return v.front() <= ca; | |
975 | } | |
976 | ||
977 | template <typename RealType, size_t Order> | |
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> | |
984 | 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> | |
989 | bool fvar<RealType, Order>::operator>=(root_type const& ca) const { | |
990 | return v.front() >= ca; | |
991 | } | |
992 | ||
993 | template <typename RealType, size_t Order> | |
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> | |
1000 | 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> | |
1005 | bool fvar<RealType, Order>::operator<(root_type const& ca) const { | |
1006 | return v.front() < ca; | |
1007 | } | |
1008 | ||
1009 | template <typename RealType, size_t Order> | |
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> | |
1016 | 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> | |
1021 | bool fvar<RealType, Order>::operator>(root_type const& ca) const { | |
1022 | return v.front() > ca; | |
1023 | } | |
1024 | ||
1025 | template <typename RealType, size_t Order> | |
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> | |
1037 | 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> | |
1057 | 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> | |
1074 | 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> | |
1105 | 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> | |
1126 | 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> | |
1149 | 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> | |
1166 | 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> | |
1198 | 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> | |
1219 | 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> | |
1231 | 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> | |
1241 | const RealType& fvar<RealType, Order>::operator[](size_t i) const { | |
1242 | return v[i]; | |
1243 | } | |
1244 | ||
1245 | template <typename RealType, size_t Order> | |
1246 | 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> | |
1264 | 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> | |
1291 | 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> | |
1308 | 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> | |
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> | |
1326 | 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> | |
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> | |
1356 | 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> | |
1362 | 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> | |
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> | |
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> | |
1388 | fvar<RealType, Order> abs(fvar<RealType, Order> const& cr) { | |
1389 | return fabs(cr); | |
1390 | } | |
1391 | ||
1392 | template <typename RealType, size_t Order> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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 |