Алгоритми множини Мандельброта: час виходу, оцінка відстані та GPU-рендеринг

1. Означення: одне рівняння, нескінченна складність

Множина Мандельброта M — це множина всіх комплексних чисел c, для яких орбіта 0 при повторному застосуванні f(z) = z² + c назавжди лишається обмеженою:

f_c(z) = z² + c Початкове значення: z₀ = 0 Орбіта: z₁ = 0² + c = c z₂ = c² + c z₃ = (c²+c)² + c ... c \in M ⟺ |z_n| не \to \infty при n \to \infty Ключовий факт: якщо |z_n| коли-небудь перевищує 2, орбіта виходить у нескінченність. Доведення: якщо |z| > 2 і |c| < 2, то |z²+c| \geq |z|² - |c| > 2|z| - |z| = |z|. Модуль зростає щонайменше на |z| за крок \to розбігається.

Межа M — нескінченно заплутана берегова лінія, видима у відео фрактального зуму, — це місце, де живуть цікаві алгоритмічні виклики. Точки строго всередині M ітерують вічно; точки ззовні виходять на якомусь скінченному кроці. Сама межа має розмірність Гаусдорфа 2 (вона має нульову площу, але настільки складна, що в топологічному сенсі заповнює площину).

Множина лежить у комплексній площині з дійсною частиною (Re c) на осі x та уявною частиною (Im c) на осі y. Уся множина вміщується в колі радіуса 2 із центром у початку координат, тож нам потрібно вибирати лише значення c з |c| ≤ 2.

2. Алгоритм часу виходу

Класичний алгоритм ітерує відповідне значення c кожного пікселя, доки |z|² > 4 (еквівалентно |z| > 2, але без обчислення квадратного кореня) або доки не досягнуто максимальної кількості ітерацій MAX_ITER:

// Псевдокод — наївний час виходу function mandelbrot(cx, cy, max_iter): zx = 0.0, zy = 0.0 for n in 0 .. max_iter: if (zx*zx + zy*zy) > 4.0: return n // вийшов — повернути кількість ітерацій zx_new = zx*zx - zy*zy + cx zy = 2.0*zx*zy + cy zx = zx_new return max_iter // не вийшов \to «всередині» множини

Повернене значення n потім відображається в колір. Найпростіше відображення призначає відтінок, пропорційний до n / max_iter. Але це дає різкі видимі смуги на кожній цілій кількості ітерацій — що приводить до прийому гладкого розфарбування, описаного в наступному розділі.

Вибір MAX_ITER

За типового вигляду (зум ≈ 1×) 256–512 ітерацій достатньо для впізнаваної деталізації. Зображення глибокого зуму вимагають мільйонів ітерацій для точок поблизу межі. Емпіричне правило: MAX_ITER ≈ 200 × (рівень зуму)^0.4. Замало ітерацій роблять множину «товстою» — зовнішні точки поблизу межі хибно класифікуються як внутрішні.

Зауваження щодо продуктивності: за роздільності 1080p (1920×1080 ≈ 2 мільйони пікселів) із 1000 ітерацій кожен наївний однопотоковий CPU-рендерер виконує ~2 мільярди множень з рухомою комою на кадр. Сучасний CPU витрачає ~2–5 секунд на кадр. GPU-шейдер виконує це менш ніж за 16 мс (60 кадрів/с), бо кожен піксель — незалежне обчислення — означення «ганебно паралельних» (embarrassingly parallel) задач.

Оптимізація: пропустити квадратний корінь

Перевірка zx*zx + zy*zy > 4.0 уникає дорогого виклику sqrt(). Інша рання оптимізація — періодично перевіряти, чи z збігся до циклу — що вказує, що точка всередині M — заощаджуючи час на «точно внутрішніх» точках.

3. Гладке розфарбування та нормалізована кількість ітерацій

Цілочислові лічильники виходу породжують кольорові смуги — концентричні кільця, чітко розділені на кожному цілому n. Нормалізована кількість ітерацій (також звана формулою Büldt або Linas) інтерполює між цілими, використовуючи логарифм модуля в момент виходу:

// Після виходу на ітерації n, z має модуль |z| = r: nu = n - log2(log2(|z|)) = n - log2( log(|z|) / log(2) ) Інтерпретація: «справжній» дробовий час виходу — це значення n, при якому |z| дорівнювало б точно 2, якби ми дозволили нецілі кроки. Використання log₂(log₂(|z|)) коригує квадратичне зростання |z|. Гладкий колір: colour(nu) — інтерполювати кольорову палітру за nu mod 1

Це усуває всі видимі смуги, даючи гладке градієнтне розфарбування, яке видно в професійних фрактальних рендерах. Дизайн палітри потім визначає кінцевий вигляд: циклічна зміна відтінку, використання синусоїдальних RGB-функцій або відображення в ручно створені точки градієнта.

Орбітальні пастки

Більш художній метод розфарбування записує мінімальну відстань від орбіти до геометричної фігури ( точки, лінії, кола чи хреста). Колір залежить від цієї мінімальної відстані. Орбітальні пастки дають разюче різні зображення — повторювані візерунки форми пастки вздовж фрактальної межі — й можуть генерувати текстури, що нагадують реальні об'єкти.

