diff --git a/chapters/FFT/code/c++/fft.cpp b/chapters/FFT/code/c++/fft.cpp index 56e2554bb..632776e10 100644 --- a/chapters/FFT/code/c++/fft.cpp +++ b/chapters/FFT/code/c++/fft.cpp @@ -4,6 +4,7 @@ #include #include #include +#include #include // These headers are for presentation not for the algorithm. @@ -23,71 +24,84 @@ constexpr T pi() { return 3.14159265358979323846264338327950288419716; } -// `cooley_tukey` does the cooley-tukey algorithm, recursively +namespace ct { + +// `recursive` does the cooley-tukey algorithm, recursively template -void cooley_tukey(Iter first, Iter last) { - auto size = last - first; - if (size >= 2) { +void recursive(Iter first, Iter last) { + std::size_t size = last - first; + auto half_size = size / 2; + if (half_size >= 1) { // 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(half_size); + for (std::size_t i = 0; i < half_size; ++i) { temp[i] = first[i * 2 + 1]; first[i] = first[i * 2]; } - for (size_t i = 0; i < size / 2; ++i) { + for (std::size_t i = 0; i < half_size; ++i) { first[i + size / 2] = temp[i]; } // recurse the splits and butterflies in each half of the range - auto split = first + size / 2; - cooley_tukey(first, split); - cooley_tukey(split, last); + auto split = first + half_size; + recursive(first, split); + recursive(split, last); // now combine each of those halves with the butterflies - for (size_t k = 0; k < size / 2; ++k) { + for (std::size_t k = 0; k < half_size; ++k) { auto w = std::exp(c64(0, -2.0 * pi() * k / size)); auto& bottom = first[k]; - auto& top = first[k + size / 2]; + auto& top = first[k + half_size]; top = bottom - w * top; bottom -= top - bottom; } } } -// note: (last - first) must be less than 2**32 - 1 template void sort_by_bit_reverse(Iter first, Iter last) { // sorts the range [first, last) in bit-reversed order, // by the method suggested by the FFT - auto size = last - first; - - for (std::uint32_t i = 0; i < size; ++i) { - auto b = i; - b = (((b & 0xaaaaaaaa) >> 1) | ((b & 0x55555555) << 1)); - b = (((b & 0xcccccccc) >> 2) | ((b & 0x33333333) << 2)); - b = (((b & 0xf0f0f0f0) >> 4) | ((b & 0x0f0f0f0f) << 4)); - b = (((b & 0xff00ff00) >> 8) | ((b & 0x00ff00ff) << 8)); - b = ((b >> 16) | (b << 16)) >> (32 - std::uint32_t(log2(size))); - if (b > i) { - swap(first[b], first[i]); + + std::size_t size = last - first; + auto log = static_cast(std::log2(size)); + auto count = log - 1; + auto bitmask = (1 << log) - 1; + + for (std::size_t i = 0; i < size; ++i) { + auto n = i; + auto a = i; + auto current_count = count; + + // this left shifts i and right shifts the shifted bit to n by bitwise or + n >>= 1; + while (n > 0) { + a = (a << 1) | (n & 1); + --current_count; + n >>= 1; + } + n = (a << current_count) & bitmask; + + if (n > i) { + swap(first[i], first[n]); } } } -// `iterative_cooley_tukey` does the cooley-tukey algorithm iteratively +// `iterative` does the cooley-tukey algorithm iteratively template -void iterative_cooley_tukey(Iter first, Iter last) { +void iterative(Iter first, Iter last) { sort_by_bit_reverse(first, last); // perform the butterfly on the range - auto size = last - first; - for (size_t stride = 2; stride <= size; stride *= 2) { + std::size_t size = last - first; + for (std::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) { + for (std::size_t j = 0; j < size; j += stride) { auto v = c64(1.0); - for (size_t k = 0; k < stride / 2; k++) { + for (std::size_t 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]); @@ -97,6 +111,41 @@ void iterative_cooley_tukey(Iter first, Iter last) { } } +template +auto sum(It const first, It const last, F f) { + using return_type = decltype(f(0, *first)); + auto ret = return_type(0); + + std::size_t i = 0; + for (auto it = first; it < last; ++it, ++i) { + ret += f(i, *it); + } + + return ret; +} + +template +void discrete(Iter const first, Iter const last) { + using namespace std::literals; + + auto const original = std::vector(first, last); + auto const size = original.size(); + + auto const i2pi = -2.0i * pi(); + + for (std::size_t i = 0; i < size; ++i) { + for (std::size_t i = 0; i < size; ++i) { + first[i] = + sum(begin(original), end(original), [&](auto const k, auto const el) { + return el * std::exp(i2pi * double(i * k) / double(size)); + }); + } + + } +} + +} // namespace ct + int main() { // initalize the FFT inputs std::random_device random_device; @@ -109,19 +158,30 @@ int main() { auto recursive = initial; auto iterative = initial; + auto discrete = initial; // Preform an FFT on the arrays. - cooley_tukey(begin(recursive), end(recursive)); - iterative_cooley_tukey(begin(iterative), end(iterative)); + ct::recursive(begin(recursive), end(recursive)); + ct::iterative(begin(iterative), end(iterative)); + ct::discrete(begin(discrete), end(discrete)); // 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) { + std::cout << std::right << std::setw(12) << "idx"; + std::cout << std::setw(12) << "|recursive|"; + std::cout << std::setw(12) << "|iterative|"; + std::cout << std::setw(12) << "|discrete|"; + std::cout << std::setw(12) << "error" << '\n'; + for (std::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) - << std::setw(16) << std::abs(it) << std::setw(16) - << (std::abs(rec) - std::abs(it)) << '\n'; + auto disc = discrete[i]; + std::cout << std::setw(12) << i; + std::cout << std::setw(12) << std::abs(rec); + std::cout << std::setw(12) << std::abs(it); + std::cout << std::setw(12) << std::abs(disc); + + auto err = std::max( + {std::abs(rec - it), std::abs(it - disc), std::abs(disc - rec)}); + std::cout << std::setw(12) << err << '\n'; } } diff --git a/chapters/tree_traversal/code/c++/tree_example.cpp b/chapters/tree_traversal/code/c++/tree_example.cpp index 519ff374e..46f408913 100644 --- a/chapters/tree_traversal/code/c++/tree_example.cpp +++ b/chapters/tree_traversal/code/c++/tree_example.cpp @@ -26,6 +26,30 @@ void dfs_recursive(node const& n) { } } +void dfs_recursive_postorder(node const& n) { + for (auto const& child : n.children) { + dfs_recursive_postorder(child); + } + std::cout << n.value << '\n'; +} + +void dfs_recursive_inorder_btree(node const& n) { + auto size = n.children.size(); + + if (size > 2) { + std::cerr << "This is not a binary tree.\n"; + std::abort(); + } + + if (size > 0) { + dfs_recursive_inorder_btree(n.children[0]); + } + std::cout << n.value << '\n'; + if (size > 1) { + dfs_recursive_inorder_btree(n.children[1]); + } +} + // Simple non-recursive scheme for DFS void dfs_stack(node const& n) { // this stack holds pointers into n's `children` vector, @@ -67,7 +91,7 @@ node create_tree(size_t num_row, size_t num_child) { std::vector vec; std::generate_n(std::back_inserter(vec), num_child, [&] { - return create_node(num_row - 1, num_child); + return create_tree(num_row - 1, num_child); }); return node{std::move(vec), num_row}; @@ -75,8 +99,8 @@ node create_tree(size_t num_row, size_t num_child) { int main() { // Creating Tree in main - auto root = create_node(3, 3); + auto root = create_tree(3, 2); dfs_recursive(root); - dfs_stack(root); - bfs_queue(root); + //dfs_stack(root); + //bfs_queue(root); }