mirror of
https://github.com/saitohirga/WSJT-X.git
synced 2025-10-31 04:50:34 -04:00
448 lines
18 KiB
C++
448 lines
18 KiB
C++
// Copyright Christopher Kormanyos 2013.
|
|
// Copyright Paul A. Bristow 2013.
|
|
// Copyright John Maddock 2013.
|
|
|
|
// Distributed under the Boost Software License, Version 1.0.
|
|
// (See accompanying file LICENSE_1_0.txt or
|
|
// copy at http://www.boost.org/LICENSE_1_0.txt).
|
|
|
|
#ifdef _MSC_VER
|
|
# pragma warning (disable : 4512) // assignment operator could not be generated.
|
|
# pragma warning (disable : 4996) // assignment operator could not be generated.
|
|
#endif
|
|
|
|
#include <iostream>
|
|
#include <limits>
|
|
#include <vector>
|
|
#include <algorithm>
|
|
#include <iomanip>
|
|
#include <iterator>
|
|
|
|
// Weisstein, Eric W. "Bessel Function Zeros." From MathWorld--A Wolfram Web Resource.
|
|
// http://mathworld.wolfram.com/BesselFunctionZeros.html
|
|
// Test values can be calculated using [@wolframalpha.com WolframAplha]
|
|
// See also http://dlmf.nist.gov/10.21
|
|
|
|
//[bessel_zero_example_1
|
|
|
|
/*`This example demonstrates calculating zeros of the Bessel, Neumann and Airy functions.
|
|
It also shows how Boost.Math and Boost.Multiprecision can be combined to provide
|
|
a many decimal digit precision. For 50 decimal digit precision we need to include
|
|
*/
|
|
|
|
#include <boost/multiprecision/cpp_dec_float.hpp>
|
|
|
|
/*`and a `typedef` for `float_type` may be convenient
|
|
(allowing a quick switch to re-compute at built-in `double` or other precision)
|
|
*/
|
|
typedef boost::multiprecision::cpp_dec_float_50 float_type;
|
|
|
|
//`To use the functions for finding zeros of the functions we need
|
|
|
|
#include <boost/math/special_functions/bessel.hpp>
|
|
|
|
//`This file includes the forward declaration signatures for the zero-finding functions:
|
|
|
|
// #include <boost/math/special_functions/math_fwd.hpp>
|
|
|
|
/*`but more details are in the full documentation, for example at
|
|
[@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]
|
|
*/
|
|
|
|
/*`This example shows obtaining both a single zero of the Bessel function,
|
|
and then placing multiple zeros into a container like `std::vector` by providing an iterator.
|
|
The signature of the single value function is:
|
|
|
|
template <class T>
|
|
inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type
|
|
cyl_bessel_j_zero(T v, // Floating-point value for Jv.
|
|
int m); // start index.
|
|
|
|
The result type is controlled by the floating-point type of parameter `v`
|
|
(but subject to the usual __precision_policy and __promotion_policy).
|
|
|
|
The signature of multiple zeros function is:
|
|
|
|
template <class T, class OutputIterator>
|
|
inline OutputIterator cyl_bessel_j_zero(T v, // Floating-point value for Jv.
|
|
int start_index, // 1-based start index.
|
|
unsigned number_of_zeros,
|
|
OutputIterator out_it); // iterator into container for zeros.
|
|
|
|
There is also a version which allows control of the __policy_section for error handling and precision.
|
|
|
|
template <class T, class OutputIterator, class Policy>
|
|
inline OutputIterator cyl_bessel_j_zero(T v, // Floating-point value for Jv.
|
|
int start_index, // 1-based start index.
|
|
unsigned number_of_zeros,
|
|
OutputIterator out_it,
|
|
const Policy& pol); // iterator into container for zeros.
|
|
|
|
*/
|
|
//] [/bessel_zero_example_1]
|
|
|
|
//[bessel_zero_example_iterator_1]
|
|
/*`We use the `cyl_bessel_j_zero` output iterator parameter `out_it`
|
|
to create a sum of 1/zeros[super 2] by defining a custom output iterator:
|
|
*/
|
|
|
|
template <class T>
|
|
struct output_summation_iterator
|
|
{
|
|
output_summation_iterator(T* p) : p_sum(p)
|
|
{}
|
|
output_summation_iterator& operator*()
|
|
{ return *this; }
|
|
output_summation_iterator& operator++()
|
|
{ return *this; }
|
|
output_summation_iterator& operator++(int)
|
|
{ return *this; }
|
|
output_summation_iterator& operator = (T const& val)
|
|
{
|
|
*p_sum += 1./ (val * val); // Summing 1/zero^2.
|
|
return *this;
|
|
}
|
|
private:
|
|
T* p_sum;
|
|
};
|
|
|
|
|
|
//] [/bessel_zero_example_iterator_1]
|
|
|
|
int main()
|
|
{
|
|
try
|
|
{
|
|
//[bessel_zero_example_2]
|
|
|
|
/*`[tip It is always wise to place code using Boost.Math inside try'n'catch blocks;
|
|
this will ensure that helpful error messages can be shown when exceptional conditions arise.]
|
|
|
|
First, evaluate a single Bessel zero.
|
|
|
|
The precision is controlled by the float-point type of template parameter `T` of `v`
|
|
so this example has `double` precision, at least 15 but up to 17 decimal digits (for the common 64-bit double).
|
|
*/
|
|
double root = boost::math::cyl_bessel_j_zero(0.0, 1);
|
|
// Displaying with default precision of 6 decimal digits:
|
|
std::cout << "boost::math::cyl_bessel_j_zero(0.0, 1) " << root << std::endl; // 2.40483
|
|
// And with all the guaranteed (15) digits:
|
|
std::cout.precision(std::numeric_limits<double>::digits10);
|
|
std::cout << "boost::math::cyl_bessel_j_zero(0.0, 1) " << root << std::endl; // 2.40482555769577
|
|
/*`But note that because the parameter `v` controls the precision of the result,
|
|
`v` [*must be a floating-point type].
|
|
So if you provide an integer type, say 0, rather than 0.0, then it will fail to compile thus:
|
|
``
|
|
root = boost::math::cyl_bessel_j_zero(0, 1);
|
|
``
|
|
with this error message
|
|
``
|
|
error C2338: Order must be a floating-point type.
|
|
``
|
|
|
|
Optionally, we can use a policy to ignore errors, C-style, returning some value
|
|
perhaps infinity or NaN, or the best that can be done. (See __user_error_handling).
|
|
|
|
To create a (possibly unwise!) policy that ignores all errors:
|
|
*/
|
|
|
|
typedef boost::math::policies::policy
|
|
<
|
|
boost::math::policies::domain_error<boost::math::policies::ignore_error>,
|
|
boost::math::policies::overflow_error<boost::math::policies::ignore_error>,
|
|
boost::math::policies::underflow_error<boost::math::policies::ignore_error>,
|
|
boost::math::policies::denorm_error<boost::math::policies::ignore_error>,
|
|
boost::math::policies::pole_error<boost::math::policies::ignore_error>,
|
|
boost::math::policies::evaluation_error<boost::math::policies::ignore_error>
|
|
> ignore_all_policy;
|
|
|
|
double inf = std::numeric_limits<double>::infinity();
|
|
double nan = std::numeric_limits<double>::quiet_NaN();
|
|
|
|
std::cout << "boost::math::cyl_bessel_j_zero(-1.0, 0) " << std::endl;
|
|
double dodgy_root = boost::math::cyl_bessel_j_zero(-1.0, 0, ignore_all_policy());
|
|
std::cout << "boost::math::cyl_bessel_j_zero(-1.0, 1) " << dodgy_root << std::endl; // 1.#QNAN
|
|
double inf_root = boost::math::cyl_bessel_j_zero(inf, 1, ignore_all_policy());
|
|
std::cout << "boost::math::cyl_bessel_j_zero(inf, 1) " << inf_root << std::endl; // 1.#QNAN
|
|
double nan_root = boost::math::cyl_bessel_j_zero(nan, 1, ignore_all_policy());
|
|
std::cout << "boost::math::cyl_bessel_j_zero(nan, 1) " << nan_root << std::endl; // 1.#QNAN
|
|
|
|
/*`Another version of `cyl_bessel_j_zero` allows calculation of multiple zeros with one call,
|
|
placing the results in a container, often `std::vector`.
|
|
For example, generate five `double` roots of J[sub v] for integral order 2.
|
|
|
|
showing the same results as column J[sub 2](x) in table 1 of
|
|
[@ http://mathworld.wolfram.com/BesselFunctionZeros.html Wolfram Bessel Function Zeros].
|
|
|
|
*/
|
|
unsigned int n_roots = 5U;
|
|
std::vector<double> roots;
|
|
boost::math::cyl_bessel_j_zero(2.0, 1, n_roots, std::back_inserter(roots));
|
|
std::copy(roots.begin(),
|
|
roots.end(),
|
|
std::ostream_iterator<double>(std::cout, "\n"));
|
|
|
|
/*`Or generate 50 decimal digit roots of J[sub v] for non-integral order `v = 71/19`.
|
|
|
|
We set the precision of the output stream and show trailing zeros to display a fixed 50 decimal digits.
|
|
*/
|
|
std::cout.precision(std::numeric_limits<float_type>::digits10); // 50 decimal digits.
|
|
std::cout << std::showpoint << std::endl; // Show trailing zeros.
|
|
|
|
float_type x = float_type(71) / 19;
|
|
float_type r = boost::math::cyl_bessel_j_zero(x, 1); // 1st root.
|
|
std::cout << "x = " << x << ", r = " << r << std::endl;
|
|
|
|
r = boost::math::cyl_bessel_j_zero(x, 20U); // 20th root.
|
|
std::cout << "x = " << x << ", r = " << r << std::endl;
|
|
|
|
std::vector<float_type> zeros;
|
|
boost::math::cyl_bessel_j_zero(x, 1, 3, std::back_inserter(zeros));
|
|
|
|
std::cout << "cyl_bessel_j_zeros" << std::endl;
|
|
// Print the roots to the output stream.
|
|
std::copy(zeros.begin(), zeros.end(),
|
|
std::ostream_iterator<float_type>(std::cout, "\n"));
|
|
|
|
/*`The Neumann function zeros are evaluated very similarly:
|
|
*/
|
|
using boost::math::cyl_neumann_zero;
|
|
|
|
double zn = cyl_neumann_zero(2., 1);
|
|
|
|
std::cout << "cyl_neumann_zero(2., 1) = " << std::endl;
|
|
//double zn0 = zn;
|
|
// std::cout << "zn0 = " << std::endl;
|
|
// std::cout << zn0 << std::endl;
|
|
//
|
|
std::cout << zn << std::endl;
|
|
// std::cout << cyl_neumann_zero(2., 1) << std::endl;
|
|
|
|
std::vector<float> nzeros(3); // Space for 3 zeros.
|
|
cyl_neumann_zero<float>(2.F, 1, nzeros.size(), nzeros.begin());
|
|
|
|
std::cout << "cyl_neumann_zero<float>(2.F, 1, " << std::endl;
|
|
// Print the zeros to the output stream.
|
|
std::copy(nzeros.begin(), nzeros.end(),
|
|
std::ostream_iterator<float>(std::cout, "\n"));
|
|
|
|
std::cout << cyl_neumann_zero(static_cast<float_type>(220)/100, 1) << std::endl;
|
|
// 3.6154383428745996706772556069431792744372398748422
|
|
|
|
/*`Finally we show how the output iterator can be used to compute a sum of zeros.
|
|
|
|
(See [@http://dx.doi.org/10.1017/S2040618500034067 Ian N. Sneddon, Infinite Sums of Bessel Zeros],
|
|
page 150 equation 40).
|
|
*/
|
|
//] [/bessel_zero_example_2]
|
|
|
|
{
|
|
//[bessel_zero_example_iterator_2]
|
|
/*`The sum is calculated for many values, converging on the analytical exact value of `1/8`.
|
|
*/
|
|
using boost::math::cyl_bessel_j_zero;
|
|
double nu = 1.;
|
|
double sum = 0;
|
|
output_summation_iterator<double> it(&sum); // sum of 1/zeros^2
|
|
cyl_bessel_j_zero(nu, 1, 10000, it);
|
|
|
|
double s = 1/(4 * (nu + 1)); // 0.125 = 1/8 is exact analytical solution.
|
|
std::cout << std::setprecision(6) << "nu = " << nu << ", sum = " << sum
|
|
<< ", exact = " << s << std::endl;
|
|
// nu = 1.00000, sum = 0.124990, exact = 0.125000
|
|
//] [/bessel_zero_example_iterator_2]
|
|
}
|
|
}
|
|
catch (std::exception& ex)
|
|
{
|
|
std::cout << "Thrown exception " << ex.what() << std::endl;
|
|
}
|
|
|
|
//[bessel_zero_example_iterator_3]
|
|
|
|
/*`Examples below show effect of 'bad' arguments that throw a `domain_error` exception.
|
|
*/
|
|
try
|
|
{ // Try a negative rank m.
|
|
std::cout << "boost::math::cyl_bessel_j_zero(-1.F, -1) " << std::endl;
|
|
float dodgy_root = boost::math::cyl_bessel_j_zero(-1.F, -1);
|
|
std::cout << "boost::math::cyl_bessel_j_zero(-1.F, -1) " << dodgy_root << std::endl;
|
|
// Throw exception Error in function boost::math::cyl_bessel_j_zero<double>(double, int):
|
|
// Order argument is -1, but must be >= 0 !
|
|
}
|
|
catch (std::exception& ex)
|
|
{
|
|
std::cout << "Throw exception " << ex.what() << std::endl;
|
|
}
|
|
|
|
/*`[note The type shown is the type [*after promotion],
|
|
using __precision_policy and __promotion_policy, from `float` to `double` in this case.]
|
|
|
|
In this example the promotion goes:
|
|
|
|
# Arguments are `float` and `int`.
|
|
# Treat `int` "as if" it were a `double`, so arguments are `float` and `double`.
|
|
# Common type is `double` - so that's the precision we want (and the type that will be returned).
|
|
# Evaluate internally as `long double` for full `double` precision.
|
|
|
|
See full code for other examples that promote from `double` to `long double`.
|
|
|
|
*/
|
|
|
|
//] [/bessel_zero_example_iterator_3]
|
|
try
|
|
{ // order v = inf
|
|
std::cout << "boost::math::cyl_bessel_j_zero(infF, 1) " << std::endl;
|
|
float infF = std::numeric_limits<float>::infinity();
|
|
float inf_root = boost::math::cyl_bessel_j_zero(infF, 1);
|
|
std::cout << "boost::math::cyl_bessel_j_zero(infF, 1) " << inf_root << std::endl;
|
|
// boost::math::cyl_bessel_j_zero(-1.F, -1)
|
|
//Thrown exception Error in function boost::math::cyl_bessel_j_zero<double>(double, int):
|
|
// Requested the -1'th zero, but the rank must be positive !
|
|
}
|
|
catch (std::exception& ex)
|
|
{
|
|
std::cout << "Thrown exception " << ex.what() << std::endl;
|
|
}
|
|
try
|
|
{ // order v = inf
|
|
double inf = std::numeric_limits<double>::infinity();
|
|
double inf_root = boost::math::cyl_bessel_j_zero(inf, 1);
|
|
std::cout << "boost::math::cyl_bessel_j_zero(inf, 1) " << inf_root << std::endl;
|
|
// Throw exception Error in function boost::math::cyl_bessel_j_zero<long double>(long double, unsigned):
|
|
// Order argument is 1.#INF, but must be finite >= 0 !
|
|
}
|
|
catch (std::exception& ex)
|
|
{
|
|
std::cout << "Thrown exception " << ex.what() << std::endl;
|
|
}
|
|
|
|
try
|
|
{ // order v = NaN
|
|
double nan = std::numeric_limits<double>::quiet_NaN();
|
|
double nan_root = boost::math::cyl_bessel_j_zero(nan, 1);
|
|
std::cout << "boost::math::cyl_bessel_j_zero(nan, 1) " << nan_root << std::endl;
|
|
// Throw exception Error in function boost::math::cyl_bessel_j_zero<long double>(long double, unsigned):
|
|
// Order argument is 1.#QNAN, but must be finite >= 0 !
|
|
}
|
|
catch (std::exception& ex)
|
|
{
|
|
std::cout << "Thrown exception " << ex.what() << std::endl;
|
|
}
|
|
|
|
try
|
|
{ // Try a negative m.
|
|
double dodgy_root = boost::math::cyl_bessel_j_zero(0.0, -1);
|
|
// warning C4146: unary minus operator applied to unsigned type, result still unsigned.
|
|
std::cout << "boost::math::cyl_bessel_j_zero(0.0, -1) " << dodgy_root << std::endl;
|
|
// boost::math::cyl_bessel_j_zero(0.0, -1) 6.74652e+009
|
|
// This *should* fail because m is unreasonably large.
|
|
|
|
}
|
|
catch (std::exception& ex)
|
|
{
|
|
std::cout << "Thrown exception " << ex.what() << std::endl;
|
|
}
|
|
|
|
try
|
|
{ // m = inf
|
|
double inf = std::numeric_limits<double>::infinity();
|
|
double inf_root = boost::math::cyl_bessel_j_zero(0.0, inf);
|
|
// warning C4244: 'argument' : conversion from 'double' to 'int', possible loss of data.
|
|
std::cout << "boost::math::cyl_bessel_j_zero(0.0, inf) " << inf_root << std::endl;
|
|
// Throw exception Error in function boost::math::cyl_bessel_j_zero<long double>(long double, int):
|
|
// Requested the 0'th zero, but must be > 0 !
|
|
|
|
}
|
|
catch (std::exception& ex)
|
|
{
|
|
std::cout << "Thrown exception " << ex.what() << std::endl;
|
|
}
|
|
|
|
try
|
|
{ // m = NaN
|
|
std::cout << "boost::math::cyl_bessel_j_zero(0.0, nan) " << std::endl ;
|
|
double nan = std::numeric_limits<double>::quiet_NaN();
|
|
double nan_root = boost::math::cyl_bessel_j_zero(0.0, nan);
|
|
// warning C4244: 'argument' : conversion from 'double' to 'int', possible loss of data.
|
|
std::cout << nan_root << std::endl;
|
|
// Throw exception Error in function boost::math::cyl_bessel_j_zero<long double>(long double, int):
|
|
// Requested the 0'th zero, but must be > 0 !
|
|
}
|
|
catch (std::exception& ex)
|
|
{
|
|
std::cout << "Thrown exception " << ex.what() << std::endl;
|
|
}
|
|
|
|
} // int main()
|
|
|
|
/*
|
|
Mathematica: Table[N[BesselJZero[71/19, n], 50], {n, 1, 20, 1}]
|
|
|
|
7.2731751938316489503185694262290765588963196701623
|
|
10.724858308883141732536172745851416647110749599085
|
|
14.018504599452388106120459558042660282427471931581
|
|
17.25249845917041718216248716654977734919590383861
|
|
20.456678874044517595180234083894285885460502077814
|
|
23.64363089714234522494551422714731959985405172504
|
|
26.819671140255087745421311470965019261522390519297
|
|
29.988343117423674742679141796661432043878868194142
|
|
33.151796897690520871250862469973445265444791966114
|
|
36.3114160002162074157243540350393860813165201842
|
|
39.468132467505236587945197808083337887765967032029
|
|
42.622597801391236474855034831297954018844433480227
|
|
45.775281464536847753390206207806726581495950012439
|
|
48.926530489173566198367766817478553992471739894799
|
|
52.076607045343002794279746041878924876873478063472
|
|
55.225712944912571393594224327817265689059002890192
|
|
58.374006101538886436775188150439025201735151418932
|
|
61.521611873000965273726742659353136266390944103571
|
|
64.66863105379093036834648221487366079456596628716
|
|
67.815145619696290925556791375555951165111460585458
|
|
|
|
Mathematica: Table[N[BesselKZero[2, n], 50], {n, 1, 5, 1}]
|
|
n |
|
|
1 | 3.3842417671495934727014260185379031127323883259329
|
|
2 | 6.7938075132682675382911671098369487124493222183854
|
|
3 | 10.023477979360037978505391792081418280789658279097
|
|
|
|
|
|
*/
|
|
|
|
/*
|
|
[bessel_zero_output]
|
|
|
|
boost::math::cyl_bessel_j_zero(0.0, 1) 2.40483
|
|
boost::math::cyl_bessel_j_zero(0.0, 1) 2.40482555769577
|
|
boost::math::cyl_bessel_j_zero(-1.0, 1) 1.#QNAN
|
|
boost::math::cyl_bessel_j_zero(inf, 1) 1.#QNAN
|
|
boost::math::cyl_bessel_j_zero(nan, 1) 1.#QNAN
|
|
5.13562230184068
|
|
8.41724414039986
|
|
11.6198411721491
|
|
14.7959517823513
|
|
17.9598194949878
|
|
|
|
x = 3.7368421052631578947368421052631578947368421052632, r = 7.2731751938316489503185694262290765588963196701623
|
|
x = 3.7368421052631578947368421052631578947368421052632, r = 67.815145619696290925556791375555951165111460585458
|
|
7.2731751938316489503185694262290765588963196701623
|
|
10.724858308883141732536172745851416647110749599085
|
|
14.018504599452388106120459558042660282427471931581
|
|
cyl_neumann_zero(2., 1) = 3.3842417671495935000000000000000000000000000000000
|
|
3.3842418193817139000000000000000000000000000000000
|
|
6.7938075065612793000000000000000000000000000000000
|
|
10.023477554321289000000000000000000000000000000000
|
|
3.6154383428745996706772556069431792744372398748422
|
|
nu = 1.00000, sum = 0.124990, exact = 0.125000
|
|
Throw exception Error in function boost::math::cyl_bessel_j_zero<double>(double, int): Order argument is -1, but must be >= 0 !
|
|
Throw exception Error in function boost::math::cyl_bessel_j_zero<long double>(long double, int): Order argument is 1.#INF, but must be finite >= 0 !
|
|
Throw exception Error in function boost::math::cyl_bessel_j_zero<long double>(long double, int): Order argument is 1.#QNAN, but must be finite >= 0 !
|
|
Throw exception Error in function boost::math::cyl_bessel_j_zero<long double>(long double, int): Requested the -1'th zero, but must be > 0 !
|
|
Throw exception Error in function boost::math::cyl_bessel_j_zero<long double>(long double, int): Requested the -2147483648'th zero, but must be > 0 !
|
|
Throw exception Error in function boost::math::cyl_bessel_j_zero<long double>(long double, int): Requested the -2147483648'th zero, but must be > 0 !
|
|
|
|
|
|
] [/bessel_zero_output]
|
|
*/
|
|
|