123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213 |
- // Copyright (c) 2014 Anton Bikineev
- // Use, modification and distribution are subject to the
- // Boost Software License, Version 1.0. (See accompanying file
- // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
- //
- // Computes test data for the derivatives of the
- // various bessel functions. v and x parameters are taken
- // from bessel_*_data.ipp files. Results of derivatives
- // are generated by the relations between the derivatives
- // and Bessel functions, which actual implementation
- // doesn't use. Results are printed to ~ 50 digits.
- //
- #include <fstream>
- #include <utility>
- #include <map>
- #include <iterator>
- #include <algorithm>
- #include <boost/multiprecision/mpfr.hpp>
- #include <boost/math/special_functions/bessel.hpp>
- template <class T>
- T bessel_j_derivative_bare(T v, T x)
- {
- return (v / x) * boost::math::cyl_bessel_j(v, x) - boost::math::cyl_bessel_j(v+1, x);
- }
- template <class T>
- T bessel_y_derivative_bare(T v, T x)
- {
- return (v / x) * boost::math::cyl_neumann(v, x) - boost::math::cyl_neumann(v+1, x);
- }
- template <class T>
- T bessel_i_derivative_bare(T v, T x)
- {
- return (v / x) * boost::math::cyl_bessel_i(v, x) + boost::math::cyl_bessel_i(v+1, x);
- }
- template <class T>
- T bessel_k_derivative_bare(T v, T x)
- {
- return (v / x) * boost::math::cyl_bessel_k(v, x) - boost::math::cyl_bessel_k(v+1, x);
- }
- template <class T>
- T sph_bessel_j_derivative_bare(T v, T x)
- {
- if((v < 0) || (floor(v) != v))
- throw std::domain_error("");
- if(v == 0)
- return -boost::math::sph_bessel(1, x);
- return boost::math::sph_bessel(itrunc(v-1), x) - ((v + 1) / x) * boost::math::sph_bessel(itrunc(v), x);
- }
- template <class T>
- T sph_bessel_y_derivative_bare(T v, T x)
- {
- if((v < 0) || (floor(v) != v))
- throw std::domain_error("");
- if(v == 0)
- return -boost::math::sph_neumann(1, x);
- return boost::math::sph_neumann(itrunc(v-1), x) - ((v + 1) / x) * boost::math::sph_neumann(itrunc(v), x);
- }
- using FloatType = boost::multiprecision::number<boost::multiprecision::mpfr_float_backend<200u> >;
- enum class BesselFamily: char
- {
- J = 0,
- Y,
- I,
- K,
- j,
- y
- };
- namespace
- {
- const unsigned kSignificand = 50u;
- const std::map<BesselFamily, std::vector<std::string> > kSourceFiles = {
- {BesselFamily::J, {"bessel_j_data.ipp", "bessel_j_int_data.ipp", "bessel_j_large_data.ipp"}},
- {BesselFamily::Y, {"bessel_y01_data.ipp", "bessel_yn_data.ipp", "bessel_yv_data.ipp"}},
- {BesselFamily::I, {"bessel_i_data.ipp", "bessel_i_int_data.ipp"}},
- {BesselFamily::K, {"bessel_k_data.ipp", "bessel_k_int_data.ipp"}},
- {BesselFamily::j, {"sph_bessel_data.ipp"}},
- {BesselFamily::y, {"sph_neumann_data.ipp"}}
- };
- FloatType (*fp)(FloatType, FloatType) = ::bessel_j_derivative_bare;
- std::string parseValue(std::string::iterator& iter)
- {
- using std::isdigit;
- auto value = std::string{};
- while (!isdigit(*iter) && *iter != '-')
- ++iter;
- while (isdigit(*iter) || *iter == '.' || *iter == 'e' || *iter == '-' || *iter == '+')
- {
- value.push_back(*iter);
- ++iter;
- }
- return value;
- }
- void replaceResultInLine(std::string& line)
- {
- using std::isdigit;
- auto iter = line.begin();
- // parse v and x values from line and convert them to FloatType
- auto v = FloatType{::parseValue(iter)};
- auto x = FloatType{::parseValue(iter)};
- auto result = fp(v, x).str(kSignificand);
- while (!isdigit(*iter) && *iter != '-')
- ++iter;
- const auto where_to_write = iter;
- while (isdigit(*iter) || *iter == '.' || *iter == 'e' || *iter == '-' || *iter == '+')
- line.erase(iter);
- line.insert(where_to_write, result.begin(), result.end());
- }
- void generateResultFile(const std::string& i_file, const std::string& o_file)
- {
- std::ifstream in{i_file.c_str()};
- std::ofstream out{o_file.c_str()};
- auto line = std::string{};
- while (!in.eof())
- {
- std::getline(in, line);
- if (__builtin_expect(line.find("SC_") != std::string::npos, 1))
- ::replaceResultInLine(line);
- out << line << std::endl;
- }
- }
- void processFiles(BesselFamily family)
- {
- const auto& family_files = kSourceFiles.find(family)->second;
- std::for_each(std::begin(family_files), std::end(family_files),
- [&](const std::string& src){
- auto new_file = src;
- const auto int_pos = new_file.find("int");
- const auto large_pos = new_file.find("large");
- const auto data_pos = new_file.find("data");
- const auto derivative_pos = (int_pos == std::string::npos ?
- (large_pos == std::string::npos ? data_pos : large_pos) :
- int_pos);
- new_file.insert(derivative_pos, "derivative_");
- ::generateResultFile(src, new_file);
- });
- }
- } // namespace
- int main(int argc, char*argv [])
- {
- auto functype = BesselFamily::J;
- auto letter = std::string{"J"};
- if(argc == 2)
- {
- if(std::strcmp(argv[1], "--Y") == 0)
- {
- functype = BesselFamily::Y;
- fp = ::bessel_y_derivative_bare;
- letter = "Y";
- }
- else if(std::strcmp(argv[1], "--I") == 0)
- {
- functype = BesselFamily::I;
- fp = ::bessel_i_derivative_bare;
- letter = "I";
- }
- else if(std::strcmp(argv[1], "--K") == 0)
- {
- functype = BesselFamily::K;
- fp = ::bessel_k_derivative_bare;
- letter = "K";
- }
- else if(std::strcmp(argv[1], "--j") == 0)
- {
- functype = BesselFamily::j;
- fp = ::sph_bessel_j_derivative_bare;
- letter = "j";
- }
- else if(std::strcmp(argv[1], "--y") == 0)
- {
- functype = BesselFamily::y;
- fp = ::sph_bessel_y_derivative_bare;
- letter = "y";
- }
- else
- BOOST_ASSERT(0);
- }
- ::processFiles(functype);
- return 0;
- }
|