diff --git a/CLAUDE.md b/CLAUDE.md index 4017846..dea2690 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -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 diff --git a/src/app.rs b/src/app.rs index ce125d3..c39b07a 100644 --- a/src/app.rs +++ b/src/app.rs @@ -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, /// Smoothed frames-per-second, recomputed each ~0.5 s window. Only advances /// while the app is actually repainting (interaction / animation / export); @@ -327,6 +366,10 @@ pub struct FractalApp { ref_center_re: Big, ref_center_im: Big, ref_half_height: f64, + /// Kind-switch morph the current `reference` was computed with. The shader + /// blends with this (not the live morph) so its delta formula always + /// matches the orbit, even while the worker lags a frame behind. + ref_morph: Option<(FractalKind, f32)>, /// Parameters of the most recent reference request (drift baseline / dedupe). last_request: Option, @@ -466,6 +509,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 +518,7 @@ impl FractalApp { ref_center_re, ref_center_im, ref_half_height, + ref_morph: None, last_request: None, #[cfg(not(target_arch = "wasm32"))] worker: crate::worker::RefWorker::spawn(), @@ -667,6 +712,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 +752,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 +825,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 +857,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 +887,16 @@ 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, + morph: Option<(FractalKind, f32)>, + ) { self.reference = Arc::new(points); + self.ref_morph = morph; self.ref_center_re = cre; self.ref_center_im = cim; self.ref_half_height = hh; @@ -865,7 +927,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 +947,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 +968,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 +981,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 +989,7 @@ impl FractalApp { key.center_re.clone(), key.center_im.clone(), key.half_height, + key.morph, ); } @@ -932,7 +998,13 @@ 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.morph, + ); self.pending = false; } } @@ -954,7 +1026,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 +1046,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 +1059,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 +1067,7 @@ impl FractalApp { key.center_re.clone(), key.center_im.clone(), key.half_height, + key.morph, ); self.last_request = Some(key); } @@ -1052,6 +1127,7 @@ impl FractalApp { aa_level: if self.antialias { 2 } else { 1 }, kind: 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 +1135,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 +1144,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 +1662,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 +1793,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 +1820,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 +1831,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 +1990,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 +2187,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 +2426,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 +2504,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 { diff --git a/src/fractal/reference.rs b/src/fractal/reference.rs index 7626761..1b66c44 100644 --- a/src/fractal/reference.rs +++ b/src/fractal/reference.rs @@ -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; + } + } } diff --git a/src/fractal/renderer.rs b/src/fractal/renderer.rs index 7698089..c107d5d 100644 --- a/src/fractal/renderer.rs +++ b/src/fractal/renderer.rs @@ -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. diff --git a/src/shaders/iterate_uniforms.wgsl b/src/shaders/iterate_uniforms.wgsl index 9f809c5..6ae3b53 100644 --- a/src/shaders/iterate_uniforms.wgsl +++ b/src/shaders/iterate_uniforms.wgsl @@ -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, // 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, // Screen dimensions screen_dim: vec2, + // 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. diff --git a/src/shaders/mandelbrot.wgsl b/src/shaders/mandelbrot.wgsl index a1537df..62cf08e 100644 --- a/src/shaders/mandelbrot.wgsl +++ b/src/shaders/mandelbrot.wgsl @@ -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, @@ -137,11 +141,11 @@ fn complex_multibrot_delta(z: vec2, e: vec2, p: vec2) -> vec2 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, e: vec2) -> vec2 { - 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, e: vec2) -> vec2 { + 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, e: vec2) -> vec2 { 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(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(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(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(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(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, e: vec2) -> vec2 { // 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) -> vec2 { - if KIND == KIND_MULTIBROT { +fn fprime_kind(kind: u32, z: vec2) -> vec2 { + 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(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(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, e: vec2) -> vec2 { + 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) -> vec2 { + 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) -> vec2 { // 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, 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, px: f32) -> Sample { // y_{n-1}, and its scaled derivative for DE). Both start at 0 (y_{-1} = 0). var e_prev = vec2(0.0, 0.0); var dzs_prev = vec2(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, 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, 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, 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; diff --git a/src/worker.rs b/src/worker.rs index be66ae0..0bc3ef6 100644 --- a/src/worker.rs +++ b/src/worker.rs @@ -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,8 @@ pub struct RefResult { pub center_im: Big, pub half_height: f64, pub points: Vec<[f32; 2]>, + /// The morph `points` was computed with (echoed from the request). + pub morph: Option<(FractalKind, f32)>, } pub struct RefWorker { @@ -68,6 +72,7 @@ impl RefWorker { center_im: req.center_im, half_height: req.half_height, points, + morph: req.morph, }) .is_err() { @@ -110,6 +115,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 +128,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)), ) } } diff --git a/tests/shader_valid.rs b/tests/shader_valid.rs index ee9b7b7..b73371d 100644 --- a/tests/shader_valid.rs +++ b/tests/shader_valid.rs @@ -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, + ); + } } } }