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)
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.
14 #include <boost/multiprecision/mpfr.hpp>
15 #include <boost/math/tools/test_data.hpp>
16 #include <boost/test/included/prg_exec_monitor.hpp>
18 #include <boost/math/special_functions/bessel.hpp>
20 using namespace boost::math::tools
;
21 using namespace boost::math
;
23 using namespace boost::multiprecision
;
26 T
bessel_j_derivative_bare(T v
, T x
)
28 return (v
/ x
) * boost::math::cyl_bessel_j(v
, x
) - boost::math::cyl_bessel_j(v
+1, x
);
32 T
bessel_y_derivative_bare(T v
, T x
)
34 return (v
/ x
) * boost::math::cyl_neumann(v
, x
) - boost::math::cyl_neumann(v
+1, x
);
38 T
bessel_i_derivative_bare(T v
, T x
)
40 return (v
/ x
) * boost::math::cyl_bessel_i(v
, x
) + boost::math::cyl_bessel_i(v
+1, x
);
44 T
bessel_k_derivative_bare(T v
, T x
)
46 return (v
/ x
) * boost::math::cyl_bessel_k(v
, x
) - boost::math::cyl_bessel_k(v
+1, x
);
50 T
sph_bessel_j_derivative_bare(T v
, T x
)
52 if((v
< 0) || (floor(v
) != v
))
53 throw std::domain_error("");
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
);
60 T
sph_bessel_y_derivative_bare(T v
, T x
)
62 if((v
< 0) || (floor(v
) != v
))
63 throw std::domain_error("");
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
);
79 int cpp_main(int argc
, char*argv
[])
81 typedef number
<mpfr_float_backend
<200> > bignum
;
83 parameter_info
<bignum
> arg1
, arg2
;
84 test_data
<bignum
> data
;
87 std::string letter
= "J";
91 if(std::strcmp(argv
[1], "--Y") == 0)
96 else if(std::strcmp(argv
[1], "--I") == 0)
101 else if(std::strcmp(argv
[1], "--K") == 0)
106 else if(std::strcmp(argv
[1], "--j") == 0)
111 else if(std::strcmp(argv
[1], "--y") == 0)
123 std::cout
<< "Welcome.\n"
124 "This program will generate spot tests for the Bessel " << letter
<< " function derivative\n\n";
126 if(0 == get_user_parameter_info(arg1
, "a"))
128 if(0 == get_user_parameter_info(arg2
, "b"))
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
;
147 data
.insert(fp
, arg2
, arg1
);
149 std::cout
<< "Any more data [y/n]?";
150 std::getline(std::cin
, line
);
151 boost::algorithm::trim(line
);
152 cont
= (line
== "y");
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
);
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());