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 #include <boost/utility/identity_type.hpp>
8 
9 BOOST_AUTO_TEST_SUITE(test_autodiff_3)
10 
BOOST_AUTO_TEST_CASE_TEMPLATE(atanh_test,T,all_float_types)11 BOOST_AUTO_TEST_CASE_TEMPLATE(atanh_test, T, all_float_types) {
12   const T eps = 3000 * test_constants_t<T>::pct_epsilon(); // percent
13   constexpr unsigned m = 5;
14   const T cx = 0.5;
15   auto x = make_fvar<T, m>(cx);
16   auto y = atanh(x);
17   // BOOST_CHECK_EQUAL(y.derivative(0) , atanh(cx)); // fails due to overload
18   BOOST_CHECK_CLOSE(y.derivative(0u), atanh(static_cast<T>(x)), eps);
19   BOOST_CHECK_CLOSE(y.derivative(1u), static_cast<T>(4) / 3, eps);
20   BOOST_CHECK_CLOSE(y.derivative(2u), static_cast<T>(16) / 9, eps);
21   BOOST_CHECK_CLOSE(y.derivative(3u), static_cast<T>(224) / 27, eps);
22   BOOST_CHECK_CLOSE(y.derivative(4u), static_cast<T>(1280) / 27, eps);
23   BOOST_CHECK_CLOSE(y.derivative(5u), static_cast<T>(31232) / 81, eps);
24 }
25 
BOOST_AUTO_TEST_CASE_TEMPLATE(atan_test,T,all_float_types)26 BOOST_AUTO_TEST_CASE_TEMPLATE(atan_test, T, all_float_types) {
27   BOOST_MATH_STD_USING
28   using namespace boost;
29 
30   const T cx = 1.0;
31   constexpr unsigned m = 5;
32   const auto x = make_fvar<T, m>(cx);
33   auto y = atan(x);
34   const auto eps = boost::math::tools::epsilon<T>() * 200; // 2eps as a percentage
35   BOOST_CHECK_CLOSE(y.derivative(0u), boost::math::constants::pi<T>() / 4, eps);
36   BOOST_CHECK_CLOSE(y.derivative(1u), T(0.5), eps);
37   BOOST_CHECK_CLOSE(y.derivative(2u), T(-0.5), eps);
38   BOOST_CHECK_CLOSE(y.derivative(3u), T(0.5), eps);
39   BOOST_CHECK_CLOSE(y.derivative(4u), T(0), eps);
40   BOOST_CHECK_CLOSE(y.derivative(5u), T(-3), eps);
41 }
42 
BOOST_AUTO_TEST_CASE_TEMPLATE(erf_test,T,all_float_types)43 BOOST_AUTO_TEST_CASE_TEMPLATE(erf_test, T, all_float_types) {
44   BOOST_MATH_STD_USING
45   using namespace boost;
46 
47   const T eps = 300 * 100 * boost::math::tools::epsilon<T>(); // percent
48   const T cx = 1.0;
49   constexpr unsigned m = 5;
50   const auto x = make_fvar<T, m>(cx);
51   auto y = erf(x);
52   BOOST_CHECK_CLOSE(y.derivative(0u), erf(static_cast<T>(x)), eps);
53   BOOST_CHECK_CLOSE(
54       y.derivative(1u),
55       T(2) / (math::constants::e<T>() * math::constants::root_pi<T>()), eps);
56   BOOST_CHECK_CLOSE(
57       y.derivative(2u),
58       T(-4) / (math::constants::e<T>() * math::constants::root_pi<T>()), eps);
59   BOOST_CHECK_CLOSE(
60       y.derivative(3u),
61       T(4) / (math::constants::e<T>() * math::constants::root_pi<T>()), eps);
62   BOOST_CHECK_CLOSE(
63       y.derivative(4u),
64       T(8) / (math::constants::e<T>() * math::constants::root_pi<T>()), eps);
65   BOOST_CHECK_CLOSE(
66       y.derivative(5u),
67       T(-40) / (math::constants::e<T>() * math::constants::root_pi<T>()), eps);
68 }
69 
BOOST_AUTO_TEST_CASE_TEMPLATE(sinc_test,T,bin_float_types)70 BOOST_AUTO_TEST_CASE_TEMPLATE(sinc_test, T, bin_float_types) {
71   BOOST_MATH_STD_USING
72   const T eps = 20000 * boost::math::tools::epsilon<T>(); // percent
73   const T cx = 1;
74   constexpr unsigned m = 5;
75   auto x = make_fvar<T, m>(cx);
76   auto y = sinc(x);
77   BOOST_CHECK_CLOSE(y.derivative(0u), sin(cx), eps);
78   BOOST_CHECK_CLOSE(y.derivative(1u), cos(cx) - sin(cx), eps);
79   BOOST_CHECK_CLOSE(y.derivative(2u), sin(cx) - 2 * cos(cx), eps);
80   BOOST_CHECK_CLOSE(y.derivative(3u), T(5) * cos(cx) - T(3) * sin(cx), eps);
81   BOOST_CHECK_CLOSE(y.derivative(4u), T(13) * sin(cx) - T(20) * cos(cx), eps);
82   BOOST_CHECK_CLOSE(y.derivative(5u), T(101) * cos(cx) - T(65) * sin(cx), eps);
83   // Test at x = 0
84   auto y2 = sinc(make_fvar<T, 10>(0));
85   BOOST_CHECK_CLOSE(y2.derivative(0u), T(1), eps);
86   BOOST_CHECK_CLOSE(y2.derivative(1u), T(0), eps);
87   BOOST_CHECK_CLOSE(y2.derivative(2u), -cx / T(3), eps);
88   BOOST_CHECK_CLOSE(y2.derivative(3u), T(0), eps);
89   BOOST_CHECK_CLOSE(y2.derivative(4u), cx / T(5), eps);
90   BOOST_CHECK_CLOSE(y2.derivative(5u), T(0), eps);
91   BOOST_CHECK_CLOSE(y2.derivative(6u), -cx / T(7), eps);
92   BOOST_CHECK_CLOSE(y2.derivative(7u), T(0), eps);
93   BOOST_CHECK_CLOSE(y2.derivative(8u), cx / T(9), eps);
94   BOOST_CHECK_CLOSE(y2.derivative(9u), T(0), eps);
95   BOOST_CHECK_CLOSE(y2.derivative(10u), -cx / T(11), eps);
96 }
97 
BOOST_AUTO_TEST_CASE_TEMPLATE(sinh_and_cosh,T,bin_float_types)98 BOOST_AUTO_TEST_CASE_TEMPLATE(sinh_and_cosh, T, bin_float_types) {
99   BOOST_MATH_STD_USING
100   const T eps = 300 * boost::math::tools::epsilon<T>(); // percent
101   const T cx = 1;
102   constexpr unsigned m = 5;
103   auto x = make_fvar<T, m>(cx);
104   auto s = sinh(x);
105   auto c = cosh(x);
106   BOOST_CHECK_CLOSE(s.derivative(0u), sinh(static_cast<T>(x)), eps);
107   BOOST_CHECK_CLOSE(c.derivative(0u), cosh(static_cast<T>(x)), eps);
108   for (auto i : boost::irange(m + 1)) {
109     BOOST_CHECK_CLOSE(s.derivative(i), static_cast<T>(i % 2 == 1 ? c : s), eps);
110     BOOST_CHECK_CLOSE(c.derivative(i), static_cast<T>(i % 2 == 1 ? s : c), eps);
111   }
112 }
113 
114 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
BOOST_AUTO_TEST_CASE_TEMPLATE(tanh_test,T,all_float_types)115 BOOST_AUTO_TEST_CASE_TEMPLATE(tanh_test, T, all_float_types) {
116   using bmp::fabs;
117   using bmp::tanh;
118   using detail::fabs;
119   using detail::tanh;
120   using std::fabs;
121   using std::tanh;
122   constexpr std::array<const char *, 6> tanh_derivatives{
123       {"0."
124        "76159415595576488811945828260479359041276859725793655159681050012195324"
125        "457663848345894752167367671442190275970155",
126        "0."
127        "41997434161402606939449673904170144491718672823077095471331144024458989"
128        "95240483056156940088623187260",
129        "-0."
130        "63970000844922450018849176930384395321921136306079914494299856318702069"
131        "34885434644440069533372017992",
132        "0."
133        "62162668077129626310653042872222339967572411755445418563968706335816206"
134        "22188951465548376863495698762",
135        "0."
136        "66509104475050167773507148092106234992757132833203125448814929383096463"
137        "47626843278089998045994094537",
138        "-5."
139        "55689355847371979760458290231697200987383372116293456019531342394708989"
140        "7942786231796317250984197038"}};
141   const T cx = 1;
142   constexpr std::size_t m = 5;
143   auto x = make_fvar<T, m>(cx);
144   auto t = tanh(x);
145   for (auto i : boost::irange(tanh_derivatives.size())) {
146     BOOST_TEST_WARN(isNearZero(t.derivative(i) -
147                                boost::lexical_cast<T>(tanh_derivatives[i])));
148   }
149 }
150 #endif
151 
BOOST_AUTO_TEST_CASE_TEMPLATE(tan_test,T,bin_float_types)152 BOOST_AUTO_TEST_CASE_TEMPLATE(tan_test, T, bin_float_types) {
153   BOOST_MATH_STD_USING
154   const T eps = 800 * boost::math::tools::epsilon<T>(); // percent
155   const T cx = boost::math::constants::third_pi<T>();
156   const T root_three = boost::math::constants::root_three<T>();
157   constexpr unsigned m = 5;
158   const auto x = make_fvar<T, m>(cx);
159   auto y = tan(x);
160   BOOST_CHECK_CLOSE(y.derivative(0u), root_three, eps);
161   BOOST_CHECK_CLOSE(y.derivative(1u), T(4), eps);
162   BOOST_CHECK_CLOSE(y.derivative(2u), T(8) * root_three, eps);
163   BOOST_CHECK_CLOSE(y.derivative(3u), T(80), eps);
164   BOOST_CHECK_CLOSE(y.derivative(4u), T(352) * root_three, eps);
165   BOOST_CHECK_CLOSE(y.derivative(5u), T(5824), eps);
166 }
167 
BOOST_AUTO_TEST_CASE_TEMPLATE(fmod_test,T,bin_float_types)168 BOOST_AUTO_TEST_CASE_TEMPLATE(fmod_test, T, bin_float_types) {
169   BOOST_MATH_STD_USING
170   constexpr unsigned m = 3;
171   const T cx = 3.25;
172   const T cy = 0.5;
173   auto x = make_fvar<T, m>(cx);
174   auto y = fmod(x, autodiff_fvar<T, m>(cy));
175   BOOST_CHECK_EQUAL(y.derivative(0u), T(0.25));
176   BOOST_CHECK_EQUAL(y.derivative(1u), T(1));
177   BOOST_CHECK_EQUAL(y.derivative(2u), T(0));
178   BOOST_CHECK_EQUAL(y.derivative(3u), T(0));
179 }
180 
BOOST_AUTO_TEST_CASE_TEMPLATE(round_and_trunc,T,all_float_types)181 BOOST_AUTO_TEST_CASE_TEMPLATE(round_and_trunc, T, all_float_types) {
182   BOOST_MATH_STD_USING
183   constexpr unsigned m = 3;
184   const T cx = 3.25;
185   auto x = make_fvar<T, m>(cx);
186   auto y = round(x);
187   BOOST_CHECK_EQUAL(y.derivative(0u), round(cx));
188   BOOST_CHECK_EQUAL(y.derivative(1u), T(0));
189   BOOST_CHECK_EQUAL(y.derivative(2u), T(0));
190   BOOST_CHECK_EQUAL(y.derivative(3u), T(0));
191   y = trunc(x);
192   BOOST_CHECK_EQUAL(y.derivative(0u), trunc(cx));
193   BOOST_CHECK_EQUAL(y.derivative(1u), T(0));
194   BOOST_CHECK_EQUAL(y.derivative(2u), T(0));
195   BOOST_CHECK_EQUAL(y.derivative(3u), T(0));
196 }
197 
BOOST_AUTO_TEST_CASE_TEMPLATE(iround_and_itrunc,T,all_float_types)198 BOOST_AUTO_TEST_CASE_TEMPLATE(iround_and_itrunc, T, all_float_types) {
199   BOOST_MATH_STD_USING
200   using namespace boost::math;
201   constexpr unsigned m = 3;
202   const T cx = 3.25;
203   auto x = make_fvar<T, m>(cx);
204   int y = iround(x);
205   BOOST_CHECK_EQUAL(y, iround(cx));
206   y = itrunc(x);
207   BOOST_CHECK_EQUAL(y, itrunc(cx));
208 }
209 
BOOST_AUTO_TEST_CASE_TEMPLATE(lambert_w0_test,T,all_float_types)210 BOOST_AUTO_TEST_CASE_TEMPLATE(lambert_w0_test, T, all_float_types) {
211   const T eps = 1000 * boost::math::tools::epsilon<T>(); // percent
212   constexpr unsigned m = 10;
213   const T cx = 3;
214   // Mathematica: N[Table[D[ProductLog[x], {x, n}], {n, 0, 10}] /. x -> 3, 52]
215   constexpr std::array<const char *, m + 1> answers{
216       {"1.049908894964039959988697070552897904589466943706341",
217        "0.1707244807388472968312949774415522047470762509741737",
218        "-0.04336545501146252734105411312976167858858970875797718",
219        "0.02321456264324789334313200360870492961288748451791104",
220        "-0.01909049778427783072663170526188353869136655225133878",
221        "0.02122935002563637629500975949987796094687564718834156",
222        "-0.02979093848448877259041971538394953658978044986784643",
223        "0.05051290266216717699803334605370337985567016837482099",
224        "-0.1004503154972645060971099914384090562800544486549660",
225        "0.2292464437392250211967939182075930820454464472006425",
226        "-0.5905839053125614593682763387470654123192290838719517"}};
227   auto x = make_fvar<T, m>(cx);
228   auto y = lambert_w0(x);
229   for (auto i : boost::irange(m + 1)) {
230     const T answer = boost::lexical_cast<T>(answers[i]);
231     BOOST_CHECK_CLOSE(y.derivative(i), answer, eps);
232   }
233   // const T cx0 = -1 / boost::math::constants::e<T>();
234   // auto edge = lambert_w0(make_fvar<T,m>(cx0));
235   // std::cout << "edge = " << edge << std::endl;
236   // edge = depth(1)(-1,inf,-inf,inf,-inf,inf,-inf,inf,-inf,inf,-inf)
237   // edge = depth(1)(-1,inf,-inf,inf,-inf,inf,-inf,inf,-inf,inf,-inf)
238   // edge =
239   // depth(1)(-1,3.68935e+19,-9.23687e+57,4.62519e+96,-2.89497e+135,2.02945e+174,-1.52431e+213,1.19943e+252,-9.75959e+290,8.14489e+329,-6.93329e+368)
240 }
241 
BOOST_AUTO_TEST_CASE_TEMPLATE(digamma_test,T,all_float_types)242 BOOST_AUTO_TEST_CASE_TEMPLATE(digamma_test, T, all_float_types) {
243   const T eps = 1000 * boost::math::tools::epsilon<T>(); // percent
244   constexpr unsigned m = 10;
245   const T cx = 3;
246   // Mathematica: N[Table[PolyGamma[n, 3], {n, 0, 10}], 52]
247   constexpr std::array<const char *, m + 1> answers{
248     {"0.9227843350984671393934879099175975689578406640600764"
249     ,"0.3949340668482264364724151666460251892189499012067984"
250     ,"-0.1541138063191885707994763230228999815299725846809978"
251     ,"0.1189394022668291490960221792470074166485057115123614"
252     ,"-0.1362661234408782319527716749688200333699420680459075"
253     ,"0.2061674381338967657421515749104633482180988039424274"
254     ,"-0.3864797149844353246542358918536669119017636069718686"
255     ,"0.8623752376394704685736020836084249051623848752441025"
256     ,"-2.228398747634885327823655450854278779627928241914664"
257     ,"6.536422382626807143525565747764891144367614117601463"
258     ,"-21.4366066287129906188428320541054572790340793874298"}};
259   auto x = make_fvar<T, m>(cx);
260   auto y = digamma(x);
261   for (auto i : boost::irange(m + 1)) {
262     const T answer = boost::lexical_cast<T>(answers[i]);
263     BOOST_CHECK_CLOSE(y.derivative(i), answer, eps);
264   }
265 }
266 
BOOST_AUTO_TEST_CASE_TEMPLATE(lgamma_test,T,all_float_types)267 BOOST_AUTO_TEST_CASE_TEMPLATE(lgamma_test, T, all_float_types) {
268   const T eps = 1000 * boost::math::tools::epsilon<T>(); // percent
269   constexpr unsigned m = 10;
270   const T cx = 3;
271   // Mathematica: N[Table[D[LogGamma[x],{x,n}] /. x->3, {n, 0, 10}], 52]
272   constexpr std::array<const char *, m + 1> answers{
273     {"0.6931471805599453094172321214581765680755001343602553"
274     ,"0.9227843350984671393934879099175975689578406640600764"
275     ,"0.3949340668482264364724151666460251892189499012067984"
276     ,"-0.1541138063191885707994763230228999815299725846809978"
277     ,"0.1189394022668291490960221792470074166485057115123614"
278     ,"-0.1362661234408782319527716749688200333699420680459075"
279     ,"0.2061674381338967657421515749104633482180988039424274"
280     ,"-0.3864797149844353246542358918536669119017636069718686"
281     ,"0.8623752376394704685736020836084249051623848752441025"
282     ,"-2.228398747634885327823655450854278779627928241914664"
283     ,"6.536422382626807143525565747764891144367614117601463"}};
284   auto x = make_fvar<T, m>(cx);
285   auto y = lgamma(x);
286   for (auto i : boost::irange(m + 1)) {
287     const T answer = boost::lexical_cast<T>(answers[i]);
288     BOOST_CHECK_CLOSE(y.derivative(i), answer, eps);
289   }
290 }
291 
BOOST_AUTO_TEST_CASE_TEMPLATE(tgamma_test,T,all_float_types)292 BOOST_AUTO_TEST_CASE_TEMPLATE(tgamma_test, T, all_float_types) {
293   const T eps = 1000 * boost::math::tools::epsilon<T>(); // percent
294   constexpr unsigned m = 10;
295   const T cx = 3;
296   // Mathematica: N[Table[D[Gamma[x],{x,n}] /. x->3, {n, 0, 10}], 52]
297   constexpr std::array<const char *, m + 1> answers{
298     {"2.0"
299     ,"1.845568670196934278786975819835195137915681328120153"
300     ,"2.492929991902693057942510065508124245503778067273315"
301     ,"3.449965013523673365279327178241708777509009968597547"
302     ,"5.521798578098737512443417699412265532987916790978887"
303     ,"8.845805593922864253981346455183370214190789096412155"
304     ,"15.86959874461221647760760269963155031595848150772695"
305     ,"27.46172054213435946038727460195592342721862288816812"
306     ,"54.64250508485402729556251663145824730270508661240771"
307     ,"96.08542140594972502872131946513104238293824803599579"
308     ,"222.0936743583156040996433943128676567542497584689499"}};
309   auto x = make_fvar<T, m>(cx);
310   auto y = tgamma(x);
311   for (auto i : boost::irange(m + 1)) {
312     const T answer = boost::lexical_cast<T>(answers[i]);
313     BOOST_CHECK_CLOSE(y.derivative(i), answer, eps);
314   }
315 }
316 
BOOST_AUTO_TEST_CASE_TEMPLATE(tgamma2_test,T,all_float_types)317 BOOST_AUTO_TEST_CASE_TEMPLATE(tgamma2_test, T, all_float_types) {
318   //const T eps = 5000 * boost::math::tools::epsilon<T>(); // ok for non-multiprecision
319   const T eps = 500000 * boost::math::tools::epsilon<T>(); // percent
320   constexpr unsigned m = 10;
321   const T cx = -1.5;
322   // Mathematica: N[Table[D[Gamma[x],{x,n}] /. x->-3/2, {n, 0, 10}], 52]
323   constexpr std::array<const char *, m + 1> answers{
324     {"2.363271801207354703064223311121526910396732608163183"
325     ,"1.661750260668596505586468565464938761014714509096807"
326     ,"23.33417984355457252918927856618603412638766668207679"
327     ,"47.02130025080143055642555842335081335790754507072526"
328     ,"1148.336052788822231948472800239024335856568111484074"
329     ,"3831.214710988836934569706027888431190714054814541186"
330     ,"138190.9008816865362698874238213771413807566436072179"
331     ,"644956.0066517306036921195893233874126907491308967028"
332     ,"3.096453684470713902448094810299787572782887316764214e7"
333     ,"1.857893143852025058151037296906468662709947415219451e8"
334     ,"1.114762466163487983067783853825224537320312784955935e10"}};
335   auto x = make_fvar<T, m>(cx);
336   auto y = tgamma(x);
337   for (auto i : boost::irange(m + 1)) {
338     const T answer = boost::lexical_cast<T>(answers[i]);
339     BOOST_CHECK_CLOSE(y.derivative(i), answer, eps);
340   }
341 }
342 
343 BOOST_AUTO_TEST_SUITE_END()
344