@@ -17,11 +17,8 @@ using std::swap;
17
17
18
18
using std::size_t ;
19
19
20
- using c64 = std::complex<double >;
21
- template <typename T>
22
- constexpr T pi () {
23
- return 3.14159265358979323846264338327950288419716 ;
24
- }
20
+ using complex = std::complex<double >;
21
+ static const double pi = 3.14159265358979323846264338327950288419716 ;
25
22
26
23
// `cooley_tukey` does the cooley-tukey algorithm, recursively
27
24
template <typename Iter>
@@ -30,12 +27,12 @@ void cooley_tukey(Iter first, Iter last) {
30
27
if (size >= 2 ) {
31
28
// split the range, with even indices going in the first half,
32
29
// and odd indices going in the last half.
33
- auto temp = std::vector<c64 >(size / 2 );
34
- for (size_t i = 0 ; i < size / 2 ; ++i) {
30
+ auto temp = std::vector<complex >(size / 2 );
31
+ for (int i = 0 ; i < size / 2 ; ++i) {
35
32
temp[i] = first[i * 2 + 1 ];
36
33
first[i] = first[i * 2 ];
37
34
}
38
- for (size_t i = 0 ; i < size / 2 ; ++i) {
35
+ for (int i = 0 ; i < size / 2 ; ++i) {
39
36
first[i + size / 2 ] = temp[i];
40
37
}
41
38
@@ -45,8 +42,8 @@ void cooley_tukey(Iter first, Iter last) {
45
42
cooley_tukey (split, last);
46
43
47
44
// now combine each of those halves with the butterflies
48
- for (size_t k = 0 ; k < size / 2 ; ++k) {
49
- auto w = std::exp (c64 (0 , -2.0 * pi< double >() * k / size));
45
+ for (int k = 0 ; k < size / 2 ; ++k) {
46
+ auto w = std::exp (complex (0 , -2.0 * pi * k / size));
50
47
51
48
auto & bottom = first[k];
52
49
auto & top = first[k + size / 2 ];
@@ -83,11 +80,11 @@ void iterative_cooley_tukey(Iter first, Iter last) {
83
80
84
81
// perform the butterfly on the range
85
82
auto size = last - first;
86
- for (size_t stride = 2 ; stride <= size; stride *= 2 ) {
87
- auto w = exp (c64 (0 , -2.0 * pi< double >() / stride));
88
- for (size_t j = 0 ; j < size; j += stride) {
89
- auto v = c64 (1.0 );
90
- for (size_t k = 0 ; k < stride / 2 ; k++) {
83
+ for (int stride = 2 ; stride <= size; stride *= 2 ) {
84
+ auto w = exp (complex (0 , -2.0 * pi / stride));
85
+ for (int j = 0 ; j < size; j += stride) {
86
+ auto v = complex (1.0 );
87
+ for (int k = 0 ; k < stride / 2 ; k++) {
91
88
first[k + j + stride / 2 ] =
92
89
first[k + j] - v * first[k + j + stride / 2 ];
93
90
first[k + j] -= (first[k + j + stride / 2 ] - first[k + j]);
@@ -103,7 +100,7 @@ int main() {
103
100
std::mt19937 rng (random_device ());
104
101
std::uniform_real_distribution<double > distribution (0.0 , 1.0 );
105
102
106
- std::array<c64 , 64 > initial;
103
+ std::array<complex , 64 > initial;
107
104
std::generate (
108
105
begin (initial), end (initial), [&] { return distribution (rng); });
109
106
@@ -117,7 +114,7 @@ int main() {
117
114
// Check if the arrays are approximately equivalent
118
115
std::cout << std::right << std::setw (16 ) << " idx" << std::setw (16 ) << " rec"
119
116
<< std::setw (16 ) << " it" << std::setw (16 ) << " subtracted" << ' \n ' ;
120
- for (int i = 0 ; i < initial.size (); ++i) {
117
+ for (size_t i = 0 ; i < initial.size (); ++i) {
121
118
auto rec = recursive[i];
122
119
auto it = iterative[i];
123
120
std::cout << std::setw (16 ) << i << std::setw (16 ) << std::abs (rec)
0 commit comments