Skip to content

Adding missing examples for Tree Traversal and FFT in C++. #114

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
wants to merge 7 commits into from
Closed

Adding missing examples for Tree Traversal and FFT in C++. #114

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
c354208
Adding dft to fft.cpp and improving the bit reversal function
Gathros May 16, 2018
513084f
Adding the two missing examples to tree_example.cpp
Gathros May 16, 2018
93a0e5e
work on modernizing Gathros' updates to fft
ubsan May 16, 2018
f49639e
minor fix in tree_example.cpp
ubsan May 16, 2018
a157047
Merge pull request #3 from ubsan/gathros-cppmissingexamples
Gathros May 20, 2018
13cb9ce
Adding comment to fft.cpp
Gathros May 20, 2018
46de9c1
updating fft.cpp
Gathros May 20, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
136 changes: 98 additions & 38 deletions chapters/FFT/code/c++/fft.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
#include <array>
#include <complex>
#include <cstdint>
#include <numeric>
#include <vector>

// These headers are for presentation not for the algorithm.
Expand All @@ -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 <typename Iter>
void cooley_tukey(Iter first, Iter last) {
auto size = last - first;
if (size >= 2) {
void recursive(Iter first, Iter last) {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As I said in another comment, I'd still argue for valarray here, it makes the implementation a lot more readable:

void fft(std::valarray<c64>& v) {
    auto size = v.size();
    if (size < 2) {
        return;
    }

    auto even = std::valarray<c64>(v[std::slice(0, size / 2, 2)]);
    auto odd  = std::valarray<c64>(v[std::slice(1, size / 2, 2)]);

    fft(even);
    fft(odd);

    for (size_t k = 0; k < size / 2; k++) {
        auto w = std::polar(1.0, -2.0 * pi<double> * k / size) * odd[k];
        v[k] = even[k] + w;
        v[k + size / 2] = even[k] - w;
    }
}

Copy link

@strega-nil strega-nil May 20, 2018 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It seems like that's more of an issue with the current implementation; we could do

template <typename Iter>
void recursive(Iter first, Iter last) {
  std::size_t size = last - first;
  if (size < 2) {
    return;
  } else if (size % 2 != 0) {
    std::cerr << "The range must be a power of two size\n";
    std::abort();
  }

  auto odds = std::vector<c64>(size / 2);
  auto evens = std::vector<c64>(size / 2);

  for (std::size_t i = 0; i < size / 2; ++i) {
    odds[i] = first[i * 2 + 1];
    evens[i] = first[i * 2];
  }

  // recurse the splits and butterflies in each half of the range
  recursive(begin(odds), end(odds));
  recursive(begin(evens), end(evens));

  // now combine each of those halves with the butterflies
  for (std::size_t k = 0; k < size / 2; ++k) {
    auto w = std::polar(1.0, -2.0 * pi<double>() * k / size) * odds[k];
    first[k] = evens[k] + w;
    first[k + size / 2] = evens[k] - w;
  }
}

// versus

void recursive(std::valarray<c64>& v) {
  auto size = v.size();
  if (size < 2) {
    return;
  } else if (size % 2 != 0) {
    std::cerr << "The range must be a power of two size\n";
    std::abort();
  }

  auto even = std::valarray<c64>(v[std::slice(0, size / 2, 2)]);
  auto odd = std::valarray<c64>(v[std::slice(1, size / 2, 2)]);

  recursive(even);
  recursive(odd);

  for (size_t k = 0; k < size / 2; k++) {
    auto w = std::polar(1.0, -2.0 * pi<double>() * k / size) * odd[k];
    v[k] = even[k] + w;
    v[k + size / 2] = even[k] - w;
  }
}

which is only 4 lines of useful extra code; I don't think the latter is significantly more readable, and I actually prefer oddsₖ = first₂ₖ₊₁ , evensₖ = first₂ₖ

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this is even nicer; moving out the split into it's own function, to make it really clear what's happening:

template <typename In, typename Out0, typename Out1>
void split(In initial, In const last, Out0 out0, Out1 out1) {
  if ((last - initial) % 2 != 0) {
    std::cerr << "The range must be a power of two size\n";
    std::abort();
  }

  while (initial != last) {
    *out0 = *initial;
    ++out0;
    ++initial;
    *out1 = *initial;
    ++out1;
    ++initial;
  }
}

// `recursive` does the cooley-tukey algorithm, recursively
template <typename Iter>
void recursive(Iter const first, Iter const last) {
  std::size_t size = last - first;
  if (size < 2) {
    return;
  }

  auto evens = std::vector<c64>(size / 2);
  auto odds = std::vector<c64>(size / 2);
  split(first, last, begin(evens), begin(odds));

  // recurse the splits and butterflies in each half of the range
  recursive(begin(odds), end(odds));
  recursive(begin(evens), end(evens));

  // now combine each of those halves with the butterflies
  for (std::size_t k = 0; k < size / 2; ++k) {
    auto w = std::polar(1.0, -2.0 * pi<double>() * k / size) * odds[k];
    first[k] = evens[k] + w;
    first[k + size / 2] = evens[k] - w;
  }
}

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I disagree with the separate split version being nicer, see my comment in the main thread. We're not really writing a standard library grade implementation and split just obfuscates a very significant step in the algorithm. I do like the proposed odds/evens version however, although I would still go with the iterator-free version with a proper return value:

void fft_calculate_recursive(std::vector<c64>& values) {
    std::size_t size = values.size();
    if (size < 2) {
        return;
    } else if (size % 2 != 0) {
        std::cerr << "The range must be a power of two size\n";
        std::abort();
    }

    auto odds = std::vector<c64>(size / 2);
    auto evens = std::vector<c64>(size / 2);

    for (std::size_t i = 0; i < size / 2; ++i) {
        odds[i] = values[i * 2 + 1];
        evens[i] = values[i * 2];
    }

    // recurse the splits and butterflies in each half of the range
    fft_calculate_recursive(odds);
    fft_calculate_recursive(evens);

    // now combine each of those halves with the butterflies
    for (std::size_t k = 0; k < size / 2; ++k) {
        auto w = std::polar(1.0, -2.0 * pi<double>() * k / size) * odds[k];
        values[k] = evens[k] + w;
        values[k + size / 2] = evens[k] - w;
    }
}

std::vector<c64> recursive(const std::vector<c64>& values) {
  std::vector<c64> result(values.begin(), values.end());
  fft_calculate_recursive(result);
  return result;
}

I think that teaching people to return actual results rather than using out parameters is way more important than any perceived genericity one gets from iterators (I'm not arguing that iterators aren't more general - I'm arguing that it isn't needed here).

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<c64>(size / 2);
for (size_t i = 0; i < size / 2; ++i) {
auto temp = std::vector<c64>(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<double>() * 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 <typename Iter>
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::size_t>(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 <typename Iter>
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<double>() / 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];
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This indentation looks incorrect

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

what is incorrect about it?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh sorry, I actually misread that line

first[k + j] -= (first[k + j + stride / 2] - first[k + j]);
Expand All @@ -97,6 +111,41 @@ void iterative_cooley_tukey(Iter first, Iter last) {
}
}

template <typename It, typename F>
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 <typename Iter>
void discrete(Iter const first, Iter const last) {
using namespace std::literals;

auto const original = std::vector<c64>(first, last);
auto const size = original.size();

auto const i2pi = -2.0i * pi<double>();

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));
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we need to cast to double here?

});
}

}
}

} // namespace ct

int main() {
// initalize the FFT inputs
std::random_device random_device;
Expand All @@ -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';
}
}
32 changes: 28 additions & 4 deletions chapters/tree_traversal/code/c++/tree_example.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -67,16 +91,16 @@ node create_tree(size_t num_row, size_t num_child) {

std::vector<node> 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};
}

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);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Commented out code.

//bfs_queue(root);
}