boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

improvement to multiprecision sines_table example

Browse Source
This commit is contained in:
pabristow
2019-05-08 11:47:50 +01:00
parent c759981a2a
commit 129daa914f
4 changed files with 107 additions and 91 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
example/Jamfile.v2
View File

@@ -61,6 +61,8 @@ test-suite examples :
[ run error_policies_example.cpp ]
[ run error_policy_example.cpp : : : <exception-handling>off:<build>no ]
[ run f_test.cpp ]
[ run fft_sines_table.cpp : : : [ requires cxx11_numeric_limits ] ]
[ run find_location_example.cpp : : : <exception-handling>off:<build>no ]
[ run find_mean_and_sd_normal.cpp : : : <exception-handling>off:<build>no ]
[ run find_root_example.cpp : : : <exception-handling>off:<build>no ]
@@ -76,7 +78,6 @@ test-suite examples :
[ run lambert_w_simple_examples.cpp : : : [ requires cxx11_numeric_limits ] ]
[ run lambert_w_precision_example.cpp : : : [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : <linkflags>-lquadmath ] [ requires cxx11_numeric_limits cxx11_explicit_conversion_operators ] ]
[ run inverse_gamma_example.cpp ]
[ run inverse_gamma_distribution_example.cpp : : : <exception-handling>off:<build>no ]
[ run laplace_example.cpp : : : <exception-handling>off:<build>no ]

48
example/big_seventh.cpp
View File

@@ -40,38 +40,47 @@ To use these floating-point types and constants, we need some includes:
#include <iostream>
#include <limits>
//` So now we can demonstrate with some trivial calculations:
/*` So now we can demonstrate with some trivial calculations:
*/
//] //[big_seventh_example_1]
int main()
{
/*`Using `typedef cpp_dec_float_50` hides the complexity of multiprecision to allow us
to define variables with 50 decimal digit precision just like built-in `double`.
//[big_seventh_example_2
/*`Using `typedef cpp_dec_float_50` hides the complexity of multiprecision,
allows us to define variables with 50 decimal digit precision just like built-in `double`.
*/
using boost::multiprecision::cpp_dec_float_50;
using boost::multiprecision::cpp_dec_float_50;
cpp_dec_float_50 seventh = cpp_dec_float_50(1) / 7;
cpp_dec_float_50 seventh = cpp_dec_float_50(1) / 7; // 1 / 7
/*`By default, output would only show the standard 6 decimal digits,
so set precision to show all 50 significant digits.
*/
std::cout.precision(std::numeric_limits<cpp_dec_float_50>::digits10);
std::cout << seventh << std::endl;
/*`By default, output would only show the standard 6 decimal digits,
so set precision to show all 50 significant digits, including any trailing zeros.
*/
std::cout.precision(std::numeric_limits<cpp_dec_float_50>::digits10);
std::cout << std::showpoint << std::endl; // Append any trailing zeros.
std::cout << seventh << std::endl;
/*`which outputs:
0.14285714285714285714285714285714285714285714285714
We can also use constants, guaranteed to be initialized with the very last bit of precision.
We can also use Boost.Math __constants like [pi],
guaranteed to be initialized with the very last bit of precision for the floating-point type.
*/
cpp_dec_float_50 circumference = boost::math::constants::pi<cpp_dec_float_50>() * 2 * seventh;
std::cout << circumference << std::endl;
std::cout << "pi = " << boost::math::constants::pi<cpp_dec_float_50>() << std::endl;
cpp_dec_float_50 circumference = boost::math::constants::pi<cpp_dec_float_50>() * 2 * seventh;
std::cout << "c = "<< circumference << std::endl;
/*`which outputs
0.89759790102565521098932668093700082405633411410717
pi = 3.1415926535897932384626433832795028841971693993751
c = 0.89759790102565521098932668093700082405633411410717
*/
//] [/big_seventh_example_1]
//] [/big_seventh_example_2]
return 0;
} // int main()
@@ -80,10 +89,11 @@ We can also use constants, guaranteed to be initialized with the very last bit o
/*
//[big_seventh_example_output
0.14285714285714285714285714285714285714285714285714
0.89759790102565521098932668093700082405633411410717
0.14285714285714285714285714285714285714285714285714
pi = 3.1415926535897932384626433832795028841971693993751
c = 0.89759790102565521098932668093700082405633411410717
//]
//] //[big_seventh_example_output]
*/

