Консольний рендерер Мандельброта з методом обурень та межі 10−308 (мовою програмування C++)
Вступ Даний практикум містить оригінальний вихідний код на C++ для високоточного рендерингу фрагментів множини Мандельброта з використанням алгоритмів оптимізації (включаючи теорію збурень) та глибокого згладжування (8x8 SSAA). Автор алгоритму та зображень: User:Aokoroko.
Ключові особливості Розрахунок опорної траєкторії на 1000 біт усього один раз. Реактивний розрахунок мільярда пікселів на апаратному double. Розрахунок можна виконати суттєво швидше, якщо використовувати білінійну апроксимацію. При використанні чисел із плаваючою комою подвійної точності (порядку 10⁻¹⁵) теорія збурень дозволяє наблизитися до рівня 10⁻³⁰⁸ - не далі. Революційний алгоритм Reference Reset to Zero. Справжній SSAA 8x8 для ідеально згладженого зображення без аліасингу. Паралелізм OpenMP для високошвидкісного багатопотокового рендерингу.
Вихідний код C++ Ниже представлений вихідний код програми, що використовувався для генерації 100-мегапіксельних вибраних зображень на Вікісховищі.
#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; for (int i = 0; i < 255; ++i) { palette[i] = static_cast<uint8_t>(std::lround(127.0 + 127.0 * std::cos(2.0 * PI * i / 255.0))); palette[i] = static_cast<uint8_t>(std::lround(127.0 + 127.0 * std::sin(2.0 * PI * i / 255.0))); palette[i] = palette[i]; } palette = 255; palette = 255; palette = 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]; greenSum += palette[paletteIndex]; redSum += palette[paletteIndex];
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; }
Приклади зображень
Посилання Офіційний репозиторій проєкту Mandelbrot CLI на GitHub — вихідний код, документація та готові релізи програми рендерингу. Реалізація теорії збурень на C++ у міжнародній бібліотеці алгоритмів Rosetta Code.