From cac558a1fb673bacdf7801b396d21eb79059c413 Mon Sep 17 00:00:00 2001 From: supersurviveur Date: Sun, 20 Sep 2026 12:25:36 +0200 Subject: [PATCH] feat: Add complex multibrot fractal --- src/app.rs | 66 ++++++++++++++++++-- src/cli.rs | 7 +++ src/fractal/buddhabrot.rs | 8 ++- src/fractal/reference.rs | 118 +++++++++++++++++++++++++++++++++++- src/fractal/renderer.rs | 4 ++ src/fractal/share.rs | 18 +++++- src/shaders/buddhabrot.wgsl | 33 +++++++--- src/shaders/colorize.wgsl | 1 + src/shaders/mandelbrot.wgsl | 62 +++++++++++++++++++ src/worker.rs | 4 ++ 10 files changed, 302 insertions(+), 19 deletions(-) diff --git a/src/app.rs b/src/app.rs index 591a4f2..d733aa2 100644 --- a/src/app.rs +++ b/src/app.rs @@ -56,6 +56,7 @@ const KINDS: &[(FractalKind, &str)] = &[ (FractalKind::Buffalo, "Buffalo"), (FractalKind::Phoenix, "Phoenix"), (FractalKind::Lambda, "Lambda"), + (FractalKind::ComplexMultibrot, "Complex Multibrot"), ]; /// UI label for a fractal kind. @@ -69,8 +70,8 @@ fn kind_label(kind: FractalKind) -> &'static str { /// The iteration formula for a kind, in human-readable notation (mirrors the /// doc comments on `FractalKind`'s variants). `power` is only used by -/// Multibrot. -fn kind_formula(kind: FractalKind, power: u32) -> String { +/// Multibrot; `complex_power` only by Complex Multibrot. +fn kind_formula(kind: FractalKind, power: u32, complex_power: (f64, f64)) -> String { match kind { FractalKind::Mandelbrot => "z = z² + c".to_string(), FractalKind::BurningShip => "z = (|Re(z)| + i|Im(z)|)² + c".to_string(), @@ -81,13 +82,16 @@ fn kind_formula(kind: FractalKind, power: u32) -> String { FractalKind::Buffalo => "z = |Re(z²)| − i|Im(z²)| + c".to_string(), FractalKind::Phoenix => "z = z² + c + p·z_prev".to_string(), FractalKind::Lambda => "z = λ·z(1 − z)".to_string(), + FractalKind::ComplexMultibrot => { + format!("z = z^({:.3}{:+.3}i) + c", complex_power.0, complex_power.1) + } } } type JuliaPreset = (&'static str, f64, f64, u32, Option<(f64, f64)>); /// Nice-looking Julia constants offered as presets. -const JULIA_PRESETS: [&[JuliaPreset]; FractalKind::Lambda as usize + 1] = [ +const JULIA_PRESETS: [&[JuliaPreset]; FractalKind::ComplexMultibrot as usize + 1] = [ &[ ("dendrite", -0.8, 0.156, 400, None), ("rabbit", -0.123, 0.745, 400, None), @@ -106,6 +110,7 @@ const JULIA_PRESETS: [&[JuliaPreset]; FractalKind::Lambda as usize + 1] = [ ("archipelago 2", -0.556, 0.253, 500, Some((-0.415, -0.267))), ], &[], + &[], ]; type SetPreset = ( @@ -120,7 +125,7 @@ type SetPreset = ( /// Curated beautiful locations offered as one-click presets. /// Each is `(name, center_re, center_im, half_height, iterations)`; the centers /// are decimals parsed at full precision so deep places stay sharp. -const SET_PRESETS: [&[SetPreset]; FractalKind::Lambda as usize + 1] = [ +const SET_PRESETS: [&[SetPreset]; FractalKind::ComplexMultibrot as usize + 1] = [ &[ ( "Seahorse Valley", @@ -171,6 +176,7 @@ const SET_PRESETS: [&[SetPreset]; FractalKind::Lambda as usize + 1] = [ Some((-0.9, -0.49)), )], &[], + &[], ]; /// Parameters a reference orbit was (or will be) computed for. Used to decide @@ -186,6 +192,7 @@ struct RequestKey { iter: u32, kind: FractalKind, power: u32, + complex_power: (f64, f64), } /// Shared state for an in-progress PNG export. The worker (a background thread @@ -273,6 +280,8 @@ pub struct FractalApp { kind: FractalKind, /// Exponent for the Multibrot kind. power: u32, + /// Complex exponent for the Complex Multibrot kind (`z^power + c`). + complex_power: (f64, f64), julia_c: (f64, f64), /// Distortion constant `p` for the Phoenix kind (`z^2 + c + p·z_{n-1}`). phoenix_p: (f64, f64), @@ -452,6 +461,7 @@ impl FractalApp { mode: FractalMode::Mandelbrot, kind: FractalKind::Mandelbrot, power: 3, + complex_power: (2.0, 0.5), julia_c: (-0.8, 0.156), phoenix_p: (-0.5, 0.0), lambda_l: (-0.5, 0.0), @@ -513,6 +523,15 @@ impl FractalApp { if let Some(p) = cli.power { self.power = p.clamp(2, 8); } + if let Some(cp) = &cli.complex_power { + let p: Vec<&str> = cp.split(',').collect(); + if let (Some(Ok(re)), Some(Ok(im))) = ( + p.first().map(|s| s.trim().parse::()), + p.get(1).map(|s| s.trim().parse::()), + ) { + self.complex_power = (re, im); + } + } self.view = Self::default_view_for(self.mode, self.kind); } if let Some(jc) = cli.julia { @@ -618,6 +637,7 @@ impl FractalApp { julia_c: self.julia_c, phoenix_p: self.phoenix_p, lambda_l: self.lambda_l, + complex_power: self.complex_power, color_scale: self.color_scale, color_offset: self.color_offset, palette: self.palette, @@ -637,6 +657,7 @@ impl FractalApp { self.julia_c = s.julia_c; self.phoenix_p = s.phoenix_p; self.lambda_l = s.lambda_l; + self.complex_power = s.complex_power; self.color_scale = s.color_scale; self.color_offset = s.color_offset; self.palette = (s.palette as usize).min(PALETTE_NAMES.len() - 1) as u32; @@ -688,6 +709,7 @@ impl FractalApp { FractalKind::Buffalo => (-0.5, -0.5, 1.5), FractalKind::Phoenix => (0.0, 0.0, 1.6), FractalKind::Lambda => (0.0, 0.0, 1.6), + FractalKind::ComplexMultibrot => (0.0, 0.0, 1.5), }; ViewState::with_center(big_from_f64(cr, 53), big_from_f64(ci, 53), hh) } @@ -704,6 +726,7 @@ impl FractalApp { iter: self.max_iterations, kind: self.kind, power: self.power, + complex_power: self.complex_power, } } @@ -729,6 +752,7 @@ impl FractalApp { || key.iter != self.max_iterations || key.kind != self.kind || key.power != self.power + || key.complex_power != self.complex_power { return true; } @@ -797,6 +821,7 @@ impl FractalApp { power: key.power, phoenix_p: key.phoenix_p, lambda_l: key.lambda_l, + complex_power: key.complex_power, }); self.pending = true; } @@ -816,6 +841,7 @@ impl FractalApp { key.power, key.phoenix_p, key.lambda_l, + key.complex_power, ) } else { compute_set_reference( @@ -827,6 +853,7 @@ impl FractalApp { key.power, key.phoenix_p, key.lambda_l, + key.complex_power, ) }; self.apply_reference( @@ -881,6 +908,7 @@ impl FractalApp { key.power, key.phoenix_p, key.lambda_l, + key.complex_power, ) } else { compute_set_reference( @@ -892,6 +920,7 @@ impl FractalApp { key.power, key.phoenix_p, key.lambda_l, + key.complex_power, ) }; self.apply_reference( @@ -921,6 +950,7 @@ impl FractalApp { 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], + complex_power: [self.complex_power.0 as f32, self.complex_power.1 as f32], de_coloring: (self.de_coloring | self.shadow) as u32, shadow: self.shadow as u32, _pad: [0; _], @@ -941,6 +971,7 @@ impl FractalApp { aspect: aspect as f32, 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], + complex_power: [self.complex_power.0 as f32, self.complex_power.1 as f32], bailout_sq: BAILOUT_SQ, kind: self.kind as u32, power: self.power, @@ -954,7 +985,7 @@ impl FractalApp { height: 0, // set by the callback from size_px total_samples: 0.0, // tracked by the renderer across frames palette: self.buddha_palette, - _pad: [0; 3], + _pad0: 0, } } @@ -1246,6 +1277,12 @@ impl FractalApp { self.lambda_l.0, self.lambda_l.1 )); } + if self.kind == FractalKind::ComplexMultibrot { + ui.label(format!( + "power = {:.6} {:+.6}i", + self.complex_power.0, self.complex_power.1 + )); + } ui.separator(); ui.label(self.kind.description()); @@ -1270,7 +1307,8 @@ impl FractalApp { ui.label( "A deep-zoom fractal explorer. It renders the Mandelbrot set \ and several related fractals (Burning Ship, Tricorn, \ - Multibrot, Celtic, Perpendicular, Buffalo, Phoenix, Lambda).", + Multibrot, Complex Multibrot, Celtic, Perpendicular, Buffalo, \ + Phoenix, Lambda).", ); ui.add_space(4.0); ui.label( @@ -1511,6 +1549,22 @@ impl FractalApp { ui.label("i"); }); } + if self.kind == FractalKind::ComplexMultibrot { + ui.horizontal(|ui| { + ui.label("power ="); + ui.add( + egui::DragValue::new(&mut self.complex_power.0) + .speed(0.01) + .range(-8.0..=8.0), + ); + ui.add( + egui::DragValue::new(&mut self.complex_power.1) + .speed(0.01) + .range(-8.0..=8.0), + ); + ui.label("i"); + }); + } if self.kind != prev_kind { self.view = Self::default_view_for(self.mode, self.kind); } diff --git a/src/cli.rs b/src/cli.rs index 3030481..e2563ef 100644 --- a/src/cli.rs +++ b/src/cli.rs @@ -17,6 +17,10 @@ pub struct Cli { #[arg(long)] pub power: Option, + /// Complex exponent for the Complex Multibrot kind (z -> z^power + c). + #[arg(long, value_name = "RE,IM")] + pub complex_power: Option, + /// Start in Julia mode with this seed constant. #[arg(long, value_name = "RE,IM")] pub julia: Option, @@ -79,6 +83,8 @@ pub enum KindArg { Buffalo, Phoenix, Lambda, + #[value(alias = "cmulti")] + ComplexMultibrot, } impl From for FractalKind { @@ -93,6 +99,7 @@ impl From for FractalKind { KindArg::Buffalo => FractalKind::Buffalo, KindArg::Phoenix => FractalKind::Phoenix, KindArg::Lambda => FractalKind::Lambda, + KindArg::ComplexMultibrot => FractalKind::ComplexMultibrot, } } } diff --git a/src/fractal/buddhabrot.rs b/src/fractal/buddhabrot.rs index 0ef3589..ca3f5e9 100644 --- a/src/fractal/buddhabrot.rs +++ b/src/fractal/buddhabrot.rs @@ -46,7 +46,11 @@ pub struct BuddhabrotUniforms { /// (yellow core, blue halo), 2 = grayscale. Display-only, like `exposure` /// — excluded from `ContentKey` so changing it doesn't reset accumulation. pub palette: u32, - pub _pad: [u32; 3], + /// Padding so `complex_power` (a vec2, 8-byte aligned in the shader) + /// starts on an 8-byte boundary. + pub _pad0: u32, + /// Complex exponent for the Complex Multibrot kind; ignored by other kinds. + pub complex_power: [f32; 2], } /// The subset of `BuddhabrotUniforms` that determines the *content* of the @@ -62,6 +66,7 @@ struct ContentKey { bailout_sq: f32, kind: u32, power: u32, + complex_power: [f32; 2], r_cap: u32, g_cap: u32, b_cap: u32, @@ -78,6 +83,7 @@ impl From<&BuddhabrotUniforms> for ContentKey { bailout_sq: u.bailout_sq, kind: u.kind, power: u.power, + complex_power: u.complex_power, r_cap: u.r_cap, g_cap: u.g_cap, b_cap: u.b_cap, diff --git a/src/fractal/reference.rs b/src/fractal/reference.rs index d0fbd57..8a1189a 100644 --- a/src/fractal/reference.rs +++ b/src/fractal/reference.rs @@ -35,6 +35,9 @@ pub enum FractalKind { Phoenix = 7, /// `z -> lambda·z(1 - z)` (logistic map). Lambda = 8, + /// `z -> z^power + c`, where `power` is a complex constant (the + /// `complex_power` argument), via the principal branch `z^p = exp(p·ln z)`. + ComplexMultibrot = 9, } impl FractalKind { @@ -53,6 +56,7 @@ impl FractalKind { FractalKind::Buffalo => "", FractalKind::Phoenix => "", FractalKind::Lambda => "", + FractalKind::ComplexMultibrot => "Like Multibrot, but the exponent itself is a complex number instead of a plain integer, via z^p = exp(p·ln z).", } } } @@ -77,6 +81,7 @@ pub fn compute_reference( power: u32, phoenix_p: (f64, f64), lambda_l: (f64, f64), + complex_power: (f64, f64), ) -> Vec<[f32; 2]> { let cr = c_re.clone().with_precision(precision).value(); let ci = c_im.clone().with_precision(precision).value(); @@ -92,6 +97,9 @@ pub fn compute_reference( // 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); @@ -162,6 +170,10 @@ pub fn compute_reference( 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) + } }; // Shift the previous iterate (only the Phoenix arm reads it). @@ -198,6 +210,32 @@ fn complex_pow(zr: &Big, zi: &Big, power: u32, precision: usize) -> (Big, Big) { (rr, ri) } +/// `true` if `x` is (numerically) zero. The f64 check is exact for a true +/// zero; only matters here to special-case `ln(0)`. +fn is_big_zero(x: &Big) -> bool { + x.to_f64().value() == 0.0 +} + +/// `(zr + i zi)^(pr + i pi)` for a complex exponent, via the principal branch +/// `z^p = exp(p·ln z)` where `ln z = ln|z| + i·arg(z)`. Used by +/// `ComplexMultibrot`; must be kept in sync with the shader's `cpow`. +/// `z = 0` is special-cased to `0` (the formula's `ln(0)` would otherwise +/// panic; this is the correct limit for the `Re(p) > 0` region the UI +/// exposes). +fn complex_pow_complex(zr: &Big, zi: &Big, pr: &Big, pi: &Big, precision: usize) -> (Big, Big) { + if is_big_zero(zr) && is_big_zero(zi) { + return (big_zero(precision), big_zero(precision)); + } + let r2 = &zr.sqr() + &zi.sqr(); + let ln_r = r2.ln() >> 1; // 0.5 * ln(r2) = ln(sqrt(r2)); exact halving. + let theta = zi.atan2(zr); + let exp_re = (pr * &ln_r - pi * &theta).with_precision(precision).value(); + let exp_im = (pr * &theta + pi * &ln_r).with_precision(precision).value(); + let mag = exp_re.exp(); + let (sin_a, cos_a) = exp_im.sin_cos(); + (&mag * &cos_a, &mag * &sin_a) +} + /// Convenience: parameter-plane ("Mandelbrot-set") reference (`z0 = 0`, /// `c = center`) for any `kind`. #[allow(clippy::too_many_arguments)] @@ -210,10 +248,21 @@ pub fn compute_set_reference( power: u32, phoenix_p: (f64, f64), lambda_l: (f64, f64), + complex_power: (f64, f64), ) -> Vec<[f32; 2]> { let zero = big_zero(precision); compute_reference( - &zero, &zero, center_re, center_im, max_iter, precision, kind, power, phoenix_p, lambda_l, + &zero, + &zero, + center_re, + center_im, + max_iter, + precision, + kind, + power, + phoenix_p, + lambda_l, + complex_power, ) } @@ -236,6 +285,7 @@ mod tests { 2, (0.0, 0.0), (0.0, 0.0), + (0.0, 0.0), ); // Independent naive f64 orbit. @@ -275,6 +325,7 @@ mod tests { 2, (0.0, 0.0), (0.0, 0.0), + (0.0, 0.0), ); assert_eq!(points.len(), 501, "interior orbit should not escape"); } @@ -293,6 +344,7 @@ mod tests { 2, (0.0, 0.0), (0.0, 0.0), + (0.0, 0.0), ); let (c_re, c_im) = (-1.75_f64, -0.03_f64); @@ -322,6 +374,7 @@ mod tests { 3, (0.0, 0.0), (0.0, 0.0), + (0.0, 0.0), ); let (c_re, c_im) = (0.3_f64, 0.2_f64); @@ -357,6 +410,7 @@ mod tests { 2, (0.0, 0.0), (0.0, 0.0), + (0.0, 0.0), ); let (mut zr, mut zi) = (0.15_f64, -0.1_f64); @@ -386,6 +440,7 @@ mod tests { 2, (0.0, 0.0), (0.0, 0.0), + (0.0, 0.0), ); let (c_re, c_im) = (-0.6_f64, 0.4_f64); @@ -416,6 +471,7 @@ mod tests { 2, (0.0, 0.0), (0.0, 0.0), + (0.0, 0.0), ); let (c_re, c_im) = (-0.7_f64, -0.2_f64); @@ -446,6 +502,7 @@ mod tests { 2, (0.0, 0.0), (0.0, 0.0), + (0.0, 0.0), ); let (c_re, c_im) = (-1.2_f64, -0.35_f64); @@ -468,8 +525,17 @@ mod tests { let cr = Big::try_from(0.5667_f64).unwrap(); let ci = Big::try_from(0.0_f64).unwrap(); let p = (-0.5_f64, 0.0_f64); - let points = - compute_set_reference(&cr, &ci, 60, 200, FractalKind::Phoenix, 2, p, (0.0, 0.0)); + let points = compute_set_reference( + &cr, + &ci, + 60, + 200, + FractalKind::Phoenix, + 2, + p, + (0.0, 0.0), + (0.0, 0.0), + ); let (c_re, c_im) = (0.5667_f64, 0.0_f64); let (mut zr, mut zi) = (0.0_f64, 0.0_f64); @@ -489,4 +555,50 @@ mod tests { zi = nzi; } } + + /// Complex Multibrot (power 2.5 + 0.3i) reference matches a naive f64 + /// iteration of `z^p = exp(p·ln z)`. + #[test] + fn complex_multibrot_reference_matches_naive_f64() { + let cr = Big::try_from(0.1_f64).unwrap(); + let ci = Big::try_from(-0.2_f64).unwrap(); + let power = (2.5_f64, 0.3_f64); + let points = compute_set_reference( + &cr, + &ci, + 60, + 200, + FractalKind::ComplexMultibrot, + 2, + (0.0, 0.0), + (0.0, 0.0), + power, + ); + + // Naive f64 complex power via z^p = exp(p * ln z), ln z = ln|z| + i*arg(z). + fn naive_cpow(zr: f64, zi: f64, pr: f64, pi: f64) -> (f64, f64) { + if zr == 0.0 && zi == 0.0 { + return (0.0, 0.0); + } + let ln_r = 0.5 * (zr * zr + zi * zi).ln(); + let theta = zi.atan2(zr); + let exp_re = pr * ln_r - pi * theta; + let exp_im = pr * theta + pi * ln_r; + let mag = exp_re.exp(); + (mag * exp_im.cos(), mag * exp_im.sin()) + } + + let (c_re, c_im) = (0.1_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}"); + let (pr, pi) = naive_cpow(zr, zi, power.0, power.1); + let nzr = pr + c_re; + let nzi = pi + c_im; + zr = nzr; + zi = nzi; + } + } } diff --git a/src/fractal/renderer.rs b/src/fractal/renderer.rs index c8213ed..30b1232 100644 --- a/src/fractal/renderer.rs +++ b/src/fractal/renderer.rs @@ -37,6 +37,7 @@ fn geom_differs(a: &Uniforms, b: &Uniforms) -> bool { || a.aa_level != b.aa_level || a.kind != b.kind || a.power != b.power + || a.complex_power != b.complex_power || a.dc_offset != b.dc_offset || a.phoenix_p != b.phoenix_p || a.de_coloring != b.de_coloring @@ -87,6 +88,9 @@ pub struct Uniforms { /// Distortion constant `l` for the Lambda map (`l·z(1 - z)`); /// ignored by other kinds. pub lambda_l: [f32; 2], + /// Complex exponent for the Complex Multibrot kind (`z^power + c`); + /// ignored by other kinds. + pub complex_power: [f32; 2], /// 0 = escape-time coloring, 1 = distance-estimation shading. pub de_coloring: u32, // 0 = classic colors, 1 = shadows diff --git a/src/fractal/share.rs b/src/fractal/share.rs index 4550b29..49a17b1 100644 --- a/src/fractal/share.rs +++ b/src/fractal/share.rs @@ -25,6 +25,8 @@ pub struct ShareState { pub phoenix_p: (f64, f64), /// Distortion constant for the Lambda kind (ignored by others). pub lambda_l: (f64, f64), + /// Complex exponent for the Complex Multibrot kind (ignored by others). + pub complex_power: (f64, f64), pub color_scale: f32, pub color_offset: f32, /// Palette index (`palette_id` in the shader). @@ -49,6 +51,7 @@ impl ShareState { FractalKind::Buffalo => "buffalo", FractalKind::Phoenix => "phoenix", FractalKind::Lambda => "lambda", + FractalKind::ComplexMultibrot => "cmulti", } )); s.push_str(&format!("&pw={}", self.power)); @@ -60,8 +63,12 @@ impl ShareState { s.push_str(&format!("&px={}&py={}", self.phoenix_p.0, self.phoenix_p.1)); s.push_str(&format!("&lx={}&ly={}", self.lambda_l.0, self.lambda_l.1)); s.push_str(&format!( - "&cs={}&co={}&pal={}", - self.color_scale, self.color_offset, self.palette + "&cpr={}&cpi={}", + self.complex_power.0, self.complex_power.1 + )); + s.push_str(&format!( + "&cs={}&co={}&pal={}&spal={}", + self.color_scale, self.color_offset, self.palette, self.shadow_palette )); s } @@ -89,6 +96,7 @@ impl ShareState { "buffalo" => FractalKind::Buffalo, "phoenix" => FractalKind::Phoenix, "lambda" => FractalKind::Lambda, + "cmulti" => FractalKind::ComplexMultibrot, _ => FractalKind::Mandelbrot, }) .unwrap_or(FractalKind::Mandelbrot), @@ -109,6 +117,10 @@ impl ShareState { map.get("lx").and_then(|s| s.parse().ok()).unwrap_or(-0.5), map.get("ly").and_then(|s| s.parse().ok()).unwrap_or(0.0), ), + complex_power: ( + map.get("cpr").and_then(|s| s.parse().ok()).unwrap_or(2.0), + map.get("cpi").and_then(|s| s.parse().ok()).unwrap_or(0.0), + ), color_scale: map.get("cs").and_then(|s| s.parse().ok()).unwrap_or(0.02), color_offset: map.get("co").and_then(|s| s.parse().ok()).unwrap_or(0.0), palette: map.get("pal").and_then(|s| s.parse().ok()).unwrap_or(0), @@ -134,6 +146,7 @@ mod tests { julia_c: (-0.123, 0.745), phoenix_p: (-0.5, 0.1), lambda_l: (-0.5, 0.0), + complex_power: (2.5, 0.3), color_scale: 0.02, color_offset: 0.25, palette: 3, @@ -149,6 +162,7 @@ mod tests { assert_eq!(d.iterations, s.iterations); assert_eq!(d.julia_c, s.julia_c); assert_eq!(d.phoenix_p, s.phoenix_p); + assert_eq!(d.complex_power, s.complex_power); assert_eq!(d.palette, s.palette); assert_eq!(d.shadow_palette, s.shadow_palette); } diff --git a/src/shaders/buddhabrot.wgsl b/src/shaders/buddhabrot.wgsl index 58eb08c..131b8b9 100644 --- a/src/shaders/buddhabrot.wgsl +++ b/src/shaders/buddhabrot.wgsl @@ -47,14 +47,15 @@ struct Uniforms { // Tonemap colour style: 0 = classic (R/G/B = raw caps), 1 = nebula // (yellow core, blue halo), 2 = grayscale. palette: u32, - // Padding to a 16-byte multiple. NOT vec3 — that type aligns to 16 - // bytes in WGSL (unlike Rust's `[u32; 3]`, which aligns to 4), which - // silently added 32 bytes instead of 16 and mismatched the Rust struct's - // size (a wgpu validation error at dispatch time: "size 96 where the - // shader expects 112"). + // Padding so `complex_power` (a vec2, 8-byte aligned) starts on an + // 8-byte boundary. NOT vec3 — that type aligns to 16 bytes in WGSL + // (unlike Rust's `[u32; 3]`, which aligns to 4), which silently added 32 + // bytes instead of 16 and mismatched the Rust struct's size (a wgpu + // validation error at dispatch time: "size 96 where the shader expects + // 112"). _pad0: u32, - _pad1: u32, - _pad2: u32, + // Complex exponent for the Complex Multibrot kind; unused by other kinds. + complex_power: vec2, }; const PALETTE_NEBULA: u32 = 0u; @@ -70,6 +71,7 @@ const KIND_PERPENDICULAR: u32 = 5u; const KIND_BUFFALO: u32 = 6u; const KIND_PHOENIX: u32 = 7u; const KIND_LAMBDA: u32 = 8u; +const KIND_COMPLEX_MULTIBROT: u32 = 9u; @group(0) @binding(0) var u: Uniforms; // Compute pass: read-write atomic histogram (3 planes of width*height, R/G/B). @@ -103,6 +105,21 @@ fn complex_pow(z: vec2, p: u32) -> vec2 { return r; } +// z^p for a complex exponent p, via the principal branch z^p = exp(p * ln z), +// ln z = ln|z| + i*arg(z). z = 0 maps to 0 (the correct limit for the +// Re(p) > 0 region the UI exposes; ln(0) would otherwise be -inf). +fn cpow(z: vec2, p: vec2) -> vec2 { + let r2 = dot(z, z); + if r2 < 1e-30 { + return vec2(0.0, 0.0); + } + let ln_r = 0.5 * log(r2); + let theta = atan2(z.y, z.x); + let mag = exp(p.x * ln_r - p.y * theta); + let ang = p.x * theta + p.y * ln_r; + return mag * vec2(cos(ang), sin(ang)); +} + // One iteration step z_n -> z_{n+1} for the current kind. `zp` is the // previous iterate (z_{n-1}), used only by the Phoenix two-term recurrence. // Must match `FractalKind` in reference.rs (the direct, non-perturbative form @@ -126,6 +143,8 @@ fn advance(z: vec2, zp: vec2, c: vec2) -> vec2 { } else if u.kind == KIND_LAMBDA { // l * z * (1 - z); c is unused (see file doc comment above). return cmul(u.lambda_l, cmul(z, vec2(1.0 - z.x, -z.y))); + } else if u.kind == KIND_COMPLEX_MULTIBROT { + return cpow(z, u.complex_power) + c; } return vec2(z.x * z.x - z.y * z.y, 2.0 * z.x * z.y) + c; // Mandelbrot } diff --git a/src/shaders/colorize.wgsl b/src/shaders/colorize.wgsl index eaaceb2..7f194d1 100644 --- a/src/shaders/colorize.wgsl +++ b/src/shaders/colorize.wgsl @@ -27,6 +27,7 @@ struct Uniforms { dc_offset: vec2, phoenix_p: vec2, lambda_l: vec2, + complex_power: vec2, de_coloring: u32, shadow: u32, }; diff --git a/src/shaders/mandelbrot.wgsl b/src/shaders/mandelbrot.wgsl index 01a2547..a2c2381 100644 --- a/src/shaders/mandelbrot.wgsl +++ b/src/shaders/mandelbrot.wgsl @@ -34,6 +34,9 @@ struct Uniforms { // Distortion constant l for the Lambda map (l*z(1 - z_{n-1})); unused // by other kinds. lambda_l: vec2, + // Complex exponent for the Complex Multibrot kind (z^power + c); unused + // by other kinds. + complex_power: vec2, // 0 = escape-time coloring, 1 = distance-estimation shading. de_coloring: u32, // 0 = classic colors, 1 = shadows @@ -49,6 +52,7 @@ const KIND_PERPENDICULAR: u32 = 5u; const KIND_BUFFALO: u32 = 6u; const KIND_PHOENIX: u32 = 7u; const KIND_LAMBDA: u32 = 8u; +const KIND_COMPLEX_MULTIBROT: u32 = 9u; @group(0) @binding(0) var u: Uniforms; @group(0) @binding(1) var ref_orbit: array>; @@ -84,6 +88,27 @@ fn conj(a: vec2) -> vec2 { return vec2(a.x, -a.y); } +// Complex division a / b. +fn cdiv(a: vec2, b: vec2) -> vec2 { + let d = dot(b, b); + return vec2(a.x * b.x + a.y * b.y, a.y * b.x - a.x * b.y) / d; +} + +// z^p for a complex exponent p, via the principal branch z^p = exp(p * ln z), +// ln z = ln|z| + i*arg(z). z = 0 maps to 0 (the correct limit for the +// Re(p) > 0 region the UI exposes; ln(0) would otherwise be -inf). +fn cpow(z: vec2, p: vec2) -> vec2 { + let r2 = dot(z, z); + if r2 < 1e-30 { + return vec2(0.0, 0.0); + } + let ln_r = 0.5 * log(r2); + let theta = atan2(z.y, z.x); + let mag = exp(p.x * ln_r - p.y * theta); + let ang = p.x * theta + p.y * ln_r; + return mag * vec2(cos(ang), sin(ang)); +} + // |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. @@ -123,6 +148,38 @@ fn multibrot_delta(z: vec2, e: vec2, p: u32) -> vec2 { return acc; } +// Number of terms kept in `complex_multibrot_delta`'s series. Truncation, not +// exactness: unlike `multibrot_delta` (a finite binomial sum for an integer +// power), a complex power has no finite expansion, so this converges rather +// than terminates. Fine as long as perturbation's usual invariant (|e| << |z|, +// kept true by rebasing) holds, since each extra term is O(w^k) smaller. +const COMPLEX_MULTIBROT_TERMS: u32 = 16u; + +// Perturbation delta for z -> z^p with a complex p: (Z+e)^p - Z^p. +// = Z^p * ((1+w)^p - 1), w = e/Z, expanded as a Taylor series in w (never +// forming 1+w, which would round tiny w away in f32 — the same reason +// `multibrot_delta` never forms Z+e directly). Series: (1+w)^p - 1 = +// sum_{k=1}^N C(p,k) w^k, with the complex binomial coefficient built up +// incrementally: C(p,k) = C(p,k-1) * (p-(k-1)) / k. +// +// Z ~ 0 (the reference start, X_0 = 0 for Mandelbrot) makes w singular; there +// (0+e)^p - 0^p = e^p exactly, so that case is handled directly via `cpow`. +fn complex_multibrot_delta(z: vec2, e: vec2, p: vec2) -> vec2 { + if dot(z, z) < 1e-20 { + return cpow(e, p); + } + let w = cdiv(e, z); + var wk = vec2(1.0, 0.0); // w^0 + var coef = vec2(1.0, 0.0); // C(p,0) + var acc = vec2(0.0, 0.0); + for (var k: u32 = 1u; k <= COMPLEX_MULTIBROT_TERMS; k = k + 1u) { + coef = cdiv(cmul(coef, p - vec2(f32(k - 1u), 0.0)), vec2(f32(k), 0.0)); + wk = cmul(wk, w); + acc = acc + cmul(coef, wk); + } + return cmul(cpow(z, p), 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. @@ -163,6 +220,8 @@ fn advance_delta(z: vec2, e: vec2) -> vec2 { // 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 u.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) } @@ -183,6 +242,9 @@ fn fprime(z: vec2) -> vec2 { } else if u.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 u.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; } diff --git a/src/worker.rs b/src/worker.rs index beb7a4b..be66ae0 100644 --- a/src/worker.rs +++ b/src/worker.rs @@ -26,6 +26,8 @@ pub struct RefRequest { pub phoenix_p: (f64, f64), /// Distortion constant for the Lambda map (ignored by other kinds). pub lambda_l: (f64, f64), + /// Complex exponent for the Complex Multibrot kind (ignored by other kinds). + pub complex_power: (f64, f64), } pub struct RefResult { @@ -107,6 +109,7 @@ fn compute(req: &RefRequest) -> Vec<[f32; 2]> { req.power, req.phoenix_p, req.lambda_l, + req.complex_power, ) } else { compute_set_reference( @@ -118,6 +121,7 @@ fn compute(req: &RefRequest) -> Vec<[f32; 2]> { req.power, req.phoenix_p, req.lambda_l, + req.complex_power, ) } }