1 // Copyright John Maddock 2012
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 #include <pch_light.hpp>
8 #define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
9 #define BOOST_TEST_MAIN
10 #include <boost/test/unit_test.hpp>
11 #include <boost/test/floating_point_comparison.hpp>
12 #include <boost/math/special_functions/math_fwd.hpp>
13 #include <boost/math/concepts/real_concept.hpp>
14 #include <boost/array.hpp>
15 #include <boost/math/special_functions/hankel.hpp>
20 # pragma warning(disable : 4756 4127) // overflow in constant arithmetic
21 // Constants are too big for float case, but this doesn't matter for test.
28 // This file tests the Hankel H1 and H2 functions.
29 // These are basically just a few spot tests, since the underlying implementation
30 // is the same as for the Bessel functions.
33 void check_close(const std::complex<T
>& a
, const std::complex<T
>& b
)
35 T tol
= boost::math::tools::epsilon
<T
>() * 3000;
36 BOOST_CHECK_CLOSE_FRACTION(a
.real(), b
.real(), tol
);
37 BOOST_CHECK_CLOSE_FRACTION(a
.imag(), b
.imag(), tol
);
41 void test_hankel(T
, const char* name
)
43 std::cout
<< "Testing type " << name
<< std::endl
;
45 static const boost::array
<boost::array
<std::complex<T
>, 4>, 16> data
=
47 // Values are v, z, J, and Y.
48 // H1 and H2 are calculated from functions.wolfram.com.
49 {{ 0, 1, static_cast<T
>(0.765197686557966551449717526102663220909274289755325241861548L), static_cast<T
>(0.0882569642156769579829267660235151628278175230906755467110438L) }},
50 {{ 20, 15.5, static_cast<T
>(0.0114685669049540880132573889173945340441319848974792913849131L), static_cast<T
>(-2.23357703561803241011539677726233479433825678625142168545338L) }},
51 {{ 202, 65, static_cast<T
>(4.03123907894335908039069177339845754619082182624766366629474e-77L), static_cast<T
>(-4.12853948338543234456930697819285129000267780751856135651604e73L
) }},
52 {{ 1.25, 2.25, static_cast<T
>(0.548918751190427983177518806565451279306141559548962787802891L), static_cast<T
>(-0.125900744882628421758627761371862848797979705547465002154794L) }},
53 {{ -20, 15.5, static_cast<T
>(0.0114685669049540880132573889173945340441319848974792913849131L), static_cast<T
>(-2.23357703561803241011539677726233479433825678625142168545338L) }},
54 {{ -20, -15.5, static_cast<T
>(0.0114685669049540880132573889173945340441319848974792913849131L), std::complex<T
>(static_cast<T
>(-2.23357703561803241011539677726233479433825678625142168545338L), static_cast<T
>(0.02293713380990817602651477783478906808826396979495858276983L))}},
55 {{ 1.25, -1.5, std::complex<T
>(static_cast<T
>(-0.335713500965919366139805990226845134897000581426417882156072L), static_cast<T
>(-0.335713500965919366139805990226845134897000581426417882156072L)), std::complex<T
>(static_cast<T
>(0.392747687664481978664420363126684555843062241632017696636768L), static_cast<T
>(-1.064174689596320710944032343580374825637063404484853460948911L)) }},
56 {{ -1.25, -1.5, std::complex<T
>(static_cast<T
>(0.515099846311769525023685454389374962206960751452481078650112L), static_cast<T
>(-0.515099846311769525023685454389374962206960751452481078650112L)), std::complex<T
>(static_cast<T
>(-0.040329260174013212604062774398331485097569895404765001721544L), static_cast<T
>(0.989870432449525837443308134380418439316351607500197155578680L)) }},
57 {{ 4, 4, static_cast<T
>(0.281129064961360106322277160229942806897088617059328870629222L), static_cast<T
>(-0.488936768533842510615657398339913206218740182079627974737267L) }},
58 {{ 4, -4, static_cast<T
>(0.281129064961360106322277160229942806897088617059328870629222L), std::complex<T
>(static_cast<T
>(-0.488936768533842510615657398339913206218740182079627974737267L), static_cast<T
>(0.562258129922720212644554320459885613794177234118657741258443L)) }},
59 {{ -4, 4, static_cast<T
>(0.281129064961360106322277160229942806897088617059328870629222L), static_cast<T
>(-0.488936768533842510615657398339913206218740182079627974737267L) }},
60 {{ -4, -4, static_cast<T
>(0.281129064961360106322277160229942806897088617059328870629222), std::complex<T
>(static_cast<T
>(-0.488936768533842510615657398339913206218740182079627974737267), static_cast<T
>(0.562258129922720212644554320459885613794177234118657741258443L)) }},
61 {{ 3, 3, static_cast<T
>(0.309062722255251643618260194946833149429135935993056794354475L), static_cast<T
>(-0.538541616105031618004703905338594463807957863604859251481262L) }},
62 {{ 3, -3, static_cast<T
>(-0.309062722255251643618260194946833149429135935993056794354475L), std::complex<T
>(static_cast<T
>(0.538541616105031618004703905338594463807957863604859251481262L), static_cast<T
>(-0.618125444510503287236520389893666298858271871986113588708949L)) }},
63 {{ -3, 3, static_cast<T
>(-0.309062722255251643618260194946833149429135935993056794354475L), static_cast<T
>(0.538541616105031618004703905338594463807957863604859251481262L) }},
64 {{ -3, -3, static_cast<T
>(0.309062722255251643618260194946833149429135935993056794354475L), std::complex<T
>(static_cast<T
>(-0.538541616105031618004703905338594463807957863604859251481262L), static_cast<T
>(0.618125444510503287236520389893666298858271871986113588708949L)) }},
67 std::complex<T
> im(0, 1);
68 for(unsigned i
= 0; i
< data
.size(); ++i
)
70 if((i
!= 2) || (std::numeric_limits
<T
>::max_exponent10
> 80))
72 check_close(boost::math::cyl_hankel_1(data
[i
][0].real(), data
[i
][1].real()), data
[i
][2] + im
* data
[i
][3]);
73 check_close(boost::math::cyl_hankel_2(data
[i
][0].real(), data
[i
][1].real()), data
[i
][2] - im
* data
[i
][3]);
76 boost::math::cyl_hankel_1(data
[i
][0].real() + 0.5f
, data
[i
][1].real()) * boost::math::constants::root_half_pi
<T
>() / sqrt(data
[i
][1]),
77 boost::math::sph_hankel_1(data
[i
][0].real(), data
[i
][1].real()));
79 boost::math::cyl_hankel_2(data
[i
][0].real() + 0.5f
, data
[i
][1].real()) * boost::math::constants::root_half_pi
<T
>() / sqrt(data
[i
][1]),
80 boost::math::sph_hankel_2(data
[i
][0].real(), data
[i
][1].real()));
86 // Instantiate a few instances to check our error handling code can cope with std::complex:
88 typedef boost::math::policies::policy
<
89 boost::math::policies::overflow_error
<boost::math::policies::throw_on_error
>,
90 boost::math::policies::denorm_error
<boost::math::policies::throw_on_error
>,
91 boost::math::policies::underflow_error
<boost::math::policies::throw_on_error
>,
92 boost::math::policies::domain_error
<boost::math::policies::throw_on_error
>,
93 boost::math::policies::pole_error
<boost::math::policies::throw_on_error
>,
94 boost::math::policies::rounding_error
<boost::math::policies::throw_on_error
>,
95 boost::math::policies::evaluation_error
<boost::math::policies::throw_on_error
>,
96 boost::math::policies::indeterminate_result_error
<boost::math::policies::throw_on_error
> > pol1
;
98 template std::complex<double> boost::math::cyl_hankel_1
<double, double, pol1
>(double, double, const pol1
&);
100 typedef boost::math::policies::policy
<
101 boost::math::policies::overflow_error
<boost::math::policies::errno_on_error
>,
102 boost::math::policies::denorm_error
<boost::math::policies::errno_on_error
>,
103 boost::math::policies::underflow_error
<boost::math::policies::errno_on_error
>,
104 boost::math::policies::domain_error
<boost::math::policies::errno_on_error
>,
105 boost::math::policies::pole_error
<boost::math::policies::errno_on_error
>,
106 boost::math::policies::rounding_error
<boost::math::policies::errno_on_error
>,
107 boost::math::policies::evaluation_error
<boost::math::policies::errno_on_error
>,
108 boost::math::policies::indeterminate_result_error
<boost::math::policies::errno_on_error
> > pol2
;
110 template std::complex<double> boost::math::cyl_hankel_1
<double, double, pol2
>(double, double, const pol2
&);
112 typedef boost::math::policies::policy
<
113 boost::math::policies::overflow_error
<boost::math::policies::ignore_error
>,
114 boost::math::policies::denorm_error
<boost::math::policies::ignore_error
>,
115 boost::math::policies::underflow_error
<boost::math::policies::ignore_error
>,
116 boost::math::policies::domain_error
<boost::math::policies::ignore_error
>,
117 boost::math::policies::pole_error
<boost::math::policies::ignore_error
>,
118 boost::math::policies::rounding_error
<boost::math::policies::ignore_error
>,
119 boost::math::policies::evaluation_error
<boost::math::policies::ignore_error
>,
120 boost::math::policies::indeterminate_result_error
<boost::math::policies::ignore_error
> > pol3
;
122 template std::complex<double> boost::math::cyl_hankel_1
<double, double, pol3
>(double, double, const pol3
&);
125 BOOST_AUTO_TEST_CASE( test_main
)
128 gsl_set_error_handler_off();
130 BOOST_MATH_CONTROL_FP
;
132 #ifndef BOOST_MATH_BUGGY_LARGE_FLOAT_CONSTANTS
133 test_hankel(0.1F
, "float");
135 test_hankel(0.1, "double");
136 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
137 test_hankel(0.1L, "long double");
139 std::cout
<< "<note>The long double tests have been disabled on this platform "
140 "either because the long double overloads of the usual math functions are "
141 "not available at all, or because they are too inaccurate for these tests "
142 "to pass.</note>" << std::endl
;