]>
Commit | Line | Data |
---|---|---|
f67539c2 TL |
1 | /////////////////////////////////////////////////////////////////////////////// |
2 | // Copyright Christopher Kormanyos 2012 - 2015, 2020. | |
3 | // Distributed under the Boost Software License, Version 1.0. | |
4 | // (See accompanying file LICENSE_1_0.txt or copy at | |
5 | // http://www.boost.org/LICENSE_1_0.txt) | |
6 | // | |
7 | ||
8 | // This example uses Boost.Multiprecision to implement | |
9 | // a high-precision Gauss-Laguerre quadrature integration. | |
10 | // The quadrature is used to calculate the airy_ai(x) function | |
11 | // for real-valued x on the positive axis with x.ge.1. | |
12 | ||
13 | // In this way, the integral representation could be seen | |
14 | // as part of a scheme to calculate real-valued Airy functions | |
15 | // on the positive axis for medium to large argument. | |
16 | // A Taylor series or hypergeometric function (not part | |
17 | // of this example) could be used for smaller arguments. | |
18 | ||
19 | // This example has been tested with decimal digits counts | |
20 | // ranging from 21...301, by adjusting the parameter | |
21 | // local::my_digits10 at compile time. | |
22 | ||
23 | // The quadrature integral representaion of airy_ai(x) used | |
24 | // in this example can be found in: | |
25 | ||
26 | // A. Gil, J. Segura, N.M. Temme, "Numerical Methods for Special | |
27 | // Functions" (SIAM Press 2007), Sect. 5.3.3, in particular Eq. 5.110, | |
28 | // page 145. Subsequently, Gil et al's book cites the another work: | |
29 | // W. Gautschi, "Computation of Bessel and Airy functions and of | |
30 | // related Gaussian quadrature formulae", BIT, 42 (2002), pp. 110-118. | |
92f5a8d4 | 31 | |
7c673cae FG |
32 | #include <cmath> |
33 | #include <cstdint> | |
34 | #include <functional> | |
35 | #include <iomanip> | |
36 | #include <iostream> | |
37 | #include <numeric> | |
f67539c2 TL |
38 | #include <sstream> |
39 | #include <tuple> | |
40 | #include <vector> | |
41 | ||
42 | #include <boost/cstdfloat.hpp> | |
7c673cae FG |
43 | #include <boost/math/constants/constants.hpp> |
44 | #include <boost/math/special_functions/cbrt.hpp> | |
f67539c2 | 45 | #include <boost/math/special_functions/bessel.hpp> |
7c673cae FG |
46 | #include <boost/math/special_functions/factorials.hpp> |
47 | #include <boost/math/special_functions/gamma.hpp> | |
48 | #include <boost/math/tools/roots.hpp> | |
f67539c2 | 49 | #include <boost/multiprecision/cpp_dec_float.hpp> |
7c673cae FG |
50 | #include <boost/noncopyable.hpp> |
51 | ||
f67539c2 | 52 | namespace gauss { namespace laguerre { |
7c673cae | 53 | |
f67539c2 | 54 | namespace util { |
7c673cae | 55 | |
f67539c2 | 56 | void progress_bar(std::ostream& os, const float percent); |
7c673cae | 57 | |
f67539c2 | 58 | void progress_bar(std::ostream& os, const float percent) |
7c673cae | 59 | { |
f67539c2 | 60 | std::stringstream strstrm; |
7c673cae | 61 | |
f67539c2 | 62 | strstrm.precision(1); |
7c673cae | 63 | |
f67539c2 | 64 | strstrm << std::fixed << percent << "%"; |
7c673cae | 65 | |
f67539c2 | 66 | os << strstrm.str() << "\r"; |
7c673cae | 67 | |
f67539c2 TL |
68 | os.flush(); |
69 | } | |
7c673cae | 70 | |
f67539c2 | 71 | } |
7c673cae | 72 | |
f67539c2 TL |
73 | namespace detail |
74 | { | |
75 | template<typename T> | |
76 | class laguerre_l_object BOOST_FINAL | |
77 | { | |
78 | public: | |
1e59de90 | 79 | laguerre_l_object(const int n, const T a) noexcept |
f67539c2 TL |
80 | : order(n), |
81 | alpha(a), | |
82 | p1 (0), | |
83 | d2 (0) { } | |
7c673cae | 84 | |
f67539c2 TL |
85 | laguerre_l_object& operator=(const laguerre_l_object& other) |
86 | { | |
87 | if(this != other) | |
88 | { | |
89 | order = other.order; | |
90 | alpha = other.alpha; | |
91 | p1 = other.p1; | |
92 | d2 = other.d2; | |
93 | } | |
7c673cae | 94 | |
f67539c2 TL |
95 | return *this; |
96 | } | |
7c673cae | 97 | |
1e59de90 | 98 | T operator()(const T& x) const noexcept |
7c673cae | 99 | { |
f67539c2 TL |
100 | // Calculate (via forward recursion): |
101 | // * the value of the Laguerre function L(n, alpha, x), called (p2), | |
102 | // * the value of the derivative of the Laguerre function (d2), | |
103 | // * and the value of the corresponding Laguerre function of | |
104 | // previous order (p1). | |
7c673cae | 105 | |
f67539c2 TL |
106 | // Return the value of the function (p2) in order to be used as a |
107 | // function object with Boost.Math root-finding. Store the values | |
108 | // of the Laguerre function derivative (d2) and the Laguerre function | |
109 | // of previous order (p1) in class members for later use. | |
7c673cae | 110 | |
f67539c2 TL |
111 | p1 = T(0); |
112 | T p2 = T(1); | |
113 | d2 = T(0); | |
7c673cae | 114 | |
f67539c2 TL |
115 | T j_plus_alpha = alpha; |
116 | T two_j_plus_one_plus_alpha_minus_x = (1 + alpha) - x; | |
7c673cae | 117 | |
f67539c2 | 118 | const T my_two = 2; |
7c673cae | 119 | |
f67539c2 TL |
120 | for(int j = 0; j < order; ++j) |
121 | { | |
122 | const T p0(p1); | |
7c673cae | 123 | |
f67539c2 TL |
124 | // Set the value of the previous Laguerre function. |
125 | p1 = p2; | |
7c673cae | 126 | |
f67539c2 TL |
127 | // Use a recurrence relation to compute the value of the Laguerre function. |
128 | p2 = ((two_j_plus_one_plus_alpha_minus_x * p1) - (j_plus_alpha * p0)) / (j + 1); | |
7c673cae | 129 | |
f67539c2 TL |
130 | ++j_plus_alpha; |
131 | two_j_plus_one_plus_alpha_minus_x += my_two; | |
132 | } | |
7c673cae | 133 | |
f67539c2 TL |
134 | // Set the value of the derivative of the Laguerre function. |
135 | d2 = ((p2 * order) - (j_plus_alpha * p1)) / x; | |
7c673cae | 136 | |
f67539c2 TL |
137 | // Return the value of the Laguerre function. |
138 | return p2; | |
139 | } | |
7c673cae | 140 | |
1e59de90 TL |
141 | const T previous () const noexcept { return p1; } |
142 | const T derivative() const noexcept { return d2; } | |
7c673cae | 143 | |
1e59de90 | 144 | static bool root_tolerance(const T& a, const T& b) noexcept |
7c673cae | 145 | { |
f67539c2 TL |
146 | using std::fabs; |
147 | ||
148 | // The relative tolerance here is: ((a - b) * 2) / (a + b). | |
149 | return ((fabs(a - b) * 2) < (fabs(a + b) * std::numeric_limits<T>::epsilon())); | |
7c673cae | 150 | } |
7c673cae | 151 | |
f67539c2 TL |
152 | private: |
153 | const int order; | |
154 | const T alpha; | |
155 | mutable T p1; | |
156 | mutable T d2; | |
157 | }; | |
7c673cae | 158 | |
f67539c2 TL |
159 | template<typename T> |
160 | class abscissas_and_weights : private boost::noncopyable | |
161 | { | |
162 | public: | |
163 | abscissas_and_weights(const int n, const T a) : order(n), | |
164 | alpha(a), | |
165 | xi (), | |
166 | wi () | |
167 | { | |
168 | if(alpha < -20.0F) | |
169 | { | |
170 | // Users can decide to perform different range checking. | |
171 | std::cout << "Range error: the order of the Laguerre function must exceed -20.0." | |
172 | << std::endl; | |
173 | } | |
174 | else | |
175 | { | |
176 | calculate(); | |
177 | } | |
178 | } | |
7c673cae | 179 | |
1e59de90 TL |
180 | const std::vector<T>& abscissa_n() const noexcept { return xi; } |
181 | const std::vector<T>& weight_n () const noexcept { return wi; } | |
7c673cae | 182 | |
f67539c2 TL |
183 | private: |
184 | const int order; | |
185 | const T alpha; | |
7c673cae | 186 | |
f67539c2 TL |
187 | std::vector<T> xi; |
188 | std::vector<T> wi; | |
7c673cae | 189 | |
f67539c2 TL |
190 | abscissas_and_weights() : order(), |
191 | alpha(), | |
192 | xi (), | |
193 | wi () { } | |
7c673cae | 194 | |
f67539c2 TL |
195 | void calculate() |
196 | { | |
197 | using std::abs; | |
7c673cae | 198 | |
f67539c2 | 199 | std::cout << "Finding the approximate roots..." << std::endl; |
7c673cae | 200 | |
f67539c2 | 201 | std::vector<std::tuple<T, T>> root_estimates; |
7c673cae | 202 | |
f67539c2 | 203 | root_estimates.reserve(static_cast<typename std::vector<std::tuple<T, T>>::size_type>(order)); |
7c673cae | 204 | |
f67539c2 | 205 | const laguerre_l_object<T> laguerre_root_object(order, alpha); |
7c673cae | 206 | |
f67539c2 TL |
207 | // Set the initial values of the step size and the running step |
208 | // to be used for finding the estimate of the first root. | |
209 | T step_size = 0.01F; | |
210 | T step = step_size; | |
7c673cae | 211 | |
f67539c2 | 212 | T first_laguerre_root = 0.0F; |
7c673cae | 213 | |
f67539c2 TL |
214 | if(alpha < -1.0F) |
215 | { | |
216 | // Iteratively step through the Laguerre function using a | |
217 | // small step-size in order to find a rough estimate of | |
218 | // the first zero. | |
7c673cae | 219 | |
f67539c2 | 220 | const bool this_laguerre_value_is_negative = (laguerre_root_object(T(0)) < 0); |
7c673cae | 221 | |
1e59de90 | 222 | constexpr int j_max = 10000; |
7c673cae | 223 | |
f67539c2 | 224 | int j = 0; |
7c673cae | 225 | |
f67539c2 TL |
226 | while((j < j_max) && (this_laguerre_value_is_negative != (laguerre_root_object(step) < 0))) |
227 | { | |
228 | // Increment the step size until the sign of the Laguerre function | |
229 | // switches. This indicates a zero-crossing, signalling the next root. | |
230 | step += step_size; | |
7c673cae | 231 | |
f67539c2 TL |
232 | ++j; |
233 | } | |
7c673cae FG |
234 | } |
235 | else | |
236 | { | |
f67539c2 TL |
237 | // Calculate an estimate of the 1st root of a generalized Laguerre |
238 | // function using either a Taylor series or an expansion in Bessel | |
239 | // function zeros. The Bessel function zeros expansion is from Tricomi. | |
7c673cae | 240 | |
f67539c2 | 241 | // Here, we obtain an estimate of the first zero of cyl_bessel_j(alpha, x). |
7c673cae | 242 | |
f67539c2 | 243 | T j_alpha_m1; |
7c673cae | 244 | |
f67539c2 TL |
245 | if(alpha < 1.4F) |
246 | { | |
247 | // For small alpha, use a short series obtained from Mathematica(R). | |
248 | // Series[BesselJZero[v, 1], {v, 0, 3}] | |
249 | // N[%, 12] | |
250 | j_alpha_m1 = ((( T(0.09748661784476F) | |
251 | * alpha - T(0.17549359276115F)) | |
252 | * alpha + T(1.54288974259931F)) | |
253 | * alpha + T(2.40482555769577F)); | |
254 | } | |
255 | else | |
256 | { | |
257 | // For larger alpha, use the first line of Eqs. 10.21.40 in the NIST Handbook. | |
258 | const T alpha_pow_third(boost::math::cbrt(alpha)); | |
259 | const T alpha_pow_minus_two_thirds(T(1) / (alpha_pow_third * alpha_pow_third)); | |
260 | ||
261 | j_alpha_m1 = alpha * ((((( + T(0.043F) | |
262 | * alpha_pow_minus_two_thirds - T(0.0908F)) | |
263 | * alpha_pow_minus_two_thirds - T(0.00397F)) | |
264 | * alpha_pow_minus_two_thirds + T(1.033150F)) | |
265 | * alpha_pow_minus_two_thirds + T(1.8557571F)) | |
266 | * alpha_pow_minus_two_thirds + T(1.0F)); | |
267 | } | |
7c673cae | 268 | |
f67539c2 TL |
269 | const T vf = ( T(order * 4U) |
270 | + T(alpha * 2U) | |
271 | + T(2U)); | |
272 | const T vf2 = vf * vf; | |
273 | const T j_alpha_m1_sqr = j_alpha_m1 * j_alpha_m1; | |
7c673cae | 274 | |
f67539c2 TL |
275 | first_laguerre_root = (j_alpha_m1_sqr * ( -T(0.6666666666667F) |
276 | + ((T(0.6666666666667F) * alpha) * alpha) | |
277 | + (T(0.3333333333333F) * j_alpha_m1_sqr) + vf2)) / (vf2 * vf); | |
7c673cae FG |
278 | } |
279 | ||
f67539c2 | 280 | bool this_laguerre_value_is_negative = (laguerre_root_object(T(0)) < 0); |
7c673cae FG |
281 | |
282 | // Re-set the initial value of the step-size based on the | |
283 | // estimate of the first root. | |
284 | step_size = first_laguerre_root / 2; | |
285 | step = step_size; | |
286 | ||
287 | // Step through the Laguerre function using a step-size | |
288 | // of dynamic width in order to find the zero crossings | |
289 | // of the Laguerre function, providing rough estimates | |
290 | // of the roots. Refine the brackets with a few bisection | |
291 | // steps, and store the results as bracketed root estimates. | |
292 | ||
f67539c2 | 293 | while(root_estimates.size() < static_cast<std::size_t>(order)) |
7c673cae FG |
294 | { |
295 | // Increment the step size until the sign of the Laguerre function | |
296 | // switches. This indicates a zero-crossing, signalling the next root. | |
297 | step += step_size; | |
298 | ||
f67539c2 | 299 | if(this_laguerre_value_is_negative != (laguerre_root_object(step) < 0)) |
7c673cae FG |
300 | { |
301 | // We have found the next zero-crossing. | |
302 | ||
303 | // Change the running sign of the Laguerre function. | |
304 | this_laguerre_value_is_negative = (!this_laguerre_value_is_negative); | |
305 | ||
306 | // We have found the first zero-crossing. Put a loose bracket around | |
307 | // the root using a window. Here, we know that the first root lies | |
308 | // between (x - step_size) < root < x. | |
309 | ||
310 | // Before storing the approximate root, perform a couple of | |
311 | // bisection steps in order to tighten up the root bracket. | |
1e59de90 | 312 | std::uintmax_t a_couple_of_iterations = 4U; |
f67539c2 | 313 | |
7c673cae | 314 | const std::pair<T, T> |
f67539c2 | 315 | root_estimate_bracket = boost::math::tools::bisect(laguerre_root_object, |
7c673cae FG |
316 | step - step_size, |
317 | step, | |
f67539c2 | 318 | laguerre_l_object<T>::root_tolerance, |
7c673cae FG |
319 | a_couple_of_iterations); |
320 | ||
321 | static_cast<void>(a_couple_of_iterations); | |
322 | ||
323 | // Store the refined root estimate as a bracketed range in a tuple. | |
f67539c2 TL |
324 | root_estimates.push_back(std::tuple<T, T>(std::get<0>(root_estimate_bracket), |
325 | std::get<1>(root_estimate_bracket))); | |
326 | ||
327 | if( (root_estimates.size() == 1U) | |
328 | || ((root_estimates.size() % 8U) == 0U) | |
329 | || (root_estimates.size() == static_cast<std::size_t>(order))) | |
330 | { | |
331 | const float progress = (100.0F * static_cast<float>(root_estimates.size())) / static_cast<float>(order); | |
332 | ||
333 | std::cout << "root_estimates.size(): " | |
334 | << root_estimates.size() | |
335 | << ", " | |
336 | ; | |
337 | ||
338 | util::progress_bar(std::cout, progress); | |
339 | } | |
7c673cae FG |
340 | |
341 | if(root_estimates.size() >= static_cast<std::size_t>(2U)) | |
342 | { | |
343 | // Determine the next step size. This is based on the distance between | |
344 | // the previous two roots, whereby the estimates of the previous roots | |
345 | // are computed by taking the average of the lower and upper range of | |
346 | // the root-estimate bracket. | |
347 | ||
f67539c2 TL |
348 | const T r0 = ( std::get<0>(*(root_estimates.crbegin() + 1U)) |
349 | + std::get<1>(*(root_estimates.crbegin() + 1U))) / 2; | |
7c673cae | 350 | |
f67539c2 TL |
351 | const T r1 = ( std::get<0>(*root_estimates.crbegin()) |
352 | + std::get<1>(*root_estimates.crbegin())) / 2; | |
7c673cae FG |
353 | |
354 | const T distance_between_previous_roots = r1 - r0; | |
355 | ||
356 | step_size = distance_between_previous_roots / 3; | |
357 | } | |
358 | } | |
359 | } | |
360 | ||
361 | const T norm_g = | |
362 | ((alpha == 0) ? T(-1) | |
363 | : -boost::math::tgamma(alpha + order) / boost::math::factorial<T>(order - 1)); | |
364 | ||
365 | xi.reserve(root_estimates.size()); | |
366 | wi.reserve(root_estimates.size()); | |
367 | ||
f67539c2 TL |
368 | std::cout << std::endl; |
369 | ||
7c673cae FG |
370 | // Calculate the abscissas and weights to full precision. |
371 | for(std::size_t i = static_cast<std::size_t>(0U); i < root_estimates.size(); ++i) | |
372 | { | |
f67539c2 TL |
373 | if( ((i % 8U) == 0U) |
374 | || ( i == root_estimates.size() - 1U)) | |
375 | { | |
376 | const float progress = (100.0F * static_cast<float>(i + 1U)) / static_cast<float>(root_estimates.size()); | |
377 | ||
378 | std::cout << "Calculating abscissas and weights. Processed " | |
379 | << (i + 1U) | |
380 | << ", " | |
381 | ; | |
382 | ||
383 | util::progress_bar(std::cout, progress); | |
384 | } | |
7c673cae FG |
385 | |
386 | // Calculate the abscissas using iterative root-finding. | |
387 | ||
388 | // Select the maximum allowed iterations, being at least 20. | |
389 | // The determination of the maximum allowed iterations is | |
390 | // based on the number of decimal digits in the numerical | |
391 | // type T. | |
1e59de90 | 392 | constexpr int local_math_tools_digits10 = |
f67539c2 TL |
393 | static_cast<int>(static_cast<boost::float_least32_t>(boost::math::tools::digits<T>()) * BOOST_FLOAT32_C(0.301)); |
394 | ||
1e59de90 | 395 | const std::uintmax_t number_of_iterations_allowed = (std::max)(20, local_math_tools_digits10 / 2); |
7c673cae | 396 | |
1e59de90 | 397 | std::uintmax_t number_of_iterations_used = number_of_iterations_allowed; |
7c673cae FG |
398 | |
399 | // Perform the root-finding using ACM TOMS 748 from Boost.Math. | |
400 | const std::pair<T, T> | |
f67539c2 TL |
401 | laguerre_root_bracket = boost::math::tools::toms748_solve(laguerre_root_object, |
402 | std::get<0>(root_estimates.at(i)), | |
403 | std::get<1>(root_estimates.at(i)), | |
404 | laguerre_l_object<T>::root_tolerance, | |
7c673cae FG |
405 | number_of_iterations_used); |
406 | ||
f67539c2 | 407 | static_cast<void>(number_of_iterations_used); |
7c673cae FG |
408 | |
409 | // Compute the Laguerre root as the average of the values from | |
410 | // the solved root bracket. | |
f67539c2 TL |
411 | const T laguerre_root = ( std::get<0>(laguerre_root_bracket) |
412 | + std::get<1>(laguerre_root_bracket)) / 2; | |
7c673cae FG |
413 | |
414 | // Calculate the weight for this Laguerre root. Here, we calculate | |
415 | // the derivative of the Laguerre function and the value of the | |
416 | // previous Laguerre function on the x-axis at the value of this | |
417 | // Laguerre root. | |
f67539c2 | 418 | static_cast<void>(laguerre_root_object(laguerre_root)); |
7c673cae FG |
419 | |
420 | // Store the abscissa and weight for this index. | |
421 | xi.push_back(laguerre_root); | |
f67539c2 | 422 | wi.push_back(norm_g / ((laguerre_root_object.derivative() * order) * laguerre_root_object.previous())); |
7c673cae | 423 | } |
f67539c2 TL |
424 | |
425 | std::cout << std::endl; | |
7c673cae | 426 | } |
f67539c2 | 427 | }; |
7c673cae | 428 | |
7c673cae | 429 | template<typename T> |
f67539c2 | 430 | struct airy_ai_object BOOST_FINAL |
7c673cae FG |
431 | { |
432 | public: | |
f67539c2 TL |
433 | airy_ai_object(const T& x) |
434 | : my_x (x), | |
435 | my_zeta (((sqrt(x) * x) * 2) / 3), | |
436 | my_factor(make_factor(my_zeta)) { } | |
437 | ||
438 | T operator()(const T& t) const | |
7c673cae | 439 | { |
7c673cae FG |
440 | using std::sqrt; |
441 | ||
f67539c2 | 442 | return my_factor / sqrt(boost::math::cbrt(2 + (t / my_zeta))); |
7c673cae FG |
443 | } |
444 | ||
f67539c2 TL |
445 | private: |
446 | const T my_x; | |
447 | const T my_zeta; | |
448 | const T my_factor; | |
7c673cae | 449 | |
f67539c2 TL |
450 | airy_ai_object() : my_x (), |
451 | my_zeta (), | |
452 | my_factor() { } | |
453 | ||
454 | static T make_factor(const T& z) | |
7c673cae | 455 | { |
f67539c2 | 456 | using std::exp; |
7c673cae FG |
457 | using std::sqrt; |
458 | ||
f67539c2 | 459 | return 1 / ((sqrt(boost::math::constants::pi<T>()) * sqrt(boost::math::cbrt(z * 48))) * (exp(z) * boost::math::tgamma(T(5) / 6))); |
7c673cae | 460 | } |
f67539c2 TL |
461 | }; |
462 | } // namespace detail | |
7c673cae | 463 | |
f67539c2 | 464 | } } // namespace gauss::laguerre |
7c673cae | 465 | |
7c673cae | 466 | |
f67539c2 TL |
467 | // A float_type is created to handle the desired number of decimal digits from `cpp_dec_float` without using __expression_templates. |
468 | struct local | |
469 | { | |
1e59de90 | 470 | static constexpr unsigned int my_digits10 = 101U; |
7c673cae | 471 | |
f67539c2 TL |
472 | typedef boost::multiprecision::number<boost::multiprecision::cpp_dec_float<my_digits10>, |
473 | boost::multiprecision::et_off> | |
474 | float_type; | |
475 | }; | |
7c673cae | 476 | |
1e59de90 | 477 | static_assert(local::my_digits10 > 20U, |
f67539c2 TL |
478 | "Error: This example is intended to have more than 20 decimal digits"); |
479 | ||
480 | int main() | |
481 | { | |
482 | // Use Gauss-Laguerre quadrature integration to compute airy_ai(x / 7) | |
483 | // with 7 <= x <= 120 and where x is incremented in steps of 1. | |
484 | ||
485 | // During development of this example, we have empirically found | |
486 | // the numbers of Gauss-Laguerre coefficients needed for convergence | |
487 | // when using various counts of base-10 digits. | |
488 | ||
489 | // Let's calibrate, for instance, the number of coefficients needed | |
490 | // at the point x = 1. | |
491 | ||
492 | // Empirical data lead to: | |
493 | // Fit[{{21.0, 3.5}, {51.0, 11.1}, {101.0, 22.5}, {201.0, 46.8}}, {1, d, d^2}, d] | |
494 | // FullSimplify[%] | |
495 | // -1.28301 + (0.235487 + 0.0000178915 d) d | |
496 | ||
497 | // We need significantly more coefficients at smaller real values than are needed | |
498 | // at larger real values because the slope derivative of airy_ai(x) gets more | |
499 | // steep as x approaches zero. | |
500 | ||
501 | // This Gauss-Laguerre quadrature is designed for airy_ai(x) with real-valued x >= 1. | |
502 | ||
1e59de90 | 503 | constexpr boost::float_least32_t d = static_cast<boost::float_least32_t>(std::numeric_limits<local::float_type>::digits10); |
f67539c2 | 504 | |
1e59de90 | 505 | constexpr boost::float_least32_t laguerre_order_factor = -1.28301F + ((0.235487F + (0.0000178915F * d)) * d); |
f67539c2 | 506 | |
1e59de90 | 507 | constexpr int laguerre_order = static_cast<int>(laguerre_order_factor * d); |
f67539c2 TL |
508 | |
509 | std::cout << "std::numeric_limits<local::float_type>::digits10: " << std::numeric_limits<local::float_type>::digits10 << std::endl; | |
510 | std::cout << "laguerre_order: " << laguerre_order << std::endl; | |
511 | ||
512 | typedef gauss::laguerre::detail::abscissas_and_weights<local::float_type> abscissas_and_weights_type; | |
513 | ||
514 | const abscissas_and_weights_type the_abscissas_and_weights(laguerre_order, local::float_type(-1) / 6); | |
515 | ||
516 | bool result_is_ok = true; | |
517 | ||
518 | for(std::uint32_t u = 7U; u < 121U; ++u) | |
7c673cae | 519 | { |
f67539c2 | 520 | const local::float_type x = local::float_type(u) / 7; |
7c673cae | 521 | |
f67539c2 | 522 | typedef gauss::laguerre::detail::airy_ai_object<local::float_type> airy_ai_object_type; |
7c673cae | 523 | |
f67539c2 | 524 | const airy_ai_object_type the_airy_ai_object(x); |
7c673cae | 525 | |
f67539c2 TL |
526 | const local::float_type airy_ai_value = |
527 | std::inner_product(the_abscissas_and_weights.abscissa_n().cbegin(), | |
528 | the_abscissas_and_weights.abscissa_n().cend(), | |
529 | the_abscissas_and_weights.weight_n().cbegin(), | |
530 | local::float_type(0U), | |
531 | std::plus<local::float_type>(), | |
532 | [&the_airy_ai_object](const local::float_type& this_abscissa, | |
533 | const local::float_type& this_weight) -> local::float_type | |
534 | { | |
535 | return the_airy_ai_object(this_abscissa) * this_weight; | |
536 | }); | |
537 | ||
538 | static const local::float_type one_third = 1.0F / local::float_type(3U); | |
539 | ||
540 | static const local::float_type one_over_pi_times_one_over_sqrt_three = | |
541 | sqrt(one_third) * boost::math::constants::one_div_pi<local::float_type>(); | |
542 | ||
543 | const local::float_type sqrt_x = sqrt(x); | |
544 | ||
545 | const local::float_type airy_ai_control = | |
546 | (sqrt_x * one_over_pi_times_one_over_sqrt_three) | |
547 | * boost::math::cyl_bessel_k(one_third, ((2.0F * x) * sqrt_x) * one_third); | |
548 | ||
549 | std::cout << std::setprecision(std::numeric_limits<local::float_type>::digits10) | |
550 | << "airy_ai_value : " | |
551 | << airy_ai_value | |
552 | << std::endl; | |
553 | ||
554 | std::cout << std::setprecision(std::numeric_limits<local::float_type>::digits10) | |
555 | << "airy_ai_control: " | |
556 | << airy_ai_control | |
557 | << std::endl; | |
558 | ||
559 | const local::float_type delta = fabs(1.0F - (airy_ai_control / airy_ai_value)); | |
560 | ||
561 | static const local::float_type tol("1E-" + boost::lexical_cast<std::string>(std::numeric_limits<local::float_type>::digits10 - 7U)); | |
7c673cae | 562 | |
f67539c2 | 563 | result_is_ok &= (delta < tol); |
7c673cae | 564 | } |
7c673cae | 565 | |
f67539c2 TL |
566 | std::cout << std::endl |
567 | << "Total... result_is_ok: " | |
568 | << std::boolalpha | |
569 | << result_is_ok | |
7c673cae | 570 | << std::endl; |
f67539c2 TL |
571 | } // int main() |
572 | ||
573 | /* | |
574 | ||
575 | ||
576 | Partial output: | |
577 | ||
578 | //[gauss_laguerre_quadrature_output_1 | |
579 | ||
580 | std::numeric_limits<local::float_type>::digits10: 101 | |
581 | laguerre_order: 2291 | |
582 | ||
583 | Finding the approximate roots... | |
584 | root_estimates.size(): 1, 0.0% | |
585 | root_estimates.size(): 8, 0.3% | |
586 | root_estimates.size(): 16, 0.7% | |
587 | ... | |
588 | root_estimates.size(): 2288, 99.9% | |
589 | root_estimates.size(): 2291, 100.0% | |
590 | ||
591 | ||
592 | Calculating abscissas and weights. Processed 1, 0.0% | |
593 | Calculating abscissas and weights. Processed 9, 0.4% | |
594 | ... | |
595 | Calculating abscissas and weights. Processed 2289, 99.9% | |
596 | Calculating abscissas and weights. Processed 2291, 100.0% | |
597 | //] [/gauss_laguerre_quadrature_output_1] | |
598 | ||
599 | //[gauss_laguerre_quadrature_output_2 | |
600 | ||
601 | airy_ai_value : 0.13529241631288141552414742351546630617494414298833070600910205475763353480226572366348710990874867334 | |
602 | airy_ai_control: 0.13529241631288141552414742351546630617494414298833070600910205475763353480226572366348710990874868323 | |
603 | airy_ai_value : 0.11392302126009621102904231059693500086750049240884734708541630001378825889924647699516200868335286103 | |
604 | airy_ai_control: 0.1139230212600962110290423105969350008675004924088473470854163000137882588992464769951620086833528582 | |
605 | ... | |
606 | airy_ai_value : 3.8990420982303275013276114626640705170145070824317976771461533035231088620152288641360519429331427451e-22 | |
607 | airy_ai_control: 3.8990420982303275013276114626640705170145070824317976771461533035231088620152288641360519429331426481e-22 | |
608 | ||
609 | Total... result_is_ok: true | |
610 | ||
611 | //] [/gauss_laguerre_quadrature_output_2] | |
612 | ||
613 | ||
614 | */ |