Zoom through the Mandelbrot set

mandelbrot.cpp

// g++ -Wall -pedantic --std=c++17 -g -O3 -fopenmp -march=native -o mandelbrot{,.cpp} -lpng -pthread
// rm -rf *.png *.gif && ./mandelbrot
// convert -delay 5 movied_*.png -loop 0 movied.gif

#include <algorithm>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <vector>

#include <png++/png.hpp>
#include <immintrin.h>

std::vector<std::uint32_t> mandelbrot_set(double center_x, double center_y, double width,
                                          std::size_t cols, std::size_t max_iter) {
    const double x_min = center_x - width / 2, x_max = center_x + width / 2,
                 y_min = center_y - width / 2, y_max = center_y + width / 2;
    const std::size_t rows = static_cast<std::size_t>(cols * (y_max - y_min) / (x_max - x_min));
    const std::size_t size = rows * cols;

    std::vector<std::uint32_t> iter_count(size, max_iter);

    const __m256d threshold = _mm256_set1_pd(4.0);
    const __m256d zero = _mm256_set1_pd(0.0);
    const __m256d one = _mm256_set1_pd(1.0);
    const __m256d two = _mm256_set1_pd(2.0);
    const __m256d kx = _mm256_set1_pd((x_max - x_min) / cols);
    const double ky = (y_max - y_min) / rows;
    const __m256d x_min_256 = _mm256_set1_pd(x_min);
    const __m256d minus_one = _mm256_set1_pd(-1);

    #pragma omp parallel for schedule(dynamic, 1)
    for (std::size_t i = 0; i < rows; ++i) {
        std::uint32_t* p_iter = &iter_count[i * cols];

        __m256d y0 = _mm256_set1_pd(y_max - ky * i);

        for (std::size_t j = 0; j < cols; j += 4) {
            __m256d x0 =
                _mm256_add_pd(x_min_256, _mm256_mul_pd(kx, _mm256_set_pd(j + 3, j + 2, j + 1, j)));

            __m256d x = zero, y = zero, x2 = zero, y2 = zero, mk = zero;

            for (std::size_t k = 0; k < max_iter; ++k) {
                y = _mm256_add_pd(_mm256_mul_pd(two, _mm256_mul_pd(x, y)), y0);
                x = _mm256_add_pd(_mm256_sub_pd(x2, y2), x0);
                x2 = _mm256_mul_pd(x, x);
                y2 = _mm256_mul_pd(y, y);

                __m256d z_norm = _mm256_add_pd(x2, y2);
                __m256d mask = _mm256_cmp_pd(z_norm, threshold, _CMP_LT_OS);
                mk = _mm256_add_pd(_mm256_and_pd(mask, one), mk);

                // Early exit
                if (_mm256_testz_pd(mask, minus_one)) {
                    break;
                }
            }

            __m128i mki = _mm256_cvtpd_epi32(mk);
            p_iter = std::copy_n(reinterpret_cast<std::uint32_t*>(&mki), 4, p_iter);
        }
    }

    return iter_count;
}

inline png::byte interpolate(double c0, double c1, std::int8_t j, std::int8_t s) {
    return (png::byte) (c0 + j * (c1 - c0) / s);
}

std::vector<png::rgb_pixel> wiki_palette() {
    std::vector<png::rgb_pixel> palette;

    palette.emplace_back(66, 30, 15);     // brown 3
    palette.emplace_back(25, 7, 26);      // dark violett
    palette.emplace_back(9, 1, 47);       // darkest blue
    palette.emplace_back(4, 4, 73);       // blue 5
    palette.emplace_back(0, 7, 100);      // blue 4
    palette.emplace_back(12, 44, 138);    // blue 3
    palette.emplace_back(24, 82, 177);    // blue 2
    palette.emplace_back(57, 125, 209);   // blue 1
    palette.emplace_back(134, 181, 229);  // blue 0
    palette.emplace_back(211, 236, 248);  // lightest blue
    palette.emplace_back(241, 233, 191);  // lightest yellow
    palette.emplace_back(248, 201, 95);   // light yellow
    palette.emplace_back(255, 170, 0);    // dirty yellow
    palette.emplace_back(204, 128, 0);    // brown 0
    palette.emplace_back(153, 87, 0);     // brown 1
    palette.emplace_back(106, 52, 3);     // brown 2

    // linear interpolation (16 to 64)
    std::vector<png::rgb_pixel> large_palette;

    for (std::size_t i = 1; i < 16; ++i) {
        const auto &c0 = palette[i - 1], &c1 = palette[i];

        for (std::int8_t j = 0; j < 4; ++j) {
            large_palette.emplace_back(interpolate(c0.red, c1.red, j, 4),
                                       interpolate(c0.green, c1.green, j, 4),
                                       interpolate(c0.blue, c1.blue, j, 4));
        }
    }

    return large_palette;
}

int main() {
    const std::size_t cols = 512;
    const std::size_t max_iter = 3000;

    double src_x = -0.21503361460851339, src_y = 0.67999116792639069;
    std::size_t n_images = 630;

    // double src_x = -0.7746806106269039, src_y = -0.1374168856037867;
    // std::size_t n_images = 630;

    double width = 4;

    auto palette = wiki_palette();

    for (std::size_t i = 0; i < n_images; ++i) {
        std::cout << src_x << ',' << src_y << ',' << width << '\n';

        auto iter_count = mandelbrot_set(src_x, src_y, width, cols, max_iter);

        png::image<png::rgb_pixel> image(cols, iter_count.size() / cols);

        for (png::uint_32 y = 0; y < image.get_height(); ++y) {
            for (png::uint_32 x = 0; x < image.get_width(); ++x) {
                auto idx = iter_count[y * cols + x] % palette.size();
                image[y][x] = palette[idx];
            }
        }

        std::stringstream ss;
        ss << "movied_" << std::setw(3) << std::setfill('0') << i << ".png";
        image.write(ss.str());

        width -= width / 20;
    }

    return 0;
}