math/tools/bessel_derivative_data_from_bessel_ipps.cpp

//  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;
}