feat: BLA iteration-skipping for deep-zoom Mandelbrot

src/app.rs

@@ -5,7 +5,7 @@ use eframe::egui_wgpu;
 use eframe::egui_wgpu::wgpu;
 
 use crate::fractal::{
-    ExportRender, FractalCallback, FractalKind, FractalRenderer, MAX_REF_POINTS, ShareState,
+    Bla, ExportRender, FractalCallback, FractalKind, FractalRenderer, MAX_REF_POINTS, ShareState,
     Uniforms,
 };
 #[cfg(target_arch = "wasm32")]
@@ -139,6 +139,11 @@ pub struct FractalApp {
 
     /// Reference orbit (`Z_n` as f32 pairs) for the current view.
     reference: Arc<Vec<[f32; 2]>>,
+    /// BLA iteration-skip table for `reference` (empty unless BLA is on and the
+    /// kind supports it). Uploaded to the GPU alongside the orbit.
+    bla: Arc<Vec<Bla>>,
+    /// Enable BLA iteration-skipping (square map only). Big speedup at deep zoom.
+    use_bla: bool,
     /// Bumped whenever `reference` is replaced, so the GPU re-uploads it.
     generation: u64,
     /// Center + zoom the current `reference` was computed at (may differ
@@ -236,6 +241,8 @@ impl FractalApp {
             antialias: false,
             de_coloring: false,
             reference: Arc::new(Vec::new()),
+            bla: Arc::new(Vec::new()),
+            use_bla: true,
             generation: 0,
             ref_center_re,
             ref_center_im,
@@ -300,6 +307,9 @@ impl FractalApp {
             if let Ok(spec) = std::env::var("MANDEL_VIEW") {
                 app.apply_view_spec(&spec);
             }
+            if std::env::var("MANDEL_NO_BLA").is_ok() {
+                app.use_bla = false; // for BLA-on vs BLA-off comparison
+            }
             if std::env::var("MANDEL_DE").is_ok() {
                 app.de_coloring = true;
             }
@@ -491,14 +501,34 @@ impl FractalApp {
         [dre, dim]
     }
 
-    fn apply_reference(&mut self, points: Vec<[f32; 2]>, cre: Big, cim: Big, hh: f64) {
+    fn apply_reference(
+        &mut self,
+        points: Vec<[f32; 2]>,
+        bla: Vec<Bla>,
+        cre: Big,
+        cim: Big,
+        hh: f64,
+    ) {
         self.reference = Arc::new(points);
+        self.bla = Arc::new(bla);
         self.ref_center_re = cre;
         self.ref_center_im = cim;
         self.ref_half_height = hh;
         self.generation = self.generation.wrapping_add(1);
     }
 
+    /// Upper bound on any pixel's `|dc|` for the current view (center→corner
+    /// distance), with margin so a reference reused after slight drift/zoom stays
+    /// valid. Used to size BLA merge radii.
+    fn bla_dc_max(&self, half_height: f64) -> f64 {
+        let aspect = if self.last_size_px.y > 0.0 {
+            (self.last_size_px.x / self.last_size_px.y) as f64
+        } else {
+            1.0
+        };
+        half_height * (1.0 + aspect * aspect).sqrt() * 2.0
+    }
+
     /// Recompute the reference orbit when needed. Native: dispatch to a worker
     /// thread and pick up completed results. Web: compute inline.
     fn ensure_reference(&mut self) {
@@ -506,6 +536,9 @@ impl FractalApp {
             let key = self.current_key();
             let precision = self.view.precision_bits();
             let max_iter = key.iter.min(MAX_REF_POINTS as u32 - 1);
+            // BLA applies to the holomorphic square map (Mandelbrot set + Julia).
+            let want_bla = self.use_bla && key.kind == FractalKind::Mandelbrot;
+            let dc_max = self.bla_dc_max(key.half_height);
 
             #[cfg(not(target_arch = "wasm32"))]
             {
@@ -519,6 +552,8 @@ impl FractalApp {
                     precision,
                     kind: key.kind,
                     power: key.power,
+                    build_bla: want_bla,
+                    dc_max,
                 });
                 self.pending = true;
             }
@@ -547,8 +582,14 @@ impl FractalApp {
                         key.power,
                     )
                 };
+                let bla = if want_bla {
+                    crate::fractal::build_bla_table(&points, dc_max)
+                } else {
+                    Vec::new()
+                };
                 self.apply_reference(
                     points,
+                    bla,
                     key.center_re.clone(),
                     key.center_im.clone(),
                     key.half_height,
@@ -560,7 +601,13 @@ impl FractalApp {
 
         #[cfg(not(target_arch = "wasm32"))]
         if let Some(res) = self.worker.try_take_latest() {
-            self.apply_reference(res.points, res.center_re, res.center_im, res.half_height);
+            self.apply_reference(
+                res.points,
+                res.bla,
+                res.center_re,
+                res.center_im,
+                res.half_height,
+            );
             self.pending = false;
         }
     }
@@ -581,7 +628,11 @@ impl FractalApp {
             power: self.power,
             dc_offset: self.dc_offset(),
             de_coloring: self.de_coloring as u32,
-            _pad: 0,
+            // Only when a table for the current reference is actually present, so
+            // a stale/empty buffer is never traversed during a kind switch.
+            use_bla: (self.use_bla
+                && self.kind == FractalKind::Mandelbrot
+                && !self.bla.is_empty()) as u32,
         }
     }
 
@@ -618,6 +669,7 @@ impl FractalApp {
             renderer.export_handles()
         };
         let reference = Arc::clone(&self.reference);
+        let bla = Arc::clone(&self.bla);
 
         let shared = Arc::new(Mutex::new(ExportShared {
             fraction: 0.0,
@@ -645,6 +697,7 @@ impl FractalApp {
                     h,
                     uniforms,
                     reference.as_slice(),
+                    bla.as_slice(),
                 );
 
                 // Render the image tile by tile, waiting for each so progress
@@ -710,6 +763,7 @@ impl FractalApp {
                     h,
                     uniforms,
                     reference.as_slice(),
+                    bla.as_slice(),
                 );
 
                 // Render tile by tile, awaiting each submission so the browser
@@ -871,6 +925,14 @@ impl FractalApp {
                  for crisp filaments at deep zoom. Exact for Mandelbrot/Multibrot, \
                  approximate for Burning Ship/Tricorn.",
             );
+        ui.add_enabled(
+            self.kind == FractalKind::Mandelbrot,
+            egui::Checkbox::new(&mut self.use_bla, "BLA skipping"),
+        )
+        .on_hover_text(
+            "Bivariate Linear Approximation: skip runs of iterations at deep zoom \
+             for a large speedup. Mandelbrot/Julia only.",
+        );
 
         ui.separator();
         // Editable center coordinates. Shown at full precision; parsed
@@ -1104,6 +1166,7 @@ impl FractalApp {
             FractalCallback {
                 uniforms,
                 reference: Arc::clone(&self.reference),
+                bla: Arc::clone(&self.bla),
                 generation: self.generation,
                 size_px,
             },

src/fractal/mod.rs

@@ -5,7 +5,7 @@ pub mod reference;
 pub mod renderer;
 pub mod share;
 
-pub use reference::{FractalKind, compute_reference, compute_set_reference};
+pub use reference::{Bla, FractalKind, build_bla_table, compute_reference, compute_set_reference};
 pub use renderer::{
     ExportRender, FractalCallback, FractalRenderer, MAX_REF_POINTS, Uniforms,
     encode_png_with_progress,

src/fractal/reference.rs

@@ -146,6 +146,150 @@ pub fn compute_set_reference(
     )
 }
 
+// ---------------------------------------------------------------------------
+// Bivariate Linear Approximation (BLA)
+//
+// Deep in a zoom the per-pixel delta stays far smaller than the reference, so
+// the nonlinear `e^2` term of the perturbation step is negligible and the step
+// is effectively linear: `e -> A e + B dc`. BLA precomputes, for runs of
+// iterations, the composed linear coefficients `(A, B)` plus a validity radius
+// `r` (the largest `|e|` for which dropping `e^2` stays within tolerance). A
+// pixel can then skip a whole run in one multiply whenever `|e| < r`.
+//
+// Runs are merged pairwise into a binary tree of levels: level `k` holds BLAs of
+// length `2^k` starting at multiples of `2^k`. The GPU walks levels high→low to
+// take the longest valid skip at the current index. Both sides recompute the
+// per-level counts from `ref_len` (count[0] = ref_len-1, count[k] = count[k-1]/2)
+// so no offset table needs to travel to the GPU — only this flat array does.
+//
+// Only the holomorphic square map (Mandelbrot/Julia, `A = 2 Z`, `B = 1`) is
+// supported; other kinds fall back to per-iteration stepping on the GPU.
+// ---------------------------------------------------------------------------
+
+/// One merged linear step: `e_{n+l} = A e_n + B dc`, valid while `|e_n| < r`.
+/// Laid out to match the WGSL `Bla` struct (two `vec2<f32>`, then `f32`, `u32`;
+/// 24-byte std430 stride).
+#[repr(C)]
+#[derive(Clone, Copy, bytemuck::Pod, bytemuck::Zeroable)]
+pub struct Bla {
+    pub a: [f32; 2],
+    pub b: [f32; 2],
+    pub r: f32,
+    pub l: u32,
+}
+
+/// Relative budget for the dropped nonlinear term, chosen near the f32 orbit
+/// storage noise floor so BLA adds no visible error over plain perturbation.
+const BLA_EPS: f64 = 1.0e-6;
+
+/// f64 working form of a BLA (merges accumulate in f64, stored as f32).
+#[derive(Clone, Copy)]
+struct BlaF {
+    ar: f64,
+    ai: f64,
+    br: f64,
+    bi: f64,
+    r: f64,
+    l: u32,
+}
+
+impl BlaF {
+    fn to_bla(self) -> Bla {
+        Bla {
+            a: [self.ar as f32, self.ai as f32],
+            b: [self.br as f32, self.bi as f32],
+            r: self.r as f32,
+            l: self.l,
+        }
+    }
+}
+
+/// Merge two consecutive BLAs (`x` then `y`) into one covering both runs.
+fn merge_bla(x: BlaF, y: BlaF, dc_max: f64) -> BlaF {
+    // A = Ay Ax ; B = Ay Bx + By  (complex).
+    let ar = y.ar * x.ar - y.ai * x.ai;
+    let ai = y.ar * x.ai + y.ai * x.ar;
+    let br = y.ar * x.br - y.ai * x.bi + y.br;
+    let bi = y.ar * x.bi + y.ai * x.br + y.bi;
+    // Valid if |e| < rx (so x holds) and |Ax e + Bx dc| < ry (so y holds):
+    //   |e| < (ry - |Bx| dc_max) / |Ax|.
+    let ax_mag = (x.ar * x.ar + x.ai * x.ai).sqrt();
+    let bx_mag = (x.br * x.br + x.bi * x.bi).sqrt();
+    let r_y = if ax_mag > 0.0 {
+        ((y.r - bx_mag * dc_max) / ax_mag).max(0.0)
+    } else {
+        x.r
+    };
+    BlaF {
+        ar,
+        ai,
+        br,
+        bi,
+        r: x.r.min(r_y),
+        l: x.l + y.l,
+    }
+}
+
+/// Build the BLA table for a square-map reference orbit. `dc_max` is an upper
+/// bound on any pixel's `|dc|` in the current view (used to size merge radii).
+/// Returns a flat array with levels concatenated (level 0 first). Empty if the
+/// orbit is too short to skip.
+pub fn build_bla_table(points: &[[f32; 2]], dc_max: f64) -> Vec<Bla> {
+    let m = points.len();
+    if m < 2 {
+        return Vec::new();
+    }
+
+    // Level 0: one single step from each index i (uses Z_i). A = 2 Z, B = 1.
+    // Dropping e^2 is within tolerance while |e| < BLA_EPS |Z| (since the kept
+    // linear term is |2 Z e|).
+    let mut level0: Vec<BlaF> = Vec::with_capacity(m - 1);
+    for z in &points[..m - 1] {
+        let (zr, zi) = (z[0] as f64, z[1] as f64);
+        let zmag = (zr * zr + zi * zi).sqrt();
+        level0.push(BlaF {
+            ar: 2.0 * zr,
+            ai: 2.0 * zi,
+            br: 1.0,
+            bi: 0.0,
+            r: BLA_EPS * zmag,
+            l: 1,
+        });
+    }
+
+    let mut levels: Vec<Vec<BlaF>> = vec![level0];
+    while levels.last().unwrap().len() >= 2 {
+        let prev = levels.last().unwrap();
+        let mut next = Vec::with_capacity(prev.len() / 2);
+        let mut i = 0;
+        while i + 1 < prev.len() {
+            next.push(merge_bla(prev[i], prev[i + 1], dc_max));
+            i += 2;
+        }
+        levels.push(next);
+    }
+
+    let mut flat = Vec::with_capacity(levels.iter().map(|l| l.len()).sum());
+    for level in &levels {
+        flat.extend(level.iter().map(|b| b.to_bla()));
+    }
+    flat
+}
+
+/// Start offset of BLA level `lv` within the flat table, computed from the orbit
+/// length exactly as the shader does (`count[0] = ref_len-1`, halving each
+/// level). Kept here so the traversal test mirrors the GPU indexing.
+#[cfg(test)]
+fn bla_level_start(ref_len: usize, lv: u32) -> (usize, usize) {
+    let mut start = 0usize;
+    let mut count = ref_len - 1;
+    for _ in 0..lv {
+        start += count;
+        count /= 2;
+    }
+    (start, count)
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -256,4 +400,126 @@ mod tests {
             zi = nzi;
         }
     }
+
+    /// Step-by-step perturbation (drops nothing): `e_{n+1} = 2 Z_n e_n + e_n^2 + dc`,
+    /// using the stored f32 orbit as `Z_n`. Returns `e` after `target_n` steps.
+    fn advance_naive(points: &[[f32; 2]], dc: (f64, f64), target_n: usize) -> (f64, f64) {
+        let (mut er, mut ei) = (0.0f64, 0.0f64);
+        for z in &points[..target_n] {
+            let (zr, zi) = (z[0] as f64, z[1] as f64);
+            let tr = 2.0 * (zr * er - zi * ei);
+            let ti = 2.0 * (zr * ei + zi * er);
+            let sr = er * er - ei * ei;
+            let si = 2.0 * er * ei;
+            er = tr + sr + dc.0;
+            ei = ti + si + dc.1;
+        }
+        (er, ei)
+    }
+
+    /// Advance `e` to exactly `target_n` steps using the BLA table — the same
+    /// walk the shader performs (longest valid skip first, else a full step),
+    /// but never skipping past `target_n`. Returns `(e, took_a_multi_step_skip)`.
+    fn advance_bla(
+        points: &[[f32; 2]],
+        table: &[Bla],
+        dc: (f64, f64),
+        target_n: usize,
+    ) -> ((f64, f64), bool) {
+        let ref_len = points.len();
+        let (mut er, mut ei) = (0.0f64, 0.0f64);
+        let mut n = 0usize;
+        let mut skipped = false;
+
+        while n < target_n {
+            let emag2 = er * er + ei * ei;
+            let mut applied = false;
+
+            // Highest level worth trying is bounded by how far we may advance.
+            let span = target_n - n;
+            let max_lv = (usize::BITS - 1 - span.leading_zeros()) as i64; // floor(log2(span))
+            let mut lv = max_lv;
+            while lv >= 0 {
+                let lvu = lv as u32;
+                let step = 1usize << lvu;
+                if n % step == 0 && n + step <= target_n {
+                    let (start, count) = bla_level_start(ref_len, lvu);
+                    let idx = n >> lvu;
+                    if idx < count {
+                        let b = table[start + idx];
+                        let r = b.r as f64;
+                        if emag2 < r * r {
+                            let (ar, ai) = (b.a[0] as f64, b.a[1] as f64);
+                            let (br, bi) = (b.b[0] as f64, b.b[1] as f64);
+                            let ner = ar * er - ai * ei + br * dc.0 - bi * dc.1;
+                            let nei = ar * ei + ai * er + br * dc.1 + bi * dc.0;
+                            er = ner;
+                            ei = nei;
+                            n += step;
+                            skipped |= step > 1;
+                            applied = true;
+                            break;
+                        }
+                    }
+                }
+                lv -= 1;
+            }
+
+            if !applied {
+                let (zr, zi) = (points[n][0] as f64, points[n][1] as f64);
+                let tr = 2.0 * (zr * er - zi * ei);
+                let ti = 2.0 * (zr * ei + zi * er);
+                let sr = er * er - ei * ei;
+                let si = 2.0 * er * ei;
+                er = tr + sr + dc.0;
+                ei = ti + si + dc.1;
+                n += 1;
+            }
+        }
+        ((er, ei), skipped)
+    }
+
+    /// The BLA walk must reproduce step-by-step perturbation at deep zoom (where
+    /// the delta is tiny and skips actually fire).
+    #[test]
+    fn bla_matches_step_by_step() {
+        // An interior center (never escapes), so `e` stays bounded and we can
+        // iterate the full orbit; zoomed so |dc| ~ 1e-9 (deep enough for big
+        // skips). Correctness of the walk is independent of which orbit we pick.
+        let cr = Big::try_from(-0.5_f64).unwrap();
+        let ci = Big::try_from(0.0_f64).unwrap();
+        let points = compute_set_reference(&cr, &ci, 800, 160, FractalKind::Mandelbrot, 2);
+        assert!(points.len() > 64, "need a long orbit to exercise BLA levels");
+
+        let half_height = 1.0e-9_f64;
+        let aspect = 1.5_f64;
+        let dc_max = half_height * (1.0 + aspect * aspect).sqrt();
+        let table = build_bla_table(&points, dc_max);
+        assert!(!table.is_empty());
+
+        let target_n = points.len() - 1;
+        // A few pixel offsets across the view (all within dc_max).
+        let offsets = [
+            (0.0, 0.0),
+            (0.6e-9, -0.4e-9),
+            (-0.9e-9, 0.3e-9),
+            (0.2e-9, 0.8e-9),
+        ];
+        let mut any_skip = false;
+        for dc in offsets {
+            let naive = advance_naive(&points, dc, target_n);
+            let (bla, skipped) = advance_bla(&points, &table, dc, target_n);
+            any_skip |= skipped;
+
+            let ref_mag = (naive.0 * naive.0 + naive.1 * naive.1).sqrt();
+            let err = ((bla.0 - naive.0).powi(2) + (bla.1 - naive.1).powi(2)).sqrt();
+            // Only the dropped e^2 differs; must stay near the BLA_EPS budget.
+            let tol = 1e-4 * ref_mag + 1e-15;
+            assert!(
+                err <= tol,
+                "dc={dc:?}: BLA {bla:?} vs naive {naive:?} (err {err:e} > tol {tol:e})"
+            );
+        }
+        assert!(any_skip, "BLA never took a multi-step skip — test is not exercising it");
+    }
 }

src/fractal/renderer.rs

@@ -13,10 +13,81 @@ use std::sync::Arc;
 
 use eframe::egui_wgpu::{self, wgpu};
 
+use super::reference::Bla;
+
 /// Maximum reference-orbit length (points) the storage buffer can hold. Also
 /// bounds the iteration count. 128k points * 8 bytes = 1 MiB.
 pub const MAX_REF_POINTS: usize = 1 << 17;
 
+/// Maximum BLA-table entries the storage buffer can hold. A full table has
+/// ~2× the orbit length (all levels summed); ~256k entries * 24 bytes = 6 MiB.
+pub const MAX_BLA_ENTRIES: usize = 2 * MAX_REF_POINTS;
+
+/// The fractal pipeline's bind group layout: uniforms (0), reference orbit (1),
+/// BLA table (2). Shared by the live renderer and [`ExportRender`].
+fn fractal_bind_group_layout_desc<'a>() -> wgpu::BindGroupLayoutDescriptor<'a> {
+    const STORAGE: wgpu::BindingType = wgpu::BindingType::Buffer {
+        ty: wgpu::BufferBindingType::Storage { read_only: true },
+        has_dynamic_offset: false,
+        min_binding_size: None,
+    };
+    wgpu::BindGroupLayoutDescriptor {
+        label: Some("fractal bind group layout"),
+        entries: &[
+            wgpu::BindGroupLayoutEntry {
+                binding: 0,
+                visibility: wgpu::ShaderStages::FRAGMENT,
+                ty: wgpu::BindingType::Buffer {
+                    ty: wgpu::BufferBindingType::Uniform,
+                    has_dynamic_offset: false,
+                    min_binding_size: None,
+                },
+                count: None,
+            },
+            wgpu::BindGroupLayoutEntry {
+                binding: 1,
+                visibility: wgpu::ShaderStages::FRAGMENT,
+                ty: STORAGE,
+                count: None,
+            },
+            wgpu::BindGroupLayoutEntry {
+                binding: 2,
+                visibility: wgpu::ShaderStages::FRAGMENT,
+                ty: STORAGE,
+                count: None,
+            },
+        ],
+    }
+}
+
+/// Build the fractal bind group from its three buffers.
+fn fractal_bind_group(
+    device: &wgpu::Device,
+    layout: &wgpu::BindGroupLayout,
+    uniform_buffer: &wgpu::Buffer,
+    ref_buffer: &wgpu::Buffer,
+    bla_buffer: &wgpu::Buffer,
+) -> wgpu::BindGroup {
+    device.create_bind_group(&wgpu::BindGroupDescriptor {
+        label: Some("fractal bind group"),
+        layout,
+        entries: &[
+            wgpu::BindGroupEntry {
+                binding: 0,
+                resource: uniform_buffer.as_entire_binding(),
+            },
+            wgpu::BindGroupEntry {
+                binding: 1,
+                resource: ref_buffer.as_entire_binding(),
+            },
+            wgpu::BindGroupEntry {
+                binding: 2,
+                resource: bla_buffer.as_entire_binding(),
+            },
+        ],
+    })
+}
+
 /// GPU-side view + coloring parameters. Layout must match `Uniforms` in the
 /// WGSL shader; total size is a multiple of 16 bytes for uniform-buffer rules.
 #[repr(C)]
@@ -45,8 +116,9 @@ pub struct Uniforms {
     pub dc_offset: [f32; 2],
     /// 0 = escape-time coloring, 1 = distance-estimation shading.
     pub de_coloring: u32,
-    /// Padding to a 16-byte multiple (uniform buffer requirement).
-    pub _pad: u32,
+    /// 0 = per-iteration stepping, 1 = BLA iteration-skipping (square map only).
+    /// This also fills the struct out to a 16-byte multiple (uniform requirement).
+    pub use_bla: u32,
 }
 
 /// Offscreen texture the fractal is rendered into, plus the bind group used to
@@ -72,9 +144,10 @@ pub struct FractalRenderer {
     bind_group_layout: wgpu::BindGroupLayout,
     uniform_buffer: wgpu::Buffer,
     ref_buffer: wgpu::Buffer,
+    bla_buffer: wgpu::Buffer,
     bind_group: wgpu::BindGroup,
     target_format: wgpu::TextureFormat,
-    /// Generation of the reference orbit currently uploaded to `ref_buffer`.
+    /// Generation of the reference orbit + BLA table currently uploaded.
     uploaded_generation: u64,
 
     /// Blit pipeline + resources that copy the cache texture to egui's surface.
@@ -108,46 +181,23 @@ impl FractalRenderer {
             mapped_at_creation: false,
         });
 
-        let bind_group_layout = device.create_bind_group_layout(&wgpu::BindGroupLayoutDescriptor {
-            label: Some("fractal bind group layout"),
-            entries: &[
-                wgpu::BindGroupLayoutEntry {
-                    binding: 0,
-                    visibility: wgpu::ShaderStages::FRAGMENT,
-                    ty: wgpu::BindingType::Buffer {
-                        ty: wgpu::BufferBindingType::Uniform,
-                        has_dynamic_offset: false,
-                        min_binding_size: None,
-                    },
-                    count: None,
-                },
-                wgpu::BindGroupLayoutEntry {
-                    binding: 1,
-                    visibility: wgpu::ShaderStages::FRAGMENT,
-                    ty: wgpu::BindingType::Buffer {
-                        ty: wgpu::BufferBindingType::Storage { read_only: true },
-                        has_dynamic_offset: false,
-                        min_binding_size: None,
-                    },
-                    count: None,
-                },
-            ],
+        let bla_buffer = device.create_buffer(&wgpu::BufferDescriptor {
+            label: Some("bla table"),
+            size: (MAX_BLA_ENTRIES * std::mem::size_of::<Bla>()) as u64,
+            usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_DST,
+            mapped_at_creation: false,
         });
 
-        let bind_group = device.create_bind_group(&wgpu::BindGroupDescriptor {
-            label: Some("fractal bind group"),
-            layout: &bind_group_layout,
-            entries: &[
-                wgpu::BindGroupEntry {
-                    binding: 0,
-                    resource: uniform_buffer.as_entire_binding(),
-                },
-                wgpu::BindGroupEntry {
-                    binding: 1,
-                    resource: ref_buffer.as_entire_binding(),
-                },
-            ],
-        });
+        let bind_group_layout =
+            device.create_bind_group_layout(&fractal_bind_group_layout_desc());
+
+        let bind_group = fractal_bind_group(
+            device,
+            &bind_group_layout,
+            &uniform_buffer,
+            &ref_buffer,
+            &bla_buffer,
+        );
 
         let pipeline_layout = device.create_pipeline_layout(&wgpu::PipelineLayoutDescriptor {
             label: Some("fractal pipeline layout"),
@@ -254,6 +304,7 @@ impl FractalRenderer {
             bind_group_layout,
             uniform_buffer,
             ref_buffer,
+            bla_buffer,
             bind_group,
             target_format,
             uploaded_generation: u64::MAX,
@@ -361,6 +412,7 @@ impl ExportRender {
         height: u32,
         uniforms: Uniforms,
         reference: &[[f32; 2]],
+        bla: &[Bla],
     ) -> Self {
         let uniform_buffer = device.create_buffer(&wgpu::BufferDescriptor {
             label: Some("export uniforms"),
@@ -381,20 +433,24 @@ impl ExportRender {
             queue.write_buffer(&ref_buffer, 0, bytemuck::cast_slice(&reference[..count]));
         }
 
-        let bind_group = device.create_bind_group(&wgpu::BindGroupDescriptor {
-            label: Some("export bind group"),
-            layout: bind_group_layout,
-            entries: &[
-                wgpu::BindGroupEntry {
-                    binding: 0,
-                    resource: uniform_buffer.as_entire_binding(),
-                },
-                wgpu::BindGroupEntry {
-                    binding: 1,
-                    resource: ref_buffer.as_entire_binding(),
-                },
-            ],
+        let bla_count = bla.len().min(MAX_BLA_ENTRIES);
+        let bla_buffer = device.create_buffer(&wgpu::BufferDescriptor {
+            label: Some("export bla table"),
+            size: (bla_count.max(1) * std::mem::size_of::<Bla>()) as u64,
+            usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_DST,
+            mapped_at_creation: false,
         });
+        if bla_count > 0 {
+            queue.write_buffer(&bla_buffer, 0, bytemuck::cast_slice(&bla[..bla_count]));
+        }
+
+        let bind_group = fractal_bind_group(
+            device,
+            bind_group_layout,
+            &uniform_buffer,
+            &ref_buffer,
+            &bla_buffer,
+        );
 
         let texture = device.create_texture(&wgpu::TextureDescriptor {
             label: Some("export target"),
@@ -583,6 +639,8 @@ pub fn encode_png_with_progress(
 pub struct FractalCallback {
     pub uniforms: Uniforms,
     pub reference: Arc<Vec<[f32; 2]>>,
+    /// BLA table for the current reference (empty when BLA is off/unsupported).
+    pub bla: Arc<Vec<Bla>>,
     pub generation: u64,
     /// Widget size in physical pixels — the cache texture resolution.
     pub size_px: [u32; 2],
@@ -612,6 +670,14 @@ impl egui_wgpu::CallbackTrait for FractalCallback {
                 0,
                 bytemuck::cast_slice(&self.reference[..count]),
             );
+            let bla_count = self.bla.len().min(MAX_BLA_ENTRIES);
+            if bla_count > 0 {
+                queue.write_buffer(
+                    &renderer.bla_buffer,
+                    0,
+                    bytemuck::cast_slice(&self.bla[..bla_count]),
+                );
+            }
             renderer.uploaded_generation = self.generation;
         }
 

src/shaders/mandelbrot.wgsl

@@ -29,6 +29,17 @@ struct Uniforms {
     dc_offset: vec2<f32>,
     // 0 = escape-time coloring, 1 = distance-estimation shading.
     de_coloring: u32,
+    // 0 = per-iteration stepping, 1 = BLA iteration-skipping (square map only).
+    use_bla: u32,
+};
+
+// One merged linear step from the BLA table: e_{n+l} = A e + B dc, valid while
+// |e| < r. Matches the Rust `Bla` struct (24-byte std430 stride).
+struct Bla {
+    a: vec2<f32>,
+    b: vec2<f32>,
+    r: f32,
+    l: u32,
 };
 
 const KIND_MANDELBROT: u32 = 0u;
@@ -38,6 +49,9 @@ const KIND_MULTIBROT: u32 = 3u;
 
 @group(0) @binding(0) var<uniform> u: Uniforms;
 @group(0) @binding(1) var<storage, read> ref_orbit: array<vec2<f32>>;
+// BLA table, levels concatenated (level 0 first). Per-level counts are recomputed
+// from `ref_len` exactly as the CPU builder laid them out.
+@group(0) @binding(2) var<storage, read> bla_table: array<Bla>;
 
 struct VsOut {
     @builtin(position) pos: vec4<f32>,
@@ -169,6 +183,45 @@ fn palette(id: u32, t: f32) -> vec3<f32> {
     return a + b * cos(6.28318530718 * (c * t + d));
 }
 
+// Result of a BLA lookup at an orbit index.
+struct Hop {
+    found: bool,
+    a: vec2<f32>,
+    b: vec2<f32>,
+    l: u32,
+};
+
+// Longest valid BLA skip starting at orbit index `n` for a delta of squared
+// magnitude `emag2`. Walks levels low->high, recomputing each level's flat-array
+// start and count from `ref_len` (count[0] = ref_len-1, halving each level).
+// Radii shrink with level and alignment is monotonic, so the valid levels form a
+// prefix: we keep the last valid one and stop at the first that fails.
+fn bla_find(n: u32, emag2: f32) -> Hop {
+    var res: Hop;
+    res.found = false;
+
+    var start: u32 = 0u;
+    var count: u32 = u.ref_len - 1u;
+    var lv: u32 = 0u;
+    loop {
+        if (count == 0u) { break; }
+        let step = 1u << lv;
+        if ((n & (step - 1u)) != 0u) { break; } // n not aligned to this level
+        let idx = n >> lv;
+        if (idx >= count) { break; }            // run would exceed the orbit
+        let b = bla_table[start + idx];
+        if (emag2 >= b.r * b.r) { break; }       // delta too large: not valid
+        res.found = true;
+        res.a = b.a;
+        res.b = b.b;
+        res.l = b.l;
+        start = start + count;
+        count = count / 2u;
+        lv = lv + 1u;
+    }
+    return res;
+}
+
 // Perturbation iterate + color a single sample. `offset` is the per-pixel
 // offset in complex units. For Mandelbrot it is the c-plane offset added every
 // step (delta starts at 0); for Julia it is the z-plane offset that seeds the
@@ -208,16 +261,34 @@ fn shade(offset: vec2<f32>, px: f32) -> vec3<f32> {
             break; // interior
         }
 
-        // Propagate the derivative of the full orbit (unaffected by rebasing,
-        // which only re-expresses the same value). Only when DE is enabled.
-        if (u.de_coloring != 0u) {
-            dz = cmul(fprime(z), dz) + dz_seed;
+        // Advance one step, or skip a whole run via BLA when the delta is small
+        // enough (square map only; other kinds keep `use_bla == 0`). The orbit
+        // derivative for DE follows the same linear map: over a run it advances
+        // by the run's own (A, B) coefficients, matching the per-step recurrence.
+        var did_skip = false;
+        if (u.use_bla != 0u) {
+            let hop = bla_find(m, dot(e, e));
+            if (hop.found) {
+                if (u.de_coloring != 0u) {
+                    dz = cmul(hop.a, dz) + cmul(hop.b, dz_seed);
+                }
+                e = cmul(hop.a, e) + cmul(hop.b, step_add);
+                m = m + hop.l;
+                n = n + hop.l;
+                did_skip = true;
+            }
+        }
+        if (!did_skip) {
+            // Propagate the derivative of the full orbit (unaffected by rebasing,
+            // which only re-expresses the same value). Only when DE is enabled.
+            if (u.de_coloring != 0u) {
+                dz = cmul(fprime(z), dz) + dz_seed;
+            }
+            // 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;
         }
-
-        // 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;
 
         // Keep the reference index valid and the delta small.
         if (m >= u.ref_len) {

src/worker.rs

@@ -9,7 +9,7 @@
 use std::sync::mpsc::{Receiver, Sender, TryRecvError, channel};
 use std::thread;
 
-use crate::fractal::{FractalKind, compute_reference, compute_set_reference};
+use crate::fractal::{Bla, FractalKind, build_bla_table, compute_reference, compute_set_reference};
 use crate::view::{Big, big_from_f64};
 
 pub struct RefRequest {
@@ -22,6 +22,10 @@ pub struct RefRequest {
     pub precision: usize,
     pub kind: FractalKind,
     pub power: u32,
+    /// Build the BLA iteration-skip table for this orbit (square map only).
+    pub build_bla: bool,
+    /// Upper bound on any pixel's `|dc|`, sizing the BLA merge radii.
+    pub dc_max: f64,
 }
 
 pub struct RefResult {
@@ -29,6 +33,7 @@ pub struct RefResult {
     pub center_im: Big,
     pub half_height: f64,
     pub points: Vec<[f32; 2]>,
+    pub bla: Vec<Bla>,
 }
 
 pub struct RefWorker {
@@ -56,12 +61,18 @@ impl RefWorker {
                     }
 
                     let points = compute(&req);
+                    let bla = if req.build_bla {
+                        build_bla_table(&points, req.dc_max)
+                    } else {
+                        Vec::new()
+                    };
                     if res_tx
                         .send(RefResult {
                             center_re: req.center_re,
                             center_im: req.center_im,
                             half_height: req.half_height,
                             points,
+                            bla,
                         })
                         .is_err()
                     {