// Орбітальна пастка: точкова пастка в (tx, ty) min_dist = +\infty для кожної ітерації (zx, zy): d = sqrt((zx-tx)² + (zy-ty)²) min_dist = min(min_dist, d) // Колір на основі min_dist (наприклад, log(min_dist) \to відтінок)

4. Виявлення внутрішньої області: перевірка періоду та тест кардіоїди

Точки всередині M ітерують аж до MAX_ITER, що робить внутрішні пікселі найдорожчим попіксельним обчисленням. Два прийоми разюче зменшують вартість обробки внутрішньої області.

Тест кардіоїди та бульбашки періоду 2

Головна кардіоїда та найбільша кругла бульбашка (з центром у c = −1) разом містять велику частку всіх внутрішніх точок. Їх можна перевірити аналітично перед ітеруванням:

// Тест кардіоїди q = (cx - 0.25)² + cy² in_cardioid = q * (q + (cx - 0.25)) < 0.25 * cy² // Тест бульбашки періоду 2 in_p2_bulb = (cx + 1)² + cy² < 0.0625 if (in_cardioid || in_p2_bulb): return max_iter одразу

Перевірка періодичності / циклів

Внутрішні орбіти зрештою входять у періодичний цикл — орбіта повторюється з певним періодом p. Виявлення циклів за Брентом або Флойдом може це ідентифікувати: зберігайте «повільну» копію z кожні k ітерацій, порівнюйте з поточним z. Якщо |z_current - z_saved| < ε, точка в M, і ітерація припиняється достроково.

// Перевірка періоду в стилі Брента (спрощено) zx_old = 0, zy_old = 0 period = 0, check = 0 for n in 0 .. max_iter: ітерувати z if |z - z_old| < 1e-10: return max_iter // цикл виявлено \to всередині M check++ if check == period: zx_old = zx; zy_old = zy check = 0; period += period // подвоїти інтервал перевірки

У поєднанні тест кардіоїди + перевірка періоду зменшують внутрішні обчислення на 40–80% за типових рівнів зуму, приблизно подвоюючи середню швидкість рендерингу.

5. Оцінка відстані для деталізації межі

Обчислення зовнішнього оцінювача відстані (DEM) дає приблизну відстань від точки c до найближчої межі M, не знаючи точно, де ця межа. Це досягається одночасним відстеженням похідної функції ітерації:

// Оцінка відстані: одночасно відстежувати z і dz/dc z = 0 + 0i, dz = 1 + 0i // dz — похідна z по c for n in 0 .. max_iter: dz = 2 * z * dz + 1 // ланцюгове правило: d/dc [z²+c] = 2z\cdot(dz/dc) + 1 z = z*z + c if |z|² > bailout²: break // Оцінка зовнішньої відстані: d \approx |z| * log(|z|) / |dz| If d < pixel_size: піксель лежить на межі або поблизу неї \to розфарбувати його.

Застосування DEM:

Згладжування (anti-aliasing): пікселі, ближчі до межі, ніж ширина одного пікселя, отримують субпіксельну вибірку (суперсемплінг або Монте- Карло).
Рендеринг країв: призначити окремий колір усім пікселям, де DEM < поріг — малює контур множини без ітерування до великих глибин.
3D-рельєфний рендеринг: використовувати DEM як поле висот; затінювати віртуальним спрямованим світлом, щоб отримати «тиснену» фрактальну рельєфну карту.

Точність: DEM — лише наближення, дійсне для зовнішньої області M. Поблизу точок Мізюревича (де критична орбіта є передперіодичною, а не періодичною) формула може переоцінювати відстані. На практиці оцінка достатньо точна для рендерингу.

6. Глибокий зум та теорія збурень

Зум до глибин понад ~10¹⁴ вичерпує точність 64-бітних чисел подвійної точності (які мають ~15–16 десяткових цифр мантиси). На таких глибинах помилки округлення накопичуються, і зображення вироджується в неправильні однокольорові області, звані глітчами.

Арифметика довільної точності

Прямолінійне рішення замінює double-числа на числа довільної точності (наприклад, 128-бітні, 256-бітні чи навіть бібліотеки на 1000+ цифр). Але вартість зростає приблизно квадратично з точністю — у 10× більше цифр → у 100× повільніше. Це робить інтерактивні рендери глибокого зуму неможливими без додаткової кмітливості.

Теорія збурень

Ключова ідея: більшість пікселів у вигляді глибокого зуму дуже близькі одне до одного в комплексній площині. Обчисліть одну опорну орбіту в центрі вікна перегляду за допомогою високоточної арифметики. Потім представте кожен інший піксель як мале збурення δ відносно опорної орбіти:

Z = опорна орбіта (обчислена з високою точністю) z = Z + δ (орбіта пікселя, мале збурення) z_{n+1} = z_n² + c (Z + δ)_{n+1} = (Z + δ)_n² + (C + Δc) ≈ Z_n² + 2·Z_n·δ_n + δ_n² + C + Δc δ_{n+1} = 2·Z_n·δ_n + δ_n² + Δc Оскільки |δ| ≪ 1, обчислюйте ітерації δ з подвійною точністю. Лише опорна Z потребує високої точності.

