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 <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);
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement