libs/math/tools/bessel_derivative_data_from_bessel_ipps.cpp

   1 //  Copyright (c) 2014 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 //  Computes test data for the derivatives of the
   7 //  various bessel functions. v and x parameters are taken
   8 //  from bessel_*_data.ipp files. Results of derivatives
   9 //  are generated by the relations between the derivatives
  10 //  and Bessel functions, which actual implementation
  11 //  doesn't use. Results are printed to ~ 50 digits.
  12 //
  13 #include <fstream>
  14 #include <utility>
  15 #include <map>
  16 #include <iterator>
  17 #include <algorithm>
  18
  19 #include <boost/multiprecision/mpfr.hpp>
  20
  21 #include <boost/math/special_functions/bessel.hpp>
  22
  23 template <class T>
  24 T bessel_j_derivative_bare(T v, T x)
  25 {
  26    return (v / x) * boost::math::cyl_bessel_j(v, x) - boost::math::cyl_bessel_j(v+1, x);
  27 }
  28
  29 template <class T>
  30 T bessel_y_derivative_bare(T v, T x)
  31 {
  32    return (v / x) * boost::math::cyl_neumann(v, x) - boost::math::cyl_neumann(v+1, x);
  33 }
  34
  35 template <class T>
  36 T bessel_i_derivative_bare(T v, T x)
  37 {
  38    return (v / x) * boost::math::cyl_bessel_i(v, x) + boost::math::cyl_bessel_i(v+1, x);
  39 }
  40
  41 template <class T>
  42 T bessel_k_derivative_bare(T v, T x)
  43 {
  44    return (v / x) * boost::math::cyl_bessel_k(v, x) - boost::math::cyl_bessel_k(v+1, x);
  45 }
  46
  47 template <class T>
  48 T sph_bessel_j_derivative_bare(T v, T x)
  49 {
  50    if((v < 0) || (floor(v) != v))
  51       throw std::domain_error("");
  52    if(v == 0)
  53       return -boost::math::sph_bessel(1, x);
  54    return boost::math::sph_bessel(itrunc(v-1), x) - ((v + 1) / x) * boost::math::sph_bessel(itrunc(v), x);
  55 }
  56
  57 template <class T>
  58 T sph_bessel_y_derivative_bare(T v, T x)
  59 {
  60    if((v < 0) || (floor(v) != v))
  61       throw std::domain_error("");
  62    if(v == 0)
  63       return -boost::math::sph_neumann(1, x);
  64    return boost::math::sph_neumann(itrunc(v-1), x) - ((v + 1) / x) * boost::math::sph_neumann(itrunc(v), x);
  65 }
  66
  67 using FloatType = boost::multiprecision::number<boost::multiprecision::mpfr_float_backend<200u> >;
  68
  69 enum class BesselFamily: char
  70 {
  71    J = 0,
  72    Y,
  73    I,
  74    K,
  75    j,
  76    y
  77 };
  78
  79 namespace
  80 {
  81
  82 const unsigned kSignificand = 50u;
  83
  84 const std::map<BesselFamily, std::vector<std::string> > kSourceFiles = {
  85    {BesselFamily::J, {"bessel_j_data.ipp", "bessel_j_int_data.ipp", "bessel_j_large_data.ipp"}},
  86    {BesselFamily::Y, {"bessel_y01_data.ipp", "bessel_yn_data.ipp", "bessel_yv_data.ipp"}},
  87    {BesselFamily::I, {"bessel_i_data.ipp", "bessel_i_int_data.ipp"}},
  88    {BesselFamily::K, {"bessel_k_data.ipp", "bessel_k_int_data.ipp"}},
  89    {BesselFamily::j, {"sph_bessel_data.ipp"}},
  90    {BesselFamily::y, {"sph_neumann_data.ipp"}}
  91 };
  92
  93 FloatType (*fp)(FloatType, FloatType) = ::bessel_j_derivative_bare;
  94
  95 std::string parseValue(std::string::iterator& iter)
  96 {
  97    using std::isdigit;
  98
  99    auto value = std::string{};
 100
 101    while (!isdigit(*iter) && *iter != '-')
 102       ++iter;
 103    while (isdigit(*iter) || *iter == '.' || *iter == 'e' || *iter == '-' || *iter == '+')
 104    {
 105       value.push_back(*iter);
 106       ++iter;
 107    }
 108    return value;
 109 }
 110
 111 void replaceResultInLine(std::string& line)
 112 {
 113    using std::isdigit;
 114
 115    auto iter = line.begin();
 116
 117    // parse v and x values from line and convert them to FloatType
 118    auto v = FloatType{::parseValue(iter)};
 119    auto x = FloatType{::parseValue(iter)};
 120    auto result = fp(v, x).str(kSignificand);
 121
 122    while (!isdigit(*iter) && *iter != '-')
 123       ++iter;
 124    const auto where_to_write = iter;
 125    while (isdigit(*iter) || *iter == '.' || *iter == 'e' || *iter == '-' || *iter == '+')
 126       line.erase(iter);
 127
 128    line.insert(where_to_write, result.begin(), result.end());
 129 }
 130
 131 void generateResultFile(const std::string& i_file, const std::string& o_file)
 132 {
 133    std::ifstream in{i_file.c_str()};
 134    std::ofstream out{o_file.c_str()};
 135
 136    auto line = std::string{};
 137    while (!in.eof())
 138    {
 139       std::getline(in, line);
 140       if (__builtin_expect(line.find("SC_") != std::string::npos, 1))
 141          ::replaceResultInLine(line);
 142       out << line << std::endl;
 143    }
 144 }
 145
 146 void processFiles(BesselFamily family)
 147 {
 148    const auto& family_files = kSourceFiles.find(family)->second;
 149
 150    std::for_each(std::begin(family_files), std::end(family_files),
 151       [&](const std::string& src){
 152          auto new_file = src;
 153
 154          const auto int_pos = new_file.find("int");
 155          const auto large_pos = new_file.find("large");
 156          const auto data_pos = new_file.find("data");
 157          const auto derivative_pos = (int_pos == std::string::npos ?
 158             (large_pos == std::string::npos ? data_pos : large_pos) :
 159             int_pos);
 160
 161          new_file.insert(derivative_pos, "derivative_");
 162
 163          ::generateResultFile(src, new_file);
 164       });
 165 }
 166
 167 } // namespace
 168
 169 int main(int argc, char*argv [])
 170 {
 171    auto functype = BesselFamily::J;
 172    auto letter = std::string{"J"};
 173
 174    if(argc == 2)
 175    {
 176       if(std::strcmp(argv[1], "--Y") == 0)
 177       {
 178          functype = BesselFamily::Y;
 179          fp = ::bessel_y_derivative_bare;
 180          letter = "Y";
 181       }
 182       else if(std::strcmp(argv[1], "--I") == 0)
 183       {
 184          functype = BesselFamily::I;
 185          fp = ::bessel_i_derivative_bare;
 186          letter = "I";
 187       }
 188       else if(std::strcmp(argv[1], "--K") == 0)
 189       {
 190          functype = BesselFamily::K;
 191          fp = ::bessel_k_derivative_bare;
 192          letter = "K";
 193       }
 194       else if(std::strcmp(argv[1], "--j") == 0)
 195       {
 196          functype = BesselFamily::j;
 197          fp = ::sph_bessel_j_derivative_bare;
 198          letter = "j";
 199       }
 200       else if(std::strcmp(argv[1], "--y") == 0)
 201       {
 202          functype = BesselFamily::y;
 203          fp = ::sph_bessel_y_derivative_bare;
 204          letter = "y";
 205       }
 206       else
 207          BOOST_ASSERT(0);
 208    }
 209
 210    ::processFiles(functype);
 211
 212    return 0;
 213 }