Це зводить високоточну арифметику до однієї орбіти на кадр (опорної), тоді як усі інші пікселі використовують дешеву double-арифметику. Глітчі через збій наближення виявляються й виправляються вибором додаткових опорних орбіт поблизу проблемних пікселів.

Найсучасніші рендерери Мандельброта (як-от Kalles Fraktaler та Mandel Machine) поєднують теорію збурень із наближенням рядами (Series Approximation, SA), яке аналітично пропускає цілі блоки ранніх ітерацій — уможливлюючи інтерактивний зум до глибин 10¹⁰⁰⁰ і далі.

7. GPU-рендеринг із шейдерами GLSL

Сучасні GPU виконують тисячі викликів шейдера паралельно. Фрагментний шейдер, що обробляє один піксель за виклик, ідеально лягає на цикл ітерації Мандельброта. Кожному пікселю не потрібні дані від жодного іншого пікселя — синхронізація не потрібна.

// Каркас фрагментного шейдера GLSL uniform vec2 u_center; // центр вигляду в комплексній площині uniform float u_scale; // одиниць на піксель uniform int u_max_iter; void main() { vec2 c = u_center + (gl_FragCoord.xy - u_resolution*0.5) * u_scale; vec2 z = vec2(0.0); float n = 0.0; for (int i = 0; i < MAX_ITER_CONST; i++) { if (dot(z, z) > 4.0) break; z = vec2(z.x*z.x - z.y*z.y + c.x, 2.0*z.x*z.y + c.y); n++; } // Гладке розфарбування float nu = n - log2(log2(length(z))); float t = nu / float(MAX_ITER_CONST); gl_FragColor = palette(t); }

Головне обмеження GPU-шейдерів — у старіших версіях GLSL межа циклу має бути константою часу компіляції. Запихання надто багатьох ітерацій переповнює кеш інструкцій GPU або перевищує тайм-аути сторожового таймера драйвера, спричиняючи скидання екрана. Практичні межі GPU — 10 000–50 000 ітерацій на піксель за інтерактивної частоти кадрів.

Подвійна точність на GPU

Споживчі відеокарти пропонують апаратну подвійну точність (64-бітну) арифметику з пропускною здатністю приблизно 1/32 від одинарної точності. Для рендерингу глибокого зуму, де потрібні double-числа, перевага GPU зменшується. Поширений обхідний шлях — double-float емуляція — представлення кожного double як пари float за алгоритмом double-float Деккера–Кнута. Це повертає точність подвійної точності за ~3–4× вартості одинарної точності, усе ще набагато швидше, ніж довільна точність на боці CPU.

Особливості реалізації у WebGL

WebGL 1 потребує розширення OES_element_index_uint для 32-бітних індексів; немає вбудованої підтримки double.
WebGL 2 підтримує EXT_color_buffer_float для 32-бітних float-цілей рендерингу (корисно для проходів накопичення).
Для інтерактивних демо рендерте на половинній роздільності й масштабуйте вгору — це вдвічі зменшує загальну кількість викликів фрагментів за помірної втрати якості.
Прогресивне уточнення: спершу відрендерте грубий прохід, потім уточніть. Користувачі сприймають грубий попередній перегляд як миттєвий відгук.

8. Множини Жюліа: сімейство родичів

Для кожного комплексного числа c відповідна множина Жюліа J(c) обчислюється тією самою ітерацією — але тепер c фіксоване, а z змінюється по всьому зображенню:

// Множина Жюліа: c — стала, z — змінна function julia(zx, zy, cx, cy, max_iter): for n in 0 .. max_iter: if (zx² + zy²) > 4.0: return n zx_new = zx² - zy² + cx zy = 2\cdotzx\cdotzy + cy zx = zx_new return max_iter

Множини Мандельброта та Жюліа тісно пов'язані: якщо c всередині M, J(c) — зв'язний фрактал; якщо c поза M, J(c) — незв'язний «пил» (топологія канторової множини). Множину Мандельброта буквально можна вважати картою всіх можливих форм Жюліа — кожна точка в M відповідає зв'язній множині Жюліа, а форма M у малих масштабах віддзеркалює сусідні множини Жюліа.

Цей зв'язок використовується в інтерактивних дослідниках: клацання по будь-якій точці c у переглядачі Мандельброта миттєво оновлює бічну панель, що показує J(c). Для естетичного дослідження значення c уздовж межі M — філаменти та антени — дають найзаплутаніші форми Жюліа.

Мандельброт як простір параметрів: технічно множина Мандельброта — це локус зв'язності квадратичного сімейства f_c(z) = z²+c. Сімейства вищого степеня (z³+c, z⁴+c, …) породжують множини Multibrot — схожі за духом, з (d−1)-кратною обертовою симетрією та відповідно складнішою структурою. Ті самі алгоритми часу виходу та оцінки відстані застосовуються з незначними модифікаціями.