]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 | |
92f5a8d4 | 117 | BOOST_ASSERT(0); |
7c673cae FG |
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 | |
92f5a8d4 | 145 | BOOST_ASSERT(0); |
7c673cae FG |
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 | } |