File:Mandelbrot Set Image 113.png

#include <atomic>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <vector>
#include <mpfr.h>
#include <omp.h>

using std::vector;

const char * CENTER_RE = "-1.028457781372052547774824579171003431968291180173";
const char * CENTER_IM = "0.361463930751890432883338309925137437795829478445";
const char * VIEW_SIZE = "1.54e-44";

const int WIDTH = 10000;
const int HEIGHT = 10000;
const int AA = 8;
const int MAX_ITER = 50000;
const double ESCAPE_RADIUS_SQUARED = 40000.0;
const int PALETTE_FRAME = 155;
const char * OUTPUT_FILE = "Mandelbrot Set Image 113.bmp";
const mpfr_prec_t PRECISION_BITS = 350;

struct Complex {
  double re;
  double im;
};

#pragma pack(push, 1)
struct BMPHeader {
  uint16_t type {
    0x4D42
  };
  uint32_t size {
    0
  };
  uint32_t reserved {
    0
  };
  uint32_t offBits {
    54
  };
  uint32_t structSize {
    40
  };
  int32_t width {
    0
  };
  int32_t height {
    0
  };
  uint16_t planes {
    1
  };
  uint16_t bitCount {
    24
  };
  uint32_t compression {
    0
  };
  uint32_t sizeImage {
    0
  };
  int32_t xPixelsPerMeter {
    2834
  };
  int32_t yPixelsPerMeter {
    2834
  };
  uint32_t colorsUsed {
    0
  };
  uint32_t colorsImportant {
    0
  };
};
#pragma pack(pop)

