libs/math/test/test_bessel_k_prime.hpp

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