File:Mandelbrot Set Image 108.png
Jump to navigation
Jump to search
Size of this preview: 600 × 600 pixels. Other resolutions: 240 × 240 pixels | 480 × 480 pixels | 768 × 768 pixels | 1,024 × 1,024 pixels | 2,048 × 2,048 pixels | 10,000 × 10,000 pixels.
Original file (10,000 × 10,000 pixels, file size: 43.14 MB, MIME type: image/png)
Summary
| Description |
Русский: Фрагмент множества Мандельброта, координаты центра: -1.7489436617686633372073552153211507258063533373824414679761, -0.0000073748967541889836640985849393311615399776865199722997 ширина изображения 1.01e-55
English: Fragment of the Mandelbrot set, coordinates: -1.7489436617686633372073552153211507258063533373824414679761, -0.0000073748967541889836640985849393311615399776865199722997 width 1.01e-55
Беларуская: Фрагмент мноства Мандэльброта, каардынаты цэнтра: -1.7489436617686633372073552153211507258063533373824414679761, -0.0000073748967541889836640985849393311615399776865199722997 шырыня 1.01e-55 |
| Date | |
| Source | Own work |
| Author | Aokoroko |
| Other versions |
|
| Source code (C++) |
#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.7489436617686633372073552153211507258063533373824414679761";
const char * CENTER_IM = "-0.0000073748967541889836640985849393311615399776865199722997";
const char * VIEW_SIZE = "1.01e-55";
const int WIDTH = 10000;
const int HEIGHT = 10000;
const int AA = 8;
const int MAX_ITER = 50000;
const double ESCAPE_RADIUS_SQUARED = 50000.0;
const int PALETTE_FRAME = 1;
const char * OUTPUT_FILE = "Mandelbrot Set Image 108.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-30 * 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;
}
|
Technical details
- 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-15, perturbation theory only allows you to zoom down to the 10-308 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. 80000 x 80000 pixels downscaled to 10000 x 10000.
- OpenMP Multi-threading: High-speed parallel computing to maximize CPU utilization.
- Software: C++ (compiled with g++), GNU C++ Compiler.
Related images
-
Previous step -
Next step
Notes
- Rosetta Code: https://rosettacode.org/wiki/Mandelbrot_set#Perturbation_Theory
- github: https://github.com/Divetoxx/Mandelbrot
Licensing
I, the copyright holder of this work, hereby publish it under the following license:
 
|
This file is made available under the Creative Commons CC0 1.0 Universal Public Domain Dedication. |
| The person who associated a work with this deed has dedicated the work to the public domain by waiving all of their rights to the work worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law. You can copy, modify, distribute and perform the work, even for commercial purposes, all without asking permission.
|
|
This image has been assessed using the Quality image guidelines and is considered a Quality image.
|
File history
Click on a date/time to view the file as it appeared at that time.
| Date/Time | Thumbnail | Dimensions | User | Comment | |
|---|---|---|---|---|---|
| current | 13:43, 13 June 2026 | | 10,000 × 10,000 (43.14 MB) | wikimediacommons>Aokoroko | Uploaded own work with UploadWizard |
File usage
There are no pages that use this file.
