diff --git a/contents/cooley_tukey/code/c++/fft.cpp b/contents/cooley_tukey/code/c++/fft.cpp index 56e2554bb..5d4772f10 100644 --- a/contents/cooley_tukey/code/c++/fft.cpp +++ b/contents/cooley_tukey/code/c++/fft.cpp @@ -17,10 +17,19 @@ using std::swap; using std::size_t; -using c64 = std::complex; -template -constexpr T pi() { - return 3.14159265358979323846264338327950288419716; +using complex = std::complex; +static const double pi = 3.14159265358979323846264338327950288419716; + +template +void dft(Iter X, Iter last) { + const auto N = last - X; + std::vector tmp(N); + for (auto i = 0; i < N; ++i) { + for (auto j = 0; j < N; ++j) { + tmp[i] += X[j] * exp(complex(0, -2.0 * M_PI * i * j / N)); + } + } + std::copy(std::begin(tmp), std::end(tmp), X); } // `cooley_tukey` does the cooley-tukey algorithm, recursively @@ -30,12 +39,12 @@ void cooley_tukey(Iter first, Iter last) { if (size >= 2) { // split the range, with even indices going in the first half, // and odd indices going in the last half. - auto temp = std::vector(size / 2); - for (size_t i = 0; i < size / 2; ++i) { + auto temp = std::vector(size / 2); + for (int i = 0; i < size / 2; ++i) { temp[i] = first[i * 2 + 1]; first[i] = first[i * 2]; } - for (size_t i = 0; i < size / 2; ++i) { + for (int i = 0; i < size / 2; ++i) { first[i + size / 2] = temp[i]; } @@ -45,8 +54,8 @@ void cooley_tukey(Iter first, Iter last) { cooley_tukey(split, last); // now combine each of those halves with the butterflies - for (size_t k = 0; k < size / 2; ++k) { - auto w = std::exp(c64(0, -2.0 * pi() * k / size)); + for (int k = 0; k < size / 2; ++k) { + auto w = std::exp(complex(0, -2.0 * pi * k / size)); auto& bottom = first[k]; auto& top = first[k + size / 2]; @@ -83,11 +92,11 @@ void iterative_cooley_tukey(Iter first, Iter last) { // perform the butterfly on the range auto size = last - first; - for (size_t stride = 2; stride <= size; stride *= 2) { - auto w = exp(c64(0, -2.0 * pi() / stride)); - for (size_t j = 0; j < size; j += stride) { - auto v = c64(1.0); - for (size_t k = 0; k < stride / 2; k++) { + for (int stride = 2; stride <= size; stride *= 2) { + auto w = exp(complex(0, -2.0 * pi / stride)); + for (int j = 0; j < size; j += stride) { + auto v = complex(1.0); + for (int k = 0; k < stride / 2; k++) { first[k + j + stride / 2] = first[k + j] - v * first[k + j + stride / 2]; first[k + j] -= (first[k + j + stride / 2] - first[k + j]); @@ -103,7 +112,7 @@ int main() { std::mt19937 rng(random_device()); std::uniform_real_distribution distribution(0.0, 1.0); - std::array initial; + std::array initial; std::generate( begin(initial), end(initial), [&] { return distribution(rng); }); @@ -117,7 +126,7 @@ int main() { // Check if the arrays are approximately equivalent std::cout << std::right << std::setw(16) << "idx" << std::setw(16) << "rec" << std::setw(16) << "it" << std::setw(16) << "subtracted" << '\n'; - for (int i = 0; i < initial.size(); ++i) { + for (size_t i = 0; i < initial.size(); ++i) { auto rec = recursive[i]; auto it = iterative[i]; std::cout << std::setw(16) << i << std::setw(16) << std::abs(rec) diff --git a/contents/cooley_tukey/cooley_tukey.md b/contents/cooley_tukey/cooley_tukey.md index 41b00170a..e414e8c16 100644 --- a/contents/cooley_tukey/cooley_tukey.md +++ b/contents/cooley_tukey/cooley_tukey.md @@ -119,7 +119,7 @@ In the end, the code looks like: {% sample lang="c" %} [import:20-39, lang:"c_cpp"](code/c/fft.c) {% sample lang="cpp" %} -[import:27-57, lang:"c_cpp"](code/c++/fft.cpp) +[import:35-66, lang:"c_cpp"](code/c++/fft.cpp) {% sample lang="hs" %} [import:6-19, lang:"haskell"](code/haskell/fft.hs) {% sample lang="py" %}