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 #include "test_autodiff.hpp"
7
8 BOOST_AUTO_TEST_SUITE(test_autodiff_2)
9
BOOST_AUTO_TEST_CASE_TEMPLATE(one_over_one_plus_x_squared,T,all_float_types)10 BOOST_AUTO_TEST_CASE_TEMPLATE(one_over_one_plus_x_squared, T, all_float_types) {
11 constexpr std::size_t m = 4;
12 const T cx(1);
13 auto f = make_fvar<T, m>(cx);
14 // f = 1 / ((f *= f) += 1);
15 f *= f;
16 f += T(1);
17 f = f.inverse();
18 BOOST_CHECK_EQUAL(f.derivative(0u), 0.5);
19 BOOST_CHECK_EQUAL(f.derivative(1u), -0.5);
20 BOOST_CHECK_EQUAL(f.derivative(2u), 0.5);
21 BOOST_CHECK_EQUAL(f.derivative(3u), 0);
22 BOOST_CHECK_EQUAL(f.derivative(4u), -3);
23 }
24
BOOST_AUTO_TEST_CASE_TEMPLATE(exp_test,T,all_float_types)25 BOOST_AUTO_TEST_CASE_TEMPLATE(exp_test, T, all_float_types) {
26 using std::exp;
27 constexpr std::size_t m = 4;
28 const T cx = 2.0;
29 const auto x = make_fvar<T, m>(cx);
30 auto y = exp(x);
31 for (auto i : boost::irange(m + 1)) {
32 // std::cout.precision(100);
33 // std::cout << "y.derivative("<<i<<") = " << y.derivative(i) << ",
34 // std::exp(cx) = " << std::exp(cx) << std::endl;
35 BOOST_CHECK_CLOSE_FRACTION(y.derivative(i), exp(cx),
36 std::numeric_limits<T>::epsilon());
37 }
38 }
39
BOOST_AUTO_TEST_CASE_TEMPLATE(pow,T,bin_float_types)40 BOOST_AUTO_TEST_CASE_TEMPLATE(pow, T, bin_float_types) {
41 const T eps = 201 * std::numeric_limits<T>::epsilon(); // percent
42 using std::log;
43 using std::pow;
44 constexpr std::size_t m = 5;
45 constexpr std::size_t n = 4;
46 const T cx = 2.0;
47 const T cy = 3.0;
48 const auto x = make_fvar<T, m>(cx);
49 const auto y = make_fvar<T, m, n>(cy);
50 auto z0 = pow(x, cy);
51 BOOST_CHECK_EQUAL(z0.derivative(0u), pow(cx, cy));
52 BOOST_CHECK_EQUAL(z0.derivative(1u), cy * pow(cx, cy - 1));
53 BOOST_CHECK_EQUAL(z0.derivative(2u), cy * (cy - 1) * pow(cx, cy - 2));
54 BOOST_CHECK_EQUAL(z0.derivative(3u),
55 cy * (cy - 1) * (cy - 2) * pow(cx, cy - 3));
56 BOOST_CHECK_EQUAL(z0.derivative(4u), 0u);
57 BOOST_CHECK_EQUAL(z0.derivative(5u), 0u);
58 auto z1 = pow(cx, y);
59 BOOST_CHECK_CLOSE(z1.derivative(0u, 0u), pow(cx, cy), eps);
60 for (auto j : boost::irange(std::size_t(1), n + 1)) {
61 BOOST_CHECK_CLOSE(z1.derivative(0u, j), pow(log(cx), j) * pow(cx, cy), eps);
62 }
63
64 for (auto i : boost::irange(std::size_t(1), m + 1)) {
65 for (auto j : boost::irange(n + 1)) {
66 BOOST_CHECK_EQUAL(z1.derivative(i, j), 0);
67 }
68 }
69
70 const auto z2 = pow(x, y);
71
72 for (auto j : boost::irange(n + 1)) {
73 BOOST_CHECK_CLOSE(z2.derivative(0u, j), pow(cx, cy) * pow(log(cx), j), eps);
74 }
75 for (auto j : boost::irange(n + 1)) {
76 BOOST_CHECK_CLOSE(z2.derivative(1u, j),
77 pow(cx, cy - 1) * pow(log(cx), static_cast<int>(j) - 1) *
78 (cy * log(cx) + j),
79 eps);
80 }
81 BOOST_CHECK_CLOSE(z2.derivative(2u, 0u), pow(cx, cy - 2) * cy * (cy - 1),
82 eps);
83 BOOST_CHECK_CLOSE(z2.derivative(2u, 1u),
84 pow(cx, cy - 2) * (cy * (cy - 1) * log(cx) + 2 * cy - 1),
85 eps);
86 for (auto j : boost::irange(std::size_t(2u), n + 1)) {
87 BOOST_CHECK_CLOSE(z2.derivative(2u, j),
88 pow(cx, cy - 2) * pow(log(cx), j - 2) *
89 (j * (2 * cy - 1) * log(cx) + (j - 1) * j +
90 (cy - 1) * cy * pow(log(cx), 2)),
91 eps);
92 }
93 BOOST_CHECK_CLOSE(z2.derivative(2u, 4u),
94 pow(cx, cy - 2) * pow(log(cx), 2) *
95 (4 * (2 * cy - 1) * log(cx) + (4 - 1) * 4 +
96 (cy - 1) * cy * pow(log(cx), 2)),
97 eps);
98 }
99
100 // TODO Tests around x=0 or y=0: pow(x,y)
BOOST_AUTO_TEST_CASE_TEMPLATE(pow2,T,bin_float_types)101 BOOST_AUTO_TEST_CASE_TEMPLATE(pow2, T, bin_float_types) {
102 const T eps = 4000 * std::numeric_limits<T>::epsilon(); // percent
103 using std::pow;
104 constexpr std::size_t m = 5;
105 constexpr std::size_t n = 5;
106 const T cx = 2;
107 const T cy = 5 / 2.0;
108 const auto x = make_fvar<T, m>(cx);
109 const auto y = make_fvar<T, 0, n>(cy);
110 const auto z = pow(x, y);
111 using namespace boost::math::constants;
112 // Mathematica: Export["pow.csv", Flatten@Table[ Simplify@D[x^y,{x,i},{y,j}]
113 // /. {x->2, y->5/2},
114 // { i, 0, 5 }, { j, 0, 5 } ] ]
115 // sed -rf pow.sed < pow.csv
116 // with pow.sed script:
117 // s/Log\[2\]\^([0-9]+)/pow(ln_two<T>(),\1)/g
118 // s/Log\[2\]/ln_two<T>()/g
119 // s/Sqrt\[2\]/root_two<T>()/g
120 // s/[0-9]\/[0-9]+/\0.0/g
121 // s/^"/static_cast<T>(/
122 // s/"$/),/
123 const T mathematica[]{
124 static_cast<T>(4 * root_two<T>()),
125 static_cast<T>(4 * root_two<T>() * ln_two<T>()),
126 static_cast<T>(4 * root_two<T>() * pow(ln_two<T>(), 2)),
127 static_cast<T>(4 * root_two<T>() * pow(ln_two<T>(), 3)),
128 static_cast<T>(4 * root_two<T>() * pow(ln_two<T>(), 4)),
129 static_cast<T>(4 * root_two<T>() * pow(ln_two<T>(), 5)),
130 static_cast<T>(5 * root_two<T>()),
131 static_cast<T>(2 * root_two<T>() * (1 + (5 * ln_two<T>()) / 2)),
132 static_cast<T>(2 * root_two<T>() * ln_two<T>() *
133 (2 + (5 * ln_two<T>()) / 2)),
134 static_cast<T>(2 * root_two<T>() * pow(ln_two<T>(), 2) *
135 (3 + (5 * ln_two<T>()) / 2)),
136 static_cast<T>(2 * root_two<T>() * pow(ln_two<T>(), 3) *
137 (4 + (5 * ln_two<T>()) / 2)),
138 static_cast<T>(2 * root_two<T>() * pow(ln_two<T>(), 4) *
139 (5 + (5 * ln_two<T>()) / 2)),
140 static_cast<T>(15 / (2 * root_two<T>())),
141 static_cast<T>(root_two<T>() * (4 + (15 * ln_two<T>()) / 4)),
142 static_cast<T>(root_two<T>() *
143 (2 + 8 * ln_two<T>() + (15 * pow(ln_two<T>(), 2)) / 4)),
144 static_cast<T>(root_two<T>() * ln_two<T>() *
145 (6 + 12 * ln_two<T>() + (15 * pow(ln_two<T>(), 2)) / 4)),
146 static_cast<T>(root_two<T>() * pow(ln_two<T>(), 2) *
147 (12 + 16 * ln_two<T>() + (15 * pow(ln_two<T>(), 2)) / 4)),
148 static_cast<T>(root_two<T>() * pow(ln_two<T>(), 3) *
149 (20 + 20 * ln_two<T>() + (15 * pow(ln_two<T>(), 2)) / 4)),
150 static_cast<T>(15 / (8 * root_two<T>())),
151 static_cast<T>((23 / 4.0 + (15 * ln_two<T>()) / 8) / root_two<T>()),
152 static_cast<T>(
153 (9 + (23 * ln_two<T>()) / 2 + (15 * pow(ln_two<T>(), 2)) / 8) /
154 root_two<T>()),
155 static_cast<T>((6 + 27 * ln_two<T>() + (69 * pow(ln_two<T>(), 2)) / 4 +
156 (15 * pow(ln_two<T>(), 3)) / 8) /
157 root_two<T>()),
158 static_cast<T>(
159 (ln_two<T>() * (24 + 54 * ln_two<T>() + 23 * pow(ln_two<T>(), 2) +
160 (15 * pow(ln_two<T>(), 3)) / 8)) /
161 root_two<T>()),
162 static_cast<T>((pow(ln_two<T>(), 2) *
163 (60 + 90 * ln_two<T>() + (115 * pow(ln_two<T>(), 2)) / 4 +
164 (15 * pow(ln_two<T>(), 3)) / 8)) /
165 root_two<T>()),
166 static_cast<T>(-15 / (32 * root_two<T>())),
167 static_cast<T>((-1 - (15 * ln_two<T>()) / 16) / (2 * root_two<T>())),
168 static_cast<T>((7 - 2 * ln_two<T>() - (15 * pow(ln_two<T>(), 2)) / 16) /
169 (2 * root_two<T>())),
170 static_cast<T>((24 + 21 * ln_two<T>() - 3 * pow(ln_two<T>(), 2) -
171 (15 * pow(ln_two<T>(), 3)) / 16) /
172 (2 * root_two<T>())),
173 static_cast<T>((24 + 96 * ln_two<T>() + 42 * pow(ln_two<T>(), 2) -
174 4 * pow(ln_two<T>(), 3) -
175 (15 * pow(ln_two<T>(), 4)) / 16) /
176 (2 * root_two<T>())),
177 static_cast<T>(
178 (ln_two<T>() *
179 (120 + 240 * ln_two<T>() + 70 * pow(ln_two<T>(), 2) -
180 5 * pow(ln_two<T>(), 3) - (15 * pow(ln_two<T>(), 4)) / 16)) /
181 (2 * root_two<T>())),
182 static_cast<T>(45 / (128 * root_two<T>())),
183 static_cast<T>((9 / 16.0 + (45 * ln_two<T>()) / 32) /
184 (4 * root_two<T>())),
185 static_cast<T>((-25 / 2.0 + (9 * ln_two<T>()) / 8 +
186 (45 * pow(ln_two<T>(), 2)) / 32) /
187 (4 * root_two<T>())),
188 static_cast<T>((-15 - (75 * ln_two<T>()) / 2 +
189 (27 * pow(ln_two<T>(), 2)) / 16 +
190 (45 * pow(ln_two<T>(), 3)) / 32) /
191 (4 * root_two<T>())),
192 static_cast<T>((60 - 60 * ln_two<T>() - 75 * pow(ln_two<T>(), 2) +
193 (9 * pow(ln_two<T>(), 3)) / 4 +
194 (45 * pow(ln_two<T>(), 4)) / 32) /
195 (4 * root_two<T>())),
196 static_cast<T>((120 + 300 * ln_two<T>() - 150 * pow(ln_two<T>(), 2) -
197 125 * pow(ln_two<T>(), 3) +
198 (45 * pow(ln_two<T>(), 4)) / 16 +
199 (45 * pow(ln_two<T>(), 5)) / 32) /
200 (4 * root_two<T>()))};
201 std::size_t k = 0;
202 for (auto i : boost::irange(m + 1)) {
203 for (auto j : boost::irange(n + 1)) {
204 BOOST_CHECK_CLOSE(z.derivative(i, j), mathematica[k++], eps);
205 }
206 }
207 }
208
BOOST_AUTO_TEST_CASE_TEMPLATE(sqrt_test,T,all_float_types)209 BOOST_AUTO_TEST_CASE_TEMPLATE(sqrt_test, T, all_float_types) {
210 using std::pow;
211 using std::sqrt;
212 constexpr std::size_t m = 5;
213 const T cx = 4.0;
214 auto x = make_fvar<T, m>(cx);
215 auto y = sqrt(x);
216 BOOST_CHECK_CLOSE_FRACTION(y.derivative(0u), sqrt(cx),
217 std::numeric_limits<T>::epsilon());
218 BOOST_CHECK_CLOSE_FRACTION(y.derivative(1u), 0.5 * pow(cx, -0.5),
219 std::numeric_limits<T>::epsilon());
220 BOOST_CHECK_CLOSE_FRACTION(y.derivative(2u), -0.5 * 0.5 * pow(cx, -1.5),
221 std::numeric_limits<T>::epsilon());
222 BOOST_CHECK_CLOSE_FRACTION(y.derivative(3u), 0.5 * 0.5 * 1.5 * pow(cx, -2.5),
223 std::numeric_limits<T>::epsilon());
224 BOOST_CHECK_CLOSE_FRACTION(y.derivative(4u),
225 -0.5 * 0.5 * 1.5 * 2.5 * pow(cx, -3.5),
226 std::numeric_limits<T>::epsilon());
227 BOOST_CHECK_CLOSE_FRACTION(y.derivative(5u),
228 0.5 * 0.5 * 1.5 * 2.5 * 3.5 * pow(cx, -4.5),
229 std::numeric_limits<T>::epsilon());
230 x = make_fvar<T, m>(0);
231 y = sqrt(x);
232 // std::cout << "sqrt(0) = " << y << std::endl; // (0,inf,-inf,inf,-inf,inf)
233 BOOST_CHECK_EQUAL(y.derivative(0u), 0);
234 for (auto i : boost::irange(std::size_t(1), m + 1)) {
235 BOOST_CHECK_EQUAL(y.derivative(i), (i % 2 == 1 ? 1 : -1) *
236 std::numeric_limits<T>::infinity());
237 }
238 }
239
BOOST_AUTO_TEST_CASE_TEMPLATE(log_test,T,all_float_types)240 BOOST_AUTO_TEST_CASE_TEMPLATE(log_test, T, all_float_types) {
241 using std::log;
242 using std::pow;
243 constexpr std::size_t m = 5;
244 const T cx = 2.0;
245 auto x = make_fvar<T, m>(cx);
246 auto y = log(x);
247 BOOST_CHECK_CLOSE_FRACTION(y.derivative(0u), log(cx),
248 std::numeric_limits<T>::epsilon());
249 BOOST_CHECK_CLOSE_FRACTION(y.derivative(1u), 1 / cx,
250 std::numeric_limits<T>::epsilon());
251 BOOST_CHECK_CLOSE_FRACTION(y.derivative(2u), -1 / pow(cx, 2),
252 std::numeric_limits<T>::epsilon());
253 BOOST_CHECK_CLOSE_FRACTION(y.derivative(3u), 2 / pow(cx, 3),
254 std::numeric_limits<T>::epsilon());
255 BOOST_CHECK_CLOSE_FRACTION(y.derivative(4u), -6 / pow(cx, 4),
256 std::numeric_limits<T>::epsilon());
257 BOOST_CHECK_CLOSE_FRACTION(y.derivative(5u), 24 / pow(cx, 5),
258 std::numeric_limits<T>::epsilon());
259 x = make_fvar<T, m>(0);
260 y = log(x);
261 // std::cout << "log(0) = " << y << std::endl; // log(0) =
262 // depth(1)(-inf,inf,-inf,inf,-inf,inf)
263 for (auto i : boost::irange(m + 1)) {
264 BOOST_CHECK_EQUAL(y.derivative(i), (i % 2 == 1 ? 1 : -1) *
265 std::numeric_limits<T>::infinity());
266 }
267 }
268
BOOST_AUTO_TEST_CASE_TEMPLATE(ylogx,T,all_float_types)269 BOOST_AUTO_TEST_CASE_TEMPLATE(ylogx, T, all_float_types) {
270 using std::log;
271 using std::pow;
272 const T eps = 100 * std::numeric_limits<T>::epsilon(); // percent
273 constexpr std::size_t m = 5;
274 constexpr std::size_t n = 4;
275 const T cx = 2.0;
276 const T cy = 3.0;
277 const auto x = make_fvar<T, m>(cx);
278 const auto y = make_fvar<T, m, n>(cy);
279 auto z = y * log(x);
280 BOOST_CHECK_EQUAL(z.derivative(0u, 0u), cy * log(cx));
281 BOOST_CHECK_EQUAL(z.derivative(0u, 1u), log(cx));
282 BOOST_CHECK_EQUAL(z.derivative(0u, 2u), 0);
283 BOOST_CHECK_EQUAL(z.derivative(0u, 3u), 0);
284 BOOST_CHECK_EQUAL(z.derivative(0u, 4u), 0);
285 for (auto i : boost::irange(1u, static_cast<unsigned>(m + 1))) {
286 BOOST_CHECK_CLOSE(z.derivative(i, 0u),
287 pow(-1, i - 1) * boost::math::factorial<T>(i - 1) * cy /
288 pow(cx, i),
289 eps);
290 BOOST_CHECK_CLOSE(
291 z.derivative(i, 1u),
292 pow(-1, i - 1) * boost::math::factorial<T>(i - 1) / pow(cx, i), eps);
293 for (auto j : boost::irange(std::size_t(2), n + 1)) {
294 BOOST_CHECK_EQUAL(z.derivative(i, j), 0u);
295 }
296 }
297 auto z1 = exp(z);
298 // RHS is confirmed by
299 // https://www.wolframalpha.com/input/?i=D%5Bx%5Ey,%7Bx,2%7D,%7By,4%7D%5D+%2F.+%7Bx-%3E2.0,+y-%3E3.0%7D
300 BOOST_CHECK_CLOSE(z1.derivative(2u, 4u),
301 pow(cx, cy - 2) * pow(log(cx), 2) *
302 (4 * (2 * cy - 1) * log(cx) + (4 - 1) * 4 +
303 (cy - 1) * cy * pow(log(cx), 2)),
304 eps);
305 }
306
BOOST_AUTO_TEST_CASE_TEMPLATE(frexp_test,T,all_float_types)307 BOOST_AUTO_TEST_CASE_TEMPLATE(frexp_test, T, all_float_types) {
308 using std::exp2;
309 using std::frexp;
310 constexpr std::size_t m = 3;
311 const T cx = 3.5;
312 const auto x = make_fvar<T, m>(cx);
313 int exp, testexp;
314 auto y = frexp(x, &exp);
315 BOOST_CHECK_EQUAL(y.derivative(0u), frexp(cx, &testexp));
316 BOOST_CHECK_EQUAL(exp, testexp);
317 BOOST_CHECK_EQUAL(y.derivative(1u), exp2(-exp));
318 BOOST_CHECK_EQUAL(y.derivative(2u), 0);
319 BOOST_CHECK_EQUAL(y.derivative(3u), 0);
320 }
321
BOOST_AUTO_TEST_CASE_TEMPLATE(ldexp_test,T,all_float_types)322 BOOST_AUTO_TEST_CASE_TEMPLATE(ldexp_test, T, all_float_types) {
323 BOOST_MATH_STD_USING
324 using boost::multiprecision::ldexp;
325 constexpr auto m = 3u;
326 const T cx = 3.5;
327 const auto x = make_fvar<T, m>(cx);
328 constexpr auto exponent = 3;
329 auto y = ldexp(x, exponent);
330 BOOST_CHECK_EQUAL(y.derivative(0u), ldexp(cx, exponent));
331 BOOST_CHECK_EQUAL(y.derivative(1u), exp2(exponent));
332 BOOST_CHECK_EQUAL(y.derivative(2u), 0);
333 BOOST_CHECK_EQUAL(y.derivative(3u), 0);
334 }
335
BOOST_AUTO_TEST_CASE_TEMPLATE(cos_and_sin,T,bin_float_types)336 BOOST_AUTO_TEST_CASE_TEMPLATE(cos_and_sin, T, bin_float_types) {
337 using std::cos;
338 using std::sin;
339 const T eps = 200 * std::numeric_limits<T>::epsilon(); // percent
340 constexpr std::size_t m = 5;
341 const T cx = boost::math::constants::third_pi<T>();
342 const auto x = make_fvar<T, m>(cx);
343 auto cos5 = cos(x);
344 BOOST_CHECK_CLOSE(cos5.derivative(0u), cos(cx), eps);
345 BOOST_CHECK_CLOSE(cos5.derivative(1u), -sin(cx), eps);
346 BOOST_CHECK_CLOSE(cos5.derivative(2u), -cos(cx), eps);
347 BOOST_CHECK_CLOSE(cos5.derivative(3u), sin(cx), eps);
348 BOOST_CHECK_CLOSE(cos5.derivative(4u), cos(cx), eps);
349 BOOST_CHECK_CLOSE(cos5.derivative(5u), -sin(cx), eps);
350 auto sin5 = sin(x);
351 BOOST_CHECK_CLOSE(sin5.derivative(0u), sin(cx), eps);
352 BOOST_CHECK_CLOSE(sin5.derivative(1u), cos(cx), eps);
353 BOOST_CHECK_CLOSE(sin5.derivative(2u), -sin(cx), eps);
354 BOOST_CHECK_CLOSE(sin5.derivative(3u), -cos(cx), eps);
355 BOOST_CHECK_CLOSE(sin5.derivative(4u), sin(cx), eps);
356 BOOST_CHECK_CLOSE(sin5.derivative(5u), cos(cx), eps);
357 // Test Order = 0 for codecov
358 auto cos0 = cos(make_fvar<T, 0>(cx));
359 BOOST_CHECK_CLOSE(cos0.derivative(0u), cos(cx), eps);
360 auto sin0 = sin(make_fvar<T, 0>(cx));
361 BOOST_CHECK_CLOSE(sin0.derivative(0u), sin(cx), eps);
362 }
363
BOOST_AUTO_TEST_CASE_TEMPLATE(acos_test,T,bin_float_types)364 BOOST_AUTO_TEST_CASE_TEMPLATE(acos_test, T, bin_float_types) {
365 const T eps = 300 * std::numeric_limits<T>::epsilon(); // percent
366 using std::acos;
367 using std::pow;
368 using std::sqrt;
369 constexpr std::size_t m = 5;
370 const T cx = 0.5;
371 auto x = make_fvar<T, m>(cx);
372 auto y = acos(x);
373 BOOST_CHECK_CLOSE(y.derivative(0u), acos(cx), eps);
374 BOOST_CHECK_CLOSE(y.derivative(1u), -1 / sqrt(1 - cx * cx), eps);
375 BOOST_CHECK_CLOSE(y.derivative(2u), -cx / pow(1 - cx * cx, 1.5), eps);
376 BOOST_CHECK_CLOSE(y.derivative(3u),
377 -(2 * cx * cx + 1) / pow(1 - cx * cx, 2.5), eps);
378 BOOST_CHECK_CLOSE(y.derivative(4u),
379 -3 * cx * (2 * cx * cx + 3) / pow(1 - cx * cx, 3.5), eps);
380 BOOST_CHECK_CLOSE(y.derivative(5u),
381 -(24 * (cx * cx + 3) * cx * cx + 9) / pow(1 - cx * cx, 4.5),
382 eps);
383 }
384
BOOST_AUTO_TEST_CASE_TEMPLATE(acosh_test,T,bin_float_types)385 BOOST_AUTO_TEST_CASE_TEMPLATE(acosh_test, T, bin_float_types) {
386 const T eps = 300 * std::numeric_limits<T>::epsilon(); // percent
387 using boost::math::acosh;
388 constexpr std::size_t m = 5;
389 const T cx = 2;
390 auto x = make_fvar<T, m>(cx);
391 auto y = acosh(x);
392 // BOOST_CHECK_EQUAL(y.derivative(0) == acosh(cx)); // FAILS! acosh(2) is
393 // overloaded for integral types
394 BOOST_CHECK_CLOSE(y.derivative(0u), acosh(static_cast<T>(x)), eps);
395 BOOST_CHECK_CLOSE(y.derivative(1u),
396 1 / boost::math::constants::root_three<T>(), eps);
397 BOOST_CHECK_CLOSE(y.derivative(2u),
398 -2 / (3 * boost::math::constants::root_three<T>()), eps);
399 BOOST_CHECK_CLOSE(y.derivative(3u),
400 1 / boost::math::constants::root_three<T>(), eps);
401 BOOST_CHECK_CLOSE(y.derivative(4u),
402 -22 / (9 * boost::math::constants::root_three<T>()), eps);
403 BOOST_CHECK_CLOSE(y.derivative(5u),
404 227 / (27 * boost::math::constants::root_three<T>()),
405 2 * eps);
406 }
407
BOOST_AUTO_TEST_CASE_TEMPLATE(asin_test,T,bin_float_types)408 BOOST_AUTO_TEST_CASE_TEMPLATE(asin_test, T, bin_float_types) {
409 const T eps = 300 * std::numeric_limits<T>::epsilon(); // percent
410 using std::asin;
411 using std::pow;
412 using std::sqrt;
413 constexpr std::size_t m = 5;
414 const T cx = 0.5;
415 auto x = make_fvar<T, m>(cx);
416 auto y = asin(x);
417 BOOST_CHECK_CLOSE(y.derivative(0u), asin(static_cast<T>(x)), eps);
418 BOOST_CHECK_CLOSE(y.derivative(1u), 1 / sqrt(1 - cx * cx), eps);
419 BOOST_CHECK_CLOSE(y.derivative(2u), cx / pow(1 - cx * cx, 1.5), eps);
420 BOOST_CHECK_CLOSE(y.derivative(3u), (2 * cx * cx + 1) / pow(1 - cx * cx, 2.5),
421 eps);
422 BOOST_CHECK_CLOSE(y.derivative(4u),
423 3 * cx * (2 * cx * cx + 3) / pow(1 - cx * cx, 3.5), eps);
424 BOOST_CHECK_CLOSE(y.derivative(5u),
425 (24 * (cx * cx + 3) * cx * cx + 9) / pow(1 - cx * cx, 4.5),
426 eps);
427 }
428
BOOST_AUTO_TEST_CASE_TEMPLATE(asin_infinity,T,all_float_types)429 BOOST_AUTO_TEST_CASE_TEMPLATE(asin_infinity, T, all_float_types) {
430 const T eps = 100 * std::numeric_limits<T>::epsilon(); // percent
431 constexpr std::size_t m = 5;
432 auto x = make_fvar<T, m>(1);
433 auto y = asin(x);
434 // std::cout << "asin(1) = " << y << std::endl; //
435 // depth(1)(1.5707963267949,inf,inf,-nan,-nan,-nan)
436 BOOST_CHECK_CLOSE(y.derivative(0u), boost::math::constants::half_pi<T>(),
437 eps); // MacOS is not exact
438 BOOST_CHECK_EQUAL(y.derivative(1u), std::numeric_limits<T>::infinity());
439 }
440
BOOST_AUTO_TEST_CASE_TEMPLATE(asin_derivative,T,bin_float_types)441 BOOST_AUTO_TEST_CASE_TEMPLATE(asin_derivative, T, bin_float_types) {
442 const T eps = 300 * std::numeric_limits<T>::epsilon(); // percent
443 using std::pow;
444 using std::sqrt;
445 constexpr std::size_t m = 4;
446 const T cx(0.5);
447 auto x = make_fvar<T, m>(cx);
448 auto y = T(1) - x * x;
449 BOOST_CHECK_EQUAL(y.derivative(0u), 1 - cx * cx);
450 BOOST_CHECK_EQUAL(y.derivative(1u), -2 * cx);
451 BOOST_CHECK_EQUAL(y.derivative(2u), -2);
452 BOOST_CHECK_EQUAL(y.derivative(3u), 0);
453 BOOST_CHECK_EQUAL(y.derivative(4u), 0);
454 y = sqrt(y);
455 BOOST_CHECK_EQUAL(y.derivative(0u), sqrt(1 - cx * cx));
456 BOOST_CHECK_CLOSE(y.derivative(1u), -cx / sqrt(1 - cx * cx), eps);
457 BOOST_CHECK_CLOSE(y.derivative(2u), -1 / pow(1 - cx * cx, 1.5), eps);
458 BOOST_CHECK_CLOSE(y.derivative(3u), -3 * cx / pow(1 - cx * cx, 2.5), eps);
459 BOOST_CHECK_CLOSE(y.derivative(4u),
460 -(12 * cx * cx + 3) / pow(1 - cx * cx, 3.5), eps);
461 y = y.inverse(); // asin'(x) = 1 / sqrt(1-x*x).
462 BOOST_CHECK_CLOSE(y.derivative(0u), 1 / sqrt(1 - cx * cx), eps);
463 BOOST_CHECK_CLOSE(y.derivative(1u), cx / pow(1 - cx * cx, 1.5), eps);
464 BOOST_CHECK_CLOSE(y.derivative(2u), (2 * cx * cx + 1) / pow(1 - cx * cx, 2.5),
465 eps);
466 BOOST_CHECK_CLOSE(y.derivative(3u),
467 3 * cx * (2 * cx * cx + 3) / pow(1 - cx * cx, 3.5), eps);
468 BOOST_CHECK_CLOSE(y.derivative(4u),
469 (24 * (cx * cx + 3) * cx * cx + 9) / pow(1 - cx * cx, 4.5),
470 eps);
471 }
472
BOOST_AUTO_TEST_CASE_TEMPLATE(asinh_test,T,bin_float_types)473 BOOST_AUTO_TEST_CASE_TEMPLATE(asinh_test, T, bin_float_types) {
474 const T eps = 300 * std::numeric_limits<T>::epsilon(); // percent
475 using boost::math::asinh;
476 constexpr std::size_t m = 5;
477 const T cx = 1;
478 auto x = make_fvar<T, m>(cx);
479 auto y = asinh(x);
480 BOOST_CHECK_CLOSE(y.derivative(0u), asinh(static_cast<T>(x)), eps);
481 BOOST_CHECK_CLOSE(y.derivative(1u), 1 / boost::math::constants::root_two<T>(),
482 eps);
483 BOOST_CHECK_CLOSE(y.derivative(2u),
484 -1 / (2 * boost::math::constants::root_two<T>()), eps);
485 BOOST_CHECK_CLOSE(y.derivative(3u),
486 1 / (4 * boost::math::constants::root_two<T>()), eps);
487 BOOST_CHECK_CLOSE(y.derivative(4u),
488 3 / (8 * boost::math::constants::root_two<T>()), eps);
489 BOOST_CHECK_CLOSE(y.derivative(5u),
490 -39 / (16 * boost::math::constants::root_two<T>()), eps);
491 }
492
BOOST_AUTO_TEST_CASE_TEMPLATE(atan2_function,T,all_float_types)493 BOOST_AUTO_TEST_CASE_TEMPLATE(atan2_function, T, all_float_types) {
494 using test_constants = test_constants_t<T>;
495 static constexpr auto m = test_constants::order;
496
497 test_detail::RandomSample<T> x_sampler{-2000, 2000};
498 test_detail::RandomSample<T> y_sampler{-2000, 2000};
499
500 for (auto i : boost::irange(test_constants::n_samples)) {
501 std::ignore = i;
502 auto x = x_sampler.next();
503 auto y = y_sampler.next();
504
505 auto autodiff_v = atan2(make_fvar<T, m>(x), make_fvar<T, m>(y));
506 auto anchor_v = atan2(x, y);
507 BOOST_CHECK_CLOSE(autodiff_v, anchor_v,
508 5000 * test_constants::pct_epsilon());
509 }
510 }
511
512 BOOST_AUTO_TEST_SUITE_END()
513