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 <array>
// using std::array
#include <vector>
// using std::vector
#include <tuple>
// using std::tuple
#include <string>
// using std::stoi
#include <iostream>
// using std::cout
#include <algorithm>
// using std::transform, std::all_of
#include <functional>
// using std::bind, std::placeholders, std::multiplies, std::plus, std::minus
#include <climits>
// using CHAR_BIT


namespace {

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

  struct Vecf64
  {
    using array_type = std::array<double_t, VLEN>;
    array_type array;

    Vecf64() = default;

    explicit Vecf64(double_t s) noexcept
    { array.fill(s); }

    Vecf64& operator=(const Vecf64&) = default;

    Vecf64& operator=(double_t s) noexcept
    {
      array.fill(s);
      return *this;
    }
    double_t& operator[](int idx) noexcept
    { return array[idx]; }

    double_t operator[](int idx) const noexcept
    { return array[idx]; }

    template<class BinaryFunctor>
    Vecf64& inplace_transform(const Vecf64& o, BinaryFunctor f)
    {
      array_type::iterator first = array.begin();
      std::transform(first, array.end(), o.array.begin(), first, f);
      return *this;
    }
    template<class BinaryFunctor>
    Vecf64 transform(const Vecf64& o, BinaryFunctor f) const
    {
      Vecf64 rtrn;
      std::transform(array.begin(), array.end(), o.array.begin(),
                     rtrn.array.begin(), f);
      return rtrn;
    }
    template<class UnaryFunctor>
    Vecf64& inplace_transform(UnaryFunctor f)
    {
      array_type::iterator first = array.begin();
      std::transform(first, array.end(), first, f);
      return *this;
    }
    template<class UnaryFunctor>
    Vecf64 transform(UnaryFunctor f) const
    {
      Vecf64 rtrn;
      std::transform(array.begin(), array.end(), rtrn.array.begin(), f);
      return rtrn;
    }

    Vecf64& operator*=(const Vecf64& o) noexcept
    { return inplace_transform(o, std::multiplies<double_t>{}); }

    Vecf64 operator*(const Vecf64& o) const noexcept
    { return transform(o, std::multiplies<double_t>{}); }

    Vecf64& operator*=(double_t s) noexcept
    {
      return inplace_transform(std::bind(std::multiplies<double_t>{},
                                         std::placeholders::_1, s));
    }
    Vecf64 operator*(double_t s) const noexcept
    {
      return transform(std::bind(std::multiplies<double_t>{},
                                 std::placeholders::_1, s));
    }
    friend Vecf64 operator*(double_t left, const Vecf64& right) noexcept
    { return right * left; }

    Vecf64& operator+=(const Vecf64& o) noexcept
    { return inplace_transform(o, std::plus<double_t>{}); }

    Vecf64 operator+(const Vecf64& o) const noexcept
    { return transform(o, std::plus<double_t>{}); }

    Vecf64& operator+=(double_t s) noexcept
    {
      return inplace_transform(std::bind(std::plus<double_t>{},
                                         std::placeholders::_1, s));
    }
    Vecf64 operator+(double_t s) const noexcept
    {
      return transform(std::bind(std::plus<double_t>{},
                                 std::placeholders::_1, s));
    }
    friend Vecf64 operator+(double_t left, const Vecf64& right) noexcept
    { return right + left; }

    Vecf64& operator-=(const Vecf64& o) noexcept
    { return inplace_transform(o, std::minus<double_t>{}); }

    Vecf64 operator-(const Vecf64& o) const noexcept
    { return transform(o, std::minus<double_t>{}); }

    Vecf64& operator-=(double_t s) noexcept
    {
      return inplace_transform(std::bind(std::minus<double_t>{},
                                         std::placeholders::_1, s));
    }
    Vecf64 operator-(double_t s) const noexcept
    {
      return transform(std::bind(std::minus<double_t>{},
                                 std::placeholders::_1, s));
    }
    friend Vecf64 operator-(double_t left, const Vecf64& right) noexcept
    {
      return right.transform(std::bind(std::minus<double_t>{},
                                       left, std::placeholders::_1));
    }
    Vecf64& inplace_square() noexcept
    {
      return inplace_transform(std::bind(std::multiplies<double_t>{},
                                         std::placeholders::_1,
                                         std::placeholders::_1));
    }
    Vecf64 square() const noexcept
    {
      return transform(std::bind(std::multiplies<double_t>{},
                                 std::placeholders::_1,
                                 std::placeholders::_1));
    }
    template<class Predicate>
    bool all(Predicate p) const
    { return std::all_of(array.begin(), array.end(), p); }
  };

  struct CrCr2Chunk
  {
    Vecf64 cr, cr2;
  };

  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
    {
      auto diverged = std::bind(std::greater<double_t>{},
                                std::placeholders::_1,
                                4.);
      return (tr + ti).all(diverged);
    }
    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:
    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();
    }
  };

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