test_jn.cpp 10 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168
  1. // Copyright John Maddock 2015.
  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. #ifdef _MSC_VER
  6. # pragma warning (disable : 4224)
  7. #endif
  8. #include <boost/math/special_functions/bessel.hpp>
  9. #include <boost/array.hpp>
  10. #include <boost/lexical_cast.hpp>
  11. #include "../../test/table_type.hpp"
  12. #include "table_helper.hpp"
  13. #include "performance.hpp"
  14. #include <iostream>
  15. typedef double T;
  16. #define SC_(x) static_cast<double>(x)
  17. static const boost::array<boost::array<typename table_type<T>::type, 3>, 8> j0_data = { {
  18. { { SC_(0.0), SC_(0.0), SC_(1.0) } },
  19. { { SC_(0.0), SC_(1.0), SC_(0.7651976865579665514497175261026632209093) } },
  20. { { SC_(0.0), SC_(-2.0), SC_(0.2238907791412356680518274546499486258252) } },
  21. { { SC_(0.0), SC_(4.0), SC_(-0.3971498098638473722865907684516980419756) } },
  22. { { SC_(0.0), SC_(-8.0), SC_(0.1716508071375539060908694078519720010684) } },
  23. { { SC_(0.0), SC_(1e-05), SC_(0.999999999975000000000156249999999565972) } },
  24. { { SC_(0.0), SC_(1e-10), SC_(0.999999999999999999997500000000000000000) } },
  25. { { SC_(0.0), SC_(-1e+01), SC_(-0.2459357644513483351977608624853287538296) } },
  26. } };
  27. static const boost::array<boost::array<T, 3>, 6> j0_tricky = { {
  28. // Big numbers make the accuracy of std::sin the limiting factor:
  29. { { SC_(0.0), SC_(1e+03), SC_(0.02478668615242017456133073111569370878617) } },
  30. { { SC_(0.0), SC_(1e+05), SC_(-0.001719201116235972192570601477073201747532) } },
  31. // test at the roots:
  32. { { SC_(0.0), T(2521642.0) / (1024 * 1024), SC_(1.80208819970046790002973759410972422387259992955354630042138e-7) } },
  33. { { SC_(0.0), T(5788221.0) / (1024 * 1024), SC_(-1.37774249380686777043369399806210229535671843632174587432454e-7) } },
  34. { { SC_(0.0), T(9074091.0) / (1024 * 1024), SC_(1.03553057441100845081018471279571355857520645127532785991335e-7) } },
  35. { { SC_(0.0), T(12364320.0) / (1024 * 1024), SC_(-3.53017140778223781420794006033810387155048392363051866610931e-9) } }
  36. } };
  37. static const boost::array<boost::array<typename table_type<T>::type, 3>, 8> j1_data = { {
  38. { { SC_(1.0), SC_(0.0), SC_(0.0) } },
  39. { { SC_(1.0), SC_(1.0), SC_(0.4400505857449335159596822037189149131274) } },
  40. { { SC_(1.0), SC_(-2.0), SC_(-0.5767248077568733872024482422691370869203) } },
  41. { { SC_(1.0), SC_(4.0), SC_(-6.604332802354913614318542080327502872742e-02) } },
  42. { { SC_(1.0), SC_(-8.0), SC_(-0.2346363468539146243812766515904546115488) } },
  43. { { SC_(1.0), SC_(1e-05), SC_(4.999999999937500000000260416666666124132e-06) } },
  44. { { SC_(1.0), SC_(1e-10), SC_(4.999999999999999999993750000000000000000e-11) } },
  45. { { SC_(1.0), SC_(-1e+01), SC_(-4.347274616886143666974876802585928830627e-02) } },
  46. } };
  47. static const boost::array<boost::array<T, 3>, 5> j1_tricky = { {
  48. // Big numbers make the accuracy of std::sin the limiting factor:
  49. { { SC_(1.0), SC_(1e+03), SC_(4.728311907089523917576071901216916285418e-03) } },
  50. { { SC_(1.0), SC_(1e+05), SC_(1.846757562882567716362123967114215743694e-03) } },
  51. // test zeros:
  52. { { SC_(1.0), T(4017834) / (1024 * 1024), SC_(3.53149033321258645807835062770856949751958513973522222203044e-7) } },
  53. { { SC_(1.0), T(7356375) / (1024 * 1024), SC_(-2.31227973111067286051984021150135526024117175836722748404342e-7) } },
  54. { { SC_(1.0), T(10667654) / (1024 * 1024), SC_(1.24591331097191900488116495350277530373473085499043086981229e-7) } },
  55. } };
  56. static const boost::array<boost::array<typename table_type<T>::type, 3>, 17> jn_data = { {
  57. // This first one is a modified test case from https://svn.boost.org/trac/boost/ticket/2733
  58. { { SC_(-1.0), SC_(1.25), SC_(-0.510623260319880467069474837274910375352924050139633057168856) } },
  59. { { SC_(2.0), SC_(0.0), SC_(0.0) } },
  60. { { SC_(-2.0), SC_(0.0), SC_(0.0) } },
  61. { { SC_(2.0), SC_(1e-02), SC_(1.249989583365885362413250958437642113452e-05) } },
  62. { { SC_(5.0), SC_(10.0), SC_(-0.2340615281867936404436949416457777864635) } },
  63. { { SC_(5.0), SC_(-10.0), SC_(0.2340615281867936404436949416457777864635) } },
  64. { { SC_(-5.0), SC_(1e+06), SC_(7.259643842453285052375779970433848914846e-04) } },
  65. { { SC_(5.0), SC_(1e+06), SC_(-0.000725964384245328505237577997043384891484649290328285235308619) } },
  66. { { SC_(-5.0), SC_(-1.0), SC_(2.497577302112344313750655409880451981584e-04) } },
  67. { { SC_(10.0), SC_(10.0), SC_(0.2074861066333588576972787235187534280327) } },
  68. { { SC_(10.0), SC_(-10.0), SC_(0.2074861066333588576972787235187534280327) } },
  69. { { SC_(10.0), SC_(-5.0), SC_(1.467802647310474131107532232606627020895e-03) } },
  70. { { SC_(-10.0), SC_(1e+06), SC_(-3.310793117604488741264958559035744460210e-04) } },
  71. { { SC_(10.0), SC_(1e+06), SC_(-0.000331079311760448874126495855903574446020957243277028930713243) } },
  72. { { SC_(1e+02), SC_(8e+01), SC_(4.606553064823477354141298259169874909670e-06) } },
  73. { { SC_(1e+03), SC_(1e+05), SC_(1.283178112502480365195139312635384057363e-03) } },
  74. { { SC_(10.0), SC_(1e-100), SC_(2.69114445546737213403880070546737213403880070546737213403880e-1010) } },
  75. } };
  76. int main()
  77. {
  78. #include "bessel_j_int_data.ipp"
  79. add_data(j0_data);
  80. add_data(j0_tricky);
  81. add_data(j1_data);
  82. add_data(j1_tricky);
  83. add_data(jn_data);
  84. add_data(bessel_j_int_data);
  85. unsigned data_total = data.size();
  86. screen_data([](const std::vector<double>& v){ return boost::math::cyl_bessel_j(static_cast<int>(v[0]), v[1]); }, [](const std::vector<double>& v){ return v[2]; });
  87. #if defined(TEST_C99) && !defined(COMPILER_COMPARISON_TABLES)
  88. screen_data([](const std::vector<double>& v){ return ::jn(static_cast<int>(v[0]), v[1]); }, [](const std::vector<double>& v){ return v[2]; });
  89. #endif
  90. #if defined(TEST_LIBSTDCXX) && !defined(COMPILER_COMPARISON_TABLES)
  91. screen_data([](const std::vector<double>& v){ return std::tr1::cyl_bessel_j(static_cast<int>(v[0]), v[1]); }, [](const std::vector<double>& v){ return v[2]; });
  92. #endif
  93. #if defined(TEST_GSL) && !defined(COMPILER_COMPARISON_TABLES)
  94. screen_data([](const std::vector<double>& v){ return gsl_sf_bessel_Jn(static_cast<int>(v[0]), v[1]); }, [](const std::vector<double>& v){ return v[2]; });
  95. #endif
  96. #if defined(TEST_RMATH) && !defined(COMPILER_COMPARISON_TABLES)
  97. screen_data([](const std::vector<double>& v){ return bessel_j(v[1], static_cast<int>(v[0])); }, [](const std::vector<double>& v){ return v[2]; });
  98. #endif
  99. unsigned data_used = data.size();
  100. std::string function = "cyl_bessel_j (integer order)[br](" + boost::lexical_cast<std::string>(data_used) + "/" + boost::lexical_cast<std::string>(data_total) + " tests selected)";
  101. std::string function_short = "cyl_bessel_j (integer order)";
  102. double time;
  103. time = exec_timed_test([](const std::vector<double>& v){ return boost::math::cyl_bessel_j(static_cast<int>(v[0]), v[1]); });
  104. std::cout << time << std::endl;
  105. #if defined(COMPILER_COMPARISON_TABLES)
  106. report_execution_time(time, std::string("Compiler Option Comparison on ") + platform_name(), "boost::math::cyl_bessel_j (integer orders)", get_compiler_options_name());
  107. #else
  108. #if !defined(COMPILER_COMPARISON_TABLES) && (defined(TEST_GSL) || defined(TEST_RMATH) || defined(TEST_C99) || defined(TEST_LIBSTDCXX))
  109. report_execution_time(time, std::string("Library Comparison with ") + std::string(compiler_name()) + std::string(" on ") + platform_name(), function, boost_name());
  110. #endif
  111. report_execution_time(time, std::string("Compiler Comparison on ") + std::string(platform_name()), function_short, compiler_name() + std::string("[br]") + boost_name());
  112. #endif
  113. //
  114. // Boost again, but with promotion to long double turned off:
  115. //
  116. #if !defined(COMPILER_COMPARISON_TABLES)
  117. if(sizeof(long double) != sizeof(double))
  118. {
  119. time = exec_timed_test([](const std::vector<double>& v){ return boost::math::cyl_bessel_j(static_cast<int>(v[0]), v[1], boost::math::policies::make_policy(boost::math::policies::promote_double<false>())); });
  120. std::cout << time << std::endl;
  121. #if !defined(COMPILER_COMPARISON_TABLES) && (defined(TEST_GSL) || defined(TEST_RMATH) || defined(TEST_C99) || defined(TEST_LIBSTDCXX))
  122. report_execution_time(time, std::string("Library Comparison with ") + std::string(compiler_name()) + std::string(" on ") + platform_name(), function, boost_name() + "[br]promote_double<false>");
  123. #endif
  124. report_execution_time(time, std::string("Compiler Comparison on ") + std::string(platform_name()), function_short, compiler_name() + std::string("[br]") + boost_name() + "[br]promote_double<false>");
  125. }
  126. #endif
  127. #if defined(TEST_C99) && !defined(COMPILER_COMPARISON_TABLES)
  128. time = exec_timed_test([](const std::vector<double>& v){ return ::jn(static_cast<int>(v[0]), v[1]); });
  129. std::cout << time << std::endl;
  130. report_execution_time(time, std::string("Library Comparison with ") + std::string(compiler_name()) + std::string(" on ") + platform_name(), function, "math.h");
  131. #endif
  132. #if defined(TEST_LIBSTDCXX) && !defined(COMPILER_COMPARISON_TABLES)
  133. time = exec_timed_test([](const std::vector<double>& v){ return std::tr1::cyl_bessel_j(static_cast<int>(v[0]), v[1]); });
  134. std::cout << time << std::endl;
  135. report_execution_time(time, std::string("Library Comparison with ") + std::string(compiler_name()) + std::string(" on ") + platform_name(), function, "tr1/cmath");
  136. #endif
  137. #if defined(TEST_GSL) && !defined(COMPILER_COMPARISON_TABLES)
  138. time = exec_timed_test([](const std::vector<double>& v){ return gsl_sf_bessel_Jn(static_cast<int>(v[0]), v[1]); });
  139. std::cout << time << std::endl;
  140. report_execution_time(time, std::string("Library Comparison with ") + std::string(compiler_name()) + std::string(" on ") + platform_name(), function, "GSL " GSL_VERSION);
  141. #endif
  142. #if defined(TEST_RMATH) && !defined(COMPILER_COMPARISON_TABLES)
  143. time = exec_timed_test([](const std::vector<double>& v){ return bessel_j(v[1], static_cast<int>(v[0])); });
  144. std::cout << time << std::endl;
  145. report_execution_time(time, std::string("Library Comparison with ") + std::string(compiler_name()) + std::string(" on ") + platform_name(), function, "Rmath " R_VERSION_STRING);
  146. #endif
  147. return 0;
  148. }