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
22 /*`[h5 Using Boost.Multiprecision to generate a high-precision array of sine coefficients for use with FFT.]
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
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.
29
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:
41 */
42 #include <boost/math/constants/constants.hpp>
43 // using boost::math::constants::pi;
44
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
51
52 #include <iostream>
53 #include <limits>
54 #include <vector>
55 #include <algorithm>
56 #include <iomanip>
57 #include <iterator>
58 #include <fstream>
59
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`.
63 */
64
65 //] [fft_sines_table_example_1]
66
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
74 "// Copyright A N Other, 2019.\n\n"
75 };
76
77 //[fft_sines_table_example_2
78
79 using boost::multiprecision::cpp_bin_float_50;
80 using boost::math::constants::pi;
81
82 //] [fft_sines_table_example_2]
83
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
main()87 int main()
88 {
89 //[fft_sines_table_example_3
90 /*`A fast Fourier transform (FFT), for example, may use a table of the values of
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
99 `boost::multiprecision::cpp_bin_float_50`
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,
103 here `cpp_bin_float_50`,
104 and a C++11 lambda function combined with `std::for_each()`.
105 */
106
107 /*`define the number of values (32) in the array of sines.
108 */
109
110 std::size_t size = 32U;
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>();
113
114 std::vector <cpp_bin_float_50> sin_values (size);
115 unsigned n = 1U;
116 // Generate the sine values.
117 std::for_each
118 (
119 sin_values.begin (),
120 sin_values.end (),
121 [&n](cpp_bin_float_50& y)
122 {
123 y = sin( pi<cpp_bin_float_50>() / pow(cpp_bin_float_50 (2), n));
124 ++n;
125 }
126 );
127
128 /*`Define the floating-point type for the generated file, either built-in
129 `double, `float, or `long double`, or a user defined type like `cpp_bin_float_50`.
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.
139 For `cpp_bin_float_50` is 50. (50 decimal digits should be ample for most applications).
140
141 */
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,
147
148 `std::streamsize precision = std::numeric_limits<cpp_bin_float_quad>::max_digits10;`
149
150 128-bit 'quad' precision of 36 decimal digits would be sufficient
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.
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).
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 */
162 const char sines_name[] = "sines.hpp"; // Assuming in same directory as .exe
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
172 { // Write prolog etc as a C++ comment.
173 std::cout << "Open file " << sines_name << " for output OK." << std::endl;
174 fout << prolog
175 << "// Table of " << sin_values.size() << " values with "
176 << precision << " decimal digits precision,\n"
177 "// generated by program fft_sines_table.cpp.\n" << std::endl;
178
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.
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.
191 fout << "\n}}; // array sines\n"; // 2nd } needed for some old GCC compiler versions.
192 break;
193 }
194 else
195 {
196 fout << ",\n";
197 i++;
198 }
199 } // for
200
201 fout.close();
202 std::cout << "Closed file " << sines_name << " for output." << std::endl;
203 }
204 //`The output file generated can be seen at [@../../example/sines.hpp]
205
206 //] [/fft_sines_table_example_3]
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]
271