122
example/fft_sines_table.cpp
View File

@@ -19,21 +19,35 @@
//[fft_sines_table_example_1
/*`[h5 Using Boost.Multiprecision to generate a high-precision array of sin coefficents for use with FFT.]
/*`[h5 Using Boost.Multiprecision to generate a high-precision array of sine coefficents for use with FFT.]
The Boost.Multiprecision library can be used for computations requiring precision
exceeding that of standard built-in types such as `float`, `double`
and `long double`. For extended-precision calculations, Boost.Multiprecision
supplies a template data type called `cpp_dec_float`. The number of decimal
digits of precision is fixed at compile-time via template parameter.
supplies a template data type called `cpp_bin_float`. The number of decimal
digits of precision is fixed at compile-time via a template parameter.
To use these floating-point types and constants, we need some includes:
One often needs to compute tables of numbers in mathematical software.
To avoid the
[@https://en.wikipedia.org/wiki/Rounding#Table-maker's_dilemma Table-maker's dilemma]
it is necessary to use a higher precision type to compute the table values so that they have
the nearest representable bit-pattern for the type, say `double`, of the table value.
This example is a program `fft_since_table.cpp` that writes a header file `sines.hpp`
containing an array of sine coefficients for use with a Fast Fourier Transform (FFT),
that can be included by the FFT program.
To use Boost.Multiprecision's high-precision floating-point types and constants, we need some includes:
*/
#include <boost/math/constants/constants.hpp>
// using boost::math::constants::pi;
#include <boost/multiprecision/cpp_dec_float.hpp>
// using boost::multiprecision::cpp_dec_float
#include <boost/multiprecision/cpp_bin_float.hpp> // for
// using boost::multiprecision::cpp_bin_float and
// using boost::multiprecision::cpp_bin_float_50;
// using boost::multiprecision::cpp_bin_float_quad;
#include <boost/array.hpp> // or <array> for std::array
#include <iostream>
#include <limits>
@@ -43,12 +57,13 @@ To use these floating-point types and constants, we need some includes:
#include <iterator>
#include <fstream>
/*`Define a text string which is a C++ comment with the program licence, copyright etc.
You could of course, tailor this to your needs, including your copyright claim.
There are versions of `array` provided by Boost.Array in `boost::array` or
the C++11 std::array, but since not all platforms provide C++11 support,
this program provides the Boost version as fallback.
/*`First, this example defines a prolog text string which is a C++ comment with the program licence, copyright etc.
(You would of course, tailor this to your needs, including *your* copyright claim).
This will appear at the top of the written header file `sines.hpp`.
*/
//] [fft_sines_table_example_1]
static const char* prolog =
{
"// Use, modification and distribution are subject to the\n"
@@ -56,27 +71,23 @@ static const char* prolog =
"// (See accompanying file LICENSE_1_0.txt\n"
"// or copy at ""http://www.boost.org/LICENSE_1_0.txt)\n\n"
"// Copyright ???? 2013.\n\n"
"// Use boost/array if std::array (C++11 feature) is not available.\n"
"#ifdef BOOST_NO_CXX11_HDR_ARRAY\n"
"#include <boost/array/array.hpp>\n"
"#else\n"
"#include <array>\n"
"#endif\n\n"
"// Copyright A N Other, 2019.\n\n"
};
//[fft_sines_table_example_2
using boost::multiprecision::cpp_dec_float_50;
using boost::multiprecision::cpp_bin_float_50;
using boost::math::constants::pi;
//] [fft_sines_table_example_2]
// VS 2010 (wrongly) requires these at file scope, not local scope in `main`.
// This program also requires `-std=c++11` option to compile using Clang and GCC.
int main()
{
/*`One often needs to compute tables of numbers in mathematical software.
A fast Fourier transform (FFT), for example, may use a table of the values of
//[fft_sines_table_example_3
/*`A fast Fourier transform (FFT), for example, may use a table of the values of
sin(([pi]/2[super n]) in its implementation details. In order to maximize the precision in
the FFT implementation, the precision of the tabulated trigonometric values
should exceed that of the built-in floating-point type used in the FFT.
@@ -85,37 +96,37 @@ The sample below computes a table of the values of sin([pi]/2[super n])
in the range 1 <= n <= 31.
This program makes use of, among other program elements, the data type
`boost::multiprecision::cpp_dec_float_50`
`boost::multiprecision::cpp_bin_float_50`
for a precision of 50 decimal digits from Boost.Multiprecision,
the value of constant [pi] retrieved from Boost.Math,
guaranteed to be initialized with the very last bit of precision for the type,
here `cpp_dec_float_50`,
here `cpp_bin_float_50`,
and a C++11 lambda function combined with `std::for_each()`.
*/
/*`define the number of values in the array.
/*`define the number of values (32) in the array of sines.
*/
std::size_t size = 32U;
cpp_dec_float_50 p = pi<cpp_dec_float_50>();
cpp_dec_float_50 p2 = boost::math::constants::pi<cpp_dec_float_50>();
//cpp_bin_float_50 p = pi<cpp_bin_float_50>();
cpp_bin_float_50 p = boost::math::constants::pi<cpp_bin_float_50>();
std::vector <cpp_dec_float_50> sin_values (size);
std::vector <cpp_bin_float_50> sin_values (size);
unsigned n = 1U;
// Generate the sine values.
std::for_each
(
sin_values.begin (),
sin_values.end (),
[&n](cpp_dec_float_50& y)
[&n](cpp_bin_float_50& y)
{
y = sin( pi<cpp_dec_float_50>() / pow(cpp_dec_float_50 (2), n));
y = sin( pi<cpp_bin_float_50>() / pow(cpp_bin_float_50 (2), n));
++n;
}
);
/*`Define the floating-point type for the generated file, either built-in
`double, `float, or `long double`, or a user defined type like `cpp_dec_float_50`.
`double, `float, or `long double`, or a user defined type like `cpp_bin_float_50`.
*/
std::string fp_type = "double";
@@ -125,26 +136,30 @@ std::cout << "Generating an `std::array` or `boost::array` for floating-point ty
/*`By default, output would only show the standard 6 decimal digits,
so set precision to show enough significant digits for the chosen floating-point type.
For `cpp_dec_float_50` is 50. (50 decimal digits should be ample for most applications).
For `cpp_bin_float_50` is 50. (50 decimal digits should be ample for most applications).
*/
std::streamsize precision = std::numeric_limits<cpp_dec_float_50>::digits10;
std::streamsize precision = std::numeric_limits<cpp_bin_float_50>::digits10;
// std::cout.precision(std::numeric_limits<cpp_dec_float_50>::digits10);
std::cout << precision << " decimal digits precision. " << std::endl;
std::cout << "Sines table precision is " << precision << " decimal digits. " << std::endl;
/*`Of course, one could also choose less, for example, 36 would be sufficient
/*`Of course, one could also choose a lower precision for the table values, for example,
`std::streamsize precision = std::numeric_limits<cpp_bin_float_quad>::max_digits10;`
128-bit 'quad' precision of 36 decimal digits would be sufficient
for the most precise current `long double` implementations using 128-bit.
In general, it should be a couple of decimal digits more (guard digits) than
`std::numeric_limits<RealType>::max_digits10` for the target system floating-point type.
If the implementation does not provide `max_digits10`, the the Kahan formula
`std::numeric_limits<RealType>::digits * 3010/10000 + 2` can be used instead.
(If the implementation does not provide `max_digits10`, the the Kahan formula
`std::numeric_limits<RealType>::digits * 3010/10000 + 2` can be used instead).
The compiler will read these values as decimal digits strings and
use the nearest representation for the floating-point type.
Now output all the sine table, to a file of your chosen name.
*/
const char sines_name[] = "sines.hpp"; // In same directory as .exe
const char sines_name[] = "sines.hpp"; // Assuming in same directory as .exe
std::ofstream fout(sines_name, std::ios_base::out); // Creates if no file exists,
// & uses default overwrite/ ios::replace.
@@ -154,19 +169,18 @@ Now output all the sine table, to a file of your chosen name.
return EXIT_FAILURE;
}
else
{
{ // Write prolog etc as a C++ comment.
std::cout << "Open file " << sines_name << " for output OK." << std::endl;
fout << prolog << "// Table of " << sin_values.size() << " values with "
fout << prolog
<< "// Table of " << sin_values.size() << " values with "
<< precision << " decimal digits precision,\n"
"// generated by program fft_sines_table.cpp.\n" << std::endl;
fout <<
"#ifdef BOOST_NO_CXX11_HDR_ARRAY""\n"
" static const boost::array<double, " << size << "> sines =\n"
"#else""\n"
" static const std::array<double, " << size << "> sines =\n"
"#endif""\n"
"{{\n"; // 2nd { needed for some GCC compiler versions.
fout << "#include <array> // std::array" << std::endl;
// Write the table of sines as a C++ array.
fout << "\nstatic const std::array<double, " << size << "> sines =\n"
"{{\n"; // 2nd { needed for some old GCC compiler versions.
fout.precision(precision);
for (unsigned int i = 0U; ;)
@@ -174,7 +188,7 @@ Now output all the sine table, to a file of your chosen name.
fout << " " << sin_values[i];
if (i == sin_values.size()-1)
{ // next is last value.
fout << "\n}};\n"; // 2nd } needed for some earlier GCC compiler versions.
fout << "\n}}; // array sines\n"; // 2nd } needed for some old GCC compiler versions.
break;
}
else
@@ -182,14 +196,14 @@ Now output all the sine table, to a file of your chosen name.
fout << ",\n";
i++;
}
}
} // for
fout.close();
std::cout << "Close file " << sines_name << " for output OK." << std::endl;
std::cout << "Closed file " << sines_name << " for output." << std::endl;
}
//`The output file generated can be seen at [@../../example/sines.hpp]
//] [/fft_sines_table_example_1]
//] [/fft_sines_table_example_3]
return EXIT_SUCCESS;

