test_bessel_k_prime.hpp 9.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179
  1. // Copyright (c) 2013 Anton Bikineev
  2. // Use, modification and distribution are subject to the
  3. // Boost Software License, Version 1.0. (See accompanying file
  4. // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
  5. #define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
  6. #include <boost/math/concepts/real_concept.hpp>
  7. #define BOOST_TEST_MAIN
  8. #include <boost/test/unit_test.hpp>
  9. #include <boost/test/tools/floating_point_comparison.hpp>
  10. #include <boost/math/special_functions/math_fwd.hpp>
  11. #include <boost/type_traits/is_floating_point.hpp>
  12. #include <boost/array.hpp>
  13. #include "functor.hpp"
  14. #include "handle_test_result.hpp"
  15. #include "table_type.hpp"
  16. #ifndef SC_
  17. # define SC_(x) static_cast<typename table_type<T>::type>(BOOST_JOIN(x, L))
  18. #endif
  19. template <class T>
  20. T cyl_bessel_k_prime_int_wrapper(T v, T x)
  21. {
  22. #ifdef BESSEL_KPN_FUNCTION_TO_TEST
  23. return static_cast<T>(
  24. BESSEL_KPN_FUNCTION_TO_TEST(
  25. boost::math::itrunc(v), x));
  26. #else
  27. return static_cast<T>(
  28. boost::math::cyl_bessel_k_prime(
  29. boost::math::itrunc(v), x));
  30. #endif
  31. }
  32. template <class Real, class T>
  33. void do_test_cyl_bessel_k_prime(const T& data, const char* type_name, const char* test_name)
  34. {
  35. #if !(defined(ERROR_REPORTING_MODE) && !defined(BESSEL_KP_FUNCTION_TO_TEST))
  36. typedef Real value_type;
  37. typedef value_type (*pg)(value_type, value_type);
  38. #ifdef BESSEL_KP_FUNCTION_TO_TEST
  39. pg funcp = BESSEL_KP_FUNCTION_TO_TEST;
  40. #elif defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
  41. pg funcp = boost::math::cyl_bessel_k_prime<value_type, value_type>;
  42. #else
  43. pg funcp = boost::math::cyl_bessel_k_prime;
  44. #endif
  45. boost::math::tools::test_result<value_type> result;
  46. std::cout << "Testing " << test_name << " with type " << type_name
  47. << "\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
  48. //
  49. // test cyl_bessel_k_prime against data:
  50. //
  51. result = boost::math::tools::test_hetero<Real>(
  52. data,
  53. bind_func<Real>(funcp, 0, 1),
  54. extract_result<Real>(2));
  55. handle_test_result(result, data[result.worst()], result.worst(), type_name, "cyl_bessel_k_prime", test_name);
  56. std::cout << std::endl;
  57. #endif
  58. }
  59. template <class Real, class T>
  60. void do_test_cyl_bessel_k_prime_int(const T& data, const char* type_name, const char* test_name)
  61. {
  62. #if !(defined(ERROR_REPORTING_MODE) && !defined(BESSEL_KPN_FUNCTION_TO_TEST))
  63. typedef Real value_type;
  64. typedef value_type (*pg)(value_type, value_type);
  65. #if defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
  66. pg funcp = cyl_bessel_k_prime_int_wrapper<value_type>;
  67. #else
  68. pg funcp = cyl_bessel_k_prime_int_wrapper;
  69. #endif
  70. boost::math::tools::test_result<value_type> result;
  71. std::cout << "Testing " << test_name << " with type " << type_name
  72. << "\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
  73. //
  74. // test cyl_bessel_k_prime against data:
  75. //
  76. result = boost::math::tools::test_hetero<Real>(
  77. data,
  78. bind_func<Real>(funcp, 0, 1),
  79. extract_result<Real>(2));
  80. handle_test_result(result, data[result.worst()], result.worst(), type_name, "cyl_bessel_k_prime (integer orders)", test_name);
  81. std::cout << std::endl;
  82. #endif
  83. }
  84. template <class T>
  85. void test_bessel(T, const char* name)
  86. {
  87. // function values calculated on wolframalpha.com
  88. static const boost::array<boost::array<T, 3>, 9> k0_prime_data = {{
  89. {{ SC_(0.0), SC_(1.0), SC_(-0.60190723019723457473754000153561733926158688996810646) }},
  90. {{ SC_(0.0), SC_(2.0), SC_(-0.1398658818165224272845988070354110238872345848415155) }},
  91. {{ SC_(0.0), SC_(4.0), SC_(-0.012483498887268431470384179980806068483841584988625846) }},
  92. {{ SC_(0.0), SC_(8.0), SC_(-0.00015536921180500113391686245062247462111706512287261616) }},
  93. {{ SC_(0.0), T(std::ldexp(1.0, -15)), SC_(-32767.99983195283164326474413165391397251047283415776) }},
  94. {{ SC_(0.0), T(std::ldexp(1.0, -30)), SC_(-1.0737418239999999900300302857268733281035379954421507e9) }},
  95. {{ SC_(0.0), T(std::ldexp(1.0, -60)), SC_(-1.1529215046068469759999999999999999816966019886812660e18) }},
  96. {{ SC_(0.0), SC_(50.0), SC_(-3.44410222671755561259185303591267155099677251348256880e-23) }},
  97. {{ SC_(0.0), SC_(100.0), SC_(-4.6798537356369092865625442420243353079749435469433535e-45) }},
  98. }};
  99. static const boost::array<boost::array<T, 3>, 9> k1_prime_data = {{
  100. {{ SC_(1.0), SC_(1.0), SC_(-1.0229316684379429080731673807482263753978066381947669) }},
  101. {{ SC_(1.0), SC_(2.0), SC_(-0.1838268136577946492950189784501873449419439168095666) }},
  102. {{ SC_(1.0), SC_(4.0), SC_(-0.014280550807670132137341240975035006345969420135630802) }},
  103. {{ SC_(1.0), SC_(8.0), SC_(-0.00016589185669844052883619221502648724960693850919283604) }},
  104. {{ SC_(1.0), T(std::ldexp(1.0, -15)), SC_(-1.0737418290065696140247028419519880092107054138744140e9) }},
  105. {{ SC_(1.0), T(std::ldexp(1.0, -30)), SC_(-1.1529215046068469862051734662283858692135761720165778e18) }},
  106. {{ SC_(1.0), T(std::ldexp(1.0, -60)), SC_(-1.329227995784915872903807060280344596602381174627566e36) }},
  107. {{ SC_(1.0), SC_(50.0), SC_(-3.47904979432384662617251257307120566286496082789299947e-23) }},
  108. {{ SC_(1.0), SC_(100.0), SC_(-4.7034267665322711118046307319041297088872889209115474e-45) }},
  109. }};
  110. static const boost::array<boost::array<T, 3>, 9> kn_prime_data = {{
  111. {{ SC_(2.0), T(std::ldexp(1.0, -30)), SC_(-4.951760157141521099596496895999999995073222803776904e27) }},
  112. {{ SC_(5.0), SC_(10.0), SC_(-0.0000666323621535481236223011866087784024278980735437002384) }},
  113. {{ SC_(-5.0), SC_(100.0), SC_(-5.3060798744208349930861060378887340340201141387578377e-45) }},
  114. {{ SC_(10.0), SC_(10.0), SC_(-0.00232426413420145080508626300083871228780582972491498296) }},
  115. {{ SC_(10.0), T(std::ldexp(1.0, -30)), SC_(-4.0637928602074079595570948641288439020852370470244381e108) }},
  116. {{ SC_(-10.0), SC_(1.0), SC_(-1.8171379399979651461891429013401068319174853467388121e9) }},
  117. {{ SC_(100.0), SC_(5.0), SC_(-1.4097486373570936520327835736048715219413065916411893e117) }},
  118. {{ SC_(100.0), SC_(80.0), SC_(-1.34557011017664184003144916855685180771861680634827508e-11) }},
  119. {{ SC_(-1000.0), SC_(700.0), SC_(-1.136342773238774160870536985092768591616106526374957e-30) }},
  120. }};
  121. static const boost::array<boost::array<T, 3>, 11> kv_prime_data = {{
  122. {{ SC_(0.5), SC_(0.875), SC_(-0.8776935068732421581818610624499915196588910540138553643355820) }},
  123. {{ SC_(0.5), SC_(1.125), SC_(-0.5541192376058293458786667962590089848709748151724170966916495) }},
  124. {{ SC_(2.25), T(std::ldexp(1.0, -30)), SC_(-1.358706605110306964608847299464328015299661532e30) }},
  125. {{ SC_(5.5), T(3217)/1024, SC_(-2.6903757178739422729800670428157504611799055394319992629519699) }},
  126. {{ SC_(-5.5), SC_(10.0), SC_(-0.000086479759593318257340087317655128751755482676477180134416784728) }},
  127. {{ SC_(-5.5), SC_(100.0), SC_(-5.4478425565190604625309457442097587701859746312164196732075323e-45) }},
  128. {{ T(10240)/1024, T(1)/1024, SC_(-2.411751224440479729811903506282248205762559999997965494837863222e42) }},
  129. {{ T(10240)/1024, SC_(10.0), SC_(-0.002324264134201450805086263000838712287805829724914982961118625775) }},
  130. {{ T(144793)/1024, SC_(100.0), SC_(-2.419425330672365273534646536102117722944744737761477017402710069e-6) }},
  131. {{ T(144793)/1024, SC_(200.0), SC_(-1.1183699286601178683373775100500418982738064865504029155187086e-67) }},
  132. {{ T(-144793)/1024, SC_(50.0), SC_(-3.906473504308773541933992099338237076647113693807893258840087e42) }},
  133. }};
  134. static const boost::array<boost::array<T, 3>, 5> kv_prime_large_data = {{
  135. {{ SC_(-0.5), static_cast<T>(ldexp(0.5, -512)), SC_(-2.75176667129887692508287667455879592490037256500173136025362e231) }},
  136. {{ SC_(0.5), static_cast<T>(ldexp(0.5, -512)), SC_(-2.75176667129887692508287667455879592490037256500173136025362e231) }},
  137. #if LDBL_MAX_10_EXP > 328
  138. {{ SC_(-1.125), static_cast<T>(ldexp(0.5, -512)), SC_(-1.67123513518264734700327664054002130440723e328) }},
  139. {{ SC_(1.125), static_cast<T>(ldexp(0.5, -512)), SC_(-1.67123513518264734700327664054002130440723e328) }},
  140. {{ SC_(0.5), static_cast<T>(ldexp(0.5, -683)), SC_(-4.5061484409559214227217449664854025793393e308) }},
  141. #else
  142. { { SC_(-1.125), static_cast<T>(ldexp(0.5, -512)), std::numeric_limits<T>::has_infinity ? -std::numeric_limits<T>::infinity() : -boost::math::tools::max_value<T>() } },
  143. { { SC_(1.125), static_cast<T>(ldexp(0.5, -512)), std::numeric_limits<T>::has_infinity ? -std::numeric_limits<T>::infinity() : -boost::math::tools::max_value<T>() } },
  144. { { SC_(0.5), static_cast<T>(ldexp(0.5, -683)), std::numeric_limits<T>::has_infinity ? -std::numeric_limits<T>::infinity() : -boost::math::tools::max_value<T>() } },
  145. #endif
  146. }};
  147. do_test_cyl_bessel_k_prime<T>(k0_prime_data, name, "Bessel K'0: Mathworld Data");
  148. do_test_cyl_bessel_k_prime<T>(k1_prime_data, name, "Bessel K'1: Mathworld Data");
  149. do_test_cyl_bessel_k_prime<T>(kn_prime_data, name, "Bessel K'n: Mathworld Data");
  150. do_test_cyl_bessel_k_prime_int<T>(k0_prime_data, name, "Bessel K'0: Mathworld Data (Integer Version)");
  151. do_test_cyl_bessel_k_prime_int<T>(k1_prime_data, name, "Bessel K'1: Mathworld Data (Integer Version)");
  152. do_test_cyl_bessel_k_prime_int<T>(kn_prime_data, name, "Bessel K'n: Mathworld Data (Integer Version)");
  153. do_test_cyl_bessel_k_prime<T>(kv_prime_data, name, "Bessel K'v: Mathworld Data");
  154. if(0 != static_cast<T>(ldexp(0.5, -512)))
  155. do_test_cyl_bessel_k_prime<T>(kv_prime_large_data, name, "Bessel K'v: Mathworld Data (large values)");
  156. #include "bessel_k_prime_int_data.ipp"
  157. do_test_cyl_bessel_k_prime<T>(bessel_k_prime_int_data, name, "Bessel K'n: Random Data");
  158. #include "bessel_k_prime_data.ipp"
  159. do_test_cyl_bessel_k_prime<T>(bessel_k_prime_data, name, "Bessel K'v: Random Data");
  160. }