feat: multibrot exponent up to 20

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, DEFAULT_HALF_HEIGHT, MAX_PRECISION_BITS, Scale, ViewState, big_from_decimal_str,
-    big_to_decimal_str, deep_scale_exp, interpolate_view, needs_deep, parse_view_spec,
+    big_from_f64, big_to_decimal_str, deep_scale_exp, interpolate_view, needs_deep,
-    precision_for,
+    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 dre = ((&self.view.center_re - &key.center_re) << k)
-        let dim = ((&self.view.center_im - &key.center_im) << k).to_f64().value();
+            .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| {

src/cli.rs

@@ -25,7 +25,7 @@ pub struct Cli {
     #[arg(long)]
     pub rendering_kind: Option<RenderingKindArg>,
 
-    /// 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<u32>,
 

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,

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() {

src/shaders/buddhabrot.wgsl

@@ -106,7 +106,7 @@ fn advance(z: vec2<f32>, zp: vec2<f32>, c: vec2<f32>) -> vec2<f32> {
     } else if KIND == KIND_TRICORN {
         return vec2<f32>(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<f32>(abs(z.x * z.x - z.y * z.y), 2.0 * z.x * z.y) + c;
     } else if KIND == KIND_PERPENDICULAR {

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;

src/shaders/mandelbrot.wgsl

@@ -187,7 +187,7 @@ fn advance_delta_kind(kind: u32, z: vec2<f32>, e: vec2<f32>) -> vec2<f32> {
         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<f32>, e: vec2<f32>) -> vec2<f32> {
 // enough to de-speckle filaments.
 fn fprime_kind(kind: u32, z: vec2<f32>) -> vec2<f32> {
     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<f32>, w: vec2<f32>, 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<f32>(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<f32>, xf: Fe, w: vec2<f32>, 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<f32>, s: i32, yt: vec2<f32>) -> Fe {
         return Fe(fprime(vec2<f32>(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);

src/view.rs

@@ -19,7 +19,6 @@ pub type Big = FBig<HalfAway, 2>;
 /// 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<core::cmp::Ordering> {
         // Normalized and positive: the exponent decides, then the mantissa.
-        Some(
+        Some(self.e.cmp(&other.e).then(self.m.partial_cmp(&other.m)?))
-            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::<f64>()
             && 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(
+        let from = ViewState::with_center(big_from_f64(-0.5, 64), big_from_f64(0.0, 64), sc(1.5));
-            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(
+        let to =
-            big_from_f64(-0.7515, bits),
+            ViewState::with_center(big_from_f64(-0.7515, bits), big_from_f64(0.1013, bits), hh);
-            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));