Compare commits

..
2 Commits
Author SHA1 Message Date
surv 61d4766088 perf: disable AA when switching fractals 2026-09-24 21:56:39 +02:00
surv 51fa5ffa02 feat: add interpolation between fractals 2026-09-24 21:43:36 +02:00
8 changed files with 644 additions and 215 deletions
+10 -1
View File
@@ -129,7 +129,16 @@ pixel is a handful of `f32` complex multiplies.
exact delta). `fprime(z)` is the derivative used for distance-estimation
(DE) shading; exact for holomorphic kinds, an approximation (`~2Z`) for the
abs-based ones. A `KIND_*` constant (from `common.wgsl`) must match the
matching `FractalKind` variant's discriminant exactly.
matching `FractalKind` variant's discriminant exactly. The per-kind bodies
are `advance_delta_kind`/`fprime_kind`; `advance_delta`/`fprime` wrap them
to blend two kinds during the kind-switch morph (`u.morph_from`,
`u.morph_w`: each step is `(1-w)·f_kind + w·f_from`, mirrored on the CPU by
the `morph` argument of `compute_reference`, in both its f64 and `FBig`
paths). The blend only exists in pipelines built with the `MORPH` override
(part of `PipelineKey`, on while `morph_w > 0`); those also skip periodicity
detection and the cardioid bypass. App side: `KindMorph` in
`app.rs`; the uniforms use the morph the *current reference* was built with
(`ref_morph`), not the live one, so orbit and delta formula never disagree.
- `src/fractal/renderer.rs` — `FractalRenderer` (wgpu pipelines, uniform +
storage buffers, bind groups), `Uniforms` (repr(C) layout that must match
the WGSL `Uniforms` struct field-for-field, including padding; it includes
+156 -12
View File
@@ -20,7 +20,7 @@ use crate::view::parse_half_height_spec;
use crate::view::parse_re_im_spec;
use crate::view::{
Big, DEFAULT_HALF_HEIGHT, ViewState, big_from_decimal_str, big_from_f64, big_to_decimal_str,
parse_view_spec, precision_for,
interpolate_view, parse_view_spec, precision_for,
};
#[cfg(not(target_arch = "wasm32"))]
use clap::Parser;
@@ -161,6 +161,35 @@ struct RequestKey {
kind: FractalKind,
power: u32,
complex_power: (f64, f64),
/// Kind-switch morph `(from_kind, weight)`, if one is running.
morph: Option<(FractalKind, f32)>,
}
/// An in-progress kind-switch animation: the iteration formula is blended per
/// step from `from` to the current kind, `(1 - w)·f_kind + w·f_from`, while the
/// camera glides from `from_view` to the new kind's default view.
struct KindMorph {
from: FractalKind,
/// Linear progress in [0, 1]; eased with smoothstep.
progress: f32,
from_view: ViewState,
to_view: ViewState,
/// Whether the morph still drives the camera. Cleared as soon as the user
/// pans/zooms, so they can take over mid-morph.
camera: bool,
}
impl KindMorph {
/// Smoothstep-eased progress.
fn eased(&self) -> f32 {
let p = self.progress.clamp(0.0, 1.0);
p * p * (3.0 - 2.0 * p)
}
/// Weight of the old kind's formula: 1 at the start, 0 at the end.
fn weight(&self) -> f32 {
1.0 - self.eased()
}
}
/// Shared state for an in-progress PNG export. The worker (a background thread
@@ -213,6 +242,12 @@ struct AnimState {
/// e-folds per second; positive zooms in, negative zooms out.
zoom_speed: f32,
/// Morph the iteration formula (and camera) when switching fractal kinds,
/// instead of cutting straight to the new kind.
kind_morph: bool,
/// Kind-switch morph duration, in seconds.
kind_morph_duration: f32,
/// Linear 2D <-> 3D transition progress in [0, 1], advanced at a constant
/// rate; `camera_state` is its smoothstep-eased value.
camera_progress: f32,
@@ -242,6 +277,8 @@ impl Default for AnimState {
lambda_angle: 0.0,
zoom: false,
zoom_speed: 0.5,
kind_morph: true,
kind_morph_duration: 1.5,
camera_progress: 0.,
camera_state: 0.,
}
@@ -309,6 +346,8 @@ pub struct FractalApp {
help_open: bool,
/// Time-based animation of colours / Julia c / Phoenix p / zoom.
anim: AnimState,
/// Kind-switch morph in progress, if any.
morph: Option<KindMorph>,
/// Smoothed frames-per-second, recomputed each ~0.5 s window. Only advances
/// while the app is actually repainting (interaction / animation / export);
@@ -327,6 +366,14 @@ pub struct FractalApp {
ref_center_re: Big,
ref_center_im: Big,
ref_half_height: f64,
/// Kind and kind-switch morph the current `reference` was computed with.
/// The shader iterates with these (not the live kind/morph) so its delta
/// formula always matches the orbit, even while the worker lags a frame
/// behind — otherwise a kind switch flashes the new kind, unblended, for
/// the frame(s) before the morphed reference arrives. `None` until the
/// first reference lands.
ref_kind: Option<FractalKind>,
ref_morph: Option<(FractalKind, f32)>,
/// Parameters of the most recent reference request (drift baseline / dedupe).
last_request: Option<RequestKey>,
@@ -466,6 +513,7 @@ impl FractalApp {
info_open: false,
help_open: false,
anim: AnimState::default(),
morph: None,
fps: 0.0,
fps_frames: 0,
fps_window_start: 0.0,
@@ -474,6 +522,8 @@ impl FractalApp {
ref_center_re,
ref_center_im,
ref_half_height,
ref_kind: None,
ref_morph: None,
last_request: None,
#[cfg(not(target_arch = "wasm32"))]
worker: crate::worker::RefWorker::spawn(),
@@ -667,6 +717,7 @@ impl FractalApp {
) {
self.mode = FractalMode::Mandelbrot;
self.view = ViewState::with_center(cre, cim, half_height);
self.morph = None;
// Presets carry a hand-tuned count; don't let the auto-scaler clobber it.
self.auto_iterations = false;
self.max_iterations = iterations.clamp(32, MAX_REF_POINTS as u32 - 1);
@@ -706,6 +757,7 @@ impl FractalApp {
/// Restore a shared state into this app.
fn apply_share(&mut self, s: &ShareState) {
self.morph = None;
self.mode = if s.julia {
FractalMode::Julia
} else {
@@ -778,6 +830,7 @@ impl FractalApp {
kind: self.kind,
power: self.power,
complex_power: self.complex_power,
morph: self.morph.as_ref().map(|m| (m.from, m.weight())),
}
}
@@ -809,11 +862,17 @@ impl FractalApp {
|| key.kind != self.kind
|| key.power != self.power
|| key.complex_power != self.complex_power
|| key.morph != self.morph.as_ref().map(|m| (m.from, m.weight()))
{
return true;
}
// Lambda in Set mode is a static fractal; don't trigger recompute on center drift.
if self.kind == FractalKind::Lambda && matches!(self.mode, FractalMode::Mandelbrot) {
// Lambda in Set mode is a static fractal; don't trigger recompute on
// center drift. (Not while morphing: the other kind's formula does
// depend on the center.)
if self.kind == FractalKind::Lambda
&& matches!(self.mode, FractalMode::Mandelbrot)
&& self.morph.is_none()
{
// But still recompute on significant zoom changes for precision
let ratio = self.view.half_height / key.half_height;
return !(0.5..=2.0).contains(&ratio);
@@ -833,8 +892,18 @@ impl FractalApp {
[dre, dim]
}
fn apply_reference(&mut self, points: Vec<[f32; 2]>, cre: Big, cim: Big, hh: f64) {
fn apply_reference(
&mut self,
points: Vec<[f32; 2]>,
cre: Big,
cim: Big,
hh: f64,
kind: FractalKind,
morph: Option<(FractalKind, f32)>,
) {
self.reference = Arc::new(points);
self.ref_kind = Some(kind);
self.ref_morph = morph;
self.ref_center_re = cre;
self.ref_center_im = cim;
self.ref_half_height = hh;
@@ -865,7 +934,7 @@ impl FractalApp {
let max_iter = key.iter.min(MAX_REF_POINTS as u32 - 1);
// Lambda in Set mode has a static fractal centered at origin.
if key.kind == FractalKind::Lambda && !key.julia {
if key.kind == FractalKind::Lambda && !key.julia && key.morph.is_none() {
key.center_re = big_from_f64(0.0, precision);
key.center_im = big_from_f64(0.0, precision);
}
@@ -885,6 +954,7 @@ impl FractalApp {
phoenix_p: key.phoenix_p,
lambda_l: key.lambda_l,
complex_power: key.complex_power,
morph: key.morph,
});
self.pending = true;
}
@@ -905,6 +975,7 @@ impl FractalApp {
key.phoenix_p,
key.lambda_l,
key.complex_power,
key.morph.map(|(k, w)| (k, w as f64)),
)
} else {
compute_set_reference(
@@ -917,6 +988,7 @@ impl FractalApp {
key.phoenix_p,
key.lambda_l,
key.complex_power,
key.morph.map(|(k, w)| (k, w as f64)),
)
};
self.apply_reference(
@@ -924,6 +996,8 @@ impl FractalApp {
key.center_re.clone(),
key.center_im.clone(),
key.half_height,
key.kind,
key.morph,
);
}
@@ -932,7 +1006,14 @@ impl FractalApp {
#[cfg(not(target_arch = "wasm32"))]
if let Some(res) = self.worker.try_take_latest() {
self.apply_reference(res.points, res.center_re, res.center_im, res.half_height);
self.apply_reference(
res.points,
res.center_re,
res.center_im,
res.half_height,
res.kind,
res.morph,
);
self.pending = false;
}
}
@@ -954,7 +1035,7 @@ impl FractalApp {
let max_iter = key.iter;
// Lambda in Set mode has a static fractal centered at origin.
if key.kind == FractalKind::Lambda && !key.julia {
if key.kind == FractalKind::Lambda && !key.julia && key.morph.is_none() {
key.center_re = big_from_f64(0.0, precision);
key.center_im = big_from_f64(0.0, precision);
}
@@ -974,6 +1055,7 @@ impl FractalApp {
key.phoenix_p,
key.lambda_l,
key.complex_power,
key.morph.map(|(k, w)| (k, w as f64)),
)
} else {
compute_set_reference(
@@ -986,6 +1068,7 @@ impl FractalApp {
key.phoenix_p,
key.lambda_l,
key.complex_power,
key.morph.map(|(k, w)| (k, w as f64)),
)
};
self.apply_reference(
@@ -993,6 +1076,8 @@ impl FractalApp {
key.center_re.clone(),
key.center_im.clone(),
key.half_height,
key.kind,
key.morph,
);
self.last_request = Some(key);
}
@@ -1050,8 +1135,9 @@ impl FractalApp {
palette_id: self.palette,
shadow_palette_id: self.shadow_palette,
aa_level: if self.antialias { 2 } else { 1 },
kind: self.kind as u32,
kind: self.ref_kind.unwrap_or(self.kind) as u32,
power: self.power,
morph_from: self.ref_morph.map_or(0, |(k, _)| k as u32),
dc_offset: self.dc_offset(),
phoenix_p: [self.phoenix_p.0 as f32, self.phoenix_p.1 as f32],
lambda_l: [self.lambda_l.0 as f32, self.lambda_l.1 as f32],
@@ -1059,6 +1145,7 @@ impl FractalApp {
de_coloring: (self.de_coloring || mode > 0) as u32,
rendering_mode: mode,
camera_direction: self.camera.direction(self.anim.camera_state).to_array(),
morph_w: self.ref_morph.map_or(0.0, |(_, w)| w),
camera_inv_proj: self
.camera
.orthographic(self.anim.camera_state)
@@ -1067,7 +1154,6 @@ impl FractalApp {
screen_dim: self.screen_dim,
light_count: gpu_lights(&self.lights).1,
cm_coef: complex_binomials(self.complex_power),
_pad: [0; _],
_pad3: [0; _],
}
}
@@ -1586,6 +1672,22 @@ impl FractalApp {
let p = self.anim.camera_progress;
self.anim.camera_state = p * p * (3.0 - 2.0 * p);
// Kind-switch morph: advance the per-iteration formula blend, and glide
// the camera to the new kind's default view unless the user took over.
if let Some(m) = &mut self.morph {
m.progress += dt as f32 / self.anim.kind_morph_duration.max(0.05);
if m.camera {
self.view = interpolate_view(&m.from_view, &m.to_view, m.eased() as f64);
}
if m.progress >= 1.0 {
if m.camera {
self.view = m.to_view.clone();
}
self.morph = None;
}
ui.ctx().request_repaint();
}
if !(self.anim.color || self.anim.zoom || julia_on || phoenix_on || lambda_on) {
return;
}
@@ -1701,9 +1803,23 @@ impl FractalApp {
});
}
if self.kind != prev_kind {
self.view = Self::default_view_for(self.mode, self.kind);
let to_view = Self::default_view_for(self.mode, self.kind);
// Buddhabrot has its own pipeline without the blended formula, so
// it keeps the instant switch.
self.morph =
(self.anim.kind_morph && self.mode != FractalMode::Buddhabrot).then(|| KindMorph {
from: prev_kind,
progress: 0.0,
from_view: self.view.clone(),
to_view: to_view.clone(),
camera: true,
});
if self.morph.is_none() {
self.view = to_view;
}
}
let prev_mode = self.mode;
ui.horizontal(|ui| {
ui.radio_value(&mut self.mode, FractalMode::Mandelbrot, "Set");
ui.radio_value(&mut self.mode, FractalMode::Julia, "Julia");
@@ -1714,6 +1830,9 @@ impl FractalApp {
progressively sharpens while the view stays still.",
);
});
if self.mode != prev_mode {
self.morph = None;
}
if self.mode == FractalMode::Buddhabrot {
self.buddhabrot_ui(ui);
@@ -1722,6 +1841,7 @@ impl FractalApp {
ui.add_space(4.);
if ui.button("Reset view").clicked() {
self.view = Self::default_view_for(self.mode, self.kind);
self.morph = None;
}
ui.add_space(8.0);
ui.small("Drag to pan · scroll to zoom toward the cursor");
@@ -1880,6 +2000,19 @@ impl FractalApp {
);
}
ui.checkbox(&mut self.anim.kind_morph, "Morph kind switch")
.on_hover_text(
"When picking another fractal, blend the old and new formulas \
at every iteration step and glide to the new default view.",
);
if self.anim.kind_morph {
ui.add(
egui::Slider::new(&mut self.anim.kind_morph_duration, 0.2..=10.0)
.text("morph s")
.logarithmic(true),
);
}
// Julia c only affects Julia mode; Phoenix p only the Phoenix kind.
if self.mode == FractalMode::Julia {
if ui.checkbox(&mut self.anim.julia, "Morph c").changed() && self.anim.julia {
@@ -2064,6 +2197,7 @@ impl FractalApp {
ui.add_space(4.);
if ui.button("Reset view").clicked() {
self.view = Self::default_view_for(self.mode, self.kind);
self.morph = None;
}
ui.add_space(8.0);
ui.small("Drag to pan · scroll to zoom toward the cursor");
@@ -2302,6 +2436,7 @@ impl FractalApp {
if ui.input(|i| i.key_pressed(egui::Key::R)) {
self.view = Self::default_view_for(self.mode, self.kind);
self.morph = None;
self.camera = Camera::new();
interacted = true;
}
@@ -2379,6 +2514,11 @@ impl FractalApp {
let now = ui.input(|i| i.time);
if interacted {
self.last_interact_time = now;
// The user is steering the camera: stop the kind-switch morph from
// overriding it (the formula blend itself carries on).
if let Some(m) = &mut self.morph {
m.camera = false;
}
}
let interacting = now - self.last_interact_time < INTERACT_SETTLE;
if interacting {
@@ -2404,8 +2544,12 @@ impl FractalApp {
self.screen_dim = [rect.width(), rect.height()];
self.camera.set_aspect_ratio(aspect as f32);
let mut uniforms = self.make_uniforms(aspect);
if interacting {
uniforms.aa_level = 1; // supersampling is wasted on the low-res pass
// Supersampling is wasted on the low-res pass, and on a kind-switch
// morph (every frame re-iterates, and the blend moves on next frame).
// `ref_morph` too: the last morphed reference outlives `morph` by a
// frame or so, until the worker delivers the plain one.
if interacting || self.morph.is_some() || self.ref_morph.is_some() {
uniforms.aa_level = 1;
}
ui.painter().add(egui_wgpu::Callback::new_paint_callback(
rect,
+364 -164
View File
@@ -9,6 +9,11 @@
//! The `(z0, c)` form serves both set types:
//! * Mandelbrot-set: `z0 = 0`, `c = view center` (the c-plane point per pixel).
//! * Julia-set: `z0 = view center`, `c = fractal constant` (fixed per view).
//!
//! While switching fractal kinds, the formula is morphed *per iteration*:
//! `Z_{n+1} = (1 - w)·f_kind(Z_n) + w·f_from(Z_n)` (see `morph` below). The
//! map is linear in the two outputs, so the GPU delta is the same blend of
//! the two kinds' deltas and perturbation/rebasing keep working unchanged.
use super::kind::FractalKind;
use crate::view::{Big, big_from_f64};
@@ -37,6 +42,9 @@ const F64_MAX_PRECISION: usize = 80;
/// Compute the reference orbit `Z_0..Z_{len-1}` where `Z_0 = z0` and
/// `Z_{n+1} = f(Z_n, c)` for the given `kind` (and `power`, for Multibrot), up
/// to `max_iter` steps at `precision` bits. Each entry is `[re, im]` in f32.
///
/// `morph = Some((from, w))` blends in a second kind's formula at every step:
/// `(1 - w)·f_kind + w·f_from` (used by the kind-switch animation).
#[allow(clippy::too_many_arguments)]
pub fn compute_reference(
z0_re: &Big,
@@ -50,20 +58,9 @@ pub fn compute_reference(
phoenix_p: (f64, f64),
lambda_l: (f64, f64),
complex_power: (f64, f64),
morph: Option<(FractalKind, f64)>,
) -> Vec<[f32; 2]> {
if precision <= F64_MAX_PRECISION {
return compute_reference_f64(
(z0_re.to_f64().value(), z0_im.to_f64().value()),
(c_re.to_f64().value(), c_im.to_f64().value()),
max_iter,
kind,
power,
phoenix_p,
lambda_l,
complex_power,
);
}
compute_reference_big(
compute_reference_inner(
z0_re,
z0_im,
c_re,
@@ -75,28 +72,86 @@ pub fn compute_reference(
phoenix_p,
lambda_l,
complex_power,
morph,
false,
)
}
/// [`compute_reference`]'s fast path for shallow views (see
/// [`F64_MAX_PRECISION`]): the same per-kind formulas in plain `f64`.
/// [`compute_reference`] plus `set_plane` (see [`StepConsts::set_plane`]):
/// picks the `f64` fast path or the `FBig` path by precision.
#[allow(clippy::too_many_arguments)]
fn compute_reference_f64(
z0: (f64, f64),
c: (f64, f64),
fn compute_reference_inner(
z0_re: &Big,
z0_im: &Big,
c_re: &Big,
c_im: &Big,
max_iter: u32,
precision: usize,
kind: FractalKind,
power: u32,
phoenix_p: (f64, f64),
lambda_l: (f64, f64),
complex_power: (f64, f64),
morph: Option<(FractalKind, f64)>,
set_plane: bool,
) -> Vec<[f32; 2]> {
// A zero-weight morph is just the plain kind; skip the second formula.
let morph = morph.filter(|&(_, w)| w != 0.0);
if precision <= F64_MAX_PRECISION {
let k = StepConstsF64 {
c: (c_re.to_f64().value(), c_im.to_f64().value()),
p: phoenix_p,
l: lambda_l,
cpow: complex_power,
power,
set_plane,
};
return compute_reference_f64(
(z0_re.to_f64().value(), z0_im.to_f64().value()),
max_iter,
kind,
&k,
morph,
);
}
let k = StepConsts {
cr: c_re.clone().with_precision(precision).value(),
ci: c_im.clone().with_precision(precision).value(),
pr: big_from_f64(phoenix_p.0, precision),
pi: big_from_f64(phoenix_p.1, precision),
lr: big_from_f64(lambda_l.0, precision),
li: big_from_f64(lambda_l.1, precision),
cpow_re: big_from_f64(complex_power.0, precision),
cpow_im: big_from_f64(complex_power.1, precision),
power,
precision,
set_plane,
};
compute_reference_big(z0_re, z0_im, max_iter, kind, &k, morph)
}
/// `f64` twin of [`StepConsts`].
struct StepConstsF64 {
c: (f64, f64),
p: (f64, f64),
l: (f64, f64),
cpow: (f64, f64),
power: u32,
set_plane: bool,
}
/// [`compute_reference`]'s fast path for shallow views (see
/// [`F64_MAX_PRECISION`]): the same per-kind formulas in plain `f64`.
fn compute_reference_f64(
z0: (f64, f64),
max_iter: u32,
kind: FractalKind,
k: &StepConstsF64,
morph: Option<(FractalKind, f64)>,
) -> Vec<[f32; 2]> {
let (cr, ci) = c;
let (mut zr, mut zi) = z0;
// Previous iterate, for the Phoenix two-term recurrence (Y_{-1} = 0).
let (mut zr_prev, mut zi_prev) = (0.0f64, 0.0f64);
let (pr, pi) = phoenix_p;
let (lr, li) = lambda_l;
let mut prev = (0.0f64, 0.0f64);
let mut points: Vec<[f32; 2]> = Vec::with_capacity(max_iter as usize + 1);
for _ in 0..=max_iter {
@@ -105,41 +160,68 @@ fn compute_reference_f64(
break;
}
let (new_zr, new_zi) = match kind {
FractalKind::Mandelbrot => ((zr + zi) * (zr - zi) + cr, 2.0 * zr * zi + ci),
FractalKind::BurningShip => (zr * zr - zi * zi + cr, (2.0 * zr * zi).abs() + ci),
FractalKind::Tricorn => (zr * zr - zi * zi + cr, ci - 2.0 * zr * zi),
FractalKind::Multibrot => {
let (mut rr, mut ri) = (1.0f64, 0.0f64);
for _ in 0..power.max(2) {
(rr, ri) = (rr * zr - ri * zi, rr * zi + ri * zr);
}
(rr + cr, ri + ci)
}
FractalKind::Celtic => ((zr * zr - zi * zi).abs() + cr, 2.0 * zr * zi + ci),
FractalKind::Perpendicular => (zr * zr - zi * zi + cr, ci - 2.0 * zr * zi.abs()),
FractalKind::Buffalo => ((zr * zr - zi * zi).abs() + cr, ci - (2.0 * zr * zi).abs()),
FractalKind::Phoenix => (
zr * zr - zi * zi + cr + (pr * zr_prev - pi * zi_prev),
2.0 * zr * zi + ci + (pr * zi_prev + pi * zr_prev),
),
FractalKind::Lambda => {
// λ·z(1 - z).
let (re2, im2) = (1.0 - zr, -zi);
let (lzr, lzi) = (lr * zr - li * zi, lr * zi + li * zr);
(lzr * re2 - lzi * im2, re2 * lzi + lzr * im2)
}
FractalKind::ComplexMultibrot => {
let (pr, pi) = complex_pow_complex_f64(zr, zi, complex_power.0, complex_power.1);
(pr + cr, pi + ci)
}
};
(zr_prev, zi_prev) = (zr, zi);
let (mut new_zr, mut new_zi) = step_f64(kind, k, zr, zi, prev);
if let Some((from, w)) = morph {
let (br, bi) = step_f64(from, k, zr, zi, prev);
new_zr += w * (br - new_zr);
new_zi += w * (bi - new_zi);
}
prev = (zr, zi);
(zr, zi) = (new_zr, new_zi);
}
points
}
/// `f64` twin of [`step`]: one step `f(Z_n)` of `kind`'s formula.
fn step_f64(
kind: FractalKind,
k: &StepConstsF64,
zr: f64,
zi: f64,
prev: (f64, f64),
) -> (f64, f64) {
let (cr, ci) = k.c;
match kind {
FractalKind::Mandelbrot => ((zr + zi) * (zr - zi) + cr, 2.0 * zr * zi + ci),
FractalKind::BurningShip => (zr * zr - zi * zi + cr, (2.0 * zr * zi).abs() + ci),
FractalKind::Tricorn => (zr * zr - zi * zi + cr, ci - 2.0 * zr * zi),
FractalKind::Multibrot => {
let (mut rr, mut ri) = (1.0f64, 0.0f64);
for _ in 0..k.power.max(2) {
(rr, ri) = (rr * zr - ri * zi, rr * zi + ri * zr);
}
(rr + cr, ri + ci)
}
FractalKind::Celtic => ((zr * zr - zi * zi).abs() + cr, 2.0 * zr * zi + ci),
FractalKind::Perpendicular => (zr * zr - zi * zi + cr, ci - 2.0 * zr * zi.abs()),
FractalKind::Buffalo => ((zr * zr - zi * zi).abs() + cr, ci - (2.0 * zr * zi).abs()),
FractalKind::Phoenix => {
let (pr, pi) = k.p;
let (zr_prev, zi_prev) = prev;
(
zr * zr - zi * zi + cr + (pr * zr_prev - pi * zi_prev),
2.0 * zr * zi + ci + (pr * zi_prev + pi * zr_prev),
)
}
FractalKind::Lambda => {
// λ·z(1 - z) (+ c on the parameter plane).
let (lr, li) = k.l;
let (re2, im2) = (1.0 - zr, -zi);
let (lzr, lzi) = (lr * zr - li * zi, lr * zi + li * zr);
let (re, im) = (lzr * re2 - lzi * im2, re2 * lzi + lzr * im2);
if k.set_plane {
(re + cr, im + ci)
} else {
(re, im)
}
}
FractalKind::ComplexMultibrot => {
let (pr, pi) = complex_pow_complex_f64(zr, zi, k.cpow.0, k.cpow.1);
(pr + cr, pi + ci)
}
}
}
/// `f64` twin of [`complex_pow_complex`] (principal branch, `0^p = 0`).
fn complex_pow_complex_f64(zr: f64, zi: f64, pr: f64, pi: f64) -> (f64, f64) {
if zr == 0.0 && zi == 0.0 {
@@ -152,38 +234,44 @@ fn complex_pow_complex_f64(zr: f64, zi: f64, pr: f64, pi: f64) -> (f64, f64) {
(mag * cos_a, mag * sin_a)
}
/// Everything a single formula step needs besides the orbit state, converted
/// to `Big` once up front.
struct StepConsts {
cr: Big,
ci: Big,
/// Phoenix distortion constant `p`.
pr: Big,
pi: Big,
/// Lambda distortion constant `l`.
lr: Big,
li: Big,
/// Complex Multibrot exponent.
cpow_re: Big,
cpow_im: Big,
power: u32,
precision: usize,
/// Parameter plane: the GPU adds `dc` every step for every kind, so the
/// Lambda map (which has no `c` of its own) is `λ·z(1 - z) + c` there.
set_plane: bool,
}
/// [`compute_reference`] at arbitrary precision (`FBig`), for deep views.
#[allow(clippy::too_many_arguments)]
fn compute_reference_big(
z0_re: &Big,
z0_im: &Big,
c_re: &Big,
c_im: &Big,
max_iter: u32,
precision: usize,
kind: FractalKind,
power: u32,
phoenix_p: (f64, f64),
lambda_l: (f64, f64),
complex_power: (f64, f64),
k: &StepConsts,
morph: Option<(FractalKind, f64)>,
) -> Vec<[f32; 2]> {
let cr = c_re.clone().with_precision(precision).value();
let ci = c_im.clone().with_precision(precision).value();
let precision = k.precision;
let morph = morph.map(|(from, w)| (from, big_from_f64(w, precision)));
let mut zr = z0_re.clone().with_precision(precision).value();
let mut zi = z0_im.clone().with_precision(precision).value();
// Previous iterate, for the Phoenix two-term recurrence (Y_{-1} = 0).
let mut zr_prev = big_zero(precision);
let mut zi_prev = big_zero(precision);
// Phoenix distortion constant `p` (a small fixed complex number).
let pr = big_from_f64(phoenix_p.0, precision);
let pi = big_from_f64(phoenix_p.1, precision);
// Lambda distortion constant `l` (a small fixed complex number).
let lr = big_from_f64(lambda_l.0, precision);
let li = big_from_f64(lambda_l.1, precision);
// Complex Multibrot exponent (a fixed complex number).
let cpow_re = big_from_f64(complex_power.0, precision);
let cpow_im = big_from_f64(complex_power.1, precision);
let mut points: Vec<[f32; 2]> = Vec::with_capacity(max_iter as usize + 1);
@@ -197,73 +285,13 @@ fn compute_reference_big(
break;
}
let (new_zr, new_zi) = match kind {
FractalKind::Mandelbrot => {
// Z^2 = (zr^2 - zi^2) + (2 zr zi) i, with zr^2 - zi^2 as
// (zr + zi)(zr - zi): one multiply instead of two squares.
let re = (&zr + &zi) * (&zr - &zi) + &cr;
let im = ((&zr * &zi) << 1) + &ci; // << 1 is exact ×2 in base 2
(re, im)
}
FractalKind::BurningShip => {
// (|zr| + i|zi|)^2 = (zr^2 - zi^2) + 2|zr zi| i.
let re = &zr.sqr() - &zi.sqr() + &cr;
let im = big_abs((&zr * &zi) << 1) + &ci;
(re, im)
}
FractalKind::Tricorn => {
// conj(z)^2 = (zr^2 - zi^2) - 2 zr zi i.
let re = &zr.sqr() - &zi.sqr() + &cr;
let im = &ci - ((&zr * &zi) << 1);
(re, im)
}
FractalKind::Multibrot => {
let (pr, pi) = complex_pow(&zr, &zi, power.max(2), precision);
(pr + &cr, pi + &ci)
}
FractalKind::Celtic => {
// |Re(z^2)| + i·Im(z^2): abs the real output of the square.
let re = big_abs(&zr.sqr() - &zi.sqr()) + &cr;
let im = ((&zr * &zi) << 1) + &ci;
(re, im)
}
FractalKind::Perpendicular => {
// (x^2 - y^2) - 2·x·|y| i: abs the imaginary input.
let re = &zr.sqr() - &zi.sqr() + &cr;
let im = if zi.to_f64().value() < 0.0 {
&ci + ((&zr * &zi) << 1)
} else {
&ci - ((&zr * &zi) << 1)
};
(re, im)
}
FractalKind::Buffalo => {
// |Re(z^2)| - |Im(z^2)| i: abs both outputs.
let re = big_abs(&zr.sqr() - &zi.sqr()) + &cr;
let im = &ci - big_abs((&zr * &zi) << 1);
(re, im)
}
FractalKind::Phoenix => {
// z^2 + c + p·z_{n-1}.
let re2 = &zr.sqr() - &zi.sqr();
let im2 = (&zr * &zi) << 1;
let pzr = &pr * &zr_prev - &pi * &zi_prev;
let pzi = &pr * &zi_prev + &pi * &zr_prev;
(re2 + &cr + pzr, im2 + &ci + pzi)
}
FractalKind::Lambda => {
// λ·z(1 - z): logistic map.
let re2 = 1 - &zr;
let im2 = -&zi;
let lzr = &lr * &zr - &li * &zi;
let lzi = &lr * &zi + &li * &zr;
(&lzr * &re2 - &lzi * &im2, re2 * lzi + lzr * im2)
}
FractalKind::ComplexMultibrot => {
let (pr, pi) = complex_pow_complex(&zr, &zi, &cpow_re, &cpow_im, precision);
(pr + &cr, pi + &ci)
}
};
let (mut new_zr, mut new_zi) = step(kind, k, &zr, &zi, &zr_prev, &zi_prev);
if let Some((from, w)) = &morph {
// (1 - w)·a + w·b = a + w·(b - a).
let (br, bi) = step(*from, k, &zr, &zi, &zr_prev, &zi_prev);
new_zr = &new_zr + &(w * &(br - &new_zr));
new_zi = &new_zi + &(w * &(bi - &new_zi));
}
// Shift the previous iterate (only the Phoenix arm reads it).
zr_prev = zr;
@@ -275,6 +303,92 @@ fn compute_reference_big(
points
}
/// One step `f(Z_n)` of `kind`'s formula (including its `+ c`), given the
/// current and previous iterate.
fn step(
kind: FractalKind,
k: &StepConsts,
zr: &Big,
zi: &Big,
zr_prev: &Big,
zi_prev: &Big,
) -> (Big, Big) {
let (cr, ci) = (&k.cr, &k.ci);
match kind {
FractalKind::Mandelbrot => {
// Z^2 = (zr^2 - zi^2) + (2 zr zi) i, with zr^2 - zi^2 as
// (zr + zi)(zr - zi): one multiply instead of two squares.
let re = (zr + zi) * (zr - zi) + cr;
let im = ((zr * zi) << 1) + ci; // << 1 is exact ×2 in base 2
(re, im)
}
FractalKind::BurningShip => {
// (|zr| + i|zi|)^2 = (zr^2 - zi^2) + 2|zr zi| i.
let re = &zr.sqr() - &zi.sqr() + cr;
let im = big_abs((zr * zi) << 1) + ci;
(re, im)
}
FractalKind::Tricorn => {
// conj(z)^2 = (zr^2 - zi^2) - 2 zr zi i.
let re = &zr.sqr() - &zi.sqr() + cr;
let im = ci - ((zr * zi) << 1);
(re, im)
}
FractalKind::Multibrot => {
let (pr, pi) = complex_pow(zr, zi, k.power.max(2), k.precision);
(pr + cr, pi + ci)
}
FractalKind::Celtic => {
// |Re(z^2)| + i·Im(z^2): abs the real output of the square.
let re = big_abs(&zr.sqr() - &zi.sqr()) + cr;
let im = ((zr * zi) << 1) + ci;
(re, im)
}
FractalKind::Perpendicular => {
// (x^2 - y^2) - 2·x·|y| i: abs the imaginary input.
let re = &zr.sqr() - &zi.sqr() + cr;
let im = if zi.to_f64().value() < 0.0 {
ci + ((zr * zi) << 1)
} else {
ci - ((zr * zi) << 1)
};
(re, im)
}
FractalKind::Buffalo => {
// |Re(z^2)| - |Im(z^2)| i: abs both outputs.
let re = big_abs(&zr.sqr() - &zi.sqr()) + cr;
let im = ci - big_abs((zr * zi) << 1);
(re, im)
}
FractalKind::Phoenix => {
// z^2 + c + p·z_{n-1}.
let re2 = &zr.sqr() - &zi.sqr();
let im2 = (zr * zi) << 1;
let pzr = &k.pr * zr_prev - &k.pi * zi_prev;
let pzi = &k.pr * zi_prev + &k.pi * zr_prev;
(re2 + cr + pzr, im2 + ci + pzi)
}
FractalKind::Lambda => {
// λ·z(1 - z): logistic map (+ c on the parameter plane).
let re2 = 1 - zr;
let im2 = -zi;
let lzr = &k.lr * zr - &k.li * zi;
let lzi = &k.lr * zi + &k.li * zr;
let re = &lzr * &re2 - &lzi * &im2;
let im = re2 * lzi + lzr * im2;
if k.set_plane {
(re + cr, im + ci)
} else {
(re, im)
}
}
FractalKind::ComplexMultibrot => {
let (pr, pi) = complex_pow_complex(zr, zi, &k.cpow_re, &k.cpow_im, k.precision);
(pr + cr, pi + ci)
}
}
}
fn big_zero(precision: usize) -> Big {
Big::from(0i32).with_precision(precision).value()
}
@@ -338,9 +452,10 @@ pub fn compute_set_reference(
phoenix_p: (f64, f64),
lambda_l: (f64, f64),
complex_power: (f64, f64),
morph: Option<(FractalKind, f64)>,
) -> Vec<[f32; 2]> {
let zero = big_zero(precision);
compute_reference(
compute_reference_inner(
&zero,
&zero,
center_re,
@@ -352,6 +467,8 @@ pub fn compute_set_reference(
phoenix_p,
lambda_l,
complex_power,
morph,
true,
)
}
@@ -375,6 +492,7 @@ mod tests {
(0.0, 0.0),
(0.0, 0.0),
(0.0, 0.0),
None,
);
// Independent naive f64 orbit.
@@ -401,37 +519,42 @@ mod tests {
}
/// The f64 fast path (shallow views) must produce the same orbit as the
/// arbitrary-precision path, for every kind, in both planes.
/// arbitrary-precision path, for every kind, in both planes, with and
/// without a kind-switch morph.
#[test]
fn f64_fast_path_matches_big() {
let bits_fast = F64_MAX_PRECISION;
let bits_big = F64_MAX_PRECISION + 64;
for kind in FractalKind::ALL {
for julia in [false, true] {
let run = |bits: usize| {
let (a, b) = (big_from_f64(-0.3, bits), big_from_f64(0.2, bits));
let (jr, ji) = (big_from_f64(-0.4, bits), big_from_f64(0.55, bits));
let args = (60, bits, kind, 3, (0.1, -0.2), (0.9, 0.3), (2.3, 0.4));
if julia {
compute_reference(
&a, &b, &jr, &ji, args.0, args.1, args.2, args.3, args.4, args.5,
args.6,
)
} else {
compute_set_reference(
&a, &b, args.0, args.1, args.2, args.3, args.4, args.5, args.6,
)
}
};
let (fast, big) = (run(bits_fast), run(bits_big));
assert_eq!(fast.len(), big.len(), "{kind:?} julia={julia}: length");
for (i, (f, b)) in fast.iter().zip(&big).enumerate() {
for k in 0..2 {
let tol = 1e-5 * (1.0 + b[k].abs());
assert!(
(f[k] - b[k]).abs() <= tol,
"{kind:?} julia={julia}: point {i} {f:?} vs {b:?}"
);
for morph in [None, Some((FractalKind::Phoenix, 0.3))] {
let run = |bits: usize| {
let (a, b) = (big_from_f64(-0.3, bits), big_from_f64(0.2, bits));
let (jr, ji) = (big_from_f64(-0.4, bits), big_from_f64(0.55, bits));
let args = (60, bits, kind, 3, (0.1, -0.2), (0.9, 0.3), (2.3, 0.4));
if julia {
compute_reference(
&a, &b, &jr, &ji, args.0, args.1, args.2, args.3, args.4, args.5,
args.6, morph,
)
} else {
compute_set_reference(
&a, &b, args.0, args.1, args.2, args.3, args.4, args.5, args.6,
morph,
)
}
};
let ctx = format!("{kind:?} julia={julia} morph={morph:?}");
let (fast, big) = (run(bits_fast), run(bits_big));
assert_eq!(fast.len(), big.len(), "{ctx}: length");
for (i, (f, b)) in fast.iter().zip(&big).enumerate() {
for k in 0..2 {
let tol = 1e-5 * (1.0 + b[k].abs());
assert!(
(f[k] - b[k]).abs() <= tol,
"{ctx}: point {i} {f:?} vs {b:?}"
);
}
}
}
}
@@ -453,6 +576,7 @@ mod tests {
(0.0, 0.0),
(0.0, 0.0),
(0.0, 0.0),
None,
);
assert_eq!(points.len(), 501, "interior orbit should not escape");
}
@@ -472,6 +596,7 @@ mod tests {
(0.0, 0.0),
(0.0, 0.0),
(0.0, 0.0),
None,
);
let (c_re, c_im) = (-1.75_f64, -0.03_f64);
@@ -502,6 +627,7 @@ mod tests {
(0.0, 0.0),
(0.0, 0.0),
(0.0, 0.0),
None,
);
let (c_re, c_im) = (0.3_f64, 0.2_f64);
@@ -538,6 +664,7 @@ mod tests {
(0.0, 0.0),
(0.0, 0.0),
(0.0, 0.0),
None,
);
let (mut zr, mut zi) = (0.15_f64, -0.1_f64);
@@ -568,6 +695,7 @@ mod tests {
(0.0, 0.0),
(0.0, 0.0),
(0.0, 0.0),
None,
);
let (c_re, c_im) = (-0.6_f64, 0.4_f64);
@@ -599,6 +727,7 @@ mod tests {
(0.0, 0.0),
(0.0, 0.0),
(0.0, 0.0),
None,
);
let (c_re, c_im) = (-0.7_f64, -0.2_f64);
@@ -630,6 +759,7 @@ mod tests {
(0.0, 0.0),
(0.0, 0.0),
(0.0, 0.0),
None,
);
let (c_re, c_im) = (-1.2_f64, -0.35_f64);
@@ -662,6 +792,7 @@ mod tests {
p,
(0.0, 0.0),
(0.0, 0.0),
None,
);
let (c_re, c_im) = (0.5667_f64, 0.0_f64);
@@ -700,6 +831,7 @@ mod tests {
(0.0, 0.0),
(0.0, 0.0),
power,
None,
);
// Naive f64 complex power via z^p = exp(p * ln z), ln z = ln|z| + i*arg(z).
@@ -728,4 +860,72 @@ mod tests {
zi = nzi;
}
}
fn set_ref(
cr: f64,
ci: f64,
kind: FractalKind,
morph: Option<(FractalKind, f64)>,
) -> Vec<[f32; 2]> {
compute_set_reference(
&Big::try_from(cr).unwrap(),
&Big::try_from(ci).unwrap(),
60,
200,
kind,
2,
(0.0, 0.0),
(0.0, 0.0),
(0.0, 0.0),
morph,
)
}
/// Morph weight 0 is the plain kind; weight 1 is entirely the from-kind.
#[test]
fn morph_endpoints_match_plain_kinds() {
let (cr, ci) = (-1.75, -0.03);
let ship = set_ref(cr, ci, FractalKind::BurningShip, None);
let mandel = set_ref(cr, ci, FractalKind::Mandelbrot, None);
let w0 = set_ref(
cr,
ci,
FractalKind::BurningShip,
Some((FractalKind::Mandelbrot, 0.0)),
);
let w1 = set_ref(
cr,
ci,
FractalKind::BurningShip,
Some((FractalKind::Mandelbrot, 1.0)),
);
assert_eq!(w0, ship);
assert_eq!(w1, mandel);
}
/// A half-way Mandelbrot / Burning Ship morph matches a naive f64
/// iteration of the per-step blend.
#[test]
fn morph_blend_matches_naive_f64() {
let (c_re, c_im) = (-0.6_f64, 0.3_f64);
let w = 0.5_f64;
let points = set_ref(
c_re,
c_im,
FractalKind::Mandelbrot,
Some((FractalKind::BurningShip, w)),
);
let (mut zr, mut zi) = (0.0_f64, 0.0_f64);
for point in &points {
let tol = 1e-4 * (1.0 + zr.abs().max(zi.abs()));
assert!((point[0] as f64 - zr).abs() < tol, "re: {point:?} vs {zr}");
assert!((point[1] as f64 - zi).abs() < tol, "im: {point:?} vs {zi}");
let re = zr * zr - zi * zi + c_re; // identical for both kinds
let im_m = 2.0 * zr * zi + c_im;
let im_b = 2.0 * (zr * zi).abs() + c_im;
zr = re;
zi = (1.0 - w) * im_m + w * im_b;
}
}
}
+15 -3
View File
@@ -42,6 +42,8 @@ fn geom_differs(a: &Uniforms, b: &Uniforms) -> bool {
|| a.dc_offset != b.dc_offset
|| a.phoenix_p != b.phoenix_p
|| a.lambda_l != b.lambda_l
|| a.morph_from != b.morph_from
|| a.morph_w != b.morph_w
|| a.de_coloring != b.de_coloring
// The iterate pass's DE clamp (`max_de`) depends on whether any
// shadow-style mode is on.
@@ -72,6 +74,8 @@ pub struct PipelineKey {
kind: u32,
julia: bool,
de: bool,
/// A kind-switch morph is in progress (`morph_w > 0`).
morph: bool,
}
impl PipelineKey {
@@ -80,14 +84,16 @@ impl PipelineKey {
kind: u.kind,
julia: u.is_julia != 0,
de: u.de_coloring != 0,
morph: u.morph_w > 0.0,
}
}
fn constants(&self) -> [(&'static str, f64); 3] {
fn constants(&self) -> [(&'static str, f64); 4] {
[
("KIND", self.kind as f64),
("IS_JULIA", self.julia as u32 as f64),
("DE", self.de as u32 as f64),
("MORPH", self.morph as u32 as f64),
]
}
}
@@ -166,7 +172,9 @@ pub struct Uniforms {
pub kind: u32,
/// Exponent for the Multibrot kind.
pub power: u32,
pub _pad: [u32; 1],
/// Kind-switch morph: the kind being blended *from* (a `FractalKind`
/// discriminant); only read when `morph_w > 0`.
pub morph_from: u32,
/// Complex offset of the view center from the reference center, so a stale
/// or reused reference (computed at a slightly different center) still maps
/// correctly. Added to every pixel's per-pixel offset.
@@ -194,7 +202,11 @@ pub struct Uniforms {
pub camera_inv_proj: [f32; 16],
/// Screen dimension
pub screen_dim: [f32; 2],
pub _pad3: [u32; 2],
/// Kind-switch morph weight: each step is `(1 - w)·f_kind + w·f_from`.
/// 0 = no morph (and the iteration pipeline is then specialized without
/// the morph path, see [`PipelineKey`]).
pub morph_w: f32,
pub _pad3: [u32; 1],
/// Complex binomial coefficients `C(complex_power, k)`, k = 1..16, two per
/// row (odd k in `[0..2]`, even k in `[2..4]`), for the Complex Multibrot
/// delta series. Derived from `complex_power` alone.
+5
View File
@@ -19,6 +19,8 @@ struct Uniforms {
kind: u32,
// Exponent for the Multibrot kind.
power: u32,
// Kind-switch morph: the kind blended *from* (see morph_w).
morph_from: u32,
dc_offset: vec2<f32>,
// Distortion constant p for the Phoenix map (z^2 + c + p*z_{n-1}); unused
// by other kinds. Placed by dc_offset so both vec2s stay 8-byte aligned.
@@ -43,6 +45,9 @@ struct Uniforms {
camera_inv_proj: mat4x4<f32>,
// Screen dimensions
screen_dim: vec2<f32>,
// Kind-switch morph weight: each step is (1 - w)*f_kind + w*f_morph_from;
// 0 = no morph. Only read by MORPH pipelines (see mandelbrot.wgsl).
morph_w: f32,
// Complex binomial coefficients C(complex_power, k) for k = 1..16, two per
// vec4 (k odd in .xy, k even in .zw), for the Complex Multibrot delta
// series. Precomputed on the CPU since they only depend on the power.
+66 -24
View File
@@ -31,6 +31,10 @@
override KIND: u32 = 0u;
override IS_JULIA: bool = false;
override DE: bool = false;
// A kind-switch morph is in progress (`u.morph_w > 0`): each step blends in a
// second kind, `u.morph_from`. That one is a runtime value (it only lives for
// the length of the animation), so only MORPH pipelines pay for its branches.
override MORPH: bool = false;
struct VsOut {
@builtin(position) pos: vec4<f32>,
@@ -137,11 +141,11 @@ fn complex_multibrot_delta(z: vec2<f32>, e: vec2<f32>, p: vec2<f32>) -> vec2<f32
return cpow(z + e, p) - cpow(z, p);
}
// One perturbation step of the current fractal's delta: e -> f(Z+e) - f(Z),
// where `z` is the reference orbit value X_m. `step_add` (dc) is added by the
// caller. Must match `FractalKind` on the CPU side.
fn advance_delta(z: vec2<f32>, e: vec2<f32>) -> vec2<f32> {
if KIND == KIND_BURNING_SHIP {
// One perturbation step of `kind`'s delta: e -> f(Z+e) - f(Z), where `z` is
// the reference orbit value X_m. `step_add` (dc) is added by the caller. Must
// match `FractalKind` on the CPU side.
fn advance_delta_kind(kind: u32, z: vec2<f32>, e: vec2<f32>) -> vec2<f32> {
if kind == KIND_BURNING_SHIP {
// (|x| + i|y|)^2 has real part x^2 - y^2 (an ordinary square delta) and
// imaginary part 2|x y|. The imaginary delta is 2(|x y| - |X Y|); diffabs
// computes it exactly, even where the product x y changes sign — which the
@@ -150,34 +154,34 @@ fn advance_delta(z: vec2<f32>, e: vec2<f32>) -> vec2<f32> {
let base = 2.0 * cmul(z, e) + cmul(e, e);
let dp = z.x * e.y + z.y * e.x + e.x * e.y;
return vec2<f32>(base.x, 2.0 * diffabs(z.x * z.y, dp));
} else if KIND == KIND_TRICORN {
} else if kind == KIND_TRICORN {
let cz = conj(z);
let ce = conj(e);
return 2.0 * cmul(cz, ce) + cmul(ce, ce);
} else if KIND == KIND_MULTIBROT {
} else if kind == KIND_MULTIBROT {
return multibrot_delta(z, e, clamp(u.power, 2u, 8u));
} else if KIND == KIND_CELTIC {
} else if kind == KIND_CELTIC {
// z^2 delta split: sq.x = delta of Re(z^2), sq.y = delta of Im(z^2).
// Celtic abs the real output, so |Re(z^2)| delta = diffabs(Re(Z^2), sq.x).
let sq = 2.0 * cmul(z, e) + cmul(e, e);
return vec2<f32>(diffabs(z.x * z.x - z.y * z.y, sq.x), sq.y);
} else if KIND == KIND_BUFFALO {
} else if kind == KIND_BUFFALO {
// Abs both outputs: real |Re(z^2)|, imag -|Im(z^2)| (Im(Z^2) = 2 X Y).
let sq = 2.0 * cmul(z, e) + cmul(e, e);
return vec2<f32>(diffabs(z.x * z.x - z.y * z.y, sq.x),
-diffabs(2.0 * z.x * z.y, sq.y));
} else if KIND == KIND_PERPENDICULAR {
} else if kind == KIND_PERPENDICULAR {
// real x^2 - y^2 (ordinary square delta), imag -2 x |y|.
// d(-2 x |y|) = -2[ X·(|Y+ey|-|Y|) + ex·|Y+ey| ]; diffabs gives |Y+ey|-|Y|.
let sq = 2.0 * cmul(z, e) + cmul(e, e);
let da = diffabs(z.y, e.y); // |Y + ey| - |Y|
let abs_yf = abs(z.y) + da; // |Y + ey|
return vec2<f32>(sq.x, -2.0 * (z.x * da + e.x * abs_yf));
} else if KIND == KIND_LAMBDA {
} else if kind == KIND_LAMBDA {
// Lambda map: z^{n+1} = λ·z·(1-z). Delta: e = λ·e·(1-2z-e).
let one_minus_2z_minus_e = vec2<f32>(1.0 - 2.0 * z.x - e.x, -2.0 * z.y - e.y);
return cmul(u.lambda_l, cmul(e, one_minus_2z_minus_e));
} else if KIND == KIND_COMPLEX_MULTIBROT {
} else if kind == KIND_COMPLEX_MULTIBROT {
return complex_multibrot_delta(z, e, u.complex_power);
}
return 2.0 * cmul(z, e) + cmul(e, e); // Mandelbrot (and Phoenix square part)
@@ -188,24 +192,59 @@ fn advance_delta(z: vec2<f32>, e: vec2<f32>) -> vec2<f32> {
// holomorphic kinds (z^2 -> 2Z, z^p -> p Z^{p-1}); for the non-holomorphic
// Burning Ship / Tricorn we use |f'| ~ |2Z|, which keeps the DE magnitude close
// enough to de-speckle filaments.
fn fprime(z: vec2<f32>) -> vec2<f32> {
if KIND == KIND_MULTIBROT {
fn fprime_kind(kind: u32, z: vec2<f32>) -> vec2<f32> {
if kind == KIND_MULTIBROT {
let p = clamp(u.power, 2u, 8u);
var zk = z; // Z^1
for (var k: u32 = 2u; k < p; k = k + 1u) {
zk = cmul(zk, z); // -> Z^{p-1}
}
return f32(p) * zk;
} else if KIND == KIND_LAMBDA {
} else if kind == KIND_LAMBDA {
// Lambda: f'(z) = λ·(1-2z).
return cmul(u.lambda_l, vec2<f32>(1.0 - 2.0 * z.x, -2.0 * z.y));
} else if KIND == KIND_COMPLEX_MULTIBROT {
} else if kind == KIND_COMPLEX_MULTIBROT {
// f'(z) = p * z^(p-1).
return cmul(u.complex_power, cpow(z, u.complex_power - vec2<f32>(1.0, 0.0)));
}
return 2.0 * z;
}
// Delta step of the current map. While switching kinds (MORPH), the map is
// blended per iteration, (1 - w)*f_kind + w*f_from; that's linear in the two
// outputs, so its delta is the same blend of both kinds' deltas (the CPU
// reference in reference.rs uses the same blend, so rebasing stays exact).
fn advance_delta(z: vec2<f32>, e: vec2<f32>) -> vec2<f32> {
let d = advance_delta_kind(KIND, z, e);
if MORPH {
return mix(d, advance_delta_kind(u.morph_from, z, e), u.morph_w);
}
return d;
}
// Derivative of the current (possibly morphing) map; blended like
// `advance_delta`.
fn fprime(z: vec2<f32>) -> vec2<f32> {
let d = fprime_kind(KIND, z);
if MORPH {
return mix(d, fprime_kind(u.morph_from, z), u.morph_w);
}
return d;
}
// Weight of the Phoenix kind's p*z_{n-1} term in the current map: 1 for plain
// Phoenix, its morph share while switching to/from Phoenix, else 0.
fn phoenix_weight() -> f32 {
var w = 0.0;
if KIND == KIND_PHOENIX {
w = select(1.0, 1.0 - u.morph_w, MORPH);
}
if MORPH && u.morph_from == KIND_PHOENIX {
w = w + u.morph_w;
}
return w;
}
// Periodicity (interior) detection, Brent-style: the full orbit value is
// saved at iterations PERIOD_FIRST_CHECK, 2x that, 4x ..., and every later
// iterate is compared against the last saved one. Returning within
@@ -229,7 +268,8 @@ fn fprime(z: vec2<f32>) -> vec2<f32> {
// Every kind here except Phoenix is
// (piecewise) conformal, so |f'| from `fprime` is the exact local scale
// factor, including the abs-folding kinds, whose folds are isometries.
// Phoenix's two-term map would need a 2x2 Jacobian, so it's excluded. So is
// Phoenix's two-term map would need a 2x2 Jacobian, so it's excluded (as
// is a kind-switch morph, for the same reason). So is
// Complex Multibrot without DE, where `fprime` would add a second `cpow`
// (log/atan2/exp) per step for a check that rarely fires on its views.
const PERIOD_FIRST_CHECK: u32 = 16u;
@@ -240,7 +280,8 @@ const PERIOD_CONFIRMATIONS: u32 = 2u;
// Whether `iterate_sample` runs periodicity detection for this kind (folds to
// a constant per pipeline).
fn periodic_enabled() -> bool {
if KIND == KIND_PHOENIX {
// A blend of two maps isn't conformal, so |f'| isn't its scale factor.
if MORPH || KIND == KIND_PHOENIX {
return false;
}
if KIND == KIND_COMPLEX_MULTIBROT && !DE {
@@ -276,7 +317,7 @@ fn iterate_sample(offset: vec2<f32>, px: f32) -> Sample {
// orbit itself, since X_1 = X_0^2 + C_ref = C_ref. That's only f32-accurate,
// so skip the test once a pixel is smaller than that error (deep zoom),
// where it could misclassify pixels right at the boundary.
if KIND == KIND_MANDELBROT && !IS_JULIA && ref_len > 1u && px > 1e-6 {
if KIND == KIND_MANDELBROT && !MORPH && !IS_JULIA && ref_len > 1u && px > 1e-6 {
let c = ref_orbit[1] + offset;
let xq = c.x - 0.25;
let q = xq * xq + c.y * c.y;
@@ -308,6 +349,7 @@ fn iterate_sample(offset: vec2<f32>, px: f32) -> Sample {
// y_{n-1}, and its scaled derivative for DE). Both start at 0 (y_{-1} = 0).
var e_prev = vec2<f32>(0.0, 0.0);
var dzs_prev = vec2<f32>(0.0, 0.0);
let phoenix_w = phoenix_weight();
var m: u32 = 0u; // reference index; invariant: y_n = xm + e, xm = X[m]
var n: u32 = 0u; // total iteration count
@@ -352,8 +394,8 @@ fn iterate_sample(offset: vec2<f32>, px: f32) -> Sample {
if !IS_JULIA {
dzs_new.x = dzs_new.x + px;
}
if KIND == KIND_PHOENIX {
dzs_new = dzs_new + cmul(u.phoenix_p, dzs_prev);
if phoenix_w > 0.0 {
dzs_new = dzs_new + phoenix_w * cmul(u.phoenix_p, dzs_prev);
dzs_prev = dzs;
}
dzs = dzs_new;
@@ -364,8 +406,8 @@ fn iterate_sample(offset: vec2<f32>, px: f32) -> Sample {
let e_old = e;
let z_old = z;
e = advance_delta(xm, e) + step_add;
if KIND == KIND_PHOENIX {
e = e + cmul(u.phoenix_p, e_prev);
if phoenix_w > 0.0 {
e = e + phoenix_w * cmul(u.phoenix_p, e_prev);
e_prev = e_old;
}
m = m + 1u;
@@ -388,7 +430,7 @@ fn iterate_sample(offset: vec2<f32>, px: f32) -> Sample {
// full value `z` (and `z2`) is unchanged by the re-expression.
// Phoenix: after rebasing the implied previous reference is Y[-1]=0,
// so the previous delta becomes the full previous value y_{n-1}.
if KIND == KIND_PHOENIX {
if phoenix_w > 0.0 {
e_prev = z_old;
}
e = z - z0;
+10
View File
@@ -28,6 +28,8 @@ pub struct RefRequest {
pub lambda_l: (f64, f64),
/// Complex exponent for the Complex Multibrot kind (ignored by other kinds).
pub complex_power: (f64, f64),
/// Kind-switch morph: `(from_kind, weight)` blended into every step.
pub morph: Option<(FractalKind, f32)>,
}
pub struct RefResult {
@@ -35,6 +37,10 @@ pub struct RefResult {
pub center_im: Big,
pub half_height: f64,
pub points: Vec<[f32; 2]>,
/// The kind and morph `points` was computed with (echoed from the
/// request).
pub kind: FractalKind,
pub morph: Option<(FractalKind, f32)>,
}
pub struct RefWorker {
@@ -68,6 +74,8 @@ impl RefWorker {
center_im: req.center_im,
half_height: req.half_height,
points,
kind: req.kind,
morph: req.morph,
})
.is_err()
{
@@ -110,6 +118,7 @@ fn compute(req: &RefRequest) -> Vec<[f32; 2]> {
req.phoenix_p,
req.lambda_l,
req.complex_power,
req.morph.map(|(k, w)| (k, w as f64)),
)
} else {
compute_set_reference(
@@ -122,6 +131,7 @@ fn compute(req: &RefRequest) -> Vec<[f32; 2]> {
req.phoenix_p,
req.lambda_l,
req.complex_power,
req.morph.map(|(k, w)| (k, w as f64)),
)
}
}
+18 -11
View File
@@ -77,23 +77,30 @@ fn mandelbrot_shader_is_valid() {
}
/// Every specialization renderer.rs can build (`PipelineKey`: kind × Julia ×
/// DE), for every fragment entry point.
/// DE × morph), for every fragment entry point.
#[test]
fn mandelbrot_shader_specializations_compile() {
let (module, info) = validate("mandelbrot.wgsl", MANDELBROT_SRC);
for kind in 0..kind_count() {
for julia in [0.0, 1.0] {
for de in [0.0, 1.0] {
let constants = [("KIND", kind as f64), ("IS_JULIA", julia), ("DE", de)];
for entry in ["fs_data", "fs_refine", "fs_color"] {
specialize(
"mandelbrot.wgsl",
&module,
&info,
naga::ShaderStage::Fragment,
entry,
&constants,
);
for morph in [0.0, 1.0] {
let constants = [
("KIND", kind as f64),
("IS_JULIA", julia),
("DE", de),
("MORPH", morph),
];
for entry in ["fs_data", "fs_refine", "fs_color"] {
specialize(
"mandelbrot.wgsl",
&module,
&info,
naga::ShaderStage::Fragment,
entry,
&constants,
);
}
}
}
}