int main() {
  const double startTime = omp_get_wtime();
  const long rawWidth = static_cast < long > (WIDTH) * AA;
  const long rawHeight = static_cast < long > (HEIGHT) * AA;

  mpfr_t centerRe, centerIm, zReMp, zImMp, tmp1, tmp2, viewSizeMp;
  mpfr_inits2(
    PRECISION_BITS,
    centerRe,
    centerIm,
    zReMp,
    zImMp,
    tmp1,
    tmp2,
    viewSizeMp,
    static_cast < mpfr_ptr > (nullptr)
  );

  mpfr_set_str(centerRe, CENTER_RE, 10, MPFR_RNDN);
  mpfr_set_str(centerIm, CENTER_IM, 10, MPFR_RNDN);
  mpfr_set_str(viewSizeMp, VIEW_SIZE, 10, MPFR_RNDN);

  const double sampleStep = mpfr_get_d(viewSizeMp, MPFR_RNDN) / rawWidth;

  vector < Complex > referenceOrbit;
  referenceOrbit.reserve(MAX_ITER + 1);

  mpfr_set_ui(zReMp, 0, MPFR_RNDN);
  mpfr_set_ui(zImMp, 0, MPFR_RNDN);

  for (int iter = 0; iter <= MAX_ITER; ++iter) {
    Complex z {
      mpfr_get_d(zReMp, MPFR_RNDN),
        mpfr_get_d(zImMp, MPFR_RNDN)
    };

    referenceOrbit.push_back(z);

    if (z.re * z.re + z.im * z.im > ESCAPE_RADIUS_SQUARED) {
      break;
    }

    mpfr_mul(tmp1, zReMp, zImMp, MPFR_RNDN);
    mpfr_sqr(tmp2, zReMp, MPFR_RNDN);
    mpfr_sqr(zReMp, zImMp, MPFR_RNDN);
    mpfr_sub(zReMp, tmp2, zReMp, MPFR_RNDN);
    mpfr_add(zReMp, zReMp, centerRe, MPFR_RNDN);
    mpfr_mul_2ui(tmp1, tmp1, 1, MPFR_RNDN);
    mpfr_add(zImMp, tmp1, centerIm, MPFR_RNDN);
  }

  const int referenceLength = static_cast < int > (referenceOrbit.size());

  mpfr_clears(
    centerRe,
    centerIm,
    zReMp,
    zImMp,
    tmp1,
    tmp2,
    viewSizeMp,
    static_cast < mpfr_ptr > (nullptr)
  );

  std::fprintf(stderr, "Reference orbit: %d points\n", referenceLength);

  const double PI = 3.14159265358979323846;
  uint8_t palette[256][3];

  for (int i = 0; i < 255; ++i) {
    palette[i][0] = static_cast < uint8_t > (
      std::lround(127.0 + 127.0 * std::cos(2.0 * PI * i / 255.0))
    );
    palette[i][1] = static_cast < uint8_t > (
      std::lround(127.0 + 127.0 * std::sin(2.0 * PI * i / 255.0))
    );
    palette[i][2] = palette[i][1];
  }

  palette[255][0] = 255;
  palette[255][1] = 255;
  palette[255][2] = 255;

  const int rowBytes = (WIDTH * 3 + 3) & ~3;
  vector < uint8_t > image(static_cast < size_t > (rowBytes) * HEIGHT, 0);

  std::atomic < int > completedRows {
    0
  };
  const Complex * reference = referenceOrbit.data();

  #pragma omp parallel for schedule(dynamic)
  for (int y = 0; y < HEIGHT; ++y) {
    uint8_t * row = image.data() + static_cast < size_t > (y) * rowBytes;

    for (int x = 0; x < WIDTH; ++x) {
      unsigned blueSum = 0;
      unsigned greenSum = 0;
      unsigned redSum = 0;

      for (int sampleY = 0; sampleY < AA; ++sampleY) {
        const double deltaCIm =
          (static_cast < long > (y) * AA + sampleY - rawHeight / 2) * sampleStep;

        for (int sampleX = 0; sampleX < AA; ++sampleX) {
          const double deltaCRe =
            (static_cast < long > (x) * AA + sampleX - rawWidth / 2) * sampleStep;

          double deltaRe = 0.0;
          double deltaIm = 0.0;
          double zRe = 0.0;
          double zIm = 0.0;

          int referenceIndex = 0;
          int iter = 0;

          while (
            iter < MAX_ITER &&
            zRe * zRe + zIm * zIm < ESCAPE_RADIUS_SQUARED
          ) {
            const double a = 2.0 * reference[referenceIndex].re + deltaRe;
            const double b = 2.0 * reference[referenceIndex].im + deltaIm;

            const double nextDeltaRe =
              a * deltaRe - b * deltaIm + deltaCRe;

            deltaIm = a * deltaIm + b * deltaRe + deltaCIm;
            deltaRe = nextDeltaRe;

            ++referenceIndex;
            ++iter;

            zRe = reference[referenceIndex].re + deltaRe;
            zIm = reference[referenceIndex].im + deltaIm;

            if (
              zRe * zRe + zIm * zIm <
              deltaRe * deltaRe + deltaIm * deltaIm ||
              referenceIndex == referenceLength - 1
            ) {
              deltaRe = zRe;
              deltaIm = zIm;
              referenceIndex = 0;
            }
          }

          const int remaining = MAX_ITER - iter;
          const uint8_t colorIndex =
            (remaining == 0) ?
            255 :
            static_cast < uint8_t > (remaining % 254);

          const int paletteIndex =
            (colorIndex == 255) ?
            255 :
            (colorIndex - PALETTE_FRAME + 255) % 255;

          blueSum += palette[paletteIndex][0];
          greenSum += palette[paletteIndex][1];
          redSum += palette[paletteIndex][2];
        }
      }

      const int samples = AA * AA;
      row[x * 3 + 0] = static_cast < uint8_t > (blueSum / samples);
      row[x * 3 + 1] = static_cast < uint8_t > (greenSum / samples);
      row[x * 3 + 2] = static_cast < uint8_t > (redSum / samples);
    }

    const int done = ++completedRows;

    if (done % 50 == 0 || done == HEIGHT) {
      std::fprintf(
        stderr,
        "\rProgress: %d/%d rows (%.1f%%)",
        done,
        HEIGHT,
        100.0 * done / HEIGHT
      );
    }
  }

  BMPHeader header;
  header.width = WIDTH;
  header.height = HEIGHT;
  header.sizeImage = static_cast < uint32_t > (image.size());
  header.size = header.sizeImage + 54;

  FILE * file = std::fopen(OUTPUT_FILE, "wb");

  if (!file) {
    std::perror(OUTPUT_FILE);
    return 1;
  }

  std::fwrite( & header, sizeof header, 1, file);
  std::fwrite(image.data(), 1, image.size(), file);
  std::fclose(file);

  std::fprintf(
    stderr,
    "\nDone: \"%s\" saved in %.1f s\n",
    OUTPUT_FILE,
    omp_get_wtime() - startTime
  );

  return 0;
}