]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // Use, modification and distribution are subject to the |
2 | // Boost Software License, Version 1.0. | |
3 | // (See accompanying file LICENSE_1_0.txt | |
4 | // or copy at http://www.boost.org/LICENSE_1_0.txt) | |
5 | ||
6 | // Copyright Paul A. Bristow 2013. | |
7 | // Copyright Christopher Kormanyos 2012, 2013. | |
8 | // Copyright John Maddock 2013. | |
9 | ||
10 | // This file is written to be included from a Quickbook .qbk document. | |
11 | // It can be compiled by the C++ compiler, and run. Any output can | |
12 | // also be added here as comment or included or pasted in elsewhere. | |
13 | // Caution: this file contains Quickbook markup as well as code | |
14 | // and comments: don't change any of the special comment markups! | |
15 | ||
16 | #ifdef _MSC_VER | |
17 | # pragma warning (disable : 4996) // -D_SCL_SECURE_NO_WARNINGS. | |
18 | #endif | |
19 | ||
20 | //[fft_sines_table_example_1 | |
21 | ||
92f5a8d4 | 22 | /*`[h5 Using Boost.Multiprecision to generate a high-precision array of sine coefficents for use with FFT.] |
7c673cae FG |
23 | |
24 | The Boost.Multiprecision library can be used for computations requiring precision | |
25 | exceeding that of standard built-in types such as `float`, `double` | |
26 | and `long double`. For extended-precision calculations, Boost.Multiprecision | |
92f5a8d4 TL |
27 | supplies a template data type called `cpp_bin_float`. The number of decimal |
28 | digits of precision is fixed at compile-time via a template parameter. | |
7c673cae | 29 | |
92f5a8d4 TL |
30 | One often needs to compute tables of numbers in mathematical software. |
31 | To avoid the | |
32 | [@https://en.wikipedia.org/wiki/Rounding#Table-maker's_dilemma Table-maker's dilemma] | |
33 | it is necessary to use a higher precision type to compute the table values so that they have | |
34 | the nearest representable bit-pattern for the type, say `double`, of the table value. | |
35 | ||
36 | This example is a program `fft_since_table.cpp` that writes a header file `sines.hpp` | |
37 | containing an array of sine coefficients for use with a Fast Fourier Transform (FFT), | |
38 | that can be included by the FFT program. | |
39 | ||
40 | To use Boost.Multiprecision's high-precision floating-point types and constants, we need some includes: | |
7c673cae FG |
41 | */ |
42 | #include <boost/math/constants/constants.hpp> | |
43 | // using boost::math::constants::pi; | |
44 | ||
92f5a8d4 TL |
45 | #include <boost/multiprecision/cpp_bin_float.hpp> // for |
46 | // using boost::multiprecision::cpp_bin_float and | |
47 | // using boost::multiprecision::cpp_bin_float_50; | |
48 | // using boost::multiprecision::cpp_bin_float_quad; | |
49 | ||
50 | #include <boost/array.hpp> // or <array> for std::array | |
7c673cae FG |
51 | |
52 | #include <iostream> | |
53 | #include <limits> | |
54 | #include <vector> | |
55 | #include <algorithm> | |
56 | #include <iomanip> | |
57 | #include <iterator> | |
58 | #include <fstream> | |
59 | ||
92f5a8d4 TL |
60 | /*`First, this example defines a prolog text string which is a C++ comment with the program licence, copyright etc. |
61 | (You would of course, tailor this to your needs, including *your* copyright claim). | |
62 | This will appear at the top of the written header file `sines.hpp`. | |
7c673cae | 63 | */ |
92f5a8d4 TL |
64 | |
65 | //] [fft_sines_table_example_1] | |
66 | ||
7c673cae FG |
67 | static const char* prolog = |
68 | { | |
69 | "// Use, modification and distribution are subject to the\n" | |
70 | "// Boost Software License, Version 1.0.\n" | |
71 | "// (See accompanying file LICENSE_1_0.txt\n" | |
72 | "// or copy at ""http://www.boost.org/LICENSE_1_0.txt)\n\n" | |
73 | ||
92f5a8d4 | 74 | "// Copyright A N Other, 2019.\n\n" |
7c673cae FG |
75 | }; |
76 | ||
92f5a8d4 | 77 | //[fft_sines_table_example_2 |
7c673cae | 78 | |
92f5a8d4 | 79 | using boost::multiprecision::cpp_bin_float_50; |
7c673cae | 80 | using boost::math::constants::pi; |
92f5a8d4 TL |
81 | |
82 | //] [fft_sines_table_example_2] | |
83 | ||
7c673cae FG |
84 | // VS 2010 (wrongly) requires these at file scope, not local scope in `main`. |
85 | // This program also requires `-std=c++11` option to compile using Clang and GCC. | |
86 | ||
87 | int main() | |
88 | { | |
92f5a8d4 TL |
89 | //[fft_sines_table_example_3 |
90 | /*`A fast Fourier transform (FFT), for example, may use a table of the values of | |
7c673cae FG |
91 | sin(([pi]/2[super n]) in its implementation details. In order to maximize the precision in |
92 | the FFT implementation, the precision of the tabulated trigonometric values | |
93 | should exceed that of the built-in floating-point type used in the FFT. | |
94 | ||
95 | The sample below computes a table of the values of sin([pi]/2[super n]) | |
96 | in the range 1 <= n <= 31. | |
97 | ||
98 | This program makes use of, among other program elements, the data type | |
92f5a8d4 | 99 | `boost::multiprecision::cpp_bin_float_50` |
7c673cae FG |
100 | for a precision of 50 decimal digits from Boost.Multiprecision, |
101 | the value of constant [pi] retrieved from Boost.Math, | |
102 | guaranteed to be initialized with the very last bit of precision for the type, | |
92f5a8d4 | 103 | here `cpp_bin_float_50`, |
7c673cae FG |
104 | and a C++11 lambda function combined with `std::for_each()`. |
105 | */ | |
106 | ||
92f5a8d4 | 107 | /*`define the number of values (32) in the array of sines. |
7c673cae FG |
108 | */ |
109 | ||
110 | std::size_t size = 32U; | |
92f5a8d4 TL |
111 | //cpp_bin_float_50 p = pi<cpp_bin_float_50>(); |
112 | cpp_bin_float_50 p = boost::math::constants::pi<cpp_bin_float_50>(); | |
7c673cae | 113 | |
92f5a8d4 | 114 | std::vector <cpp_bin_float_50> sin_values (size); |
7c673cae FG |
115 | unsigned n = 1U; |
116 | // Generate the sine values. | |
117 | std::for_each | |
118 | ( | |
119 | sin_values.begin (), | |
120 | sin_values.end (), | |
92f5a8d4 | 121 | [&n](cpp_bin_float_50& y) |
7c673cae | 122 | { |
92f5a8d4 | 123 | y = sin( pi<cpp_bin_float_50>() / pow(cpp_bin_float_50 (2), n)); |
7c673cae FG |
124 | ++n; |
125 | } | |
126 | ); | |
127 | ||
128 | /*`Define the floating-point type for the generated file, either built-in | |
92f5a8d4 | 129 | `double, `float, or `long double`, or a user defined type like `cpp_bin_float_50`. |
7c673cae FG |
130 | */ |
131 | ||
132 | std::string fp_type = "double"; | |
133 | ||
134 | std::cout << "Generating an `std::array` or `boost::array` for floating-point type: " | |
135 | << fp_type << ". " << std::endl; | |
136 | ||
137 | /*`By default, output would only show the standard 6 decimal digits, | |
138 | so set precision to show enough significant digits for the chosen floating-point type. | |
92f5a8d4 TL |
139 | For `cpp_bin_float_50` is 50. (50 decimal digits should be ample for most applications). |
140 | ||
7c673cae | 141 | */ |
92f5a8d4 TL |
142 | std::streamsize precision = std::numeric_limits<cpp_bin_float_50>::digits10; |
143 | ||
144 | std::cout << "Sines table precision is " << precision << " decimal digits. " << std::endl; | |
145 | ||
146 | /*`Of course, one could also choose a lower precision for the table values, for example, | |
7c673cae | 147 | |
92f5a8d4 | 148 | `std::streamsize precision = std::numeric_limits<cpp_bin_float_quad>::max_digits10;` |
7c673cae | 149 | |
92f5a8d4 | 150 | 128-bit 'quad' precision of 36 decimal digits would be sufficient |
7c673cae FG |
151 | for the most precise current `long double` implementations using 128-bit. |
152 | In general, it should be a couple of decimal digits more (guard digits) than | |
153 | `std::numeric_limits<RealType>::max_digits10` for the target system floating-point type. | |
92f5a8d4 TL |
154 | (If the implementation does not provide `max_digits10`, the the Kahan formula |
155 | `std::numeric_limits<RealType>::digits * 3010/10000 + 2` can be used instead). | |
7c673cae FG |
156 | |
157 | The compiler will read these values as decimal digits strings and | |
158 | use the nearest representation for the floating-point type. | |
159 | ||
160 | Now output all the sine table, to a file of your chosen name. | |
161 | */ | |
92f5a8d4 | 162 | const char sines_name[] = "sines.hpp"; // Assuming in same directory as .exe |
7c673cae FG |
163 | |
164 | std::ofstream fout(sines_name, std::ios_base::out); // Creates if no file exists, | |
165 | // & uses default overwrite/ ios::replace. | |
166 | if (fout.is_open() == false) | |
167 | { // failed to open OK! | |
168 | std::cout << "Open file " << sines_name << " failed!" << std::endl; | |
169 | return EXIT_FAILURE; | |
170 | } | |
171 | else | |
92f5a8d4 | 172 | { // Write prolog etc as a C++ comment. |
7c673cae | 173 | std::cout << "Open file " << sines_name << " for output OK." << std::endl; |
92f5a8d4 TL |
174 | fout << prolog |
175 | << "// Table of " << sin_values.size() << " values with " | |
7c673cae FG |
176 | << precision << " decimal digits precision,\n" |
177 | "// generated by program fft_sines_table.cpp.\n" << std::endl; | |
178 | ||
92f5a8d4 TL |
179 | fout << "#include <array> // std::array" << std::endl; |
180 | ||
181 | // Write the table of sines as a C++ array. | |
182 | fout << "\nstatic const std::array<double, " << size << "> sines =\n" | |
183 | "{{\n"; // 2nd { needed for some old GCC compiler versions. | |
7c673cae FG |
184 | fout.precision(precision); |
185 | ||
186 | for (unsigned int i = 0U; ;) | |
187 | { | |
188 | fout << " " << sin_values[i]; | |
189 | if (i == sin_values.size()-1) | |
190 | { // next is last value. | |
92f5a8d4 | 191 | fout << "\n}}; // array sines\n"; // 2nd } needed for some old GCC compiler versions. |
7c673cae FG |
192 | break; |
193 | } | |
194 | else | |
195 | { | |
196 | fout << ",\n"; | |
197 | i++; | |
198 | } | |
92f5a8d4 | 199 | } // for |
7c673cae FG |
200 | |
201 | fout.close(); | |
92f5a8d4 | 202 | std::cout << "Closed file " << sines_name << " for output." << std::endl; |
7c673cae FG |
203 | } |
204 | //`The output file generated can be seen at [@../../example/sines.hpp] | |
92f5a8d4 TL |
205 | |
206 | //] [/fft_sines_table_example_3] | |
7c673cae FG |
207 | |
208 | return EXIT_SUCCESS; | |
209 | ||
210 | } // int main() | |
211 | ||
212 | /* | |
213 | //[fft_sines_table_example_output | |
214 | ||
215 | The printed table is: | |
216 | ||
217 | 1 | |
218 | 0.70710678118654752440084436210484903928483593768847 | |
219 | 0.38268343236508977172845998403039886676134456248563 | |
220 | 0.19509032201612826784828486847702224092769161775195 | |
221 | 0.098017140329560601994195563888641845861136673167501 | |
222 | 0.049067674327418014254954976942682658314745363025753 | |
223 | 0.024541228522912288031734529459282925065466119239451 | |
224 | 0.012271538285719926079408261951003212140372319591769 | |
225 | 0.0061358846491544753596402345903725809170578863173913 | |
226 | 0.003067956762965976270145365490919842518944610213452 | |
227 | 0.0015339801862847656123036971502640790799548645752374 | |
228 | 0.00076699031874270452693856835794857664314091945206328 | |
229 | 0.00038349518757139558907246168118138126339502603496474 | |
230 | 0.00019174759731070330743990956198900093346887403385916 | |
231 | 9.5873799095977345870517210976476351187065612851145e-05 | |
232 | 4.7936899603066884549003990494658872746866687685767e-05 | |
233 | 2.3968449808418218729186577165021820094761474895673e-05 | |
234 | 1.1984224905069706421521561596988984804731977538387e-05 | |
235 | 5.9921124526424278428797118088908617299871778780951e-06 | |
236 | 2.9960562263346607504548128083570598118251878683408e-06 | |
237 | 1.4980281131690112288542788461553611206917585861527e-06 | |
238 | 7.4901405658471572113049856673065563715595930217207e-07 | |
239 | 3.7450702829238412390316917908463317739740476297248e-07 | |
240 | 1.8725351414619534486882457659356361712045272098287e-07 | |
241 | 9.3626757073098082799067286680885620193236507169473e-08 | |
242 | 4.681337853654909269511551813854009695950362701667e-08 | |
243 | 2.3406689268274552759505493419034844037886207223779e-08 | |
244 | 1.1703344634137277181246213503238103798093456639976e-08 | |
245 | 5.8516723170686386908097901008341396943900085051757e-09 | |
246 | 2.9258361585343193579282304690689559020175857150074e-09 | |
247 | 1.4629180792671596805295321618659637103742615227834e-09 | |
248 | */ | |
249 | ||
250 | //] [/fft_sines_table_example_output] | |
251 | ||
252 | //[fft_sines_table_example_check | |
253 | ||
254 | /*` | |
255 | The output can be copied as text and readily integrated into a given source | |
256 | code. Alternatively, the output can be written to a text or even be used | |
257 | within a self-written automatic code generator as this example. | |
258 | ||
259 | A computer algebra system can be used to verify the results obtained from | |
260 | Boost.Math and Boost.Multiprecision. For example, the __Mathematica | |
261 | computer algebra system can obtain a similar table with the command: | |
262 | ||
263 | Table[N[Sin[Pi / (2^n)], 50], {n, 1, 31, 1}] | |
264 | ||
265 | The __WolframAlpha computational knowledge engine can also be used to generate | |
266 | this table. The same command can be pasted into the compute box. | |
267 | ||
268 | */ | |
269 | ||
270 | //] [/fft_sines_table_example_check] |