From c36c8b9829e9f69061eb610bc452f9273430a458 Mon Sep 17 00:00:00 2001 From: supersurviveur Date: Tue, 15 Sep 2026 11:29:44 +0200 Subject: [PATCH] feat: add other fractals --- src/app.rs | 106 +++++++++++++++++++--- src/fractal/mod.rs | 2 +- src/fractal/reference.rs | 169 +++++++++++++++++++++++++++++++----- src/fractal/renderer.rs | 6 ++ src/fractal/share.rs | 32 +++++-- src/shaders/mandelbrot.wgsl | 80 ++++++++++++++++- src/worker.rs | 15 +++- 7 files changed, 363 insertions(+), 47 deletions(-) diff --git a/src/app.rs b/src/app.rs index 11a46cd..7e6ffad 100644 --- a/src/app.rs +++ b/src/app.rs @@ -5,10 +5,11 @@ use eframe::egui_wgpu; use eframe::egui_wgpu::wgpu; use crate::fractal::{ - ExportRender, FractalCallback, FractalRenderer, MAX_REF_POINTS, ShareState, Uniforms, + ExportRender, FractalCallback, FractalKind, FractalRenderer, MAX_REF_POINTS, ShareState, + Uniforms, }; #[cfg(target_arch = "wasm32")] -use crate::fractal::{compute_mandelbrot_reference, compute_reference}; +use crate::fractal::{compute_reference, compute_set_reference}; use crate::view::{ Big, DEFAULT_HALF_HEIGHT, ViewState, big_from_decimal_str, big_from_f64, big_to_decimal_str, precision_for, @@ -26,6 +27,23 @@ pub enum FractalMode { Julia, } +/// Selectable fractal formulas, with UI labels. +const KINDS: &[(FractalKind, &str)] = &[ + (FractalKind::Mandelbrot, "Mandelbrot"), + (FractalKind::BurningShip, "Burning Ship"), + (FractalKind::Tricorn, "Tricorn"), + (FractalKind::Multibrot, "Multibrot"), +]; + +/// UI label for a fractal kind. +fn kind_label(kind: FractalKind) -> &'static str { + KINDS + .iter() + .find(|(k, _)| *k == kind) + .map(|(_, name)| *name) + .unwrap_or("Mandelbrot") +} + /// Nice-looking Julia constants offered as presets. const JULIA_PRESETS: &[(&str, f64, f64)] = &[ ("dendrite", -0.8, 0.156), @@ -73,6 +91,8 @@ struct RequestKey { julia: bool, julia_c: (f64, f64), iter: u32, + kind: FractalKind, + power: u32, } /// Shared state for an in-progress PNG export. The worker (a background thread @@ -89,6 +109,10 @@ struct ExportShared { pub struct FractalApp { view: ViewState, mode: FractalMode, + /// Iteration formula. + kind: FractalKind, + /// Exponent for the Multibrot kind. + power: u32, julia_c: (f64, f64), max_iterations: u32, color_scale: f32, @@ -182,6 +206,8 @@ impl FractalApp { let mut app = Self { view, mode: FractalMode::Mandelbrot, + kind: FractalKind::Mandelbrot, + power: 3, julia_c: (-0.8, 0.156), max_iterations: 512, color_scale: 0.15, @@ -219,6 +245,20 @@ impl FractalApp { // Debug/testing hooks. #[cfg(not(target_arch = "wasm32"))] { + if let Ok(k) = std::env::var("MANDEL_KIND") { + app.kind = match k.trim().to_ascii_lowercase().as_str() { + "burningship" | "burning_ship" | "ship" => FractalKind::BurningShip, + "tricorn" | "mandelbar" => FractalKind::Tricorn, + "multibrot" | "multi" => FractalKind::Multibrot, + _ => FractalKind::Mandelbrot, + }; + if let Ok(p) = std::env::var("MANDEL_POWER") + && let Ok(p) = p.trim().parse::() + { + app.power = p.clamp(2, 8); + } + app.view = Self::default_view_for(app.mode, app.kind); + } if let Ok(jc) = std::env::var("MANDEL_JULIA") { let p: Vec<&str> = jc.split(',').collect(); if let (Some(Ok(re)), Some(Ok(im))) = ( @@ -227,7 +267,7 @@ impl FractalApp { ) { app.mode = FractalMode::Julia; app.julia_c = (re, im); - app.view = Self::default_view_for(FractalMode::Julia); + app.view = Self::default_view_for(FractalMode::Julia, app.kind); } } if let Ok(frag) = std::env::var("MANDEL_SHARE") @@ -294,6 +334,7 @@ impl FractalApp { let sig_digits = sig_digits_for(self.view.precision_bits()); ShareState { julia: matches!(self.mode, FractalMode::Julia), + kind: self.kind, center_re: big_to_decimal_str(&self.view.center_re, sig_digits), center_im: big_to_decimal_str(&self.view.center_im, sig_digits), half_height: self.view.half_height, @@ -311,6 +352,7 @@ impl FractalApp { } else { FractalMode::Mandelbrot }; + self.kind = s.kind; self.julia_c = s.julia_c; self.color_scale = s.color_scale; self.color_offset = s.color_offset; @@ -340,14 +382,20 @@ impl FractalApp { format!("#{fragment}") } - /// Default view for a given fractal mode. - fn default_view_for(mode: FractalMode) -> ViewState { - match mode { - FractalMode::Mandelbrot => ViewState::default(), - FractalMode::Julia => { - ViewState::with_center(big_from_f64(0.0, 53), big_from_f64(0.0, 53), 1.5) - } + /// Default view for a given set type and fractal kind. The Julia (dynamical) + /// plane is centered on the origin for every kind; the parameter plane frames + /// each kind's interesting region. + fn default_view_for(mode: FractalMode, kind: FractalKind) -> ViewState { + if mode == FractalMode::Julia { + return ViewState::with_center(big_from_f64(0.0, 53), big_from_f64(0.0, 53), 1.5); } + let (cr, ci, hh) = match kind { + FractalKind::Mandelbrot => (-0.5, 0.0, 1.25), + FractalKind::BurningShip => (-0.5, -0.5, 1.3), + FractalKind::Tricorn => (-0.25, 0.0, 1.6), + FractalKind::Multibrot => (0.0, 0.0, 1.5), + }; + ViewState::with_center(big_from_f64(cr, 53), big_from_f64(ci, 53), hh) } fn current_key(&self) -> RequestKey { @@ -358,6 +406,8 @@ impl FractalApp { julia: matches!(self.mode, FractalMode::Julia), julia_c: self.julia_c, iter: self.max_iterations, + kind: self.kind, + power: self.power, } } @@ -378,6 +428,8 @@ impl FractalApp { if key.julia != matches!(self.mode, FractalMode::Julia) || key.julia_c != self.julia_c || key.iter != self.max_iterations + || key.kind != self.kind + || key.power != self.power { return true; } @@ -422,6 +474,8 @@ impl FractalApp { julia_c: key.julia_c, max_iter, precision, + kind: key.kind, + power: key.power, }); self.pending = true; } @@ -437,13 +491,17 @@ impl FractalApp { &ji, max_iter, precision, + key.kind, + key.power, ) } else { - compute_mandelbrot_reference( + compute_set_reference( &key.center_re, &key.center_im, max_iter, precision, + key.kind, + key.power, ) }; self.apply_reference( @@ -476,7 +534,10 @@ impl FractalApp { is_julia: matches!(self.mode, FractalMode::Julia) as u32, palette_id: self.palette, aa_level: if self.antialias { 2 } else { 1 }, + kind: self.kind.shader_id(), + power: self.power, dc_offset: self.dc_offset(), + _pad: [0, 0], } } @@ -676,8 +737,25 @@ impl FractalApp { ui.heading("Fractal Explorer"); ui.separator(); + // Fractal formula. Switching kinds jumps to a sensible default view, + // since interesting regions differ between fractals. + let prev_kind = self.kind; + egui::ComboBox::from_label("fractal") + .selected_text(kind_label(self.kind)) + .show_ui(ui, |ui| { + for &(kind, name) in KINDS { + ui.selectable_value(&mut self.kind, kind, name); + } + }); + if self.kind == FractalKind::Multibrot { + ui.add(egui::Slider::new(&mut self.power, 2..=8).text("power")); + } + if self.kind != prev_kind { + self.view = Self::default_view_for(self.mode, self.kind); + } + ui.horizontal(|ui| { - ui.radio_value(&mut self.mode, FractalMode::Mandelbrot, "Mandelbrot"); + ui.radio_value(&mut self.mode, FractalMode::Mandelbrot, "Set"); ui.radio_value(&mut self.mode, FractalMode::Julia, "Julia"); }); @@ -705,7 +783,7 @@ impl FractalApp { }); } - if self.mode == FractalMode::Mandelbrot { + if self.mode == FractalMode::Mandelbrot && self.kind == FractalKind::Mandelbrot { ui.label("places:"); ui.horizontal_wrapped(|ui| { for &(name, re, im, half_height, iter) in MANDEL_PLACES { @@ -854,7 +932,7 @@ impl FractalApp { ui.separator(); if ui.button("Reset view").clicked() { - self.view = Self::default_view_for(self.mode); + self.view = Self::default_view_for(self.mode, self.kind); } ui.add_space(8.0); ui.small("Drag to pan · scroll to zoom toward the cursor"); diff --git a/src/fractal/mod.rs b/src/fractal/mod.rs index ccbe6e0..709f6e5 100644 --- a/src/fractal/mod.rs +++ b/src/fractal/mod.rs @@ -5,7 +5,7 @@ pub mod reference; pub mod renderer; pub mod share; -pub use reference::{compute_mandelbrot_reference, compute_reference}; +pub use reference::{FractalKind, compute_reference, compute_set_reference}; pub use renderer::{ ExportRender, FractalCallback, FractalRenderer, MAX_REF_POINTS, Uniforms, encode_png_with_progress, diff --git a/src/fractal/reference.rs b/src/fractal/reference.rs index f4ae348..697af30 100644 --- a/src/fractal/reference.rs +++ b/src/fractal/reference.rs @@ -1,24 +1,51 @@ //! High-precision reference-orbit computation for perturbation rendering. //! -//! We iterate `Z_{n+1} = Z_n^2 + C` at high precision (`dashu-float`), storing -//! each `Z_n` as an `f32` pair. Every pixel is then rendered on the GPU as a -//! small `f32` delta from this orbit — that is what makes deep zoom cheap. See -//! `shaders/mandelbrot.wgsl` for the delta side. +//! We iterate the fractal's formula `Z_{n+1} = f(Z_n, C)` at high precision +//! (`dashu-float`), storing each `Z_n` as an `f32` pair. Every pixel is then +//! rendered on the GPU as a small `f32` delta from this orbit — that is what +//! makes deep zoom cheap. See `shaders/mandelbrot.wgsl` for the delta side; the +//! delta formula there must match the orbit formula here. //! -//! The `(z0, c)` form serves both fractals: -//! * Mandelbrot: `z0 = 0`, `c = view center` (the c-plane point per pixel). -//! * Julia: `z0 = view center`, `c = julia constant` (fixed for all pixels). +//! 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). use crate::view::Big; +/// The iteration formula. Must be kept in sync with `advance_delta` and the +/// `KIND_*` constants in the shader. +#[derive(Clone, Copy, PartialEq, Eq, Debug)] +pub enum FractalKind { + /// `z -> z^2 + c`. + Mandelbrot, + /// `z -> (|Re z| + i|Im z|)^2 + c`. + BurningShip, + /// `z -> conj(z)^2 + c` (the Mandelbar). + Tricorn, + /// `z -> z^power + c` (power >= 2). + Multibrot, +} + +impl FractalKind { + /// Integer id matching the shader's `KIND_*` constants. + pub fn shader_id(self) -> u32 { + match self { + FractalKind::Mandelbrot => 0, + FractalKind::BurningShip => 1, + FractalKind::Tricorn => 2, + FractalKind::Multibrot => 3, + } + } +} + /// Reference orbit escapes once |Z|^2 exceeds this. Kept larger than the pixel /// bailout so pixels escaping alongside the reference can still reach their /// bailout before the stored orbit runs out. const REFERENCE_ESCAPE_SQ: f64 = 1.0e10; /// Compute the reference orbit `Z_0..Z_{len-1}` where `Z_0 = z0` and -/// `Z_{n+1} = Z_n^2 + c`, up to `max_iter` steps at `precision` bits. Each entry -/// is `[re, im]` in f32. +/// `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. pub fn compute_reference( z0_re: &Big, z0_im: &Big, @@ -26,6 +53,8 @@ pub fn compute_reference( c_im: &Big, max_iter: u32, precision: usize, + kind: FractalKind, + power: u32, ) -> Vec<[f32; 2]> { let cr = c_re.clone().with_precision(precision).value(); let ci = c_im.clone().with_precision(precision).value(); @@ -45,15 +74,33 @@ pub fn compute_reference( break; } - // Z = Z^2 + C, with Z^2 = (zr^2 - zi^2) + (2 zr zi) i. - let zr2 = zr.sqr(); - let zi2 = zi.sqr(); - let new_zr = ((&zr2 - &zi2) + &cr).with_precision(precision).value(); - let two_zr_zi = (&zr * &zi) << 1; // exact multiply-by-2 in base 2 - let new_zi = (two_zr_zi + &ci).with_precision(precision).value(); + let (new_zr, new_zi) = match kind { + FractalKind::Mandelbrot => { + // Z^2 = (zr^2 - zi^2) + (2 zr zi) i. + let re = &zr.sqr() - &zi.sqr() + &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) + } + }; - zr = new_zr; - zi = new_zi; + zr = new_zr.with_precision(precision).value(); + zi = new_zi.with_precision(precision).value(); } points @@ -63,15 +110,40 @@ fn big_zero(precision: usize) -> Big { Big::from(0i32).with_precision(precision).value() } -/// Convenience: Mandelbrot reference (`z0 = 0`, `c = center`). -pub fn compute_mandelbrot_reference( +/// Absolute value of a `Big`. The sign check via f64 is exact except for values +/// so tiny that |x| ≈ x either way — negligible against the f32 orbit storage. +fn big_abs(x: Big) -> Big { + if x.to_f64().value() < 0.0 { -x } else { x } +} + +/// `(zr + i zi)^power` by repeated complex multiply at `precision` bits. +fn complex_pow(zr: &Big, zi: &Big, power: u32, precision: usize) -> (Big, Big) { + let mut rr = Big::from(1i32).with_precision(precision).value(); + let mut ri = big_zero(precision); + for _ in 0..power { + // (rr + i ri)(zr + i zi) = (rr zr - ri zi) + (rr zi + ri zr) i. + let nr = (&rr * zr - &ri * zi).with_precision(precision).value(); + let ni = (&rr * zi + &ri * zr).with_precision(precision).value(); + rr = nr; + ri = ni; + } + (rr, ri) +} + +/// Convenience: parameter-plane ("Mandelbrot-set") reference (`z0 = 0`, +/// `c = center`) for any `kind`. +pub fn compute_set_reference( center_re: &Big, center_im: &Big, max_iter: u32, precision: usize, + kind: FractalKind, + power: u32, ) -> Vec<[f32; 2]> { let zero = big_zero(precision); - compute_reference(&zero, &zero, center_re, center_im, max_iter, precision) + compute_reference( + &zero, &zero, center_re, center_im, max_iter, precision, kind, power, + ) } #[cfg(test)] @@ -84,7 +156,7 @@ mod tests { fn reference_matches_naive_f64() { let cr = Big::try_from(-0.75_f64).unwrap(); let ci = Big::try_from(0.1_f64).unwrap(); - let points = compute_mandelbrot_reference(&cr, &ci, 60, 200); + let points = compute_set_reference(&cr, &ci, 60, 200, FractalKind::Mandelbrot, 2); // Independent naive f64 orbit. let (c_re, c_im) = (-0.75_f64, 0.1_f64); @@ -108,10 +180,52 @@ mod tests { fn interior_orbit_runs_full_length() { let cr = Big::try_from(-0.2_f64).unwrap(); let ci = Big::try_from(0.0_f64).unwrap(); - let points = compute_mandelbrot_reference(&cr, &ci, 500, 120); + let points = compute_set_reference(&cr, &ci, 500, 120, FractalKind::Mandelbrot, 2); assert_eq!(points.len(), 501, "interior orbit should not escape"); } + /// Burning Ship reference matches a naive f64 iteration of the same formula. + #[test] + fn burning_ship_reference_matches_naive_f64() { + let cr = Big::try_from(-1.75_f64).unwrap(); + let ci = Big::try_from(-0.03_f64).unwrap(); + let points = compute_set_reference(&cr, &ci, 60, 200, FractalKind::BurningShip, 2); + + let (c_re, c_im) = (-1.75_f64, -0.03_f64); + 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 nzr = zr * zr - zi * zi + c_re; + let nzi = 2.0 * (zr * zi).abs() + c_im; + zr = nzr; + zi = nzi; + } + } + + /// Multibrot (power 3) reference matches a naive f64 cube iteration. + #[test] + fn multibrot3_reference_matches_naive_f64() { + let cr = Big::try_from(0.3_f64).unwrap(); + let ci = Big::try_from(0.2_f64).unwrap(); + let points = compute_set_reference(&cr, &ci, 60, 200, FractalKind::Multibrot, 3); + + let (c_re, c_im) = (0.3_f64, 0.2_f64); + 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}"); + // z^3 = z * z^2. + let (r2, i2) = (zr * zr - zi * zi, 2.0 * zr * zi); + let nzr = zr * r2 - zi * i2 + c_re; + let nzi = zr * i2 + zi * r2 + c_im; + zr = nzr; + zi = nzi; + } + } + /// Julia orbit (fixed c, z0 = center) matches a naive f64 iteration. #[test] fn julia_reference_matches_naive_f64() { @@ -119,7 +233,16 @@ mod tests { let z0_im = Big::try_from(-0.1_f64).unwrap(); let c_re = Big::try_from(-0.8_f64).unwrap(); let c_im = Big::try_from(0.156_f64).unwrap(); - let points = compute_reference(&z0_re, &z0_im, &c_re, &c_im, 60, 200); + let points = compute_reference( + &z0_re, + &z0_im, + &c_re, + &c_im, + 60, + 200, + FractalKind::Mandelbrot, + 2, + ); let (mut zr, mut zi) = (0.15_f64, -0.1_f64); let (cr, ci) = (-0.8_f64, 0.156_f64); diff --git a/src/fractal/renderer.rs b/src/fractal/renderer.rs index b245778..d2a6ba0 100644 --- a/src/fractal/renderer.rs +++ b/src/fractal/renderer.rs @@ -35,10 +35,16 @@ pub struct Uniforms { pub palette_id: u32, /// Supersampling factor per axis: 1 = off, 2 = 2×2 (4 samples). pub aa_level: u32, + /// Iteration formula (`FractalKind::shader_id`). + pub kind: u32, + /// Exponent for the Multibrot kind. + pub power: 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. pub dc_offset: [f32; 2], + /// Padding to a 16-byte multiple (uniform buffer requirement). + pub _pad: [u32; 2], } /// Offscreen texture the fractal is rendered into, plus the bind group used to diff --git a/src/fractal/share.rs b/src/fractal/share.rs index c4f20e8..2a6452b 100644 --- a/src/fractal/share.rs +++ b/src/fractal/share.rs @@ -2,15 +2,18 @@ //! iterations, Julia constant, coloring) as a compact URL fragment so deep-zoom //! locations can be shared or bookmarked. //! -//! Format: `m=m&re=&im=&hh=&it=&cs=&co=` with +//! Format: `m=m&f=&re=&im=&hh=&it=&cs=&co=` with //! `m=j&jr=&ji=` added for Julia. `re`/`im` are full-precision decimal //! strings. use std::collections::HashMap; +use crate::fractal::FractalKind; + #[derive(Clone, Debug)] pub struct ShareState { pub julia: bool, + pub kind: FractalKind, pub center_re: String, pub center_im: String, pub half_height: f64, @@ -24,14 +27,21 @@ impl ShareState { pub fn encode(&self) -> String { let mut s = String::new(); s.push_str(if self.julia { "m=j" } else { "m=m" }); + s.push_str(&format!("&f={}", match self.kind { + FractalKind::Mandelbrot => "mandel", + FractalKind::BurningShip => "burning", + FractalKind::Multibrot => "multi", + FractalKind::Tricorn => "tricorn", + })); s.push_str(&format!( "&re={}&im={}&hh={}&it={}", self.center_re, self.center_im, self.half_height, self.iterations )); - if self.julia { - s.push_str(&format!("&jr={}&ji={}", self.julia_c.0, self.julia_c.1)); - } - s.push_str(&format!("&cs={}&co={}", self.color_scale, self.color_offset)); + s.push_str(&format!("&jr={}&ji={}", self.julia_c.0, self.julia_c.1)); + s.push_str(&format!( + "&cs={}&co={}", + self.color_scale, self.color_offset + )); s } @@ -46,6 +56,16 @@ impl ShareState { Some(ShareState { julia: map.get("m").map(|m| *m == "j").unwrap_or(false), + kind: map + .get("f") + .map(|f| match *f { + "mandel" => FractalKind::Mandelbrot, + "multi" => FractalKind::Multibrot, + "burning" => FractalKind::BurningShip, + "tricorn" => FractalKind::Tricorn, + _ => FractalKind::Mandelbrot + }) + .unwrap_or(FractalKind::Mandelbrot), center_re: (*map.get("re")?).to_string(), center_im: (*map.get("im")?).to_string(), half_height: map.get("hh")?.parse().ok()?, @@ -68,6 +88,7 @@ mod tests { fn round_trip() { let s = ShareState { julia: true, + kind: FractalKind::Tricorn, center_re: "-0.743643887037158704752191506114774".into(), center_im: "0.131825904205311970493132056385139".into(), half_height: 1.5e-20, @@ -78,6 +99,7 @@ mod tests { }; let d = ShareState::decode(&s.encode()).unwrap(); assert_eq!(d.julia, s.julia); + assert_eq!(d.kind, s.kind); assert_eq!(d.center_re, s.center_re); assert_eq!(d.center_im, s.center_im); assert_eq!(d.half_height, s.half_height); diff --git a/src/shaders/mandelbrot.wgsl b/src/shaders/mandelbrot.wgsl index 419367c..67a69ae 100644 --- a/src/shaders/mandelbrot.wgsl +++ b/src/shaders/mandelbrot.wgsl @@ -22,9 +22,18 @@ struct Uniforms { is_julia: u32, palette_id: u32, aa_level: u32, + // Iteration formula: 0 Mandelbrot, 1 Burning Ship, 2 Tricorn, 3 Multibrot. + kind: u32, + // Exponent for the Multibrot kind. + power: u32, dc_offset: vec2, }; +const KIND_MANDELBROT: u32 = 0u; +const KIND_BURNING_SHIP: u32 = 1u; +const KIND_TRICORN: u32 = 2u; +const KIND_MULTIBROT: u32 = 3u; + @group(0) @binding(0) var u: Uniforms; @group(0) @binding(1) var ref_orbit: array>; @@ -54,6 +63,73 @@ fn cmul(a: vec2, b: vec2) -> vec2 { return vec2(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x); } +// Complex conjugate. +fn conj(a: vec2) -> vec2 { + return vec2(a.x, -a.y); +} + +// |c + d| - |c|, evaluated exactly (no catastrophic cancellation even when the +// sum crosses zero). This is what makes the Burning Ship delta correct through +// the sign flips that happen all along the axes, where the ship's detail lives. +fn diffabs(c: f32, d: f32) -> f32 { + let cd = c + d; + if (c >= 0.0) { + return select(-(2.0 * c + d), d, cd >= 0.0); + } + return select(-d, 2.0 * c + d, cd > 0.0); +} + +// Binomial coefficient C(n, k) as f32 (exact for the small powers we use). +fn binom(n: u32, k: u32) -> f32 { + var num = 1.0; + var den = 1.0; + for (var i: u32 = 0u; i < k; i = i + 1u) { + num = num * f32(n - i); + den = den * f32(i + 1u); + } + return num / den; +} + +// Perturbation delta for z -> z^p: sum_{k=1}^{p} C(p,k) Z^{p-k} e^k. Expanded so +// the large z^p term is never formed (that would cancel catastrophically). +fn multibrot_delta(z: vec2, e: vec2, p: u32) -> vec2 { + var zp: array, 9>; // Z^0 .. Z^8 + zp[0] = vec2(1.0, 0.0); + for (var j: u32 = 1u; j <= p; j = j + 1u) { + zp[j] = cmul(zp[j - 1u], z); + } + var acc = vec2(0.0, 0.0); + var ek = vec2(1.0, 0.0); // e^0 + for (var k: u32 = 1u; k <= p; k = k + 1u) { + ek = cmul(ek, e); // e^k + acc = acc + binom(p, k) * cmul(zp[p - k], ek); + } + return acc; +} + +// 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, e: vec2) -> vec2 { + if (u.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 + // old sign(X)sign(Y) shortcut got wrong whenever the delta was large + // enough to flip it (all the time at shallow zoom). + 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 (u.kind == KIND_TRICORN) { + let cz = conj(z); + let ce = conj(e); + return 2.0 * cmul(cz, ce) + cmul(ce, ce); + } else if (u.kind == KIND_MULTIBROT) { + return multibrot_delta(z, e, clamp(u.power, 2u, 8u)); + } + return 2.0 * cmul(z, e) + cmul(e, e); // Mandelbrot +} + // Smooth cyclic palettes (Inigo Quilez cosine palettes), selected by id. fn palette(id: u32, t: f32) -> vec3 { if (id == 4u) { @@ -107,8 +183,8 @@ fn shade(offset: vec2) -> vec3 { break; // interior } - // Advance the delta: e = 2*X_m*e + e^2 (+ dc for Mandelbrot). - e = 2.0 * cmul(xm, e) + cmul(e, e) + step_add; + // Advance the delta by this fractal's formula (+ dc for the set plane). + e = advance_delta(xm, e) + step_add; m = m + 1u; n = n + 1u; diff --git a/src/worker.rs b/src/worker.rs index cd1e63c..5ff4327 100644 --- a/src/worker.rs +++ b/src/worker.rs @@ -9,7 +9,7 @@ use std::sync::mpsc::{Receiver, Sender, TryRecvError, channel}; use std::thread; -use crate::fractal::{compute_mandelbrot_reference, compute_reference}; +use crate::fractal::{FractalKind, compute_reference, compute_set_reference}; use crate::view::{Big, big_from_f64}; pub struct RefRequest { @@ -20,6 +20,8 @@ pub struct RefRequest { pub julia_c: (f64, f64), pub max_iter: u32, pub precision: usize, + pub kind: FractalKind, + pub power: u32, } pub struct RefResult { @@ -97,8 +99,17 @@ fn compute(req: &RefRequest) -> Vec<[f32; 2]> { &ji, req.max_iter, req.precision, + req.kind, + req.power, ) } else { - compute_mandelbrot_reference(&req.center_re, &req.center_im, req.max_iter, req.precision) + compute_set_reference( + &req.center_re, + &req.center_im, + req.max_iter, + req.precision, + req.kind, + req.power, + ) } }