libs/math/tools/bessel_derivative_data.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. Results of derivatives
   8 //  are generated by the relations between the derivatives
   9 //  and Bessel functions, which actual implementation
  10 //  doesn't use. Results are printed to ~ 50 digits.
  11 //
  12 #include <fstream>
  13
  14 #include <boost/multiprecision/mpfr.hpp>
  15 #include <boost/math/tools/test_data.hpp>
  16 #include <boost/test/included/prg_exec_monitor.hpp>
  17
  18 #include <boost/math/special_functions/bessel.hpp>
  19
  20 using namespace boost::math::tools;
  21 using namespace boost::math;
  22 using namespace std;
  23 using namespace boost::multiprecision;
  24
  25 template <class T>
  26 T bessel_j_derivative_bare(T v, T x)
  27 {
  28    return (v / x) * boost::math::cyl_bessel_j(v, x) - boost::math::cyl_bessel_j(v+1, x);
  29 }
  30
  31 template <class T>
  32 T bessel_y_derivative_bare(T v, T x)
  33 {
  34    return (v / x) * boost::math::cyl_neumann(v, x) - boost::math::cyl_neumann(v+1, x);
  35 }
  36
  37 template <class T>
  38 T bessel_i_derivative_bare(T v, T x)
  39 {
  40    return (v / x) * boost::math::cyl_bessel_i(v, x) + boost::math::cyl_bessel_i(v+1, x);
  41 }
  42
  43 template <class T>
  44 T bessel_k_derivative_bare(T v, T x)
  45 {
  46    return (v / x) * boost::math::cyl_bessel_k(v, x) - boost::math::cyl_bessel_k(v+1, x);
  47 }
  48
  49 template <class T>
  50 T sph_bessel_j_derivative_bare(T v, T x)
  51 {
  52    if((v < 0) || (floor(v) != v))
  53       throw std::domain_error("");
  54    if(v == 0)
  55       return -boost::math::sph_bessel(1, x);
  56    return boost::math::sph_bessel(itrunc(v-1), x) - ((v + 1) / x) * boost::math::sph_bessel(itrunc(v), x);
  57 }
  58
  59 template <class T>
  60 T sph_bessel_y_derivative_bare(T v, T x)
  61 {
  62    if((v < 0) || (floor(v) != v))
  63       throw std::domain_error("");
  64    if(v == 0)
  65       return -boost::math::sph_neumann(1, x);
  66    return boost::math::sph_neumann(itrunc(v-1), x) - ((v + 1) / x) * boost::math::sph_neumann(itrunc(v), x);
  67 }
  68
  69 enum
  70 {
  71    func_J = 0,
  72    func_Y,
  73    func_I,
  74    func_K,
  75    func_j,
  76    func_y
  77 };
  78
  79 int cpp_main(int argc, char*argv [])
  80 {
  81    typedef number<mpfr_float_backend<200> > bignum;
  82
  83    parameter_info<bignum> arg1, arg2;
  84    test_data<bignum> data;
  85
  86    int functype = 0;
  87    std::string letter = "J";
  88
  89    if(argc == 2)
  90    {
  91       if(std::strcmp(argv[1], "--Y") == 0)
  92       {
  93          functype = func_Y;
  94          letter = "Y";
  95       }
  96       else if(std::strcmp(argv[1], "--I") == 0)
  97       {
  98          functype = func_I;
  99          letter = "I";
 100       }
 101       else if(std::strcmp(argv[1], "--K") == 0)
 102       {
 103          functype = func_K;
 104          letter = "K";
 105       }
 106       else if(std::strcmp(argv[1], "--j") == 0)
 107       {
 108          functype = func_j;
 109          letter = "j";
 110       }
 111       else if(std::strcmp(argv[1], "--y") == 0)
 112       {
 113          functype = func_y;
 114          letter = "y";
 115       }
 116       else
 117          BOOST_ASSERT(0);
 118    }
 119
 120    bool cont;
 121    std::string line;
 122
 123    std::cout << "Welcome.\n"
 124       "This program will generate spot tests for the Bessel " << letter << " function derivative\n\n";
 125    do{
 126       if(0 == get_user_parameter_info(arg1, "a"))
 127          return 1;
 128       if(0 == get_user_parameter_info(arg2, "b"))
 129          return 1;
 130
 131       bignum (*fp)(bignum, bignum) = 0;
 132       if(functype == func_J)
 133          fp = bessel_j_derivative_bare;
 134       else if(functype == func_I)
 135          fp = bessel_i_derivative_bare;
 136       else if(functype == func_K)
 137          fp = bessel_k_derivative_bare;
 138       else if(functype == func_Y)
 139          fp = bessel_y_derivative_bare;
 140       else if(functype == func_j)
 141          fp = sph_bessel_j_derivative_bare;
 142       else if(functype == func_y)
 143          fp = sph_bessel_y_derivative_bare;
 144       else
 145          BOOST_ASSERT(0);
 146
 147       data.insert(fp, arg2, arg1);
 148
 149       std::cout << "Any more data [y/n]?";
 150       std::getline(std::cin, line);
 151       boost::algorithm::trim(line);
 152       cont = (line == "y");
 153    }while(cont);
 154
 155    std::cout << "Enter name of test data file [default=bessel_j_derivative_data.ipp]";
 156    std::getline(std::cin, line);
 157    boost::algorithm::trim(line);
 158    if(line == "")
 159       line = "bessel_j_derivative_data.ipp";
 160    std::ofstream ofs(line.c_str());
 161    line.erase(line.find('.'));
 162    ofs << std::scientific << std::setprecision(50);
 163    write_code(ofs, data, line.c_str());
 164
 165    return 0;
 166 }