diff --git a/src/app.rs b/src/app.rs index a1de971..f527dcc 100644 --- a/src/app.rs +++ b/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>, + /// 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>, + /// 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, + 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, }, diff --git a/src/fractal/mod.rs b/src/fractal/mod.rs index 709f6e5..1f989d7 100644 --- a/src/fractal/mod.rs +++ b/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, diff --git a/src/fractal/reference.rs b/src/fractal/reference.rs index 697af30..fb7ad92 100644 --- a/src/fractal/reference.rs +++ b/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`, 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 { + 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 = 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![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"); + } } diff --git a/src/fractal/renderer.rs b/src/fractal/renderer.rs index 31c4bb0..ccedd5e 100644 --- a/src/fractal/renderer.rs +++ b/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::()) 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::()) 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>, + /// BLA table for the current reference (empty when BLA is off/unsupported). + pub bla: Arc>, 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; } diff --git a/src/shaders/mandelbrot.wgsl b/src/shaders/mandelbrot.wgsl index fbec524..17f08b4 100644 --- a/src/shaders/mandelbrot.wgsl +++ b/src/shaders/mandelbrot.wgsl @@ -29,6 +29,17 @@ struct Uniforms { dc_offset: vec2, // 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, + b: vec2, + r: f32, + l: u32, }; const KIND_MANDELBROT: u32 = 0u; @@ -38,6 +49,9 @@ const KIND_MULTIBROT: u32 = 3u; @group(0) @binding(0) var u: Uniforms; @group(0) @binding(1) var ref_orbit: array>; +// 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 bla_table: array; struct VsOut { @builtin(position) pos: vec4, @@ -169,6 +183,45 @@ fn palette(id: u32, t: f32) -> vec3 { return a + b * cos(6.28318530718 * (c * t + d)); } +// Result of a BLA lookup at an orbit index. +struct Hop { + found: bool, + a: vec2, + b: vec2, + 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, px: f32) -> vec3 { 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) { diff --git a/src/worker.rs b/src/worker.rs index 5ff4327..8355303 100644 --- a/src/worker.rs +++ b/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, } 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() {