libs/math/test/test_bessel_j_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/math/constants/constants.hpp>
  13 #include <boost/type_traits/is_floating_point.hpp>
  14 #include <boost/array.hpp>
  15 #include "functor.hpp"
  16
  17 #include "handle_test_result.hpp"
  18 #include "table_type.hpp"
  19
  20 #ifndef SC_
  21 #  define SC_(x) static_cast<typename table_type<T>::type>(BOOST_JOIN(x, L))
  22 #endif
  23
  24 template <class Real, class T>
  25 void do_test_cyl_bessel_j_prime(const T& data, const char* type_name, const char* test_name)
  26 {
  27 #if !(defined(ERROR_REPORTING_MODE) && !defined(BESSEL_JP_FUNCTION_TO_TEST))
  28    typedef Real                   value_type;
  29
  30    typedef value_type (*pg)(value_type, value_type);
  31 #ifdef BESSEL_JP_FUNCTION_TO_TEST
  32    pg funcp = BESSEL_JP_FUNCTION_TO_TEST;
  33 #elif defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
  34    pg funcp = boost::math::cyl_bessel_j_prime<value_type, value_type>;
  35 #else
  36    pg funcp = boost::math::cyl_bessel_j_prime;
  37 #endif
  38
  39    boost::math::tools::test_result<value_type> result;
  40
  41    std::cout << "Testing " << test_name << " with type " << type_name
  42       << "\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
  43
  44    //
  45    // test cyl_bessel_j against data:
  46    //
  47    result = boost::math::tools::test_hetero<Real>(
  48       data,
  49       bind_func<Real>(funcp, 0, 1),
  50       extract_result<Real>(2));
  51    handle_test_result(result, data[result.worst()], result.worst(), type_name, "cyl_bessel_j_prime", test_name);
  52    std::cout << std::endl;
  53 #endif
  54 }
  55
  56 template <class T>
  57 T cyl_bessel_j_prime_int_wrapper(T v, T x)
  58 {
  59 #ifdef BESSEL_JPN_FUNCTION_TO_TEST
  60    return static_cast<T>(BESSEL_JPN_FUNCTION_TO_TEST(boost::math::itrunc(v), x));
  61 #else
  62    return static_cast<T>(boost::math::cyl_bessel_j_prime(boost::math::itrunc(v), x));
  63 #endif
  64 }
  65
  66
  67 template <class Real, class T>
  68 void do_test_cyl_bessel_j_prime_int(const T& data, const char* type_name, const char* test_name)
  69 {
  70 #if !(defined(ERROR_REPORTING_MODE) && !defined(BESSEL_JPN_FUNCTION_TO_TEST))
  71    typedef Real                   value_type;
  72
  73    typedef value_type (*pg)(value_type, value_type);
  74 #if defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
  75    pg funcp = cyl_bessel_j_prime_int_wrapper<value_type>;
  76 #else
  77    pg funcp = cyl_bessel_j_prime_int_wrapper;
  78 #endif
  79
  80    boost::math::tools::test_result<value_type> result;
  81
  82    std::cout << "Testing " << test_name << " with type " << type_name
  83       << "\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
  84
  85    //
  86    // test cyl_bessel_j against data:
  87    //
  88    result = boost::math::tools::test_hetero<Real>(
  89       data,
  90       bind_func<Real>(funcp, 0, 1),
  91       extract_result<Real>(2));
  92    handle_test_result(result, data[result.worst()], result.worst(), type_name, "cyl_bessel_j_prime (integer orders)", test_name);
  93    std::cout << std::endl;
  94 #endif
  95 }
  96
  97 template <class Real, class T>
  98 void do_test_sph_bessel_j_prime(const T& data, const char* type_name, const char* test_name)
  99 {
 100 #if !(defined(ERROR_REPORTING_MODE) && !defined(BESSEL_JPS_FUNCTION_TO_TEST))
 101    typedef Real                   value_type;
 102
 103    typedef value_type (*pg)(unsigned, value_type);
 104 #ifdef BESSEL_JPS_FUNCTION_TO_TEST
 105    pg funcp = BESSEL_JPS_FUNCTION_TO_TEST;
 106 #elif defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
 107    pg funcp = boost::math::sph_bessel_prime<value_type>;
 108 #else
 109    pg funcp = boost::math::sph_bessel_prime;
 110 #endif
 111
 112    boost::math::tools::test_result<value_type> result;
 113
 114    std::cout << "Testing " << test_name << " with type " << type_name
 115       << "\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
 116
 117    //
 118    // test sph_bessel against data:
 119    //
 120    result = boost::math::tools::test_hetero<Real>(
 121       data,
 122       bind_func_int1<Real>(funcp, 0, 1),
 123       extract_result<Real>(2));
 124    handle_test_result(result, data[result.worst()], result.worst(), type_name, "sph_bessel_prime", test_name);
 125    std::cout << std::endl;
 126 #endif
 127 }
 128
 129 template <class T>
 130 void test_bessel_prime(T, const char* name)
 131 {
 132    //
 133    // The actual test data is rather verbose, so it's in a separate file
 134    //
 135    // The contents are as follows, each row of data contains
 136    // three items, input value a, input value b and J'(a, b):
 137    //
 138     // function values calculated on http://functions.wolfram.com/
 139     static const boost::array<boost::array<typename table_type<T>::type, 3>, 8> j0_data = {{
 140        {{ SC_(0.0), SC_(0.0), SC_(0.0) }},
 141         {{ SC_(0.0), SC_(1.0), SC_(-0.440050585744933515959682203718914913127) }},
 142         {{ SC_(0.0), SC_(-2.0), SC_(0.576724807756873387202448242269137086920) }},
 143         {{ SC_(0.0), SC_(4.0), SC_(0.06604332802354913614318542080327502873) }},
 144         {{ SC_(0.0), SC_(-8.0), SC_(0.2346363468539146243812766515904546115488) }},
 145         {{ SC_(0.0), SC_(1e-05), SC_(-0.499999999993750000000026041666666612413194e-5) }},
 146         {{ SC_(0.0), SC_(1e-10), SC_(-0.499999999999999999999375000000000000000000e-10) }},
 147         {{ SC_(0.0), SC_(-1e+01), SC_(0.0434727461688614366697487680258592883062724) }},
 148     }};
 149     static const boost::array<boost::array<T, 3>, 6> j0_tricky = {{
 150         // Big numbers make the accuracy of std::sin the limiting factor:
 151        {{ SC_(0.0), SC_(1e+03), SC_(-0.00472831190708952391757607190121691628542) }},
 152         {{ SC_(0.0), SC_(1e+05), SC_(-0.0018467575628825677163621239671142157437) }},
 153         // test at the regular Bessel roots:
 154         {{ SC_(0.0), T(2521642)/(1024 * 1024), SC_(-0.519147572225778564548541576612898453392794) }},
 155         {{ SC_(0.0), T(5788221)/(1024 * 1024), SC_(0.34026483151709114336072749629487476476084) }},
 156         {{ SC_(0.0), T(9074091)/(1024 * 1024), SC_(-0.271452311894657014854145327490965399410) }},
 157         {{ SC_(0.0), T(12364320)/(1024 * 1024), SC_(0.2324598316641066033541448467171088144257742) }}
 158     }};
 159
 160     static const boost::array<boost::array<typename table_type<T>::type, 3>, 8> j1_data = {{
 161         {{ SC_(1.0), SC_(0.0), SC_(0.5) }},
 162         {{ SC_(1.0), SC_(1.0), SC_(0.325147100813033035490035322383748307781902) }},
 163         {{ SC_(1.0), SC_(-2.0), SC_(-0.064471624737201025549396666484619917634997) }},
 164         {{ SC_(1.0), SC_(4.0), SC_(-0.38063897785796008825079441325087928479376) }},
 165         {{ SC_(1.0), SC_(-8.0), SC_(0.1423212637808145780432098264031651746248) }},
 166         {{ SC_(1.0), SC_(1e-05), SC_(0.499999999981250000000130208333332953559028) }},
 167         {{ SC_(1.0), SC_(1e-10), SC_(0.499999999999999999998125000000000000000001) }},
 168         {{ SC_(1.0), SC_(-1e+01), SC_(-0.250283039068234478864735739287914682660226) }},
 169     }};
 170     static const boost::array<boost::array<T, 3>, 5> j1_tricky = {{
 171         // Big numbers make the accuracy of std::sin the limiting factor:
 172         {{ SC_(1.0), SC_(1e+03), SC_(0.024781957840513085037413155043792491869881) }},
 173         {{ SC_(1.0), SC_(1e+05), SC_(-0.0017192195838116010182477650983128728897) }},
 174         // test at the regular Bessel roots:
 175         {{ SC_(1.0), T(4017834)/(1024*1024), SC_(-0.4027594878673806944036073218740057200193405151367) }},
 176         {{ SC_(1.0), T(7356375)/(1024*1024), SC_(0.3001157854852247730548242543591186404228210449219) }},
 177         {{ SC_(1.0), T(10667654)/(1024*1024), SC_(-0.2497048893045206718888096020236844196915626525879) }},
 178     }};
 179
 180     static const boost::array<boost::array<typename table_type<T>::type, 3>, 17> jn_data = {{
 181         {{ SC_(-1.0), SC_(1.25), SC_(-0.2374074770153809244011000600949046202003956) }},
 182         {{ SC_(2.0), SC_(0.0), SC_(0.0) }},
 183         {{ SC_(-2.0), SC_(0.0), SC_(0.0) }},
 184         {{ SC_(2.0), SC_(1e-02), SC_(0.00249995833352864539930612069540679799606337) }},
 185         {{ SC_(5.0), SC_(10.0), SC_(-0.102571922008611714904101858221407144485) }},
 186         {{ SC_(5.0), SC_(-10.0), SC_(-0.102571922008611714904101858221407144485) }},
 187         {{ SC_(-5.0), SC_(1e+06), SC_(-0.0003310524513007044105585859534523271988) }},
 188         {{ SC_(5.0), SC_(1e+06), SC_(0.0003310524513007044105585859534523271988) }},
 189         {{ SC_(-5.0), SC_(-1.0), SC_(-0.001227850313053782886909720690402218190791576) }},
 190         {{ SC_(10.0), SC_(10.0), SC_(0.08436957863176118824849051273337698304165) }},
 191         {{ SC_(10.0), SC_(-10.0), SC_(-0.08436957863176118824849051273337698304165) }},
 192         {{ SC_(10.0), SC_(-5.0), SC_(-0.00258467784485473925206548854676116157568106) }},
 193         {{ SC_(-10.0), SC_(1e+06), SC_(-0.0007259518037193243350387875733893635962) }},
 194         {{ SC_(10.0), SC_(1e+06), SC_(-0.0007259518037193243350387875733893635962) }},
 195         {{ SC_(1e+02), SC_(8e+01), SC_(3.5036060582489177538508950593467499997755e-06) }},
 196         {{ SC_(1e+03), SC_(1e+05), SC_(-0.0021724469777608393409850758227465776486) }},
 197         {{ SC_(10.0), SC_(1e-100), SC_(2.69114445546737213403880070546737213403880070547e-909) }},
 198     }};
 199     do_test_cyl_bessel_j_prime<T>(j0_data, name, "Bessel J0': Mathworld Data");
 200     do_test_cyl_bessel_j_prime<T>(j0_tricky, name, "Bessel J0': Mathworld Data (Tricky cases)");
 201     do_test_cyl_bessel_j_prime<T>(j1_data, name, "Bessel J1': Mathworld Data");
 202     do_test_cyl_bessel_j_prime<T>(j1_tricky, name, "Bessel J1': Mathworld Data (tricky cases)");
 203     do_test_cyl_bessel_j_prime<T>(jn_data, name, "Bessel JN': Mathworld Data");
 204
 205     do_test_cyl_bessel_j_prime_int<T>(j0_data, name, "Bessel J0': Mathworld Data (Integer Version)");
 206     do_test_cyl_bessel_j_prime_int<T>(j0_tricky, name, "Bessel J0': Mathworld Data (Tricky cases) (Integer Version)");
 207     do_test_cyl_bessel_j_prime_int<T>(j1_data, name, "Bessel J1': Mathworld Data (Integer Version)");
 208     do_test_cyl_bessel_j_prime_int<T>(j1_tricky, name, "Bessel J1': Mathworld Data (tricky cases) (Integer Version)");
 209     do_test_cyl_bessel_j_prime_int<T>(jn_data, name, "Bessel JN': Mathworld Data (Integer Version)");
 210
 211     static const boost::array<boost::array<T, 3>, 21> jv_data = {{
 212          {{ T(22.5), T(0), SC_(0.0) }},
 213          {{ T(2457)/1024, T(1)/1024, SC_(9.35477929043111040277363766198320562099360690e-6) }},
 214          {{ SC_(5.5), T(3217)/1024, SC_(0.042165579369684463582791278988393873) }},
 215          {{ SC_(-5.5), T(3217)/1024, SC_(3.361570113176257957139775812778503494) }},
 216          {{ SC_(-5.5), SC_(1e+04), SC_(0.007593311396019034252155600098309836289) }},
 217          {{ SC_(5.5), SC_(1e+04), SC_(-0.00245022241637437956702428797044365092) }},
 218          {{ SC_(5.5), SC_(1e+06), SC_(-0.000279242826717266554062248256927185394) }},
 219          {{ SC_(5.125), SC_(1e+06), SC_(0.0001830632695189459708211614700642271) }},
 220          {{ SC_(5.875), SC_(1e+06), SC_(-0.0006474276718101871487286860109203539) }},
 221          {{ SC_(0.5), SC_(101.0), SC_(0.070640819172197226936337703929857171981702865) }},
 222          {{ SC_(-5.5), SC_(1e+04), SC_(0.007593311396019034252155600098309836289) }},
 223          {{ SC_(-5.5), SC_(1e+06), SC_(-0.0007474243882060190346457525218941411076) }},
 224          {{ SC_(-0.5), SC_(101.0), SC_(-0.036238035321276062532981494694583591262302408) }},
 225          {{ T(-10486074) / (1024*1024), T(1)/512, SC_(-7.0724447469115535625316241941528969321944e35) }},
 226          {{ T(-10486074) / (1024*1024), SC_(15.0), SC_(-0.15994088796049823354364759206656917967697690) }},
 227          {{ T(10486074) / (1024*1024), SC_(1e+02), SC_(-0.05778764167290516644655950658602424434253) }},
 228          {{ T(10486074) / (1024*1024), SC_(2e+04), SC_(-0.00091101010794789360775314125410690740803) }},
 229          {{ T(-10486074) / (1024*1024), SC_(1e+02), SC_(-0.057736130385111563671838499496767877709471701) }},
 230          {{ SC_(1.5), T(8034)/1024, SC_(0.2783550354042687982259490073096357) }},
 231          {{ SC_(8.5), boost::math::constants::pi<T>() * 4, SC_(-0.194590144622675911618596506265006877277074) }},
 232          {{ SC_(-8.5), boost::math::constants::pi<T>() * 4, SC_(-0.014516314554743677558496402742690038592728) }},
 233     }};
 234     do_test_cyl_bessel_j_prime<T>(jv_data, name, "Bessel J': Mathworld Data");
 235     static const boost::array<boost::array<T, 3>, 4> jv_large_data = {{
 236 #if LDBL_MAX_10_EXP > 308
 237       {{ SC_(-0.5), static_cast<T>(std::ldexp(0.5, -683)), SC_(-2.8687031947358902542073388638943588627056993e308) }},
 238 #else
 239       {{ SC_(-0.5), static_cast<T>(std::ldexp(0.5, -450)), SC_(-1.7688953183288445554095310240218576026580197125814e203) }},
 240 #endif
 241       {{ SC_(256.0), SC_(512.0), SC_(0.032286467266411904239327492993951594201583145) }},
 242       {{ SC_(-256.0), SC_(8.0), SC_(4.6974301387555891979202431551474684165419e-352) }},
 243       {{ SC_(-2.5), SC_(4.0), SC_(-0.3580070651681080294136741901878543615958139) }},
 244     }};
 245     if(jv_large_data[0][1] != 0)
 246       do_test_cyl_bessel_j_prime<T>(jv_large_data, name, "Bessel J': Mathworld Data (large values)");
 247
 248 #include "bessel_j_prime_int_data.ipp"
 249     do_test_cyl_bessel_j_prime<T>(bessel_j_prime_int_data, name, "Bessel JN': Random Data");
 250
 251 #include "bessel_j_prime_data.ipp"
 252     do_test_cyl_bessel_j_prime<T>(bessel_j_prime_data, name, "Bessel J': Random Data");
 253
 254 #include "bessel_j_prime_large_data.ipp"
 255     do_test_cyl_bessel_j_prime<T>(bessel_j_prime_large_data, name, "Bessel J': Random Data (Tricky large values)");
 256
 257 #include "sph_bessel_prime_data.ipp"
 258     do_test_sph_bessel_j_prime<T>(sph_bessel_prime_data, name, "Bessel j': Random Data");
 259
 260     //
 261     // Some special cases:
 262     //
 263     BOOST_CHECK_EQUAL(boost::math::cyl_bessel_j_prime(T(1), T(0)), T(0.5));
 264     BOOST_CHECK_EQUAL(boost::math::cyl_bessel_j_prime(T(-1), T(0)), T(-0.5));
 265     BOOST_CHECK_EQUAL(boost::math::cyl_bessel_j_prime(T(2), T(0)), T(0));
 266
 267     //
 268     // Special cases that are errors:
 269     //
 270     BOOST_MATH_CHECK_THROW(boost::math::sph_bessel_prime(1, T(0)), std::domain_error);
 271     BOOST_MATH_CHECK_THROW(boost::math::sph_bessel_prime(100000, T(0)), std::domain_error);
 272     BOOST_MATH_CHECK_THROW(boost::math::cyl_bessel_j_prime(T(-2.5), T(0)), std::domain_error);
 273     BOOST_MATH_CHECK_THROW(boost::math::cyl_bessel_j_prime(T(-2.5), T(-2)), std::domain_error);
 274     BOOST_MATH_CHECK_THROW(boost::math::cyl_bessel_j_prime(T(2.5), T(-2)), std::domain_error);
 275 }
 276