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. 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.
19 #include <boost/multiprecision/mpfr.hpp>
21 #include <boost/math/special_functions/bessel.hpp>
24 T
bessel_j_derivative_bare(T v
, T x
)
26 return (v
/ x
) * boost::math::cyl_bessel_j(v
, x
) - boost::math::cyl_bessel_j(v
+1, x
);
30 T
bessel_y_derivative_bare(T v
, T x
)
32 return (v
/ x
) * boost::math::cyl_neumann(v
, x
) - boost::math::cyl_neumann(v
+1, x
);
36 T
bessel_i_derivative_bare(T v
, T x
)
38 return (v
/ x
) * boost::math::cyl_bessel_i(v
, x
) + boost::math::cyl_bessel_i(v
+1, x
);
42 T
bessel_k_derivative_bare(T v
, T x
)
44 return (v
/ x
) * boost::math::cyl_bessel_k(v
, x
) - boost::math::cyl_bessel_k(v
+1, x
);
48 T
sph_bessel_j_derivative_bare(T v
, T x
)
50 if((v
< 0) || (floor(v
) != v
))
51 throw std::domain_error("");
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
);
58 T
sph_bessel_y_derivative_bare(T v
, T x
)
60 if((v
< 0) || (floor(v
) != v
))
61 throw std::domain_error("");
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
);
67 using FloatType
= boost::multiprecision::number
<boost::multiprecision::mpfr_float_backend
<200u> >;
69 enum class BesselFamily
: char
82 const unsigned kSignificand
= 50u;
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"}}
93 FloatType (*fp
)(FloatType
, FloatType
) = ::bessel_j_derivative_bare
;
95 std::string
parseValue(std::string::iterator
& iter
)
99 auto value
= std::string
{};
101 while (!isdigit(*iter
) && *iter
!= '-')
103 while (isdigit(*iter
) || *iter
== '.' || *iter
== 'e' || *iter
== '-' || *iter
== '+')
105 value
.push_back(*iter
);
111 void replaceResultInLine(std::string
& line
)
115 auto iter
= line
.begin();
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
);
122 while (!isdigit(*iter
) && *iter
!= '-')
124 const auto where_to_write
= iter
;
125 while (isdigit(*iter
) || *iter
== '.' || *iter
== 'e' || *iter
== '-' || *iter
== '+')
128 line
.insert(where_to_write
, result
.begin(), result
.end());
131 void generateResultFile(const std::string
& i_file
, const std::string
& o_file
)
133 std::ifstream in
{i_file
.c_str()};
134 std::ofstream out
{o_file
.c_str()};
136 auto line
= std::string
{};
139 std::getline(in
, line
);
140 if (__builtin_expect(line
.find("SC_") != std::string::npos
, 1))
141 ::replaceResultInLine(line
);
142 out
<< line
<< std::endl
;
146 void processFiles(BesselFamily family
)
148 const auto& family_files
= kSourceFiles
.find(family
)->second
;
150 std::for_each(std::begin(family_files
), std::end(family_files
),
151 [&](const std::string
& src
){
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
) :
161 new_file
.insert(derivative_pos
, "derivative_");
163 ::generateResultFile(src
, new_file
);
169 int main(int argc
, char*argv
[])
171 auto functype
= BesselFamily::J
;
172 auto letter
= std::string
{"J"};
176 if(std::strcmp(argv
[1], "--Y") == 0)
178 functype
= BesselFamily::Y
;
179 fp
= ::bessel_y_derivative_bare
;
182 else if(std::strcmp(argv
[1], "--I") == 0)
184 functype
= BesselFamily::I
;
185 fp
= ::bessel_i_derivative_bare
;
188 else if(std::strcmp(argv
[1], "--K") == 0)
190 functype
= BesselFamily::K
;
191 fp
= ::bessel_k_derivative_bare
;
194 else if(std::strcmp(argv
[1], "--j") == 0)
196 functype
= BesselFamily::j
;
197 fp
= ::sph_bessel_j_derivative_bare
;
200 else if(std::strcmp(argv
[1], "--y") == 0)
202 functype
= BesselFamily::y
;
203 fp
= ::sph_bessel_y_derivative_bare
;
210 ::processFiles(functype
);