Eigen Mandelbrot

/**
 * Author: F. Philipp
 *
 * Translation of
 * https://benchmarksgame.alioth.debian.org/u64q/program.php?test=mandelbrot&lang=rust&id=4
 * with minor changes
 *
 * Contains a bug resulting in slightly different results
 */

#include <Eigen/Dense>


#include <vector>
// using std::vector
#include <string>
// using std::stoi
#include <iostream>
// using std::cout
#include <algorithm>
// using std::transform
#include <climits>
// using CHAR_BIT


namespace {

  constexpr int MAX_ITER = 50;
  constexpr int VLEN = CHAR_BIT;
  using double_t = double;

  using Vecf64 = Eigen::Array<double_t, VLEN, 1>;

  template<class T>
  using vector = std::vector<T, Eigen::aligned_allocator<T> >;

  struct CrCr2Chunk
  {
    Vecf64 cr, cr2;

    EIGEN_MAKE_ALIGNED_OPERATOR_NEW
  };

  class Mandelbrot8
  {
    Vecf64 zr, zi, tr, ti, cr;
    double_t ci, ci2;

    void advance(int iterations) noexcept
    {
      for(int i = 0; i < iterations; ++i) {
        zi = 2. * zr * zi + ci;
        zr = tr - ti + cr;
        tr = zr.square();
        ti = zi.square();
      }
    }
    bool all_diverged() const noexcept
    {
      return ((tr + ti) > 4).all();
    }
    unsigned char to_byte() const noexcept
    {
      const Vecf64 trti = tr + ti;
      unsigned char accu = 0;
      for(int i = 0; i < VLEN; ++i) {
        if(trti[i] <= 4.)
          accu |= 0x80 >> i;
      }
      return accu;
    }
  public:
    EIGEN_MAKE_ALIGNED_OPERATOR_NEW

    explicit Mandelbrot8(double_t ci) noexcept
      : ci(ci),
        ci2(ci * ci)
    {}
    unsigned char operator()(const CrCr2Chunk& col_invariant) noexcept
    {
      zr = col_invariant.cr;
      zi = ci;
      tr = col_invariant.cr2;
      ti = ci2;
      cr = col_invariant.cr;
      advance(4);
      constexpr int n = MAX_ITER / 5 - 1;
      for(int i = 0; i < n; ++i) {
        if(all_diverged())
          return 0;
        advance(4);
      }
      return to_byte();
    }
  };

  vector<CrCr2Chunk> make_column_invariant(int bytes_per_line)
  {
    const int pixels_per_line = bytes_per_line * VLEN;
    const double_t inv = 2. / pixels_per_line;
    vector<CrCr2Chunk> rtrn(bytes_per_line);
#   pragma omp parallel for
    for(int col_byte = 0; col_byte < bytes_per_line; ++col_byte) {
      CrCr2Chunk& chunk = rtrn[col_byte];
      for(int bit = 0; bit < VLEN; ++bit) {
        const int col = col_byte * VLEN + bit;
        chunk.cr[bit] = col * inv - 1.5;
      }
      chunk.cr2 = chunk.cr.square();
    }
    return rtrn;
  }

  std::vector<unsigned char>
  make_bitmap(const vector<CrCr2Chunk>& column_invariants)
  {
    const int bytes_per_line = static_cast<int>(column_invariants.size());
    const int pixels_per_line = bytes_per_line * VLEN;
    const double_t inv = 2. / pixels_per_line;
    std::vector<unsigned char> bitmap(pixels_per_line * bytes_per_line);
#   pragma omp parallel for
    for(int row = 0; row < pixels_per_line; ++row) {
      const int start = row * bytes_per_line;
      const double_t ci = row * inv - 1.;
      std::transform(column_invariants.begin(), column_invariants.end(),
                     bitmap.begin() + start,
                     Mandelbrot8(ci));
    }
    return bitmap;
  }
  void print_pbm(const std::vector<unsigned char>& bitmap,
                 int pixels_per_line,
                 std::ostream& binary_output)
  {
    binary_output << "P4\n" << pixels_per_line
                  << ' ' << pixels_per_line
                  << '\n';
    binary_output.write(reinterpret_cast<const char*>(bitmap.data()),
                        bitmap.size());
    binary_output.flush();
  }
}

int main(int argc, char** argv)
{
  /* number of pixels in each dimension */
  int pixels_per_line = 200;
  if(argc > 1)
    pixels_per_line = std::stoi(argv[1]);
  const int bytes_per_line = pixels_per_line / VLEN;
  pixels_per_line = bytes_per_line * VLEN;
  const vector<CrCr2Chunk> column_invariants =
    make_column_invariant(bytes_per_line);
  const std::vector<unsigned char> bitmap = make_bitmap(column_invariants);
  print_pbm(bitmap, pixels_per_line, std::cout);
}