ceph/src/boost/libs/math/test/test_autodiff_3.cpp

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