Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- /**
- * 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);
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement