]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // Copyright John Maddock 2006. |
2 | // Copyright Paul A. Bristow 2007, 2009 | |
3 | // Use, modification and distribution are subject to the | |
4 | // Boost Software License, Version 1.0. (See accompanying file | |
5 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
6 | ||
7 | #ifdef _MSC_VER | |
8 | # pragma warning (disable : 4756) // overflow in constant arithmetic | |
9 | #endif | |
10 | ||
11 | #include <boost/math/concepts/real_concept.hpp> | |
12 | #define BOOST_TEST_MAIN | |
13 | #include <boost/test/unit_test.hpp> | |
14 | #include <boost/test/floating_point_comparison.hpp> | |
15 | #include <boost/math/special_functions/math_fwd.hpp> | |
b32b8144 | 16 | #include <boost/math/special_functions/legendre.hpp> |
7c673cae | 17 | #include <boost/math/constants/constants.hpp> |
b32b8144 | 18 | #include <boost/multiprecision/cpp_bin_float.hpp> |
7c673cae FG |
19 | #include <boost/array.hpp> |
20 | #include "functor.hpp" | |
21 | ||
22 | #include "handle_test_result.hpp" | |
23 | #include "table_type.hpp" | |
24 | ||
25 | #ifndef SC_ | |
26 | #define SC_(x) static_cast<typename table_type<T>::type>(BOOST_JOIN(x, L)) | |
27 | #endif | |
28 | ||
29 | template <class Real, class T> | |
30 | void do_test_legendre_p(const T& data, const char* type_name, const char* test_name) | |
31 | { | |
32 | typedef Real value_type; | |
33 | ||
34 | typedef value_type (*pg)(int, value_type); | |
35 | pg funcp; | |
36 | ||
37 | #if !(defined(ERROR_REPORTING_MODE) && !defined(LEGENDRE_P_FUNCTION_TO_TEST)) | |
38 | #ifdef LEGENDRE_P_FUNCTION_TO_TEST | |
39 | funcp = LEGENDRE_P_FUNCTION_TO_TEST; | |
40 | #elif defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS) | |
41 | funcp = boost::math::legendre_p<value_type>; | |
42 | #else | |
43 | funcp = boost::math::legendre_p; | |
44 | #endif | |
45 | ||
46 | boost::math::tools::test_result<value_type> result; | |
47 | ||
48 | std::cout << "Testing " << test_name << " with type " << type_name | |
49 | << "\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n"; | |
50 | ||
51 | // | |
52 | // test legendre_p against data: | |
53 | // | |
54 | result = boost::math::tools::test_hetero<Real>( | |
55 | data, | |
56 | bind_func_int1<Real>(funcp, 0, 1), | |
57 | extract_result<Real>(2)); | |
58 | handle_test_result(result, data[result.worst()], result.worst(), type_name, "legendre_p", test_name); | |
59 | #endif | |
60 | ||
61 | typedef value_type (*pg2)(unsigned, value_type); | |
62 | #if !(defined(ERROR_REPORTING_MODE) && !defined(LEGENDRE_Q_FUNCTION_TO_TEST)) | |
63 | #ifdef LEGENDRE_Q_FUNCTION_TO_TEST | |
64 | pg2 funcp2 = LEGENDRE_Q_FUNCTION_TO_TEST; | |
65 | #elif defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS) | |
66 | pg2 funcp2 = boost::math::legendre_q<value_type>; | |
67 | #else | |
68 | pg2 funcp2 = boost::math::legendre_q; | |
69 | #endif | |
70 | ||
71 | // | |
72 | // test legendre_q against data: | |
73 | // | |
74 | result = boost::math::tools::test_hetero<Real>( | |
75 | data, | |
76 | bind_func_int1<Real>(funcp2, 0, 1), | |
77 | extract_result<Real>(3)); | |
78 | handle_test_result(result, data[result.worst()], result.worst(), type_name, "legendre_q", test_name); | |
79 | ||
80 | std::cout << std::endl; | |
81 | #endif | |
82 | } | |
83 | ||
84 | template <class Real, class T> | |
85 | void do_test_assoc_legendre_p(const T& data, const char* type_name, const char* test_name) | |
86 | { | |
87 | #if !(defined(ERROR_REPORTING_MODE) && !defined(LEGENDRE_PA_FUNCTION_TO_TEST)) | |
88 | typedef Real value_type; | |
89 | ||
90 | typedef value_type (*pg)(int, int, value_type); | |
91 | #ifdef LEGENDRE_PA_FUNCTION_TO_TEST | |
92 | pg funcp = LEGENDRE_PA_FUNCTION_TO_TEST; | |
93 | #elif defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS) | |
94 | pg funcp = boost::math::legendre_p<value_type>; | |
95 | #else | |
96 | pg funcp = boost::math::legendre_p; | |
97 | #endif | |
98 | ||
99 | boost::math::tools::test_result<value_type> result; | |
100 | ||
101 | std::cout << "Testing " << test_name << " with type " << type_name | |
102 | << "\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n"; | |
103 | ||
104 | // | |
105 | // test legendre_p against data: | |
106 | // | |
107 | result = boost::math::tools::test_hetero<Real>( | |
108 | data, | |
109 | bind_func_int2<Real>(funcp, 0, 1, 2), | |
110 | extract_result<Real>(3)); | |
111 | handle_test_result(result, data[result.worst()], result.worst(), type_name, "legendre_p (associated)", test_name); | |
112 | std::cout << std::endl; | |
113 | #endif | |
114 | } | |
115 | ||
116 | template <class T> | |
117 | void test_legendre_p(T, const char* name) | |
118 | { | |
119 | // | |
120 | // The actual test data is rather verbose, so it's in a separate file | |
121 | // | |
122 | // The contents are as follows, each row of data contains | |
123 | // three items, input value a, input value b and erf(a, b): | |
124 | // | |
125 | # include "legendre_p.ipp" | |
126 | ||
127 | do_test_legendre_p<T>(legendre_p, name, "Legendre Polynomials: Small Values"); | |
128 | ||
129 | # include "legendre_p_large.ipp" | |
130 | ||
131 | do_test_legendre_p<T>(legendre_p_large, name, "Legendre Polynomials: Large Values"); | |
132 | ||
133 | # include "assoc_legendre_p.ipp" | |
134 | ||
135 | do_test_assoc_legendre_p<T>(assoc_legendre_p, name, "Associated Legendre Polynomials: Small Values"); | |
136 | ||
137 | } | |
138 | ||
139 | template <class T> | |
140 | void test_spots(T, const char* t) | |
141 | { | |
142 | std::cout << "Testing basic sanity checks for type " << t << std::endl; | |
143 | // | |
144 | // basic sanity checks, tolerance is 100 epsilon: | |
145 | // | |
146 | T tolerance = boost::math::tools::epsilon<T>() * 100; | |
147 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(1, static_cast<T>(0.5L)), static_cast<T>(0.5L), tolerance); | |
148 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-1, static_cast<T>(0.5L)), static_cast<T>(1L), tolerance); | |
149 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(4, static_cast<T>(0.5L)), static_cast<T>(-0.2890625000000000000000000000000000000000L), tolerance); | |
150 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-4, static_cast<T>(0.5L)), static_cast<T>(-0.4375000000000000000000000000000000000000L), tolerance); | |
151 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(7, static_cast<T>(0.5L)), static_cast<T>(0.2231445312500000000000000000000000000000L), tolerance); | |
152 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-7, static_cast<T>(0.5L)), static_cast<T>(0.3232421875000000000000000000000000000000L), tolerance); | |
153 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(40, static_cast<T>(0.5L)), static_cast<T>(-0.09542943523261546936538467572384923220258L), tolerance); | |
154 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-40, static_cast<T>(0.5L)), static_cast<T>(-0.1316993126940266257030910566308990611306L), tolerance); | |
155 | ||
156 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(4, 2, static_cast<T>(0.5L)), static_cast<T>(4.218750000000000000000000000000000000000L), tolerance); | |
157 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-4, 2, static_cast<T>(0.5L)), static_cast<T>(5.625000000000000000000000000000000000000L), tolerance); | |
158 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(7, 5, static_cast<T>(0.5L)), static_cast<T>(-5696.789530152175143607977274672800795328L), tolerance); | |
159 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-7, 4, static_cast<T>(0.5L)), static_cast<T>(465.1171875000000000000000000000000000000L), tolerance); | |
160 | if(std::numeric_limits<T>::max_exponent > std::numeric_limits<float>::max_exponent) | |
161 | { | |
162 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(40, 30, static_cast<T>(0.5L)), static_cast<T>(-7.855722083232252643913331343916012143461e45L), tolerance); | |
163 | } | |
164 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-40, 20, static_cast<T>(0.5L)), static_cast<T>(4.966634149702370788037088925152355134665e30L), tolerance); | |
165 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(4, 2, static_cast<T>(-0.5L)), static_cast<T>(4.218750000000000000000000000000000000000L), tolerance); | |
166 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-4, 2, static_cast<T>(-0.5L)), static_cast<T>(-5.625000000000000000000000000000000000000L), tolerance); | |
167 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(7, 5, static_cast<T>(-0.5L)), static_cast<T>(-5696.789530152175143607977274672800795328L), tolerance); | |
168 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-7, 4, static_cast<T>(-0.5L)), static_cast<T>(465.1171875000000000000000000000000000000L), tolerance); | |
169 | if(std::numeric_limits<T>::max_exponent > std::numeric_limits<float>::max_exponent) | |
170 | { | |
171 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(40, 30, static_cast<T>(-0.5L)), static_cast<T>(-7.855722083232252643913331343916012143461e45L), tolerance); | |
172 | } | |
173 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-40, 20, static_cast<T>(-0.5L)), static_cast<T>(-4.966634149702370788037088925152355134665e30L), tolerance); | |
174 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(4, -2, static_cast<T>(0.5L)), static_cast<T>(0.01171875000000000000000000000000000000000L), tolerance); | |
175 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-4, -2, static_cast<T>(0.5L)), static_cast<T>(0.04687500000000000000000000000000000000000L), tolerance); | |
176 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(7, -5, static_cast<T>(0.5L)), static_cast<T>(0.00002378609812640364935569308025139290054701L), tolerance); | |
177 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-7, -4, static_cast<T>(0.5L)), static_cast<T>(0.0002563476562500000000000000000000000000000L), tolerance); | |
178 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(40, -30, static_cast<T>(0.5L)), static_cast<T>(-2.379819988646847616996471299410611801239e-48L), tolerance); | |
179 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-40, -20, static_cast<T>(0.5L)), static_cast<T>(4.356454600748202401657099008867502679122e-33L), tolerance); | |
180 | ||
181 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_q(1, static_cast<T>(0.5L)), static_cast<T>(-0.7253469278329725771511886907693685738381L), tolerance); | |
182 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_q(4, static_cast<T>(0.5L)), static_cast<T>(0.4401745259867706044988642951843745400835L), tolerance); | |
183 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_q(7, static_cast<T>(0.5L)), static_cast<T>(-0.3439152932669753451878700644212067616780L), tolerance); | |
184 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_q(40, static_cast<T>(0.5L)), static_cast<T>(0.1493671665503550095010454949479907886011L), tolerance); | |
185 | } | |
186 | ||
b32b8144 FG |
187 | template <class T> |
188 | void test_legendre_p_prime() | |
189 | { | |
190 | T tolerance = 100*boost::math::tools::epsilon<T>(); | |
191 | T x = -1; | |
192 | while (x <= 1) | |
193 | { | |
194 | // P_0'(x) = 0 | |
195 | BOOST_CHECK_SMALL(::boost::math::legendre_p_prime<T>(0, x), tolerance); | |
196 | // Reflection formula for P_{-1}(x) = P_{0}(x): | |
197 | BOOST_CHECK_SMALL(::boost::math::legendre_p_prime<T>(-1, x), tolerance); | |
198 | ||
199 | // P_1(x) = x, so P_1'(x) = 1: | |
200 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(1, x), static_cast<T>(1), tolerance); | |
201 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-2, x), static_cast<T>(1), tolerance); | |
202 | ||
203 | // P_2(x) = 3x^2/2 + k => P_2'(x) = 3x | |
204 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(2, x), 3*x, tolerance); | |
205 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-3, x), 3*x, tolerance); | |
206 | ||
207 | // P_3(x) = (5x^3 - 3x)/2 => P_3'(x) = (15x^2 - 3)/2: | |
208 | T xsq = x*x; | |
209 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(3, x), (15*xsq - 3)/2, tolerance); | |
210 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-4, x), (15*xsq -3)/2, tolerance); | |
211 | ||
212 | // P_4(x) = (35x^4 - 30x^2 +3)/8 => P_4'(x) = (5x/2)*(7x^2 - 3) | |
213 | T expected = 5*x*(7*xsq - 3)/2; | |
214 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(4, x), expected, tolerance); | |
215 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-5, x), expected, tolerance); | |
216 | ||
217 | // P_5(x) = (63x^5 - 70x^3 + 15x)/8 => P_5'(x) = (315*x^4 - 210*x^2 + 15)/8 | |
218 | T x4 = xsq*xsq; | |
219 | expected = (315*x4 - 210*xsq + 15)/8; | |
220 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(5, x), expected, tolerance); | |
221 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-6, x), expected, tolerance); | |
222 | ||
223 | // P_6(x) = (231x^6 -315*x^4 +105x^2 -5)/16 => P_6'(x) = (6*231*x^5 - 4*315*x^3 + 105x)/16 | |
224 | expected = 21*x*(33*x4 - 30*xsq + 5)/8; | |
225 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(6, x), expected, tolerance); | |
226 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-7, x), expected, tolerance); | |
227 | ||
228 | // Mathematica: D[LegendreP[7, x],x] | |
229 | T x6 = x4*xsq; | |
230 | expected = 7*(429*x6 -495*x4 + 135*xsq - 5)/16; | |
231 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(7, x), expected, tolerance); | |
232 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-8, x), expected, tolerance); | |
233 | ||
234 | // Mathematica: D[LegendreP[8, x],x] | |
235 | // The naive polynomial evaluation algorithm is going to get worse from here out, so this will be enough. | |
236 | expected = 9*x*(715*x6 - 1001*x4 + 385*xsq - 35)/16; | |
237 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(8, x), expected, tolerance); | |
238 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-9, x), expected, tolerance); | |
239 | ||
240 | x += static_cast<T>(1)/static_cast<T>(pow(T(2), T(4))); | |
241 | } | |
242 | ||
243 | int n = 0; | |
244 | while (n < 5000) | |
245 | { | |
246 | T expected = n*(n+1)*boost::math::constants::half<T>(); | |
247 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(n, (T) 1), expected, tolerance); | |
248 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-n - 1, (T) 1), expected, tolerance); | |
249 | if (n & 1) | |
250 | { | |
251 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(n, (T) -1), expected, tolerance); | |
252 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-n - 1, (T) -1), expected, tolerance); | |
253 | } | |
254 | else | |
255 | { | |
256 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(n, (T) -1), -expected, tolerance); | |
257 | BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-n - 1, (T) -1), -expected, tolerance); | |
258 | } | |
259 | ++n; | |
260 | } | |
261 | } | |
262 | ||
263 | template<class Real> | |
264 | void test_legendre_p_zeros() | |
265 | { | |
266 | std::cout << "Testing Legendre zeros on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n"; | |
267 | using std::sqrt; | |
268 | using std::abs; | |
269 | using boost::math::legendre_p_zeros; | |
270 | using boost::math::legendre_p; | |
271 | using boost::math::constants::third; | |
272 | Real tol = std::numeric_limits<Real>::epsilon(); | |
273 | ||
274 | // Check the trivial cases: | |
275 | std::vector<Real> zeros = legendre_p_zeros<Real>(1); | |
276 | BOOST_ASSERT(zeros.size() == 1); | |
277 | BOOST_CHECK_SMALL(zeros[0], tol); | |
278 | ||
279 | zeros = legendre_p_zeros<Real>(2); | |
280 | BOOST_ASSERT(zeros.size() == 1); | |
281 | BOOST_CHECK_CLOSE_FRACTION(zeros[0], (Real) 1/ sqrt(static_cast<Real>(3)), tol); | |
282 | ||
283 | zeros = legendre_p_zeros<Real>(3); | |
284 | BOOST_ASSERT(zeros.size() == 2); | |
285 | BOOST_CHECK_SMALL(zeros[0], tol); | |
286 | BOOST_CHECK_CLOSE_FRACTION(zeros[1], sqrt(static_cast<Real>(3)/static_cast<Real>(5)), tol); | |
287 | ||
288 | zeros = legendre_p_zeros<Real>(4); | |
289 | BOOST_ASSERT(zeros.size() == 2); | |
290 | BOOST_CHECK_CLOSE_FRACTION(zeros[0], sqrt( (15-2*sqrt(static_cast<Real>(30)))/static_cast<Real>(35) ), tol); | |
291 | BOOST_CHECK_CLOSE_FRACTION(zeros[1], sqrt( (15+2*sqrt(static_cast<Real>(30)))/static_cast<Real>(35) ), tol); | |
292 | ||
293 | ||
294 | zeros = legendre_p_zeros<Real>(5); | |
295 | BOOST_ASSERT(zeros.size() == 3); | |
296 | BOOST_CHECK_SMALL(zeros[0], tol); | |
297 | BOOST_CHECK_CLOSE_FRACTION(zeros[1], third<Real>()*sqrt( (35 - 2*sqrt(static_cast<Real>(70)))/static_cast<Real>(7) ), 2*tol); | |
298 | BOOST_CHECK_CLOSE_FRACTION(zeros[2], third<Real>()*sqrt( (35 + 2*sqrt(static_cast<Real>(70)))/static_cast<Real>(7) ), 2*tol); | |
299 | ||
300 | // Don't take the tolerances too seriously. | |
301 | // The other test shows that the zeros are estimated more accurately than the function! | |
302 | for (int n = 6; n < 130; ++n) | |
303 | { | |
304 | zeros = legendre_p_zeros<Real>(n); | |
305 | if (n & 1) | |
306 | { | |
307 | BOOST_CHECK(zeros.size() == (n-1)/2 +1); | |
308 | BOOST_CHECK_SMALL(zeros[0], tol); | |
309 | } | |
310 | else | |
311 | { | |
312 | // Zero is not a zero of the odd Legendre polynomials | |
313 | BOOST_CHECK(zeros.size() == n/2); | |
314 | BOOST_CHECK(zeros[0] > 0); | |
315 | BOOST_CHECK_SMALL(legendre_p(n, zeros[0]), 550*tol); | |
316 | } | |
317 | Real previous_zero = zeros[0]; | |
318 | for (int k = 1; k < zeros.size(); ++k) | |
319 | { | |
320 | Real next_zero = zeros[k]; | |
321 | BOOST_CHECK(next_zero > previous_zero); | |
322 | ||
323 | std::string err = "Tolerance failed for (n, k) = (" + boost::lexical_cast<std::string>(n) + "," + boost::lexical_cast<std::string>(k) + ")\n"; | |
324 | if (n < 40) | |
325 | { | |
326 | BOOST_CHECK_MESSAGE( abs(legendre_p(n, next_zero)) < 100*tol, | |
327 | err); | |
328 | } | |
329 | else | |
330 | { | |
331 | BOOST_CHECK_MESSAGE( abs(legendre_p(n, next_zero)) < 1000*tol, | |
332 | err); | |
333 | } | |
334 | previous_zero = next_zero; | |
335 | } | |
336 | // The zeros of orthogonal polynomials are contained strictly in (a, b). | |
337 | BOOST_CHECK(previous_zero < 1); | |
338 | } | |
339 | return; | |
340 | } | |
341 | ||
342 | int test_legendre_p_zeros_double_ulp(int min_x, int max_n) | |
343 | { | |
344 | std::cout << "Testing ULP distance for Legendre zeros.\n"; | |
345 | using std::abs; | |
346 | using boost::math::legendre_p_zeros; | |
347 | using boost::math::float_distance; | |
348 | using boost::multiprecision::cpp_bin_float_quad; | |
349 | ||
350 | double max_float_distance = 0; | |
351 | for (int n = min_x; n < max_n; ++n) | |
352 | { | |
353 | std::vector<double> double_zeros = legendre_p_zeros<double>(n); | |
354 | std::vector<cpp_bin_float_quad> quad_zeros = legendre_p_zeros<cpp_bin_float_quad>(n); | |
355 | BOOST_ASSERT(quad_zeros.size() == double_zeros.size()); | |
356 | for (int k = 0; k < (int)double_zeros.size(); ++k) | |
357 | { | |
358 | double d = abs(float_distance(double_zeros[k], quad_zeros[k].convert_to<double>())); | |
359 | if (d > max_float_distance) | |
360 | { | |
361 | max_float_distance = d; | |
362 | } | |
363 | } | |
364 | } | |
365 | ||
366 | return (int) max_float_distance; | |
367 | } |