25
example/sines.hpp
View File

@@ -3,23 +3,14 @@
// (See accompanying file LICENSE_1_0.txt
// or copy at http://www.boost.org/LICENSE_1_0.txt)
// Copyright ???? 2013.
// Use boost/array if std::array (C++11 feature) is not available.
#ifdef BOOST_NO_CXX11_HDR_ARRAY
#include <boost/array/array.hpp>
#else
#include <array>
#endif
// Copyright A N Other, 2019.
// Table of 32 values with 50 decimal digits precision,
// generated by program fft_sines_table.cpp.
#ifdef BOOST_NO_CXX11_HDR_ARRAY
static const boost::array<double, 32> sines =
#else
static const std::array<double, 32> sines =
#endif
#include <array> // std::array
static const std::array<double, 32> sines =
{{
1,
0.70710678118654752440084436210484903928483593768847,
@@ -27,7 +18,7 @@
0.19509032201612826784828486847702224092769161775195,
0.098017140329560601994195563888641845861136673167501,
0.049067674327418014254954976942682658314745363025753,
0.024541228522912288031734529459282925065466119239451,
0.024541228522912288031734529459282925065466119239452,
0.012271538285719926079408261951003212140372319591769,
0.0061358846491544753596402345903725809170578863173913,
0.003067956762965976270145365490919842518944610213452,
@@ -44,13 +35,13 @@
1.4980281131690112288542788461553611206917585861527e-06,
7.4901405658471572113049856673065563715595930217207e-07,
3.7450702829238412390316917908463317739740476297248e-07,
1.8725351414619534486882457659356361712045272098287e-07,
1.8725351414619534486882457659356361712045272098286e-07,
9.3626757073098082799067286680885620193236507169473e-08,
4.681337853654909269511551813854009695950362701667e-08,
2.3406689268274552759505493419034844037886207223779e-08,
1.1703344634137277181246213503238103798093456639976e-08,
5.8516723170686386908097901008341396943900085051757e-09,
5.8516723170686386908097901008341396943900085051756e-09,
2.9258361585343193579282304690689559020175857150074e-09,
1.4629180792671596805295321618659637103742615227834e-09,
7.3145903963357984046044319684941757518633453150407e-10
}};
}}; // array sines