Fractals/Mandelbrot CLI: Renderer with Perturbation Theory

< Back to Fractals book Introduction This page contains the original C++ source code used to render high-precision fragments of the Mandelbrot set using perturbation theory and 8x8 Super-Sampling Anti-Aliasing (SSAA). Created by User:Aokoroko. Key Features High-Precision Reference: The 1000-bit reference trajectory is computed exactly once per zoom layer. Hardware-Native Performance: Blazing-fast math for billions of pixels utilizing hardware-native double registers. Bilinear approximation: The calculation can be performed significantly faster by using bilinear approximation. When using double-precision floating-point numbers (on the order of 10⁻¹⁵), perturbation theory only allows you to zoom down to the 10⁻³⁰⁸ level—no further. Innovative Algorithm: Revolutionary *Reference Reset to Zero* implementation. True 8x8 SSAA: Pristine, anti-aliased image quality with 64 independent samples per pixel. OpenMP Multi-threading: High-speed parallel computing to maximize CPU utilization. C++ Source Code #include <atomic> #include <cmath> #include <cstdint> #include <cstdio> #include <vector> #include <algorithm> #include <mpfr.h> #include <omp.h> using std::vector; using std::min; const char * CENTER_RE = "-1.99999543561201124623198345433951143502785679245726844745821388800402678499411681518036306219179273434395557574279985918047221291197081186140687781560831995"; const char * CENTER_IM = "-0.00000000000000000000000026198152173811047783694060060607013913873144250985383083459221663448338433592617272786772587281530484110756597337683912309313885172"; const char * VIEW_SIZE = "1.15e-119"; const int WIDTH = 2160; const int HEIGHT = 2160; const int AA = 8; const int MAX_ITER = 50000; const double ESCAPE_RADIUS_SQUARED = 50000.0; const int PALETTE_FRAME = 200; const char * OUTPUT_FILE = "Mandelbrot Set Image 112.bmp"; const mpfr_prec_t PRECISION_BITS = 1000; const int REF_SIZE = MAX_ITER + 200; 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(REF_SIZE); mpfr_set_ui(zReMp, 0, MPFR_RNDN); mpfr_set_ui(zImMp, 0, MPFR_RNDN); for (int iter = 0; iter < REF_SIZE - 1; ++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); std::fprintf(stderr, "Precomputing skip100 matrices...\n"); vector<Complex> coeff_A(REF_SIZE, {1.0, 0.0}); vector<Complex> coeff_B(REF_SIZE, {0.0, 0.0}); vector<double> rad_R(REF_SIZE, 2.0); vector<double> aS_squared(REF_SIZE, 0.0); for (int i = 0; i < referenceLength; ++i) { double r2 = referenceOrbit[i].re * referenceOrbit[i].re + referenceOrbit[i].im * referenceOrbit[i].im; aS_squared[i] = (r2 < ESCAPE_RADIUS_SQUARED) ? r2 : 0.0; } const int loop_limit = min(static_cast<int>(MAX_ITER), referenceLength - 105); #pragma omp parallel for for (int i = 0; i < loop_limit; ++i) { double min_r2 = ESCAPE_RADIUS_SQUARED; for (int k = 0; k < 100; ++k) { if (i + k >= referenceLength) break; if (aS_squared[i + k] < min_r2) min_r2 = aS_squared[i + k]; } rad_R[i] = std::sqrt(min_r2); for (int k = 0; k < 100; ++k) { if (i + k >= referenceLength) break; double r_re = referenceOrbit[i + k].re; double r_im = referenceOrbit[i + k].im; double next_A_re = 2.0 * (r_re * coeff_A[i].re - r_im * coeff_A[i].im); double next_A_im = 2.0 * (r_re * coeff_A[i].im + r_im * coeff_A[i].re); double next_B_re = 2.0 * (r_re * coeff_B[i].re - r_im * coeff_B[i].im) + 1.0; double next_B_im = 2.0 * (r_re * coeff_B[i].im + r_im * coeff_B[i].re); coeff_A[i].re = next_A_re; coeff_A[i].im = next_A_im; coeff_B[i].re = next_B_re; coeff_B[i].im = next_B_im; } } 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) { if (zRe * zRe + zIm * zIm >= ESCAPE_RADIUS_SQUARED) { break; } double eps_abs2 = deltaRe * deltaRe + deltaIm * deltaIm; double limit_r2 = 1e-60 * rad_R[referenceIndex] * rad_R[referenceIndex]; if (eps_abs2 < limit_r2 && (referenceIndex + 100 < loop_limit) && (iter + 100 < MAX_ITER)) { double backup_deltaRe = deltaRe; double backup_deltaIm = deltaIm; int backup_refIdx = referenceIndex; int backup_iter = iter; double next_eps_re = (coeff_A[referenceIndex].re * deltaRe - coeff_A[referenceIndex].im * deltaIm) + (coeff_B[referenceIndex].re * deltaCRe - coeff_B[referenceIndex].im * deltaCIm); double next_eps_im = (coeff_A[referenceIndex].re * deltaIm + coeff_A[referenceIndex].im * deltaRe) + (coeff_B[referenceIndex].re * deltaCIm + coeff_B[referenceIndex].im * deltaCRe); deltaRe = next_eps_re; deltaIm = next_eps_im; referenceIndex += 100; iter += 100; zRe = reference[referenceIndex].re + deltaRe; zIm = reference[referenceIndex].im + deltaIm; if (zRe * zRe + zIm * zIm >= ESCAPE_RADIUS_SQUARED) { deltaRe = backup_deltaRe; deltaIm = backup_deltaIm; referenceIndex = backup_refIdx; iter = backup_iter; } else { continue; } } 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 >= loop_limit) { 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! Saved to %s in %.2f seconds.\n", OUTPUT_FILE, omp_get_wtime() - startTime); return 0; } Rendered Examples External Links Official Mandelbrot CLI Repository on GitHub — source code, documentation, and pre-compiled releases. C++ Perturbation Theory Optimization showcased in the global algorithm repository on Rosetta Code. Category:Book:Fractals