Comments 21
чтобы рендерить с шириной e-52, нужна спец-математика (в плане что обычные double/longdouble не вмещают такую точность). что используете, какие библиотеки? и думаю, что много кто, узнав о красоте множества Мандельброта, сам писал подобное рисовалово. Я точно писал, но исходники накрылись с тех пор, и математика была “обычной”.
float (4 байта) до e-38, double (8 байт) до e-308, long double данные разнятся, от -308 до -4932. Так что спец математика не слишком нужна.
Когда считаешь y=x^2-1 (вроде там такая итерация), а разные иксы отличаются на 1e-52, а сами имеют размер 1e-0, то есть условное (-0.5, -0.5) и окрестности радиусом 1е-52, в лоб не посчитать, потеряешь точность. Вот и интересует, какие трюки использовались при рендере ММБ с таким разрешением.
А, понял проблему. Автор ничего не говорит об этом, его обе программы используют простой вариант с умножением и обычными типами данных.
Однократный расчёт опоры: Сверхтяжелый BigFloat-радар MPFR вычисляет точную траекторию всего для одной-единственной центральной точки кадра и строго ОДИН раз в начале рендеринга.
Аппаратное ускорение на double: Весь остальной массив экрана (миллионы супер-пикселей) рассчитывается параллельно на бешеной скорости чистых, аппаратных регистров
doubleпроцессора, вычисляя лишь микроскопические отклонения (дельты) от центральной оси. Скорость генерации взлетела в 1000 раз!
Обычно под такие координаты берут GMP или собирают самописный класс длинной арифметики на AVX-инструкциях
#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.7580082869932816847911302993841322510696136102893411525626408121299664464945680592947275781494";
const char * CENTER_IM = "-0.0120395123143965484759056667813645320687762663207866292850653888764442902115988533819610775691";
const char * VIEW_SIZE = "1.27e-90";
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 = 70;
const char * OUTPUT_FILE = "Mandelbrot Set Image 115.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;
}
Да, ни одна цивилизация древности не оставила изображений множества Мандельброта, значит все разговоры об их продвинутости не слишком основательны. Никакие индейцы майя и другие, кому приписывают будто бы высокий уровень развития, множество Мандельброта не открыли.
Ну да, ни арабские и индийские математики, ни египетские астрономы, ни составитель календарей майя не изобрели микроволновок, айфонов и интернет. Ничтожества, правда?
А те, кто приручили огонь, развели огороды и применил хинин- неразвитые существа, не тем в жизни занимались
Кажется это про альтернативные исторические теории, где мегаразвитые цивилизации прошлого с помощью лазерных пил под руководством инопланетян создавали сооружения, которые сегодня не умеют создавать, а не про то, что раньше жили идиоты
вы слишком истерично реагируете на безобидное высказывание. думаю, автор хотел сказать, что по сравнению со знаниями современной цивилизации, люди прошлого (особенно изолированные группы, вроде индейцев) просто дикари необразованные. и это факт. для людей будущего таковыми дикарями будем уже мы. некоторые индивиды любят возвеличивать любую реальную или выдуманную цивилизацию, принижая заслуги современной человеческой. мотив мне не понятен, думается, источник в каких-то психологических комплексах.
@Dead_Bit, слово "истерично" неуместно.
автор пока не высказался о том, что он имел в виду.
по тексту - противопоставлялись древние цивилизации и достижения автора:
Это сделал я. Остальное - компьютер. Поиск красивых фрагментов цветных версий множества Мандельброта - вот я этим и занимаюсь.
Я написал программы на ГитХабе причем README не только "Language English" но и "Язык Русский"!
Восемь фракталов множества Мандельброта я создал сам и загрузил на Викисклад.
...
Это своего рода «виртуальная выставка» лучших работ проекта.
А пробовали поиграть с абберациями?
Отличный пример того, как детское увлечение фракталами из советских журналов дорастает до собственных оптимизированных утилит с OpenMP) Исходный код на паскале в цикле прямо вернул в школьные годы
кстати, ваши косинусоидальные волны приравнивают красный к зеленому, в результате на картинках только два цвета, синий и желтый. я для своей визуализации делал программу, которая рисует 256 переходов по заданным смещениям и амплитудам каналов, волны могли к тому же иметь больше одного периода на 256 цветов, и смотрел, нравится ли мне получившаяся гамма, если да, то впихивал её в визуализатор и рендерил картинки в ней, было весело. Советую подобрать что-то поинтереснее.
Множество Мандельброта