]>
git.proxmox.com Git - ceph.git/blob - ceph/src/boost/libs/math/tools/bessel_data.cpp
1 // Copyright (c) 2007 John Maddock
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 various bessel functions using
7 // archived - deliberately naive - version of the code.
8 // We'll rely on the high precision of mp_t to get us out of
9 // trouble and not worry about how long the calculations take.
10 // This provides a reasonably independent set of test data to
11 // compare against newly added asymptotic expansions etc.
16 #include <boost/math/tools/test_data.hpp>
17 #include <boost/math/special_functions/bessel.hpp>
19 using namespace boost::math::tools
;
20 using namespace boost::math
;
21 using namespace boost::math::detail
;
24 // Compute J(v, x) and Y(v, x) simultaneously by Steed's method, see
25 // Barnett et al, Computer Physics Communications, vol 8, 377 (1974)
27 int bessel_jy_bare(T v
, T x
, T
* J
, T
* Y
, int kind
= need_j
|need_y
)
29 // Jv1 = J_(v+1), Yv1 = Y_(v+1), fv = J_(v+1) / J_v
30 // Ju1 = J_(u+1), Yu1 = Y_(u+1), fu = J_(u+1) / J_u
31 T u
, Jv
, Ju
, Yv
, Yv1
, Yu
, Yu1
, fv
, fu
;
32 T W
, p
, q
, gamma
, current
, prev
, next
;
37 using namespace boost::math::tools
;
38 using namespace boost::math::constants
;
43 v
= -v
; // v is non-negative from here
44 kind
= need_j
|need_y
; // need both for reflection formula
46 n
= real_cast
<int>(v
+ 0.5L);
47 u
= v
- n
; // -1/2 <= u < 1/2
51 *J
= *Y
= policies::raise_domain_error
<T
>("",
52 "Real argument x=%1% must be non-negative, complex number result not supported", x
, policies::policy
<>());
57 *J
= *Y
= policies::raise_overflow_error
<T
>(
58 "", 0, policies::policy
<>());
62 // x is positive until reflection
63 W
= T(2) / (x
* pi
<T
>()); // Wronskian
64 if (x
<= 2) // x in (0, 2]
66 if(temme_jy(u
, x
, &Yu
, &Yu1
, policies::policy
<>())) // Temme series
74 for (k
= 1; k
<= n
; k
++) // forward recurrence for Y
76 next
= 2 * (u
+ k
) * current
/ x
- prev
;
82 CF1_jy(v
, x
, &fv
, &s
, policies::policy
<>()); // continued fraction CF1
83 Jv
= W
/ (Yv
* fv
- Yv1
); // Wronskian relation
85 else // x in (2, \infty)
88 CF1_jy(v
, x
, &fv
, &s
, policies::policy
<>());
89 // tiny initial value to prevent overflow
90 T init
= sqrt(tools::min_value
<T
>());
93 for (k
= n
; k
> 0; k
--) // backward recurrence for J
95 next
= 2 * (u
+ k
) * current
/ x
- prev
;
99 T ratio
= (s
* init
) / current
; // scaling ratio
100 // can also call CF1() to get fu, not much difference in precision
102 CF2_jy(u
, x
, &p
, &q
, policies::policy
<>()); // continued fraction CF2
103 T t
= u
/ x
- fu
; // t = J'/J
105 Ju
= sign(current
) * sqrt(W
/ (q
+ gamma
* (p
- t
)));
107 Jv
= Ju
* ratio
; // normalization
110 Yu1
= Yu
* (u
/x
- p
- q
/gamma
);
115 for (k
= 1; k
<= n
; k
++) // forward recurrence for Y
117 next
= 2 * (u
+ k
) * current
/ x
- prev
;
126 T z
= (u
+ n
% 2) * pi
<T
>();
127 *J
= cos(z
) * Jv
- sin(z
) * Yv
; // reflection formula
128 *Y
= sin(z
) * Jv
+ cos(z
) * Yv
;
142 T
cyl_bessel_j_bare(T v
, T x
)
145 bessel_jy_bare(v
, x
, &j
, &y
);
147 std::cout
<< progress
++ << ": J(" << v
<< ", " << x
<< ") = " << j
<< std::endl
;
150 throw std::domain_error("");
156 T
cyl_bessel_i_bare(T v
, T x
)
161 // better have integer v:
164 T r
= cyl_bessel_i_bare(v
, -x
);
165 if(tools::real_cast
<int>(v
) & 1)
170 return policies::raise_domain_error
<T
>(
172 "Got x = %1%, but we need x >= 0", x
, policies::policy
<>());
176 return (v
== 0) ? 1 : 0;
179 boost::math::detail::bessel_ik(v
, x
, &I
, &K
, 0xffff, policies::policy
<>());
181 std::cout
<< progress
++ << ": I(" << v
<< ", " << x
<< ") = " << I
<< std::endl
;
184 throw std::domain_error("");
190 T
cyl_bessel_k_bare(T v
, T x
)
195 return policies::raise_domain_error
<T
>(
197 "Got x = %1%, but we need x > 0", x
, policies::policy
<>());
201 return (v
== 0) ? policies::raise_overflow_error
<T
>("", 0, policies::policy
<>())
202 : policies::raise_domain_error
<T
>(
204 "Got x = %1%, but we need x > 0", x
, policies::policy
<>());
207 bessel_ik(v
, x
, &I
, &K
, 0xFFFF, policies::policy
<>());
209 std::cout
<< progress
++ << ": K(" << v
<< ", " << x
<< ") = " << K
<< std::endl
;
212 throw std::domain_error("");
218 T
cyl_neumann_bare(T v
, T x
)
221 bessel_jy(v
, x
, &j
, &y
, 0xFFFF, policies::policy
<>());
223 std::cout
<< progress
++ << ": Y(" << v
<< ", " << x
<< ") = " << y
<< std::endl
;
226 throw std::domain_error("");
232 T
sph_bessel_j_bare(T v
, T x
)
234 std::cout
<< progress
++ << ": j(" << v
<< ", " << x
<< ") = ";
235 if((v
< 0) || (floor(v
) != v
))
236 throw std::domain_error("");
237 T r
= sqrt(constants::pi
<T
>() / (2 * x
)) * cyl_bessel_j_bare(v
+0.5, x
);
238 std::cout
<< r
<< std::endl
;
243 T
sph_bessel_y_bare(T v
, T x
)
245 std::cout
<< progress
++ << ": y(" << v
<< ", " << x
<< ") = ";
246 if((v
< 0) || (floor(v
) != v
))
247 throw std::domain_error("");
248 T r
= sqrt(constants::pi
<T
>() / (2 * x
)) * cyl_neumann_bare(v
+0.5, x
);
249 std::cout
<< r
<< std::endl
;
263 int main(int argc
, char* argv
[])
265 std::cout
<< std::setprecision(17) << std::scientific
;
266 std::cout
<< sph_bessel_j_bare(0., 0.1185395751953125e4
) << std::endl
;
267 std::cout
<< sph_bessel_j_bare(22., 0.6540834903717041015625) << std::endl
;
269 std::cout
<< std::setprecision(40) << std::scientific
;
271 parameter_info
<mp_t
> arg1
, arg2
;
272 test_data
<mp_t
> data
;
275 std::string letter
= "J";
279 if(std::strcmp(argv
[1], "--Y") == 0)
284 else if(std::strcmp(argv
[1], "--I") == 0)
289 else if(std::strcmp(argv
[1], "--K") == 0)
294 else if(std::strcmp(argv
[1], "--j") == 0)
299 else if(std::strcmp(argv
[1], "--y") == 0)
311 std::cout
<< "Welcome.\n"
312 "This program will generate spot tests for the Bessel " << letter
<< " function\n\n";
314 get_user_parameter_info(arg1
, "v");
315 get_user_parameter_info(arg2
, "x");
316 mp_t (*fp
)(mp_t
, mp_t
);
317 if(functype
== func_J
)
318 fp
= cyl_bessel_j_bare
;
319 else if(functype
== func_I
)
320 fp
= cyl_bessel_i_bare
;
321 else if(functype
== func_K
)
322 fp
= cyl_bessel_k_bare
;
323 else if(functype
== func_Y
)
324 fp
= cyl_neumann_bare
;
325 else if(functype
== func_j
)
326 fp
= sph_bessel_j_bare
;
327 else if(functype
== func_y
)
328 fp
= sph_bessel_y_bare
;
332 data
.insert(fp
, arg1
, arg2
);
334 std::cout
<< "Any more data [y/n]?";
335 std::getline(std::cin
, line
);
336 boost::algorithm::trim(line
);
337 cont
= (line
== "y");
340 std::cout
<< "Enter name of test data file [default=bessel_j_data.ipp]";
341 std::getline(std::cin
, line
);
342 boost::algorithm::trim(line
);
344 line
= "bessel_j_data.ipp";
345 std::ofstream
ofs(line
.c_str());
346 line
.erase(line
.find('.'));
347 ofs
<< std::scientific
<< std::setprecision(40);
348 write_code(ofs
, data
, line
.c_str());