diff --git a/src/app.rs b/src/app.rs index b40d121..e66114d 100644 --- a/src/app.rs +++ b/src/app.rs @@ -19,14 +19,28 @@ use crate::lights::{Light, gpu_lights}; use crate::view::parse_half_height_spec; use crate::view::parse_re_im_spec; use crate::view::{ - Big, DEFAULT_HALF_HEIGHT, MAX_PRECISION_BITS, Scale, ViewState, big_from_decimal_str, big_from_f64, - big_to_decimal_str, deep_scale_exp, interpolate_view, needs_deep, parse_view_spec, - precision_for, + Big, DEFAULT_HALF_HEIGHT, MAX_PRECISION_BITS, Scale, ViewState, big_from_decimal_str, + big_from_f64, big_to_decimal_str, deep_scale_exp, interpolate_view, needs_deep, + parse_view_spec, precision_for, }; #[cfg(not(target_arch = "wasm32"))] use clap::Parser; const BAILOUT_SQ: f32 = 1.0e6; + +/// Pixel bailout |z|^2 for `kind`. For Multibrot z^p, one step from |z| = R +/// (with the reference within 2R, which rebasing guarantees) must stay a +/// finite f32, |z|^2 included: (2R)^(2p) <= 2^126, i.e. R^2 <= 2^(126/p - 2). +/// Otherwise inf - inf turns into NaN, which never compares above the bailout +/// and paints exterior pixels as interior. Unchanged for p <= 5; still well +/// above the escape radius (<= 2) at the maximum power. +fn bailout_sq(kind: FractalKind, power: u32) -> f32 { + if kind == FractalKind::Multibrot { + BAILOUT_SQ.min((126.0 / power.max(2) as f32 - 2.0).exp2()) + } else { + BAILOUT_SQ + } +} /// Cap on exported image dimension (px), to stay within GPU texture limits. const MAX_EXPORT_DIM: u32 = 8192 * 16; /// While the user is actively panning/zooming, the fractal is rendered into a @@ -1034,7 +1048,7 @@ impl FractalApp { FractalMode::Mandelbrot }; self.kind = s.kind; - self.power = s.power.clamp(2, 200); + self.power = s.power.clamp(2, 20); self.julia_c = s.julia_c; self.phoenix_p = s.phoenix_p; self.lambda_l = s.lambda_l; @@ -1078,10 +1092,18 @@ impl FractalApp { /// 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), Scale::from_f64(1.5)); + return ViewState::with_center( + big_from_f64(0.0, 53), + big_from_f64(0.0, 53), + Scale::from_f64(1.5), + ); } let (cr, ci, hh) = kind.default_set_view(); - ViewState::with_center(big_from_f64(cr, 53), big_from_f64(ci, 53), Scale::from_f64(hh)) + ViewState::with_center( + big_from_f64(cr, 53), + big_from_f64(ci, 53), + Scale::from_f64(hh), + ) } /// The request key for the current state. Its `iter` is the reference @@ -1110,8 +1132,12 @@ impl FractalApp { fn drift_from(&self, key: &RequestKey) -> f64 { let hh = self.view.half_height; let k = -hh.exponent() as isize; - let dre = ((&self.view.center_re - &key.center_re) << k).to_f64().value(); - let dim = ((&self.view.center_im - &key.center_im) << k).to_f64().value(); + let dre = ((&self.view.center_re - &key.center_re) << k) + .to_f64() + .value(); + let dim = ((&self.view.center_im - &key.center_im) << k) + .to_f64() + .value(); (dre * dre + dim * dim).sqrt() / hh.scaled_f64(-hh.exponent()) } @@ -1422,7 +1448,7 @@ impl FractalApp { ref_len: self.reference.len() as u32, color_offset: self.color_offset, color_scale: self.color_scale, - bailout_sq: BAILOUT_SQ, + bailout_sq: bailout_sq(self.ref_kind.unwrap_or(self.kind), self.power), is_julia: matches!(self.mode, FractalMode::Julia) as u32, palette_id: self.palette, shadow_palette_id: self.shadow_palette, @@ -1465,7 +1491,7 @@ impl FractalApp { 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, + bailout_sq: bailout_sq(self.kind, self.power), kind: self.kind as u32, power: self.power, r_cap: self.buddha_r_cap, @@ -2028,7 +2054,11 @@ impl FractalApp { if self.anim.zoom && self.anim.zoom_speed != 0.0 { let max_hh = Scale::from_f64(DEFAULT_HALF_HEIGHT * 4.0); let factor = (-(self.anim.zoom_speed as f64) * dt).exp(); - let target = self.view.half_height.mul_f64(factor).clamp(Scale::MIN, max_hh); + let target = self + .view + .half_height + .mul_f64(factor) + .clamp(Scale::MIN, max_hh); let f = target.ratio(self.view.half_height); if (f - 1.0).abs() > 1.0e-9 { self.view @@ -2074,7 +2104,7 @@ impl FractalApp { } }); if self.kind == FractalKind::Multibrot { - ui.add(egui::Slider::new(&mut self.power, 2..=8).text("power")); + ui.add(egui::Slider::new(&mut self.power, 2..=20).text("power")); } if self.kind == FractalKind::Phoenix { ui.horizontal(|ui| { diff --git a/src/cli.rs b/src/cli.rs index 6ca797f..e9532c0 100644 --- a/src/cli.rs +++ b/src/cli.rs @@ -25,7 +25,7 @@ pub struct Cli { #[arg(long)] pub rendering_kind: Option, - /// Exponent for the Multibrot kind (z -> z^power + c), clamped to [2, 8]. + /// Exponent for the Multibrot kind (z -> z^power + c), clamped to [2, 20]. #[arg(long)] pub power: Option, diff --git a/src/fractal/kind.rs b/src/fractal/kind.rs index bdf8e17..fa7fc49 100644 --- a/src/fractal/kind.rs +++ b/src/fractal/kind.rs @@ -16,7 +16,7 @@ pub enum FractalKind { BurningShip = 1, /// `z -> conj(z)^2 + c` (the Mandelbar). Tricorn = 2, - /// `z -> z^power + c` (power >= 2). + /// `z -> z^power + c` (integer power in [2, 20]). Multibrot = 3, /// `z -> |Re(z^2)| + i·Im(z^2) + c` (abs on the real output of the square). Celtic = 4, diff --git a/src/fractal/reference.rs b/src/fractal/reference.rs index 3b1afc0..cfb4357 100644 --- a/src/fractal/reference.rs +++ b/src/fractal/reference.rs @@ -564,13 +564,17 @@ mod tests { fn f64_fast_path_matches_big() { let bits_fast = F64_MAX_PRECISION; let bits_big = F64_MAX_PRECISION + 64; - for kind in FractalKind::ALL { + let cases = FractalKind::ALL + .into_iter() + .map(|kind| (kind, 3)) + .chain([(FractalKind::Multibrot, 20)]); // highest supported power + for (kind, power) in cases { for julia in [false, true] { 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)); + let args = (60, bits, kind, power, (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, @@ -583,7 +587,7 @@ mod tests { ) } }; - let ctx = format!("{kind:?} julia={julia} morph={morph:?}"); + let ctx = format!("{kind:?} power={power} 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() { diff --git a/src/shaders/buddhabrot.wgsl b/src/shaders/buddhabrot.wgsl index c7d1a23..5dee1eb 100644 --- a/src/shaders/buddhabrot.wgsl +++ b/src/shaders/buddhabrot.wgsl @@ -106,7 +106,7 @@ fn advance(z: vec2, zp: vec2, c: vec2) -> vec2 { } else if KIND == KIND_TRICORN { return vec2(z.x * z.x - z.y * z.y, -2.0 * z.x * z.y) + c; } else if KIND == KIND_MULTIBROT { - return complex_pow(z, clamp(u.power, 2u, 8u)) + c; + return complex_pow(z, clamp(u.power, 2u, MULTIBROT_MAX_POWER)) + c; } else if KIND == KIND_CELTIC { return vec2(abs(z.x * z.x - z.y * z.y), 2.0 * z.x * z.y) + c; } else if KIND == KIND_PERPENDICULAR { diff --git a/src/shaders/common.wgsl b/src/shaders/common.wgsl index ce91242..9c38379 100644 --- a/src/shaders/common.wgsl +++ b/src/shaders/common.wgsl @@ -43,6 +43,10 @@ const KIND_MANDELBROT: u32 = 0u; const KIND_BURNING_SHIP: u32 = 1u; const KIND_TRICORN: u32 = 2u; const KIND_MULTIBROT: u32 = 3u; +// Highest Multibrot power (the UI/CLI/share-link clamp in app.rs matches). +// `bailout_sq` in app.rs shrinks the bailout with the power so z^p stays a +// finite f32. +const MULTIBROT_MAX_POWER: u32 = 20u; const KIND_CELTIC: u32 = 4u; const KIND_PERPENDICULAR: u32 = 5u; const KIND_BUFFALO: u32 = 6u; diff --git a/src/shaders/mandelbrot.wgsl b/src/shaders/mandelbrot.wgsl index 3410b27..bdef869 100644 --- a/src/shaders/mandelbrot.wgsl +++ b/src/shaders/mandelbrot.wgsl @@ -187,7 +187,7 @@ fn advance_delta_kind(kind: u32, z: vec2, e: vec2) -> vec2 { let ce = conj(e); return 2.0 * cmul(cz, ce) + cmul(ce, ce); } else if kind == KIND_MULTIBROT { - return multibrot_delta(z, e, clamp(u.power, 2u, 8u)); + return multibrot_delta(z, e, clamp(u.power, 2u, MULTIBROT_MAX_POWER)); } 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). @@ -222,7 +222,7 @@ fn advance_delta_kind(kind: u32, z: vec2, e: vec2) -> vec2 { // enough to de-speckle filaments. fn fprime_kind(kind: u32, z: vec2) -> vec2 { if kind == KIND_MULTIBROT { - let p = clamp(u.power, 2u, 8u); + let p = clamp(u.power, 2u, MULTIBROT_MAX_POWER); var zk = z; // Z^1 for (var k: u32 = 2u; k < p; k = k + 1u) { zk = cmul(zk, z); // -> Z^{p-1} @@ -382,7 +382,7 @@ fn advance_delta_scaled_kind(kind: u32, x: vec2, w: vec2, sc: f32, se: let cw = conj(w); return 2.0 * cmul(cx, cw) + sc * cmul(cw, cw); } else if kind == KIND_MULTIBROT { - return cmul(w, multibrot_sum(x, x + sc * w, clamp(u.power, 2u, 8u))); + return cmul(w, multibrot_sum(x, x + sc * w, clamp(u.power, 2u, MULTIBROT_MAX_POWER))); } else if kind == KIND_CELTIC { let sq = 2.0 * cmul(x, w) + sc * cmul(w, w); return vec2(diffabs_scaled(x.x * x.x - x.y * x.y, sq.x, se), sq.y); @@ -434,7 +434,7 @@ fn deep_step_kind(kind: u32, x: vec2, xf: Fe, w: vec2, sc: f32, s: i32 } var deg = 2; if kind == KIND_MULTIBROT { - deg = i32(clamp(u.power, 2u, 8u)); + deg = i32(clamp(u.power, 2u, MULTIBROT_MAX_POWER)); } let xk = ldexp2_sat(xf.m, xf.e - ue); let se = s - ue; @@ -545,7 +545,7 @@ fn deep_fprime(xf: Fe, w: vec2, s: i32, yt: vec2) -> Fe { return Fe(fprime(vec2(0.0, 0.0)), 0); } if KIND == KIND_MULTIBROT { - let p = clamp(u.power, 2u, 8u); + let p = clamp(u.power, 2u, MULTIBROT_MAX_POWER); var ym = yf.m; // m^(p-1) for (var k: u32 = 2u; k < p; k = k + 1u) { ym = cmul(ym, yf.m); diff --git a/src/view.rs b/src/view.rs index 6fff78e..6daae08 100644 --- a/src/view.rs +++ b/src/view.rs @@ -19,7 +19,6 @@ pub type Big = FBig; /// Half-height (complex units) of the default view; also the zoom-1 reference. pub const DEFAULT_HALF_HEIGHT: f64 = 1.25; - /// Below this pixel size (complex units per pixel) the GPU renders with the /// deep pipeline, whose per-pixel deltas start out as an f32 mantissa times /// `2^scale_exp`. Plain f32 stays exact as long as the smallest per-pixel @@ -89,7 +88,10 @@ impl Scale { return Self::MIN; } if m.is_infinite() { - return Scale { m: 1.0, e: i32::MAX / 2 }; + return Scale { + m: 1.0, + e: i32::MAX / 2, + }; } // Bring m into [1, 2) through its own binary exponent (exact). let k = m.log2().floor() as i32; @@ -175,11 +177,7 @@ impl Scale { impl PartialOrd for Scale { fn partial_cmp(&self, other: &Self) -> Option { // Normalized and positive: the exponent decides, then the mantissa. - Some( - self.e - .cmp(&other.e) - .then(self.m.partial_cmp(&other.m)?), - ) + Some(self.e.cmp(&other.e).then(self.m.partial_cmp(&other.m)?)) } } @@ -197,7 +195,12 @@ impl core::fmt::Display for Scale { } // Out of f64's range: round the exact decimal expansion instead. let sig = f.precision().map_or(17, |p| p + 1); - let dec = self.to_big().to_decimal().value().with_precision(sig).value(); + let dec = self + .to_big() + .to_decimal() + .value() + .with_precision(sig) + .value(); let repr = dec.repr(); let digits = repr.significand().to_string(); let digits = digits.trim_end_matches('0'); @@ -227,7 +230,11 @@ impl FromStr for Scale { if let Ok(x) = s.parse::() && x.is_normal() { - return if x > 0.0 { Ok(Self::from_f64(x)) } else { Err(()) }; + return if x > 0.0 { + Ok(Self::from_f64(x)) + } else { + Err(()) + }; } // Too small (or large) for f64: go through an exact decimal. let dec = DBig::from_str(s).map_err(|_| ())?; @@ -490,7 +497,8 @@ mod tests { #[test] fn interpolate_view_hits_exact_endpoints() { let bits = precision_for(sc(1.0)); - let from = ViewState::with_center(big_from_f64(-0.5, bits), big_from_f64(0.0, bits), sc(1.5)); + let from = + ViewState::with_center(big_from_f64(-0.5, bits), big_from_f64(0.0, bits), sc(1.5)); let to = ViewState::with_center( big_from_f64(-0.7515, precision_for(sc(1e-20))), big_from_f64(0.1013, precision_for(sc(1e-20))), @@ -515,7 +523,8 @@ mod tests { #[test] fn interpolate_view_keeps_target_offset_bounded() { let bits = precision_for(sc(1.0)); - let from = ViewState::with_center(big_from_f64(-0.5, bits), big_from_f64(0.0, bits), sc(1.5)); + let from = + ViewState::with_center(big_from_f64(-0.5, bits), big_from_f64(0.0, bits), sc(1.5)); let to = ViewState::with_center( big_from_f64(-0.7515, precision_for(sc(1e-20))), big_from_f64(0.1013, precision_for(sc(1e-20))), @@ -539,7 +548,14 @@ mod tests { #[test] fn scale_parse_display_round_trip() { - for s in ["1.25", "1e-20", "1.5e-20", "3.7e-4000", "1e-400", "9.99999e-310"] { + for s in [ + "1.25", + "1e-20", + "1.5e-20", + "3.7e-4000", + "1e-400", + "9.99999e-310", + ] { let a: Scale = s.parse().unwrap(); let b: Scale = a.to_string().parse().unwrap(); assert_eq!(a, b, "{s} -> {a}"); @@ -565,7 +581,10 @@ mod tests { assert_eq!(sc(1.0).exponent(), 0); assert_eq!(sc(0.75).exponent(), -1); assert_eq!(sc(4.0).scaled_f64(-2), 1.0); - assert_eq!(a.scaled_f64(-a.exponent()), a.mul_f64(1.0).scaled_f64(-a.exponent())); + assert_eq!( + a.scaled_f64(-a.exponent()), + a.mul_f64(1.0).scaled_f64(-a.exponent()) + ); assert!((1.0..2.0).contains(&a.scaled_f64(-a.exponent()))); assert_eq!(a.to_f64(), 0.0); assert_eq!(Scale::MIN.mul_f64(0.5), Scale::MIN); @@ -577,18 +596,11 @@ mod tests { /// same geometric pace as the half-height. #[test] fn interpolate_view_past_f64_range() { - let from = ViewState::with_center( - big_from_f64(-0.5, 64), - big_from_f64(0.0, 64), - sc(1.5), - ); + let from = ViewState::with_center(big_from_f64(-0.5, 64), big_from_f64(0.0, 64), sc(1.5)); let hh: Scale = "1e-1000".parse().unwrap(); let bits = precision_for(hh); - let to = ViewState::with_center( - big_from_f64(-0.7515, bits), - big_from_f64(0.1013, bits), - hh, - ); + let to = + ViewState::with_center(big_from_f64(-0.7515, bits), big_from_f64(0.1013, bits), hh); let end = interpolate_view(&from, &to, 1.0); assert_eq!(end.half_height, hh); assert_eq!(re_im_f64(&end), re_im_f64(&to));