1 // Copyright Nick Thompson, 2017
2 // Copyright John Maddock 2017
3 // Use, modification and distribution are subject to the
4 // Boost Software License, Version 1.0.
5 // (See accompanying file LICENSE_1_0.txt
6 // or copy at http://www.boost.org/LICENSE_1_0.txt)
14 #include <boost/math/constants/constants.hpp>
15 #include <boost/math/special_functions/cbrt.hpp>
16 #include <boost/math/special_functions/factorials.hpp>
17 #include <boost/math/special_functions/gamma.hpp>
18 #include <boost/math/tools/roots.hpp>
19 #include <boost/noncopyable.hpp>
21 #define CPP_BIN_FLOAT 1
22 #define CPP_DEC_FLOAT 2
23 #define CPP_MPFR_FLOAT 3
25 //#define MP_TYPE CPP_BIN_FLOAT
26 #define MP_TYPE CPP_DEC_FLOAT
27 //#define MP_TYPE CPP_MPFR_FLOAT
31 struct digits_characteristics
33 static const int digits10
= 300;
34 static const int guard_digits
= 6;
38 #if (MP_TYPE == CPP_BIN_FLOAT)
39 #include <boost/multiprecision/cpp_bin_float.hpp>
40 namespace mp
= boost::multiprecision
;
41 typedef mp::number
<mp::cpp_bin_float
<digits_characteristics::digits10
+ digits_characteristics::guard_digits
>, mp::et_off
> mp_type
;
42 #elif (MP_TYPE == CPP_DEC_FLOAT)
43 #include <boost/multiprecision/cpp_dec_float.hpp>
44 namespace mp
= boost::multiprecision
;
45 typedef mp::number
<mp::cpp_dec_float
<digits_characteristics::digits10
+ digits_characteristics::guard_digits
>, mp::et_off
> mp_type
;
46 #elif (MP_TYPE == CPP_MPFR_FLOAT)
47 #include <boost/multiprecision/mpfr.hpp>
48 namespace mp
= boost::multiprecision
;
49 typedef mp::number
<mp::mpfr_float_backend
<digits_characteristics::digits10
+ digits_characteristics::guard_digits
>, mp::et_off
> mp_type
;
51 #error MP_TYPE is undefined
55 class laguerre_function_object
58 laguerre_function_object(const int n
, const T a
) : order(n
),
63 laguerre_function_object(const laguerre_function_object
& other
) : order(other
.order
),
68 ~laguerre_function_object() { }
70 T
operator()(const T
& x
) const
72 // Calculate (via forward recursion):
73 // * the value of the Laguerre function L(n, alpha, x), called (p2),
74 // * the value of the derivative of the Laguerre function (d2),
75 // * and the value of the corresponding Laguerre function of
76 // previous order (p1).
78 // Return the value of the function (p2) in order to be used as a
79 // function object with Boost.Math root-finding. Store the values
80 // of the Laguerre function derivative (d2) and the Laguerre function
81 // of previous order (p1) in class members for later use.
87 T
j_plus_alpha(alpha
);
88 T
two_j_plus_one_plus_alpha_minus_x(1 + alpha
- x
);
94 for(j
= 0; j
< order
; ++j
)
98 // Set the value of the previous Laguerre function.
101 // Use a recurrence relation to compute the value of the Laguerre function.
102 p2
= ((two_j_plus_one_plus_alpha_minus_x
* p1
) - (j_plus_alpha
* p0
)) / (j
+ 1);
105 two_j_plus_one_plus_alpha_minus_x
+= my_two
;
108 // Set the value of the derivative of the Laguerre function.
109 d2
= ((p2
* j
) - (j_plus_alpha
* p1
)) / x
;
111 // Return the value of the Laguerre function.
115 const T
& previous () const { return p1
; }
116 const T
& derivative() const { return d2
; }
118 static bool root_tolerance(const T
& a
, const T
& b
)
122 // The relative tolerance here is: ((a - b) * 2) / (a + b).
123 return (abs((a
- b
) * 2) < ((a
+ b
) * boost::math::tools::epsilon
<T
>()));
132 laguerre_function_object();
134 const laguerre_function_object
& operator=(const laguerre_function_object
&);
138 class guass_laguerre_abscissas_and_weights
: private boost::noncopyable
141 guass_laguerre_abscissas_and_weights(const int n
, const T a
) : order(n
),
149 // TBD: If we ever boostify this, throw a range error here.
150 // If so, then also document it in the docs.
151 std::cout
<< "Range error: the order of the Laguerre function must exceed -20.0." << std::endl
;
159 virtual ~guass_laguerre_abscissas_and_weights() { }
161 const std::vector
<T
>& abscissas() const { return xi
; }
162 const std::vector
<T
>& weights () const { return wi
; }
164 bool get_valid() const { return valid
; }
178 std::cout
<< "finding approximate roots..." << std::endl
;
180 std::vector
<boost::math::tuple
<T
, T
> > root_estimates
;
182 root_estimates
.reserve(static_cast<typename
std::vector
<boost::math::tuple
<T
, T
> >::size_type
>(order
));
184 const laguerre_function_object
<T
> laguerre_object(order
, alpha
);
186 // Set the initial values of the step size and the running step
187 // to be used for finding the estimate of the first root.
191 T first_laguerre_root
= 0.0F
;
193 bool first_laguerre_root_has_been_found
= true;
197 // Iteratively step through the Laguerre function using a
198 // small step-size in order to find a rough estimate of
201 bool this_laguerre_value_is_negative
= (laguerre_object(mp_type(0)) < 0);
203 static const int j_max
= 10000;
207 for(j
= 0; (j
< j_max
) && (this_laguerre_value_is_negative
!= (laguerre_object(step
) < 0)); ++j
)
209 // Increment the step size until the sign of the Laguerre function
210 // switches. This indicates a zero-crossing, signalling the next root.
216 first_laguerre_root_has_been_found
= false;
220 // We have found the first zero-crossing. Put a loose bracket around
221 // the root using a window. Here, we know that the first root lies
222 // between (x - step_size) < root < x.
224 // Before storing the approximate root, perform a couple of
225 // bisection steps in order to tighten up the root bracket.
226 boost::uintmax_t a_couple_of_iterations
= 3U;
227 const std::pair
<T
, T
>
228 first_laguerre_root
= boost::math::tools::bisect(laguerre_object
,
231 laguerre_function_object
<T
>::root_tolerance
,
232 a_couple_of_iterations
);
234 static_cast<void>(a_couple_of_iterations
);
239 // Calculate an estimate of the 1st root of a generalized Laguerre
240 // function using either a Taylor series or an expansion in Bessel
241 // function zeros. The Bessel function zeros expansion is from Tricomi.
243 // Here, we obtain an estimate of the first zero of J_alpha(x).
249 // For small alpha, use a short series obtained from Mathematica(R).
250 // Series[BesselJZero[v, 1], {v, 0, 3}]
252 j_alpha_m1
= ((( 0.09748661784476F
253 * alpha
- 0.17549359276115F
)
254 * alpha
+ 1.54288974259931F
)
255 * alpha
+ 2.40482555769577F
);
259 // For larger alpha, use the first line of Eqs. 10.21.40 in the NIST Handbook.
260 const T
alpha_pow_third(boost::math::cbrt(alpha
));
261 const T
alpha_pow_minus_two_thirds(T(1) / (alpha_pow_third
* alpha_pow_third
));
263 j_alpha_m1
= alpha
* ((((( + 0.043F
264 * alpha_pow_minus_two_thirds
- 0.0908F
)
265 * alpha_pow_minus_two_thirds
- 0.00397F
)
266 * alpha_pow_minus_two_thirds
+ 1.033150F
)
267 * alpha_pow_minus_two_thirds
+ 1.8557571F
)
268 * alpha_pow_minus_two_thirds
+ 1.0F
);
271 const T vf
= ((order
* 4.0F
) + (alpha
* 2.0F
) + 2.0F
);
272 const T vf2
= vf
* vf
;
273 const T j_alpha_m1_sqr
= j_alpha_m1
* j_alpha_m1
;
275 first_laguerre_root
= (j_alpha_m1_sqr
* (-0.6666666666667F
+ ((0.6666666666667F
* alpha
) * alpha
) + (0.3333333333333F
* j_alpha_m1_sqr
) + vf2
)) / (vf2
* vf
);
278 if(first_laguerre_root_has_been_found
)
280 bool this_laguerre_value_is_negative
= (laguerre_object(mp_type(0)) < 0);
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;
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.
293 while(static_cast<int>(root_estimates
.size()) < order
)
295 // Increment the step size until the sign of the Laguerre function
296 // switches. This indicates a zero-crossing, signalling the next root.
299 if(this_laguerre_value_is_negative
!= (laguerre_object(step
) < 0))
301 // We have found the next zero-crossing.
303 // Change the running sign of the Laguerre function.
304 this_laguerre_value_is_negative
= (!this_laguerre_value_is_negative
);
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.
310 // Before storing the approximate root, perform a couple of
311 // bisection steps in order to tighten up the root bracket.
312 boost::uintmax_t a_couple_of_iterations
= 3U;
313 const std::pair
<T
, T
>
314 root_estimate_bracket
= boost::math::tools::bisect(laguerre_object
,
317 laguerre_function_object
<T
>::root_tolerance
,
318 a_couple_of_iterations
);
320 static_cast<void>(a_couple_of_iterations
);
322 // Store the refined root estimate as a bracketed range in a tuple.
323 root_estimates
.push_back(boost::math::tuple
<T
, T
>(root_estimate_bracket
.first
,
324 root_estimate_bracket
.second
));
326 if(root_estimates
.size() >= static_cast<std::size_t>(2U))
328 // Determine the next step size. This is based on the distance between
329 // the previous two roots, whereby the estimates of the previous roots
330 // are computed by taking the average of the lower and upper range of
331 // the root-estimate bracket.
333 const T r0
= ( boost::math::get
<0>(*(root_estimates
.rbegin() + 1U))
334 + boost::math::get
<1>(*(root_estimates
.rbegin() + 1U))) / 2;
336 const T r1
= ( boost::math::get
<0>(*root_estimates
.rbegin())
337 + boost::math::get
<1>(*root_estimates
.rbegin())) / 2;
339 const T distance_between_previous_roots
= r1
- r0
;
341 step_size
= distance_between_previous_roots
/ 3;
347 ((alpha
== 0) ? T(-1)
348 : -boost::math::tgamma(alpha
+ order
) / boost::math::factorial
<T
>(order
- 1));
350 xi
.reserve(root_estimates
.size());
351 wi
.reserve(root_estimates
.size());
353 // Calculate the abscissas and weights to full precision.
354 for(std::size_t i
= static_cast<std::size_t>(0U); i
< root_estimates
.size(); ++i
)
356 std::cout
<< "calculating abscissa and weight for index: " << i
<< std::endl
;
358 // Calculate the abscissas using iterative root-finding.
360 // Select the maximum allowed iterations, being at least 20.
361 // The determination of the maximum allowed iterations is
362 // based on the number of decimal digits in the numerical
364 const int my_digits10
= static_cast<int>(static_cast<float>(boost::math::tools::digits
<T
>()) * 0.301F
);
365 const boost::uintmax_t number_of_iterations_allowed
= (std::max
)(20, my_digits10
/ 2);
367 boost::uintmax_t number_of_iterations_used
= number_of_iterations_allowed
;
369 // Perform the root-finding using ACM TOMS 748 from Boost.Math.
370 const std::pair
<T
, T
>
371 laguerre_root_bracket
= boost::math::tools::toms748_solve(laguerre_object
,
372 boost::math::get
<0>(root_estimates
[i
]),
373 boost::math::get
<1>(root_estimates
[i
]),
374 laguerre_function_object
<T
>::root_tolerance
,
375 number_of_iterations_used
);
377 // Based on the result of *each* root-finding operation, re-assess
378 // the validity of the Guass-Laguerre abscissas and weights object.
379 valid
&= (number_of_iterations_used
< number_of_iterations_allowed
);
381 // Compute the Laguerre root as the average of the values from
382 // the solved root bracket.
383 const T laguerre_root
= ( laguerre_root_bracket
.first
384 + laguerre_root_bracket
.second
) / 2;
386 // Calculate the weight for this Laguerre root. Here, we calculate
387 // the derivative of the Laguerre function and the value of the
388 // previous Laguerre function on the x-axis at the value of this
390 static_cast<void>(laguerre_object(laguerre_root
));
392 // Store the abscissa and weight for this index.
393 xi
.push_back(laguerre_root
);
394 wi
.push_back(norm_g
/ ((laguerre_object
.derivative() * order
) * laguerre_object
.previous()));
403 struct gauss_laguerre_ai
406 gauss_laguerre_ai(const T X
) : x(X
)
411 zeta
= ((sqrt(x
) * x
) * 2) / 3;
413 const T zeta_times_48_pow_sixth
= sqrt(boost::math::cbrt(zeta
* 48));
415 factor
= 1 / ((sqrt(boost::math::constants::pi
<T
>()) * zeta_times_48_pow_sixth
) * (exp(zeta
) * gamma_of_five_sixths()));
418 gauss_laguerre_ai(const gauss_laguerre_ai
& other
) : x (other
.x
),
420 factor(other
.factor
) { }
422 T
operator()(const T
& t
) const
426 return factor
/ sqrt(boost::math::cbrt(2 + (t
/ zeta
)));
434 static const T
& gamma_of_five_sixths()
436 static const T value
= boost::math::tgamma(T(5) / 6);
441 const gauss_laguerre_ai
& operator=(const gauss_laguerre_ai
&);
445 T
gauss_laguerre_airy_ai(const T x
)
447 static const float digits_factor
= static_cast<float>(std::numeric_limits
<mp_type
>::digits10
) / 300.0F
;
448 static const int laguerre_order
= static_cast<int>(600.0F
* digits_factor
);
450 static const guass_laguerre_abscissas_and_weights
<T
> abscissas_and_weights(laguerre_order
, -T(1) / 6);
454 if(abscissas_and_weights
.get_valid())
456 const gauss_laguerre_ai
<T
> this_gauss_laguerre_ai(x
);
459 std::inner_product(abscissas_and_weights
.abscissas().begin(),
460 abscissas_and_weights
.abscissas().end(),
461 abscissas_and_weights
.weights().begin(),
464 [&this_gauss_laguerre_ai
](const T
& this_abscissa
, const T
& this_weight
) -> T
466 return this_gauss_laguerre_ai(this_abscissa
) * this_weight
;
471 // TBD: Consider an error message.
472 airy_ai_result
= T(0);
475 return airy_ai_result
;
481 // Use Gauss-Laguerre integration to compute airy_ai(120 / 7).
490 // 3.8990420982303275013276114626640705170145070824318e-22
493 // 3.899042098230327501327611462664070517014507082431797677146153303523108862015228
494 // 864136051942933142648e-22
497 // 3.899042098230327501327611462664070517014507082431797677146153303523108862015228
498 // 86413605194293314264788265460938200890998546786740097437064263800719644346113699
499 // 77010905030516409847054404055843899790277e-22
502 // 3.899042098230327501327611462664070517014507082431797677146153303523108862015228
503 // 86413605194293314264788265460938200890998546786740097437064263800719644346113699
504 // 77010905030516409847054404055843899790277083960877617919088116211775232728792242
505 // 9346416823281460245814808276654088201413901972239996130752528e-22
508 // 3.899042098230327501327611462664070517014507082431797677146153303523108862015228
509 // 86413605194293314264788265460938200890998546786740097437064263800719644346113699
510 // 77010905030516409847054404055843899790277083960877617919088116211775232728792242
511 // 93464168232814602458148082766540882014139019722399961307525276722937464859521685
512 // 42826483602153339361960948844649799257455597165900957281659632186012043089610827
513 // 78871305322190941528281744734605934497977375094921646511687434038062987482900167
514 // 45127557400365419545e-22
516 // Mathematica(R) or Wolfram's Alpha:
517 // N[AiryAi[120 / 7], 300]
518 std::cout
<< std::setprecision(digits_characteristics::digits10
)
519 << gauss_laguerre_airy_ai(mp_type(120) / 7)