ceph/src/boost/libs/math/example/bessel_zeros_example_1.cpp

   1
   2 // Copyright Christopher Kormanyos 2013.
   3 // Copyright Paul A. Bristow 2013.
   4 // Copyright John Maddock 2013.
   5
   6 // Distributed under the Boost Software License, Version 1.0.
   7 // (See accompanying file LICENSE_1_0.txt or
   8 // copy at http://www.boost.org/LICENSE_1_0.txt).
   9
  10 #ifdef _MSC_VER
  11 #  pragma warning (disable : 4512) // assignment operator could not be generated.
  12 #  pragma warning (disable : 4996) // assignment operator could not be generated.
  13 #endif
  14
  15 #include <iostream>
  16 #include <limits>
  17 #include <vector>
  18 #include <algorithm>
  19 #include <iomanip>
  20 #include <iterator>
  21
  22 // Weisstein, Eric W. "Bessel Function Zeros." From MathWorld--A Wolfram Web Resource.
  23 // http://mathworld.wolfram.com/BesselFunctionZeros.html
  24 // Test values can be calculated using [@wolframalpha.com WolframAplha]
  25 // See also http://dlmf.nist.gov/10.21
  26
  27 //[bessel_zeros_example_1
  28
  29 /*`This example demonstrates calculating zeros of the Bessel and Neumann functions.
  30 It also shows how Boost.Math and Boost.Multiprecision can be combined to provide
  31 a many decimal digit precision. For 50 decimal digit precision we need to include
  32 */
  33
  34   #include <boost/multiprecision/cpp_dec_float.hpp>
  35
  36 /*`and a `typedef` for `float_type` may be convenient
  37 (allowing a quick switch to re-compute at built-in `double` or other precision)
  38 */
  39   typedef boost::multiprecision::cpp_dec_float_50 float_type;
  40
  41 //`To use the functions for finding zeros of the functions we need
  42
  43   #include <boost/math/special_functions/bessel.hpp>
  44
  45 //`This file includes the forward declaration signatures for the zero-finding functions:
  46
  47 //  #include <boost/math/special_functions/math_fwd.hpp>
  48
  49 /*`but more details are in the full documentation, for example at
  50 [@http://www.boost.org/doc/libs/1_53_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/bessel/bessel_over.html Boost.Math Bessel functions].
  51 */
  52
  53 /*`This example shows obtaining both a single zero of the Bessel function,
  54 and then placing multiple zeros into a container like `std::vector` by providing an iterator.
  55 */
  56 //] [/bessel_zeros_example_1]
  57
  58 /*The signature of the single value function is:
  59
  60   template <class T>
  61   inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type
  62     cyl_bessel_j_zero(
  63            T v,      // Floating-point value for Jv.
  64            int m);   // start index.
  65
  66 The result type is controlled by the floating-point type of parameter `v`
  67 (but subject to the usual __precision_policy and __promotion_policy).
  68
  69 The signature of multiple zeros function is:
  70
  71   template <class T, class OutputIterator>
  72   inline OutputIterator cyl_bessel_j_zero(
  73                                 T v,                      // Floating-point value for Jv.
  74                                 int start_index,          // 1-based start index.
  75                                 unsigned number_of_zeros, // How many zeros to generate
  76                                 OutputIterator out_it);   // Destination for zeros.
  77
  78 There is also a version which allows control of the __policy_section for error handling and precision.
  79
  80   template <class T, class OutputIterator, class Policy>
  81   inline OutputIterator cyl_bessel_j_zero(
  82                                 T v,                      // Floating-point value for Jv.
  83                                 int start_index,          // 1-based start index.
  84                                 unsigned number_of_zeros, // How many zeros to generate
  85                                 OutputIterator out_it,    // Destination for zeros.
  86                                 const Policy& pol);       // Policy to use.
  87 */
  88
  89 int main()
  90 {
  91   try
  92   {
  93 //[bessel_zeros_example_2
  94
  95 /*`[tip It is always wise to place code using Boost.Math inside try'n'catch blocks;
  96 this will ensure that helpful error messages are shown when exceptional conditions arise.]
  97
  98 First, evaluate a single Bessel zero.
  99
 100 The precision is controlled by the float-point type of template parameter `T` of `v`
 101 so this example has `double` precision, at least 15 but up to 17 decimal digits (for the common 64-bit double).
 102 */
 103 //    double root = boost::math::cyl_bessel_j_zero(0.0, 1);
 104 //    // Displaying with default precision of 6 decimal digits:
 105 //    std::cout << "boost::math::cyl_bessel_j_zero(0.0, 1) " << root << std::endl; // 2.40483
 106 //    // And with all the guaranteed (15) digits:
 107 //    std::cout.precision(std::numeric_limits<double>::digits10);
 108 //    std::cout << "boost::math::cyl_bessel_j_zero(0.0, 1) " << root << std::endl; // 2.40482555769577
 109 /*`But note that because the parameter `v` controls the precision of the result,
 110 `v` [*must be a floating-point type].
 111 So if you provide an integer type, say 0, rather than 0.0, then it will fail to compile thus:
 112 ``
 113     root = boost::math::cyl_bessel_j_zero(0, 1);
 114 ``
 115 with this error message
 116 ``
 117   error C2338: Order must be a floating-point type.
 118 ``
 119
 120 Optionally, we can use a policy to ignore errors, C-style, returning some value,
 121 perhaps infinity or NaN, or the best that can be done. (See __user_error_handling).
 122
 123 To create a (possibly unwise!) policy `ignore_all_policy` that ignores all errors:
 124 */
 125
 126   typedef boost::math::policies::policy<
 127     boost::math::policies::domain_error<boost::math::policies::ignore_error>,
 128     boost::math::policies::overflow_error<boost::math::policies::ignore_error>,
 129     boost::math::policies::underflow_error<boost::math::policies::ignore_error>,
 130     boost::math::policies::denorm_error<boost::math::policies::ignore_error>,
 131     boost::math::policies::pole_error<boost::math::policies::ignore_error>,
 132     boost::math::policies::evaluation_error<boost::math::policies::ignore_error>
 133               > ignore_all_policy;
 134  //`Examples of use of this `ignore_all_policy` are
 135
 136     double inf = std::numeric_limits<double>::infinity();
 137     double nan = std::numeric_limits<double>::quiet_NaN();
 138
 139     double dodgy_root = boost::math::cyl_bessel_j_zero(-1.0, 1, ignore_all_policy());
 140     std::cout << "boost::math::cyl_bessel_j_zero(-1.0, 1) " << dodgy_root << std::endl; // 1.#QNAN
 141     double inf_root = boost::math::cyl_bessel_j_zero(inf, 1, ignore_all_policy());
 142     std::cout << "boost::math::cyl_bessel_j_zero(inf, 1) " << inf_root << std::endl; // 1.#QNAN
 143     double nan_root = boost::math::cyl_bessel_j_zero(nan, 1, ignore_all_policy());
 144     std::cout << "boost::math::cyl_bessel_j_zero(nan, 1) " << nan_root << std::endl; // 1.#QNAN
 145
 146 /*`Another version of `cyl_bessel_j_zero`  allows calculation of multiple zeros with one call,
 147 placing the results in a container, often `std::vector`.
 148 For example, generate and display the first five `double` roots of J[sub v] for integral order 2,
 149 as column ['J[sub 2](x)] in table 1 of
 150 [@ http://mathworld.wolfram.com/BesselFunctionZeros.html Wolfram Bessel Function Zeros].
 151 */
 152     unsigned int n_roots = 5U;
 153     std::vector<double> roots;
 154     boost::math::cyl_bessel_j_zero(2.0, 1, n_roots, std::back_inserter(roots));
 155     std::copy(roots.begin(),
 156               roots.end(),
 157               std::ostream_iterator<double>(std::cout, "\n"));
 158
 159 /*`Or we can use Boost.Multiprecision to generate 50 decimal digit roots of ['J[sub v]]
 160 for non-integral order `v= 71/19 == 3.736842`, expressed as an exact-integer fraction
 161 to generate the most accurate value possible for all floating-point types.
 162
 163 We set the precision of the output stream, and show trailing zeros to display a fixed 50 decimal digits.
 164 */
 165     std::cout.precision(std::numeric_limits<float_type>::digits10); // 50 decimal digits.
 166     std::cout << std::showpoint << std::endl; // Show trailing zeros.
 167
 168     float_type x = float_type(71) / 19;
 169     float_type r = boost::math::cyl_bessel_j_zero(x, 1); // 1st root.
 170     std::cout << "x = " << x << ", r = " << r << std::endl;
 171
 172     r = boost::math::cyl_bessel_j_zero(x, 20U); // 20th root.
 173     std::cout << "x = " << x << ", r = " << r << std::endl;
 174
 175     std::vector<float_type> zeros;
 176     boost::math::cyl_bessel_j_zero(x, 1, 3, std::back_inserter(zeros));
 177
 178     std::cout << "cyl_bessel_j_zeros" << std::endl;
 179     // Print the roots to the output stream.
 180     std::copy(zeros.begin(), zeros.end(),
 181               std::ostream_iterator<float_type>(std::cout, "\n"));
 182 //] [/bessel_zeros_example_2]
 183   }
 184   catch (std::exception const& ex)
 185   {
 186     std::cout << "Thrown exception " << ex.what() << std::endl;
 187   }
 188
 189  } // int main()
 190
 191  /*
 192
 193  Output:
 194
 195    Description: Autorun "J:\Cpp\big_number\Debug\bessel_zeros_example_1.exe"
 196   boost::math::cyl_bessel_j_zero(-1.0, 1) 3.83171
 197   boost::math::cyl_bessel_j_zero(inf, 1) 1.#QNAN
 198   boost::math::cyl_bessel_j_zero(nan, 1) 1.#QNAN
 199   5.13562
 200   8.41724
 201   11.6198
 202   14.796
 203   17.9598
 204
 205   x = 3.7368421052631578947368421052631578947368421052632, r = 7.2731751938316489503185694262290765588963196701623
 206   x = 3.7368421052631578947368421052631578947368421052632, r = 67.815145619696290925556791375555951165111460585458
 207   cyl_bessel_j_zeros
 208   7.2731751938316489503185694262290765588963196701623
 209   10.724858308883141732536172745851416647110749599085
 210   14.018504599452388106120459558042660282427471931581
 211
 212 */
 213