surv/mandelbrot
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

feat: Add complex multibrot fractal

Browse Source
This commit is contained in:
surv
2026-09-20 13:37:33 +02:00
parent c7d687c107
commit 06b52fe954
10 changed files with 311 additions and 19 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/app.rs
+60 -6
View File
@@ -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::<f64>()),
p.get(1).map(|s| s.trim().parse::<f64>()),
) {
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);
}
src/cli.rs
+7
View File
@@ -17,6 +17,10 @@ pub struct Cli {
#[arg(long)]
pub power: Option<u32>,
/// Complex exponent for the Complex Multibrot kind (z -> z^power + c).
#[arg(long, value_name = "RE,IM")]
pub complex_power: Option<String>,
/// Start in Julia mode with this seed constant.
#[arg(long, value_name = "RE,IM")]
pub julia: Option<String>,
@@ -79,6 +83,8 @@ pub enum KindArg {
Buffalo,
Phoenix,
Lambda,
#[value(alias = "cmulti")]
ComplexMultibrot,
}
impl From<KindArg> for FractalKind {
@@ -93,6 +99,7 @@ impl From<KindArg> for FractalKind {
KindArg::Buffalo => FractalKind::Buffalo,
KindArg::Phoenix => FractalKind::Phoenix,
KindArg::Lambda => FractalKind::Lambda,
KindArg::ComplexMultibrot => FractalKind::ComplexMultibrot,
}
}
}
src/fractal/buddhabrot.rs
+7 -1
View File
@@ -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,
src/fractal/reference.rs
+115 -3
View File
@@ -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;
}
}
}
src/fractal/renderer.rs
+4
View File
@@ -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
src/fractal/share.rs
+16 -2
View File
@@ -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);
}
src/shaders/buddhabrot.wgsl
+26 -7
View File
@@ -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<u32> — 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<u32> — 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<f32>,
};
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<uniform> 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<f32>, p: u32) -> vec2<f32> {
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<f32>, p: vec2<f32>) -> vec2<f32> {
let r2 = dot(z, z);
if r2 < 1e-30 {
return vec2<f32>(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<f32>(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<f32>, zp: vec2<f32>, c: vec2<f32>) -> vec2<f32> {
} 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<f32>(1.0 - z.x, -z.y)));
} else if u.kind == KIND_COMPLEX_MULTIBROT {
return cpow(z, u.complex_power) + c;
}
return vec2<f32>(z.x * z.x - z.y * z.y, 2.0 * z.x * z.y) + c; // Mandelbrot
}
src/shaders/colorize.wgsl
+1
View File
@@ -27,6 +27,7 @@ struct Uniforms {
dc_offset: vec2<f32>,
phoenix_p: vec2<f32>,
lambda_l: vec2<f32>,
complex_power: vec2<f32>,
de_coloring: u32,
shadow: u32,
};
src/shaders/mandelbrot.wgsl
+71
View File
@@ -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<f32>,
// Complex exponent for the Complex Multibrot kind (z^power + c); unused
// by other kinds.
complex_power: vec2<f32>,
// 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<uniform> u: Uniforms;
@group(0) @binding(1) var<storage, read> ref_orbit: array<vec2<f32>>;
@@ -84,6 +88,27 @@ fn conj(a: vec2<f32>) -> vec2<f32> {
return vec2<f32>(a.x, -a.y);
}
// Complex division a / b.
fn cdiv(a: vec2<f32>, b: vec2<f32>) -> vec2<f32> {
let d = dot(b, b);
return vec2<f32>(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<f32>, p: vec2<f32>) -> vec2<f32> {
let r2 = dot(z, z);
if r2 < 1e-30 {
return vec2<f32>(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<f32>(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,47 @@ fn multibrot_delta(z: vec2<f32>, e: vec2<f32>, p: u32) -> vec2<f32> {
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.
//
// When |e| << |Z| (the common case: it's the whole reason perturbation
// works), forming Z+e directly would round e away in f32, so instead expand
// = Z^p * ((1+w)^p - 1), w = e/Z, as a Taylor series in w: (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. Unlike `multibrot_delta`
// (a finite binomial sum for an integer power), this only *converges* — and
// only for |w| < 1 — rather than terminating exactly.
//
// Right after a rebase (or near a reference point close to zero, where w is
// singular), e is *not* small relative to Z — that's normal perturbation
// dynamics, not a deep-zoom edge case — and the series above would diverge.
// But forming Z+e directly is numerically safe exactly there (e isn't many
// orders of magnitude smaller than Z), so fall back to a plain subtraction.
fn complex_multibrot_delta(z: vec2<f32>, e: vec2<f32>, p: vec2<f32>) -> vec2<f32> {
// |w|^2 = |e|^2 / |Z|^2; inf or nan (Z ~ 0, or both ~ 0) correctly fails
// the `< 0.25` test below and falls through to the direct branch.
let w2 = dot(e, e) / dot(z, z);
if w2 < 0.25 {
let w = cdiv(e, z);
var wk = vec2<f32>(1.0, 0.0); // w^0
var coef = vec2<f32>(1.0, 0.0); // C(p,0)
var acc = vec2<f32>(0.0, 0.0);
for (var k: u32 = 1u; k <= COMPLEX_MULTIBROT_TERMS; k = k + 1u) {
coef = cdiv(cmul(coef, p - vec2<f32>(f32(k - 1u), 0.0)), vec2<f32>(f32(k), 0.0));
wk = cmul(wk, w);
acc = acc + cmul(coef, wk);
}
return cmul(cpow(z, p), acc);
}
return cpow(z + e, p) - cpow(z, p);
}
// One perturbation step of the current fractal's delta: e -> f(Z+e) - f(Z),
// where `z` is the reference orbit value X_m. `step_add` (dc) is added by the
// caller. Must match `FractalKind` on the CPU side.
@@ -163,6 +229,8 @@ fn advance_delta(z: vec2<f32>, e: vec2<f32>) -> vec2<f32> {
// Lambda map: z^{n+1} = λ·z·(1-z). Delta: e = λ·e·(1-2z-e).
let one_minus_2z_minus_e = vec2<f32>(1.0 - 2.0 * z.x - e.x, -2.0 * z.y - e.y);
return cmul(u.lambda_l, cmul(e, one_minus_2z_minus_e));
} else if 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 +251,9 @@ fn fprime(z: vec2<f32>) -> vec2<f32> {
} else if u.kind == KIND_LAMBDA {
// Lambda: f'(z) = λ·(1-2z).
return cmul(u.lambda_l, vec2<f32>(1.0 - 2.0 * z.x, -2.0 * z.y));
} else if u.kind == KIND_COMPLEX_MULTIBROT {
// f'(z) = p * z^(p-1).
return cmul(u.complex_power, cpow(z, u.complex_power - vec2<f32>(1.0, 0.0)));
}
return 2.0 * z;
}
src/worker.rs
+4
View File
@@ -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,
)
}
}