From bfd3d18f3d7fe415a6dcd8a23aed0615ca186ba5 Mon Sep 17 00:00:00 2001 From: supersurviveur Date: Fri, 25 Sep 2026 15:19:33 +0200 Subject: [PATCH] feat: Add arbitrary precision numbers at deep zooms --- CLAUDE.md | 65 +++- src/app.rs | 87 +++-- src/fractal/mod.rs | 2 +- src/fractal/reference.rs | 150 +++++++- src/fractal/renderer.rs | 64 +++- src/headless.rs | 6 +- src/shaders/iterate_uniforms.wgsl | 4 + src/shaders/mandelbrot.wgsl | 596 +++++++++++++++++++++++++++++- src/view.rs | 43 ++- src/worker.rs | 6 +- tests/shader_valid.rs | 35 +- 11 files changed, 954 insertions(+), 104 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index c8e156c..1cbe8f4 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -8,7 +8,11 @@ A deep-zoom fractal explorer (Rust + wgpu + egui + WGSL). It zooms past the ~10¹³× limit of plain `f64` using **perturbation theory**: one high-precision reference orbit is computed on the CPU (arbitrary precision via `dashu-float`), and every pixel is rendered on the GPU as a cheap `f32` delta from it, with -rebasing to avoid glitches. The `f32` GPU tier reaches roughly 10³⁰×. Runs +rebasing to avoid glitches. Plain `f32` deltas run out of exponent range +once a pixel is ~2^-124 wide (~10³⁴× at 1080p), so from 2^-122 per pixel +(`view::DEEP_PIXEL_SIZE`) a `DEEP` shader variant starts each pixel with +rescaled deltas (f32 mantissa × 2^i32), reaching ~10³⁰⁰× (the `f64` limit of +`half_height`, `view::MIN_HALF_HEIGHT`). Runs natively (Vulkan/Metal/DX12) and in the browser (WebGPU only — WebGL2 can't do storage buffers, which the fragment shader needs for the reference orbit). @@ -90,8 +94,15 @@ what makes deep zoom cheap — one expensive high-precision orbit, then every pixel is a handful of `f32` complex multiplies. - `src/view.rs` — `ViewState`; center is arbitrary-precision `FBig` (`Big` - type alias), pixel scale stays `f64` (still in-range at 10³⁰×). Precision - (bits) scales with zoom depth (`precision_for`). + type alias), pixel scale stays `f64` (so zoom is clamped at + `MIN_HALF_HEIGHT` = 1e-300). Precision (bits) scales with zoom depth + (`precision_for`). `needs_deep` switches rendering to the deep pipeline + once a pixel of the full-resolution render is below `DEEP_PIXEL_SIZE` + (2^-122; the f32 path is exact down to 2^-124 with AA's quarter-pixel + offsets, measured, and the deep path is ~40% slower, so the switch is as + late as that allows). `deep_scale_exp` gives the scale exponent. + `make_uniforms(aspect, height_px)` takes that full-resolution height, the + same during the interaction-downscaled pass so the pipeline doesn't flip. - `src/fractal/kind.rs` — the `FractalKind` enum (Mandelbrot, Burning Ship, Tricorn, Multibrot, Celtic, Perpendicular, Buffalo, Phoenix, Lambda, Complex Multibrot) plus everything that only needs to switch on it: @@ -104,7 +115,12 @@ pixel is a handful of `f32` complex multiplies. precision ≤ `F64_MAX_PRECISION` (80 bits, i.e. shallow views) it takes a plain-`f64` fast path (`compute_reference_f64`), so each kind's formula exists twice in this file (f64 + `FBig`) and both must stay in sync; - `f64_fast_path_matches_big` checks they agree. Requests are made with 1.5× + `f64_fast_path_matches_big` checks they agree. The result is a `RefOrbit`: + `points` plus a parallel `exps`. A point below 2^-100 (only possible on the + `FBig` path) is stored as a normalized mantissa with its exponent in `exps` + (the true value is `points[n]·2^exps[n]`). That happens when the orbit + passes near 0 at a deep minibrot. `has_scaled()` then forces the deep + pipeline, the only one that reads `exps`. Requests are made with 1.5× iteration headroom (`reference_iterations` in `app.rs`), so auto-iterations creeping up during a zoom doesn't recompute the orbit every frame. - `src/shaders/*.wgsl` — none of these are standalone WGSL modules; WGSL has @@ -125,8 +141,8 @@ pixel is a handful of `f32` complex multiplies. pipeline creation. Read those constants in the shader, never `u.kind` / `u.is_julia` / `u.de_coloring` (they're still uploaded for layout reasons). `renderer.rs` builds one pipeline set per `PipelineKey` lazily on first - use, and `tests/shader_valid.rs` compiles every kind × Julia × DE variant to - SPIR-V. So a new kind needs no pipeline-list change, only its `KIND_*` + use, and `tests/shader_valid.rs` compiles every kind × Julia × DE × morph × + deep variant to SPIR-V. So a new kind needs no pipeline-list change, only its `KIND_*` constant. `buddhabrot.wgsl` does the same with its own `override KIND`. Interior pixels exit early through **periodicity detection**. It uses Brent-style checkpoints plus two guards: the cycle's multiplier must be @@ -157,12 +173,41 @@ pixel is a handful of `f32` complex multiplies. detection and the cardioid bypass. App side: `KindMorph` in `app.rs`; the uniforms use the morph the *current reference* was built with (`ref_morph`), not the live one, so orbit and delta formula never disagree. + **Deep views** (`DEEP` override, `u.scale_exp != 0`) handle zooms where + f32 deltas underflow. `make_uniforms` sets `scale_exp = E` (≈ log2 of the + half-height) and uploads `span`/`dc_offset` × 2^-E. The per-pixel `offset` + and `px` are therefore in units of 2^E. `iterate_sample` first runs a + **deep prologue**: + - The delta is carried as `w·2^sx` and the DE derivative as `v·2^sv` + (separate exponents, since they drift apart near the critical point). + - Each step goes through `advance_delta_scaled` → + `deep_step_kind` → `advance_delta_scaled_kind`. These return the step at + its own output scale `t`. Next to the critical point (X tiny or 0), the + linear term vanishes and the step's value is ~e^p, far below 2^sx. + `deep_step_kind` measures X and e in a common unit (the kinds are + p-homogeneous) and the loop moves `sx` there. Assuming the e² terms merely + flush when negligible was wrong exactly there: pixels near deep minibrots + lost their delta and followed the reference forever. + - Rebasing uses X at full range (`ref_fe`). + - Once `|e| > 2^DEEP_EXIT_LOG2` (and dzs is normal), the state converts to + f32 and the ordinary loop continues from the same `n`/`m`. + - Periodicity detection restarts after the prologue with a sentinel save, + because saving the hand-off `z` (an arbitrary phase) made exterior pixels + shadowing a periodic nucleus reference read as interior. + + The deep path is exact at any depth: forcing it everywhere (raise + `DEEP_PIXEL_SIZE`, raise `DEEP_EXIT_LOG2` to about -8) must reproduce the + plain f32 renders on non-chaotic views. That's the check to rerun after + changing it. Known gaps: Lambda's critical point is 1/2, so its step keeps + the input scale. Lambda set mode's reference sits at the origin, so it + never reaches deep zooms anyway. - `src/fractal/renderer.rs` — `FractalRenderer` (wgpu pipelines, uniform + storage buffers, bind groups), `Uniforms` (repr(C) layout that must match the WGSL `Uniforms` struct field-for-field, including padding; it includes CPU-precomputed data: `cm_coef`, the Complex Multibrot binomial - coefficients from `app.rs::complex_binomials`, and `light_count` for the - packed `GpuLight` buffer from `lights.rs::gpu_lights`), and + coefficients from `app.rs::complex_binomials`, `light_count` for the + packed `GpuLight` buffer from `lights.rs::gpu_lights`, and `scale_exp`, + the deep view scale), `ref_exp_buffer` (binding 3, `RefOrbit::exps`), and `FractalCallback` (the `egui_wgpu::CallbackTrait` impl: `prepare()` uploads changed buffers and decides whether to re-run the iterate pass, the cheap colourise pass, or just blit the cached texture). Also `ExportRender`, a @@ -195,7 +240,9 @@ Touches, in order: `kind.rs` (enum variant + `ALL` slot + `label`/ `description`/`formula`/`share_tag`/`from_share_tag`/`default_set_view` arms), `reference.rs` (CPU iteration formula arm, and a test comparing against a naive `f64` iteration), `common.wgsl` (matching `KIND_*` const), -`mandelbrot.wgsl` (matching `advance_delta`/`fprime` arms), `buddhabrot.wgsl` +`mandelbrot.wgsl` (matching `advance_delta`/`fprime` arms, plus the deep +path's `advance_delta_scaled_kind` arm, its degree in `deep_step_kind` and, +if not z²-like, a `deep_fprime` arm), `buddhabrot.wgsl` (matching arm in `advance()`, if the kind makes sense as a Buddhabrot), `renderer.rs` `Uniforms` (only if the kind needs a new per-kind constant, e.g. Phoenix's `phoenix_p`), `app.rs` (`JULIA_PRESETS`/`SET_PRESETS` slot, diff --git a/src/app.rs b/src/app.rs index c22103f..d6e34f7 100644 --- a/src/app.rs +++ b/src/app.rs @@ -12,15 +12,16 @@ use crate::camera::Camera; use crate::cli::Cli; use crate::fractal::{ BuddhabrotCallback, BuddhabrotRenderer, BuddhabrotUniforms, ExportRender, FractalCallback, - FractalKind, FractalRenderer, MAX_REF_POINTS, ShareState, Uniforms, compute_reference, - compute_set_reference, + FractalKind, FractalRenderer, MAX_REF_POINTS, RefOrbit, ShareState, Uniforms, + compute_reference, compute_set_reference, }; 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, ViewState, big_from_decimal_str, big_from_f64, big_to_decimal_str, - interpolate_view, parse_view_spec, precision_for, + Big, DEFAULT_HALF_HEIGHT, MIN_HALF_HEIGHT, 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; @@ -161,7 +162,7 @@ pub(crate) struct RefJob { #[cfg(not(target_arch = "wasm32"))] impl RefJob { /// Iterate the reference orbit at full precision (the expensive part). - pub(crate) fn compute(&self) -> Vec<[f32; 2]> { + pub(crate) fn compute(&self) -> RefOrbit { let key = &self.key; let precision = self.precision; let morph = key.morph.map(|(k, w)| (k, w as f64)); @@ -548,7 +549,7 @@ pub struct FractalApp { fps_window_start: f64, /// Reference orbit (`Z_n` as f32 pairs) for the current view. - reference: Arc>, + reference: Arc, /// Bumped whenever `reference` is replaced, so the GPU re-uploads it. generation: u64, /// Center + zoom the current `reference` was computed at (may differ @@ -708,7 +709,7 @@ impl FractalApp { fps: 0.0, fps_frames: 0, fps_window_start: 0.0, - reference: Arc::new(Vec::new()), + reference: Arc::new(RefOrbit::default()), generation: 0, ref_center_re, ref_center_im, @@ -1149,20 +1150,32 @@ impl FractalApp { self.drift_from(key) > 0.5 * self.view.half_height || !(0.5..=2.0).contains(&ratio) } - /// Complex offset of the live view center from the reference center, in f32. - fn dc_offset(&self) -> [f32; 2] { + /// Complex offset of the live view center from the reference center. + fn dc_offset(&self) -> (f64, f64) { let dre = (&self.view.center_re - &self.ref_center_re) .to_f64() - .value() as f32; + .value(); let dim = (&self.view.center_im - &self.ref_center_im) .to_f64() - .value() as f32; - [dre, dim] + .value(); + (dre, dim) + } + + /// Binary exponent of the deep (rescaled) view scale, or 0 for the plain + /// f32 path (see `Uniforms::scale_exp`). Deep when a pixel of a render + /// `height_px` tall is too small for f32 (`needs_deep`), or when the + /// reference orbit holds points only the deep pipeline can read. + fn scale_exp(&self, height_px: f64) -> i32 { + if needs_deep(self.view.half_height, height_px) || self.reference.has_scaled() { + deep_scale_exp(self.view.half_height) + } else { + 0 + } } fn apply_reference( &mut self, - points: Vec<[f32; 2]>, + points: RefOrbit, cre: Big, cim: Big, hh: f64, @@ -1182,7 +1195,7 @@ impl FractalApp { /// rendering to build its own `ExportRender` without going through /// `egui_wgpu`'s callback machinery. #[cfg(not(target_arch = "wasm32"))] - pub(crate) fn reference_points(&self) -> &[[f32; 2]] { + pub(crate) fn reference_points(&self) -> &RefOrbit { &self.reference } @@ -1327,7 +1340,7 @@ impl FractalApp { /// Install the orbit computed for `job` (from `reference_job`) as the /// current reference, along with the iteration count it was made for. #[cfg(not(target_arch = "wasm32"))] - pub(crate) fn finish_reference(&mut self, job: RefJob, points: Vec<[f32; 2]>) { + pub(crate) fn finish_reference(&mut self, job: RefJob, points: RefOrbit) { self.max_iterations = job.max_iterations; let key = job.key; self.apply_reference( @@ -1383,11 +1396,18 @@ impl FractalApp { .zoom_at_pixel(pos.x as f64, pos.y as f64, height_px, factor); } - pub(crate) fn make_uniforms(&self, aspect: f64) -> Uniforms { + /// `height_px` is the full-resolution render height (it decides whether + /// the view needs the deep pipeline); pass the same value during + /// interaction's downscaled pass, so the pipeline doesn't flip. + pub(crate) fn make_uniforms(&self, aspect: f64, height_px: f64) -> Uniforms { let (span_x, span_y) = self.view.span(aspect); let mode = self.effective_rendering_mode(); + // Deep views upload the geometry pre-multiplied by 2^-E (exact). + let scale_exp = self.scale_exp(height_px); + let inv_scale = 2f64.powi(-scale_exp); + let (dc_re, dc_im) = self.dc_offset(); Uniforms { - span: [span_x as f32, span_y as f32], + span: [(span_x * inv_scale) as f32, (span_y * inv_scale) as f32], max_iter: self.max_iterations.min(MAX_REF_POINTS as u32 - 1), ref_len: self.reference.len() as u32, color_offset: self.color_offset, @@ -1400,7 +1420,7 @@ impl FractalApp { kind: self.ref_kind.unwrap_or(self.kind) as u32, power: self.power, morph_from: self.ref_morph.map_or(0, |(k, _)| k as u32), - dc_offset: self.dc_offset(), + dc_offset: [(dc_re * inv_scale) as f32, (dc_im * inv_scale) 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], @@ -1416,7 +1436,7 @@ impl FractalApp { screen_dim: self.screen_dim, light_count: gpu_lights(&self.lights).1, cm_coef: complex_binomials(self.complex_power), - _pad3: [0; _], + scale_exp, } } @@ -1476,7 +1496,7 @@ impl FractalApp { let scale = self.export_scale.max(1.0); let w = ((self.last_size_px.x * scale).round() as u32).clamp(16, MAX_EXPORT_DIM); let h = ((self.last_size_px.y * scale).round() as u32).clamp(16, MAX_EXPORT_DIM); - let uniforms = self.make_uniforms(w as f64 / h as f64); + let uniforms = self.make_uniforms(w as f64 / h as f64, h as f64); let device = rs.device.clone(); let queue = rs.queue.clone(); @@ -1507,14 +1527,7 @@ impl FractalApp { .unwrap_or_else(|| format!("fractal-{}.png", unix_timestamp())); std::thread::spawn(move || { let er = ExportRender::new( - &device, - &queue, - &handles, - w, - h, - uniforms, - reference.as_slice(), - &lights, + &device, &queue, &handles, w, h, uniforms, &reference, &lights, ); let sh = Arc::clone(&shared); let png = @@ -1535,14 +1548,7 @@ impl FractalApp { const RENDER_END: f32 = 0.6; wasm_bindgen_futures::spawn_local(async move { let er = ExportRender::new( - &device, - &queue, - &handles, - w, - h, - uniforms, - reference.as_slice(), - &lights, + &device, &queue, &handles, w, h, uniforms, &reference, &lights, ); // Render tile by tile, awaiting each submission so the browser @@ -2451,7 +2457,7 @@ impl FractalApp { && hh > 0.0 && hh.is_finite() { - self.view.half_height = hh; + self.view.half_height = hh.max(MIN_HALF_HEIGHT); self.view.sync_precision(); } self.zoom_edited = false; @@ -2862,7 +2868,12 @@ impl FractalApp { self.screen_dim = [rect.width(), rect.height()]; self.camera.set_aspect_ratio(aspect as f32); - let mut uniforms = self.make_uniforms(aspect); + // Full-resolution height (3D included), not the interaction-downscaled one. + let mut full_height = (rect.height() * ppp).round() as f64; + if self.effective_rendering_mode() == 2 { + full_height *= self.render_scale_3d as f64; + } + let mut uniforms = self.make_uniforms(aspect, full_height); // Supersampling is wasted on the low-res pass, and on a kind-switch // morph (every frame re-iterates, and the blend moves on next frame). // `ref_morph` too: the last morphed reference outlives `morph` by a diff --git a/src/fractal/mod.rs b/src/fractal/mod.rs index 8037fd3..f046e58 100644 --- a/src/fractal/mod.rs +++ b/src/fractal/mod.rs @@ -9,7 +9,7 @@ pub mod share; pub use buddhabrot::{BuddhabrotCallback, BuddhabrotRenderer, BuddhabrotUniforms}; pub use kind::FractalKind; -pub use reference::{compute_reference, compute_set_reference}; +pub use reference::{RefOrbit, compute_reference, compute_set_reference}; #[cfg(not(target_arch = "wasm32"))] pub use renderer::PipelineKey; #[cfg(target_arch = "wasm32")] diff --git a/src/fractal/reference.rs b/src/fractal/reference.rs index 457e984..9557b30 100644 --- a/src/fractal/reference.rs +++ b/src/fractal/reference.rs @@ -39,6 +39,86 @@ const REFERENCE_ESCAPE_SQ: f64 = 1.0e10; /// below 0.1% of a pixel. const F64_MAX_PRECISION: usize = 80; +/// Orbit points with a magnitude below `2^TINY_LOG2` are stored normalized +/// (mantissa + exponent, see [`RefOrbit::exps`]): f32's smallest normal is +/// ~2^-126, and the GPU's deep (rescaled) phase needs these points' exact +/// value to decide rebasing. The margin keeps a few mantissa bits clear of +/// the subnormal range for the smaller component. +const TINY_LOG2: i32 = -100; + +/// A reference orbit as uploaded to the GPU. +#[derive(Clone, Debug, Default, PartialEq)] +pub struct RefOrbit { + /// `Z_n` as f32 pairs. For points with a non-zero `exps[n]`, a mantissa + /// instead: the true value is `points[n] * 2^exps[n]`. + pub points: Vec<[f32; 2]>, + /// Per-point binary exponent (same length as `points`). Non-zero only for + /// points too small for f32's exponent range (see [`TINY_LOG2`]); only + /// the deep shader pipeline reads it, so [`Self::has_scaled`] forces it. + pub exps: Vec, +} + +impl RefOrbit { + fn with_capacity(n: usize) -> Self { + Self { + points: Vec::with_capacity(n), + exps: Vec::with_capacity(n), + } + } + + fn push(&mut self, point: [f32; 2]) { + self.points.push(point); + self.exps.push(0); + } + + /// Whether any point is stored as mantissa + exponent, i.e. the orbit + /// can only be read by the deep pipeline. + pub fn has_scaled(&self) -> bool { + self.exps.iter().any(|&e| e != 0) + } +} + +impl core::ops::Deref for RefOrbit { + type Target = [[f32; 2]]; + fn deref(&self) -> &Self::Target { + &self.points + } +} + +impl<'a> IntoIterator for &'a RefOrbit { + type Item = &'a [f32; 2]; + type IntoIter = core::slice::Iter<'a, [f32; 2]>; + fn into_iter(self) -> Self::IntoIter { + self.points.iter() + } +} + +/// `floor(log2|x|)`, or `None` for zero. Exact (from the binary +/// representation), and works far below f64's range. +fn big_log2_floor(x: &Big) -> Option { + let repr = x.repr(); + let digits = repr.digits(); + (digits > 0).then(|| repr.exponent() + digits as isize - 1) +} + +/// Store `(zr, zi)` into `orbit`, as plain f32 unless its magnitude is below +/// `2^TINY_LOG2`, in which case both components share an exponent `k` and +/// the stored mantissa `Z * 2^-k` has its larger component in `[0.5, 1)`. +fn push_big_point(orbit: &mut RefOrbit, zr: &Big, zi: &Big) { + let (lr, li) = (big_log2_floor(zr), big_log2_floor(zi)); + let top = lr.max(li); + match top { + Some(top) if top < TINY_LOG2 as isize => { + let k = top + 1; + let mr = (zr.clone() << -k).to_f64().value() as f32; + let mi = (zi.clone() << -k).to_f64().value() as f32; + orbit.points.push([mr, mi]); + orbit.exps.push(k as i32); + } + _ => orbit.push([zr.to_f64().value() as f32, zi.to_f64().value() as f32]), + } +} + /// Compute the reference orbit `Z_0..Z_{len-1}` where `Z_0 = z0` and /// `Z_{n+1} = f(Z_n, c)` for the given `kind` (and `power`, for Multibrot), up /// to `max_iter` steps at `precision` bits. Each entry is `[re, im]` in f32. @@ -59,7 +139,7 @@ pub fn compute_reference( lambda_l: (f64, f64), complex_power: (f64, f64), morph: Option<(FractalKind, f64)>, -) -> Vec<[f32; 2]> { +) -> RefOrbit { // A zero-weight morph is just the plain kind; skip the second formula. let morph = morph.filter(|&(_, w)| w != 0.0); if precision <= F64_MAX_PRECISION { @@ -110,12 +190,12 @@ fn compute_reference_f64( kind: FractalKind, k: &StepConstsF64, morph: Option<(FractalKind, f64)>, -) -> Vec<[f32; 2]> { +) -> RefOrbit { let (mut zr, mut zi) = z0; // Previous iterate, for the Phoenix two-term recurrence (Y_{-1} = 0). let mut prev = (0.0f64, 0.0f64); - let mut points: Vec<[f32; 2]> = Vec::with_capacity(max_iter as usize + 1); + let mut points = RefOrbit::with_capacity(max_iter as usize + 1); for _ in 0..=max_iter { points.push([zr as f32, zi as f32]); if zr * zr + zi * zi > REFERENCE_ESCAPE_SQ { @@ -217,7 +297,7 @@ fn compute_reference_big( kind: FractalKind, k: &StepConsts, morph: Option<(FractalKind, f64)>, -) -> Vec<[f32; 2]> { +) -> RefOrbit { let precision = k.precision; let morph = morph.map(|(from, w)| (from, big_from_f64(w, precision))); @@ -227,13 +307,12 @@ fn compute_reference_big( let mut zr_prev = big_zero(precision); let mut zi_prev = big_zero(precision); - let mut points: Vec<[f32; 2]> = Vec::with_capacity(max_iter as usize + 1); + let mut points = RefOrbit::with_capacity(max_iter as usize + 1); for _ in 0..=max_iter { - let fr = zr.to_f64().value() as f32; - let fi = zi.to_f64().value() as f32; - points.push([fr, fi]); + push_big_point(&mut points, &zr, &zi); + let [fr, fi] = *points.points.last().unwrap(); let mag = (fr as f64) * (fr as f64) + (fi as f64) * (fi as f64); if mag > REFERENCE_ESCAPE_SQ { break; @@ -403,7 +482,7 @@ pub fn compute_set_reference( lambda_l: (f64, f64), complex_power: (f64, f64), morph: Option<(FractalKind, f64)>, -) -> Vec<[f32; 2]> { +) -> RefOrbit { let zero = big_zero(precision); compute_reference( &zero, @@ -510,6 +589,52 @@ mod tests { } } + /// Orbit points below f32's range are stored as a normalized mantissa + /// plus exponent (for the deep GPU phase); every other point stays a + /// plain f32 with exponent 0. + #[test] + fn tiny_points_are_stored_normalized() { + let bits = 400; + // c = -1 + δ: X_2 = c(c + 1) = -δ + δ², far below f32's range. + let delta = 1e-45_f64; + let cr = big_from_f64(-1.0, bits) + big_from_f64(delta, bits); + let ci = big_from_f64(0.0, bits); + let orbit = compute_set_reference( + &cr, + &ci, + 3, + bits, + FractalKind::Mandelbrot, + 2, + (0.0, 0.0), + (0.0, 0.0), + (0.0, 0.0), + None, + ); + assert_eq!(orbit.exps.len(), orbit.points.len()); + assert!(orbit.has_scaled()); + assert_eq!(&orbit.exps[..2], &[0, 0], "X_0 = 0 and X_1 = c are plain"); + let [mr, mi] = orbit.points[2]; + let k = orbit.exps[2]; + assert!(k < TINY_LOG2, "exponent {k}"); + assert!( + (0.5..1.0).contains(&mr.abs()), + "mantissa {mr} not normalized" + ); + assert_eq!(mi, 0.0); + let x2 = mr as f64 * 2f64.powi(k); + assert!( + (x2 + delta).abs() < 1e-6 * delta, + "X_2 = {x2}, expected {}", + -delta + ); + + // A shallow orbit stays entirely plain. + let plain = set_ref(-0.75, 0.1, FractalKind::Mandelbrot, None); + assert!(!plain.has_scaled()); + assert!(plain.exps.iter().all(|&e| e == 0)); + } + /// A point inside the main cardioid never escapes: full-length orbit. #[test] fn interior_orbit_runs_full_length() { @@ -842,12 +967,7 @@ mod tests { } } - fn set_ref( - cr: f64, - ci: f64, - kind: FractalKind, - morph: Option<(FractalKind, f64)>, - ) -> Vec<[f32; 2]> { + fn set_ref(cr: f64, ci: f64, kind: FractalKind, morph: Option<(FractalKind, f64)>) -> RefOrbit { compute_set_reference( &Big::try_from(cr).unwrap(), &Big::try_from(ci).unwrap(), diff --git a/src/fractal/renderer.rs b/src/fractal/renderer.rs index 0ac355a..f05c0d5 100644 --- a/src/fractal/renderer.rs +++ b/src/fractal/renderer.rs @@ -14,6 +14,7 @@ use std::sync::Arc; use eframe::egui_wgpu::{self, wgpu}; +use super::reference::RefOrbit; use crate::lights::{GpuLight, Light, MAX_LIGHT_COUNT, gpu_lights}; /// Maximum reference-orbit length (points) the storage buffer can hold. Also @@ -51,6 +52,7 @@ fn geom_differs(a: &Uniforms, b: &Uniforms) -> bool { || a.power != b.power || a.complex_power != b.complex_power || a.dc_offset != b.dc_offset + || a.scale_exp != b.scale_exp || a.phoenix_p != b.phoenix_p || a.lambda_l != b.lambda_l || a.morph_from != b.morph_from @@ -87,6 +89,9 @@ pub struct PipelineKey { de: bool, /// A kind-switch morph is in progress (`morph_w > 0`). morph: bool, + /// Deep view: the delta starts out in rescaled (mantissa + exponent) + /// form (`scale_exp != 0`). + deep: bool, } impl PipelineKey { @@ -96,15 +101,17 @@ impl PipelineKey { julia: u.is_julia != 0, de: u.de_coloring != 0, morph: u.morph_w > 0.0, + deep: u.scale_exp != 0, } } - fn constants(&self) -> [(&'static str, f64); 4] { + fn constants(&self) -> [(&'static str, f64); 5] { [ ("KIND", self.kind as f64), ("IS_JULIA", self.julia as u32 as f64), ("DE", self.de as u32 as f64), ("MORPH", self.morph as u32 as f64), + ("DEEP", self.deep as u32 as f64), ] } } @@ -217,7 +224,11 @@ pub struct Uniforms { /// 0 = no morph (and the iteration pipeline is then specialized without /// the morph path, see [`PipelineKey`]). pub morph_w: f32, - pub _pad3: [u32; 1], + /// Binary exponent `E` of the deep (rescaled) view scale: `span` and + /// `dc_offset` are uploaded multiplied by `2^-E`, so they stay inside + /// f32's exponent range at any depth. Non-zero exactly when the deep + /// pipeline is used (see [`PipelineKey`] and `mandelbrot.wgsl`'s `DEEP`). + pub scale_exp: i32, /// Complex binomial coefficients `C(complex_power, k)`, k = 1..16, two per /// row (odd k in `[0..2]`, even k in `[2..4]`), for the Complex Multibrot /// delta series. Derived from `complex_power` alone. @@ -279,6 +290,7 @@ pub struct FractalRenderer { bind_group_layout: wgpu::BindGroupLayout, uniform_buffer: wgpu::Buffer, ref_buffer: wgpu::Buffer, + ref_exp_buffer: wgpu::Buffer, lights_buffer: wgpu::Buffer, bind_group: wgpu::BindGroup, target_format: wgpu::TextureFormat, @@ -331,6 +343,15 @@ impl FractalRenderer { mapped_at_creation: false, }); + // Per-point exponents of the reference orbit (`RefOrbit::exps`), only + // read by deep pipelines. + let ref_exp_buffer = device.create_buffer(&wgpu::BufferDescriptor { + label: Some("reference orbit exponents"), + size: (MAX_REF_POINTS * std::mem::size_of::()) as u64, + usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_DST, + mapped_at_creation: false, + }); + let lights_buffer = device.create_buffer(&wgpu::BufferDescriptor { label: Some("lights parameters"), size: std::mem::size_of::<[GpuLight; MAX_LIGHT_COUNT]>() as u64, @@ -374,6 +395,16 @@ impl FractalRenderer { }, count: None, }, + wgpu::BindGroupLayoutEntry { + binding: 3, + 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, + }, ], }); @@ -393,6 +424,10 @@ impl FractalRenderer { binding: 2, resource: lights_buffer.as_entire_binding(), }, + wgpu::BindGroupEntry { + binding: 3, + resource: ref_exp_buffer.as_entire_binding(), + }, ], }); @@ -590,6 +625,7 @@ impl FractalRenderer { bind_group_layout, uniform_buffer, ref_buffer, + ref_exp_buffer, lights_buffer, bind_group, target_format, @@ -892,7 +928,7 @@ impl ExportRender { width: u32, height: u32, uniforms: Uniforms, - reference: &[[f32; 2]], + reference: &RefOrbit, lights: &[Light], ) -> Self { let uniform_buffer = device.create_buffer(&wgpu::BufferDescriptor { @@ -910,8 +946,19 @@ impl ExportRender { usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_DST, mapped_at_creation: false, }); + let ref_exp_buffer = device.create_buffer(&wgpu::BufferDescriptor { + label: Some("export reference orbit exponents"), + size: (count.max(1) * std::mem::size_of::()) as u64, + usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_DST, + mapped_at_creation: false, + }); if count > 0 { queue.write_buffer(&ref_buffer, 0, bytemuck::cast_slice(&reference[..count])); + queue.write_buffer( + &ref_exp_buffer, + 0, + bytemuck::cast_slice(&reference.exps[..count]), + ); } // Only read by the shadow branch's custom-lights palette; harmless @@ -942,6 +989,10 @@ impl ExportRender { binding: 2, resource: lights_buffer.as_entire_binding(), }, + wgpu::BindGroupEntry { + binding: 3, + resource: ref_exp_buffer.as_entire_binding(), + }, ], }); @@ -1391,7 +1442,7 @@ pub struct FractalCallback { /// Lights buffer contents, from [`gpu_lights`] (its count is in /// `uniforms.light_count`). pub lights: [GpuLight; MAX_LIGHT_COUNT], - pub reference: Arc>, + pub reference: Arc, pub generation: u64, /// Widget size in physical pixels — the cache texture resolution. pub size_px: [u32; 2], @@ -1457,6 +1508,11 @@ impl egui_wgpu::CallbackTrait for FractalCallback { 0, bytemuck::cast_slice(&self.reference[..count]), ); + queue.write_buffer( + &renderer.ref_exp_buffer, + 0, + bytemuck::cast_slice(&self.reference.exps[..count]), + ); renderer.uploaded_generation = self.generation; } diff --git a/src/headless.rs b/src/headless.rs index fd12d88..dd8a595 100644 --- a/src/headless.rs +++ b/src/headless.rs @@ -54,7 +54,7 @@ pub fn run(cli: Cli) -> Result<(), String> { let (device, queue) = pollster::block_on(request_device())?; let format = wgpu::TextureFormat::Bgra8Unorm; let renderer = FractalRenderer::new(&device, format); - let uniforms = app.make_uniforms(width as f64 / height as f64); + let uniforms = app.make_uniforms(width as f64 / height as f64, height as f64); let handles = renderer.export_handles(&device, &uniforms); let er = ExportRender::new( @@ -263,7 +263,7 @@ fn run_animation( let (png_tx, png_rx) = mpsc::sync_channel::<(usize, Vec, u32, bool)>(threads * 2); let png_rx = Mutex::new(png_rx); std::thread::scope(|scope| { - let (ref_tx, ref_rx) = mpsc::sync_channel::<(usize, Vec<[f32; 2]>)>(threads * 2); + let (ref_tx, ref_rx) = mpsc::sync_channel::<(usize, crate::fractal::RefOrbit)>(threads * 2); for _ in 0..threads { let ref_tx = ref_tx.clone(); let (jobs, next_job, failed) = (&jobs, &next_job, &failed); @@ -317,7 +317,7 @@ fn run_animation( apply_frame(&mut app, i as u32); app.finish_reference(jobs[i].clone(), points); - let uniforms = app.make_uniforms(aspect); + let uniforms = app.make_uniforms(aspect, height as f64); let handles = pipelines .entry(PipelineKey::from_uniforms(&uniforms)) .or_insert_with(|| renderer.export_handles(&device, &uniforms)); diff --git a/src/shaders/iterate_uniforms.wgsl b/src/shaders/iterate_uniforms.wgsl index 2accfd7..54d7047 100644 --- a/src/shaders/iterate_uniforms.wgsl +++ b/src/shaders/iterate_uniforms.wgsl @@ -48,6 +48,10 @@ struct Uniforms { // Kind-switch morph weight: each step is (1 - w)*f_kind + w*f_morph_from; // 0 = no morph. Only read by MORPH pipelines (see mandelbrot.wgsl). morph_w: f32, + // Deep views: binary exponent E of the view scale. `span` and `dc_offset` + // are uploaded multiplied by 2^-E so they stay in f32's range; 0 = not + // deep (plain f32 values). Only read by DEEP pipelines (mandelbrot.wgsl). + scale_exp: i32, // Complex binomial coefficients C(complex_power, k) for k = 1..16, two per // vec4 (k odd in .xy, k even in .zw), for the Complex Multibrot delta // series. Precomputed on the CPU since they only depend on the power. diff --git a/src/shaders/mandelbrot.wgsl b/src/shaders/mandelbrot.wgsl index 8e99528..3410b27 100644 --- a/src/shaders/mandelbrot.wgsl +++ b/src/shaders/mandelbrot.wgsl @@ -11,6 +11,11 @@ // true value |y| drops below the delta |e|, or the reference runs out, we reset // the reference index to 0 and carry the full value as the new delta (valid // because X_0 = 0). +// +// Past ~1e30 zoom the deltas themselves leave f32's exponent range (smallest +// normal ~1.2e-38), so `DEEP` pipelines start each pixel in a rescaled form, +// e = w * 2^s with an f32 mantissa `w` and an i32 exponent `s`, and hand over +// to the plain f32 loop once the delta is big enough (see `iterate_sample`). @group(0) @binding(0) var u: Uniforms; @group(0) @binding(1) var ref_orbit: array>; @@ -20,6 +25,11 @@ // Only read by the adaptive-AA refine pass (`fs_refine`): the 1-sample-per- // pixel data texture written by `fs_data`, which decides where to supersample. @group(1) @binding(0) var coarse_tex: texture_2d; +// Per-point binary exponents of the reference orbit (`RefOrbit::exps`): the +// true X[m] is ref_orbit[m] * 2^ref_exp[m]. Non-zero only for points below +// f32's range, which only deep references contain; only `DEEP` pipelines read +// it (see `ref_at`). +@group(0) @binding(3) var ref_exp: array; // Pipeline-overridable specialization constants, set per pipeline from the // uniforms' `kind` / `is_julia` / `de_coloring` (see `PipelineKey` in @@ -35,6 +45,10 @@ override DE: bool = false; // second kind, `u.morph_from`. That one is a runtime value (it only lives for // the length of the animation), so only MORPH pipelines pay for its branches. override MORPH: bool = false; +// Deep view (`u.scale_exp != 0`): the per-pixel offset `dc`, the pixel size +// and `u.span` / `u.dc_offset` are all in units of 2^scale_exp, and each pixel +// starts in the rescaled deep phase (see `iterate_sample`). +override DEEP: bool = false; struct VsOut { @builtin(position) pos: vec4, @@ -82,14 +96,18 @@ fn diffabs(c: f32, d: f32) -> f32 { // tiny, but that only perturbs `s` by a relative f32 epsilon, and the result // is `e * s`, so the delta keeps full relative precision. fn multibrot_delta(z: vec2, e: vec2, p: u32) -> vec2 { - let y = z + e; + return cmul(e, multibrot_sum(z, z + e, p)); +} + +// The sum in `multibrot_delta`, sum_{k=0}^{p-1} y^k Z^{p-1-k} with y = Z+e. +fn multibrot_sum(z: vec2, y: vec2, p: u32) -> vec2 { var s = vec2(1.0, 0.0); var zj = vec2(1.0, 0.0); for (var j: u32 = 1u; j < p; j = j + 1u) { zj = cmul(zj, z); // Z^j s = cmul(s, y) + zj; } - return cmul(e, s); + return s; } // Maximum number of terms in `complex_multibrot_delta`'s series (matches the @@ -242,6 +260,324 @@ fn fprime(z: vec2) -> vec2 { return d; } +// ---- Deep (rescaled) phase helpers ------------------------------------- +// +// A deep value is an f32 mantissa times 2^exponent. All scaling is by exact +// powers of two, so it never rounds. + +const LN2: f32 = 0.6931471805599453; +// Where `ldexp_sat` saturates. Only ever compared against small values +// (see `diffabs_scaled`), and far enough from f32's max that doubling it, +// or squaring a value of the size it's compared with, stays finite. +const LDEXP_SAT: f32 = 1.2676506e30; // 2^100 + +// x * 2^k for any k. WGSL's `ldexp` is only defined for exponents inside +// f32's range, so go through `frexp`: results below the smallest normal +// flush to 0, results above 2^100 saturate to ±2^100. +fn ldexp_sat(x: f32, k: i32) -> f32 { + if x == 0.0 { + return 0.0; + } + let f = frexp(x); + let ex = f.exp + k; + if ex > 100 { + return select(-LDEXP_SAT, LDEXP_SAT, x > 0.0); + } + if ex < -125 { + return 0.0; + } + return ldexp(f.fract, ex); +} + +fn ldexp2_sat(v: vec2, k: i32) -> vec2 { + return vec2(ldexp_sat(v.x, k), ldexp_sat(v.y, k)); +} + +// A complex number as mantissa * 2^e, with the mantissa's larger component in +// [0.5, 1) (or exactly 0, with e = 0). +struct Fe { + m: vec2, + e: i32, +}; + +fn fe_make(v: vec2, e: i32) -> Fe { + let a = max(abs(v.x), abs(v.y)); + if a == 0.0 { + return Fe(vec2(0.0, 0.0), 0); + } + let k = frexp(a).exp; + return Fe(ldexp2_sat(v, -k), e + k); +} + +// Reference point X[m] as f32 (points stored normalized flush to 0 here; +// only deep references have any, see `ref_exp`). +fn ref_at(m: u32) -> vec2 { + let x = ref_orbit[m]; + if DEEP { + let k = ref_exp[m]; + if k != 0 { + return ldexp2_sat(x, k); + } + } + return x; +} + +// X[m] with its full exponent range (deep phase only). +fn ref_fe(m: u32) -> Fe { + return fe_make(ref_orbit[m], ref_exp[m]); +} + +// Complex log of m * 2^e (m != 0). +fn clog_fe(m: vec2, e: i32) -> vec2 { + return vec2(0.5 * log(dot(m, m)) + f32(e) * LN2, atan2(m.y, m.x)); +} + +fn cexp(a: vec2) -> vec2 { + return exp(a.x) * vec2(cos(a.y), sin(a.y)); +} + +// diffabs(c, 2^s * d) / 2^s = diffabs(c / 2^s, d): diffabs is positively +// homogeneous. Saturating c / 2^s is harmless: once it dwarfs |d| the result +// is just ±d. +fn diffabs_scaled(c: f32, d: f32, s: i32) -> f32 { + return diffabs(ldexp_sat(c, -s), d); +} + +// Sum of two `Fe`s, at the larger one's exponent. +fn fe_add(a: Fe, b: Fe) -> Fe { + if a.m.x == 0.0 && a.m.y == 0.0 { + return b; + } + if b.m.x == 0.0 && b.m.y == 0.0 { + return a; + } + let e = max(a.e, b.e); + return fe_make(ldexp2_sat(a.m, a.e - e) + ldexp2_sat(b.m, b.e - e), e); +} + +// One deep step's delta, (f(X + e) - f(X)) / 2^t: the mantissa `w` and its +// exponent `t`. `t` is the input scale s except next to the critical point, +// where the linear part of the step vanishes and the result is ~e^2, far +// below the input's scale. +struct DeepStep { + w: vec2, + t: i32, +}; + +// The per-kind formula of a deep step: (f(X + e) - f(X)) / 2^t for the +// delta e = w * 2^s, given X measured in units of 2^u (`x` = X / 2^u) and +// `sc` = 2^se, se = s - u, the delta's scale in those units. The result's +// scale is t = s + (p-1)·u for a degree-p kind, since every kind here but +// Lambda is p-homogeneous in (X, e) jointly (`diffabs_scaled` rescales the +// fold-point comparisons the same way). Usually u = 0 (x = X, sc = 2^s, +// t = s); see `deep_step_kind` for when it isn't. Lambda always gets u = 0. +// Every kind needs an arm here too. +fn advance_delta_scaled_kind(kind: u32, x: vec2, w: vec2, sc: f32, se: i32) -> vec2 { + if kind == KIND_BURNING_SHIP { + let base = 2.0 * cmul(x, w) + sc * cmul(w, w); + let dp = x.x * w.y + x.y * w.x + sc * w.x * w.y; + return vec2(base.x, 2.0 * diffabs_scaled(x.x * x.y, dp, se)); + } else if kind == KIND_TRICORN { + let cx = conj(x); + 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))); + } 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); + } else if kind == KIND_BUFFALO { + 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), + -diffabs_scaled(2.0 * x.x * x.y, sq.y, se)); + } else if kind == KIND_PERPENDICULAR { + let sq = 2.0 * cmul(x, w) + sc * cmul(w, w); + let da = diffabs_scaled(x.y, w.y, se); // (|Y + ey| - |Y|) / 2^s + let abs_yf = abs(x.y) + sc * da; // |Y + ey| * 2^(s-t) + return vec2(sq.x, -2.0 * (x.x * da + w.x * abs_yf)); + } else if kind == KIND_LAMBDA { + let t = vec2(1.0 - 2.0 * x.x - sc * w.x, -2.0 * x.y - sc * w.y); + return cmul(u.lambda_l, cmul(w, t)); + } + return 2.0 * cmul(x, w) + sc * cmul(w, w); // Mandelbrot (and Phoenix square part) +} + +// Below this exponent a reference point counts as next to the critical point +// 0 (see `deep_step_kind`). Above it, the e^2 terms 2^s * w^2 can only flush +// to 0 when they are below 2^-50 of the linear ones. +const DEEP_X_NEAR_LOG2: i32 = -60; + +// One deep step of `kind` for the delta e = w * 2^s from X (`x` as f32, `xf` +// at full range). Usually the input scale is kept (t = s). But when X is +// tiny (next to the critical point, e.g. at a minibrot's period), the linear +// term vanishes and the step's value is ~e^p, which would flush to 0 at +// scale s. There every z^p-like kind is p-homogeneous in (X, e) jointly, so +// both are measured in units of 2^u (u = the larger one's exponent) and the +// result lands at t = s + (p-1)·u (`ue` below). +fn deep_step_kind(kind: u32, x: vec2, xf: Fe, w: vec2, sc: f32, s: i32) -> DeepStep { + if kind == KIND_COMPLEX_MULTIBROT { + return complex_multibrot_step(xf, w, s); + } + let x_zero = xf.m.x == 0.0 && xf.m.y == 0.0; + // Lambda's critical point is 1/2 and its step has a constant linear + // term (λ·e), so it never needs this. + if kind == KIND_LAMBDA || (!x_zero && xf.e >= DEEP_X_NEAR_LOG2) { + return DeepStep(advance_delta_scaled_kind(kind, x, w, sc, s), s); + } + let kw = deep_log2(w, vec2(0.0, 0.0)); + if kw == DEEP_ZERO { + return DeepStep(w, s); // e = 0: f(X) - f(X) + } + var ue = s + kw; + if !x_zero { + ue = max(ue, xf.e); + } + var deg = 2; + if kind == KIND_MULTIBROT { + deg = i32(clamp(u.power, 2u, 8u)); + } + let xk = ldexp2_sat(xf.m, xf.e - ue); + let se = s - ue; + let dw = advance_delta_scaled_kind(kind, xk, w, ldexp_sat(1.0, se), se); + return DeepStep(dw, s + (deg - 1) * ue); +} + +// Scaled `complex_multibrot_delta` with X as a full-range `Fe`. Same series +// as the f32 version, rewritten as X^(p-1) * w * sum_k C(p,k) r^(k-1) +// (r = e/X) so nothing is formed at the delta's true scale; X^(p-1)'s own +// exponent goes into the result's `t`. When |e/X| >= 0.5, X is itself tiny +// (|X| <= 2|e|), so both terms of the direct form are taken in log space at +// the larger one's scale. Branch-cut crossings leave the deep phase before +// stepping (`deep_cut_crossing`). +fn complex_multibrot_step(xf: Fe, w: vec2, s: i32) -> DeepStep { + let p = u.complex_power; + let wf = fe_make(w, s); + if wf.m.x == 0.0 && wf.m.y == 0.0 { + return DeepStep(vec2(0.0, 0.0), s); + } + if xf.m.x == 0.0 && xf.m.y == 0.0 { + // e^p. + let l = cmul(p, clog_fe(wf.m, wf.e)); + let k = i32(floor(l.x / LN2)); + return DeepStep(cexp(l - vec2(f32(k) * LN2, 0.0)), k); + } + let r = ldexp2_sat(cdiv(wf.m, xf.m), wf.e - xf.e); + if dot(r, r) < 0.25 { + var acc = cm_coef(1u); + var rk = r; // r^(k-1) + for (var k: u32 = 2u; k <= COMPLEX_MULTIBROT_TERMS; k = k + 1u) { + acc = acc + cmul(cm_coef(k), rk); + rk = cmul(rk, r); + if dot(rk, rk) < 1e-18 * dot(acc, acc) { + break; + } + } + // X^(p-1) = cexp(l) = cexp(l - k·ln2) * 2^k. + let l = cmul(p - vec2(1.0, 0.0), clog_fe(xf.m, xf.e)); + let k = i32(floor(l.x / LN2)); + let x_pm1 = cexp(l - vec2(f32(k) * LN2, 0.0)); + return DeepStep(cmul(cmul(w, x_pm1), acc), s + k); + } + // (X + e)^p - X^p, both in log space (X + e may be exactly 0). + let xs = ldexp2_sat(xf.m, xf.e - s); + let yf = fe_make(xs + w, s); + let lb = cmul(p, clog_fe(xf.m, xf.e)); + var k = i32(floor(lb.x / LN2)); + var la = vec2(0.0, 0.0); + let y_zero = yf.m.x == 0.0 && yf.m.y == 0.0; + if !y_zero { + la = cmul(p, clog_fe(yf.m, yf.e)); + k = max(k, i32(floor(la.x / LN2))); + } + let kl = vec2(f32(k) * LN2, 0.0); + var ya = vec2(0.0, 0.0); + if !y_zero { + ya = cexp(la - kl); + } + return DeepStep(ya - cexp(lb - kl), k); +} + +// Deep-phase twin of `advance_delta` (same morph blend, at the larger of the +// two kinds' output scales). +fn advance_delta_scaled(x: vec2, xf: Fe, w: vec2, sc: f32, s: i32) -> DeepStep { + let a = deep_step_kind(KIND, x, xf, w, sc, s); + if MORPH { + let b = deep_step_kind(u.morph_from, x, xf, w, sc, s); + let t = max(a.t, b.t); + return DeepStep(mix(ldexp2_sat(a.w, a.t - t), ldexp2_sat(b.w, b.t - t), u.morph_w), t); + } + return a; +} + +// Whether this step would take Complex Multibrot's X + e across the branch +// cut (see `complex_multibrot_delta`). The delta then jumps to the size of X, +// so the deep phase ends and the f32 loop takes the step. +fn deep_cut_crossing(xf: Fe, w: vec2, s: i32) -> bool { + let cm = KIND == KIND_COMPLEX_MULTIBROT || (MORPH && u.morph_from == KIND_COMPLEX_MULTIBROT); + if !cm { + return false; + } + let yn = xf.m + ldexp2_sat(w, s - xf.e); // (X + e) / 2^xe + return xf.m.x < 0.0 && ((xf.m.y < 0.0) != (yn.y < 0.0)); +} + +// Binary exponent of the largest component of a pair of complex mantissas +// (`frexp` convention: |x| < 2^k), or DEEP_ZERO when both are exactly 0. +const DEEP_ZERO: i32 = -100000; +fn deep_log2(a: vec2, b: vec2) -> i32 { + let m = max(max(abs(a.x), abs(a.y)), max(abs(b.x), abs(b.y))); + if m == 0.0 { + return DEEP_ZERO; + } + return frexp(m).exp; +} + +// f'(y) at the full value y = X + w * 2^s, as an `Fe`. Plain `fprime` in f32 +// unless y is below f32's comfortable range, which only happens next to the +// critical point 0 (after a rebase), where X is itself tiny: then y is +// formed in the 2^s-scaled domain and f' ~ p·y^(p-1) keeps its exponent. +fn deep_fprime(xf: Fe, w: vec2, s: i32, yt: vec2) -> Fe { + if MORPH || max(abs(yt.x), abs(yt.y)) >= DEEP_TINY { + return Fe(fprime(yt), 0); + } + let yf = fe_make(ldexp2_sat(xf.m, xf.e - s) + w, s); + if yf.m.x == 0.0 && yf.m.y == 0.0 { + return Fe(fprime(vec2(0.0, 0.0)), 0); + } + if KIND == KIND_MULTIBROT { + let p = clamp(u.power, 2u, 8u); + var ym = yf.m; // m^(p-1) + for (var k: u32 = 2u; k < p; k = k + 1u) { + ym = cmul(ym, yf.m); + } + return fe_make(f32(p) * ym, yf.e * i32(p - 1u)); + } else if KIND == KIND_LAMBDA { + return Fe(fprime(vec2(0.0, 0.0)), 0); // λ(1 - 2y) ~ λ + } else if KIND == KIND_COMPLEX_MULTIBROT { + // p·y^(p-1) = p·exp(l), l = (p-1)·ln y; keep exp(Re l)'s exponent. + let p = u.complex_power; + let l = cmul(p - vec2(1.0, 0.0), clog_fe(yf.m, yf.e)); + let k = i32(floor(l.x / LN2)); + return fe_make(cmul(p, cexp(vec2(l.x - f32(k) * LN2, l.y))), k); + } + // z^2-like kinds: 2y (|f'| = |2y| for the abs variants too, see fprime). + return Fe(2.0 * yf.m, yf.e); +} + +// The deep phase hands over to the f32 loop once the delta's magnitude +// reaches 2^DEEP_EXIT_LOG2: by then |e|^2 is still a normal f32 and the +// pixel offset dc (< 2^-99 on deep views) is below f32 rounding of e. The DE +// derivative only has to be a comfortably normal f32 (DEEP_EXIT_DZ_LOG2). +const DEEP_EXIT_LOG2: i32 = -48; +const DEEP_EXIT_DZ_LOG2: i32 = -100; +// Below this, a full value y goes through `deep_fprime`'s extended path. +const DEEP_TINY: f32 = 7.888609e-31; // 2^-100 +// The mantissas are renormalized once their exponent drifts past ±this. +const DEEP_RENORM_LOG2: i32 = 16; +// A reference point can only matter for rebasing when it's within this many +// binades above the delta's scale (|w| < 2^DEEP_RENORM_LOG2). +const DEEP_NEAR_LOG2: i32 = 24; + // Weight of the Phoenix kind's p*z_{n-1} term in the current map: 1 for plain // Phoenix, its morph share while switching to/from Phoenix, else 0. fn phoenix_weight() -> f32 { @@ -331,13 +667,14 @@ struct Sample { // Perturbation iterate 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 initial delta -// (c is fixed, so nothing is added per step). +// (c is fixed, so nothing is added per step). In `DEEP` pipelines both +// `offset` and `px` are in units of 2^u.scale_exp. fn iterate_sample(offset: vec2, px: f32) -> Sample { // Loop invariants, read once instead of on every iteration. let max_iter = u.max_iter; let bailout_sq = u.bailout_sq; let ref_len = u.ref_len; - let z0 = ref_orbit[0]; // reference start (0 for Mandelbrot, center for Julia) + let z0 = ref_at(0u); // reference start (0 for Mandelbrot, center for Julia) // Main cardioid / period-2 bulb bypass: those points never escape, so skip // iterating them (they'd otherwise all burn the full max_iter). `offset` is @@ -345,7 +682,7 @@ fn iterate_sample(offset: vec2, px: f32) -> Sample { // orbit itself, since X_1 = X_0^2 + C_ref = C_ref. That's only f32-accurate, // so skip the test once a pixel is smaller than that error (deep zoom), // where it could misclassify pixels right at the boundary. - if KIND == KIND_MANDELBROT && !MORPH && !IS_JULIA && ref_len > 1u && px > 1e-6 { + if KIND == KIND_MANDELBROT && !MORPH && !IS_JULIA && !DEEP && ref_len > 1u && px > 1e-6 { let c = ref_orbit[1] + offset; let xq = c.x - 0.25; let q = xq * xq + c.y * c.y; @@ -361,6 +698,8 @@ fn iterate_sample(offset: vec2, px: f32) -> Sample { // seeds the delta and nothing is added per step. var step_add = offset; var e = vec2(0.0, 0.0); + // Pixel size in complex units (`px` is pre-scaled in DEEP pipelines). + var px_t = px; // Orbit derivative for distance estimation, pre-multiplied by the pixel // size `px`. For the set plane it is px·d/dc (starts at 0, gains +px each // step); for Julia it is px·d/dz0 (starts at px). The raw derivative grows @@ -371,7 +710,7 @@ fn iterate_sample(offset: vec2, px: f32) -> Sample { if IS_JULIA { step_add = vec2(0.0, 0.0); e = offset; - dzs = vec2(px, 0.0); + dzs = vec2(px_t, 0.0); } // Previous-iterate state for the Phoenix two-term recurrence (delta of // y_{n-1}, and its scaled derivative for DE). Both start at 0 (y_{-1} = 0). @@ -403,8 +742,247 @@ fn iterate_sample(offset: vec2, px: f32) -> Sample { var cand_d2 = 3.0e38; var mult2_cand = 1.0; + // Deep phase (DEEP pipelines only): iterate the delta as e = w * 2^sx, + // with `w` an f32 mantissa kept near 1 by renormalizing and `sx` an i32 + // exponent, while |e| is too small for f32 (it starts at the pixel offset, + // ~2^scale_exp). The DE derivative is linear in the same way and carried + // as dzs = v * 2^sv, with its own exponent: near the critical point (after + // a rebase) f' is tiny, so dzs and e can drift far apart. Once both are + // big enough (or a rebase makes the delta large), the state is converted + // to plain f32 and the loop below carries on from the same n / m. + if DEEP { + let scale_e = u.scale_exp; + var sx = scale_e; + var sv = scale_e; + var sc = ldexp_sat(1.0, sx); + var w = vec2(0.0, 0.0); + var v = vec2(0.0, 0.0); + var w_prev = vec2(0.0, 0.0); + var v_prev = vec2(0.0, 0.0); + // Per-step additions (dc and px for the set plane), in units of 2^sx + // and 2^sv. + var d = offset; + var pd = px; + if IS_JULIA { + w = offset; + v = vec2(px, 0.0); + d = vec2(0.0, 0.0); + pd = 0.0; + } + let z0f = ref_fe(0u); + var xf = z0f; + var xt = z0; + var xf_old = xf; + var xt_old = xt; + var w_old = w; + var rebase_exit = false; + var cut_exit = false; + loop { + // Full value y = X + e as f32 (e flushes to 0 when negligible). + let yt = xt + ldexp2_sat(w, sx); + if dot(yt, yt) > bailout_sq { + escaped = true; + break; + } + if n >= max_iter { + break; + } + if deep_cut_crossing(xf, w, sx) { + cut_exit = true; + break; + } + + if DE { + let fp = deep_fprime(xf, w, sx, yt); + if fp.e == 0 { + var v_new = cmul(fp.m, v); + if !IS_JULIA { + v_new.x = v_new.x + pd; + } + if phoenix_w > 0.0 { + v_new = v_new + phoenix_w * cmul(u.phoenix_p, v_prev); + v_prev = v; + } + v = v_new; + } else { + // f' is below f32's range (next to the critical point): + // sum the terms at the largest one's scale and move + // there, like the delta below. + var acc = fe_make(cmul(fp.m, v), fp.e + sv); + if !IS_JULIA { + acc = fe_add(acc, fe_make(vec2(px, 0.0), scale_e)); + } + if phoenix_w > 0.0 { + acc = fe_add(acc, fe_make(phoenix_w * cmul(u.phoenix_p, v_prev), sv)); + } + if acc.m.x == 0.0 && acc.m.y == 0.0 { + if phoenix_w > 0.0 { + v_prev = v; + } + v = acc.m; + } else { + if phoenix_w > 0.0 { + v_prev = ldexp2_sat(v, sv - acc.e); + } + v = acc.m; + sv = acc.e; + if !IS_JULIA { + pd = ldexp_sat(px, scale_e - sv); + } + } + } + } + w_old = w; + let st = advance_delta_scaled(xt, xf, w, sc, sx); + if st.t == sx { + w = st.w + d; + if phoenix_w > 0.0 { + w = w + phoenix_w * cmul(u.phoenix_p, w_prev); + w_prev = w_old; + } + } else { + // The step's value is at another scale (next to the + // critical point it is ~e^2, far below 2^sx): add dc and the + // Phoenix term at the largest addend's scale and move there. + var acc = fe_make(st.w, st.t); + if !IS_JULIA { + acc = fe_add(acc, fe_make(offset, scale_e)); + } + if phoenix_w > 0.0 { + acc = fe_add(acc, fe_make(phoenix_w * cmul(u.phoenix_p, w_prev), sx)); + } + if acc.m.x == 0.0 && acc.m.y == 0.0 { + w = acc.m; + if phoenix_w > 0.0 { + w_prev = w_old; + } + } else { + // w_old (the pre-step delta) is still needed at the new + // scale: it's the next previous delta, and a rebase reads it. + w_old = ldexp2_sat(w_old, sx - acc.e); + if phoenix_w > 0.0 { + w_prev = w_old; + } + w = acc.m; + sx = acc.e; + sc = ldexp_sat(1.0, sx); + if !IS_JULIA { + d = ldexp2_sat(offset, scale_e - sx); + } + } + } + m = m + 1u; + n = n + 1u; + + if m >= ref_len { + escaped = true; // see the f32 loop's reference-exhausted case + break; + } + xf_old = xf; + xt_old = xt; + xf = ref_fe(m); + xt = ldexp2_sat(xf.m, xf.e); + + // Rebase test |X + e| < |e|, in units of 2^sx. Only possible when + // |X| is within a few binades of |e| (|w| < 2^DEEP_RENORM_LOG2). + if xf.e - sx < DEEP_NEAR_LOG2 { + let q = ldexp2_sat(xf.m, xf.e - sx) + w; + if dot(q, q) < dot(w, w) { + // The new delta y - X[0] stays tiny only if X[0] is + // (always, for the set plane), and for Phoenix only if the + // previous full value y_{n-1} (its new previous delta) is. + let z0_small = (z0f.m.x == 0.0 && z0f.m.y == 0.0) + || z0f.e - sx < DEEP_NEAR_LOG2; + let prev_small = phoenix_w == 0.0 || xf_old.e - sx < DEEP_NEAR_LOG2; + if !(z0_small && prev_small) { + rebase_exit = true; + break; + } + if phoenix_w > 0.0 { + w_prev = ldexp2_sat(xf_old.m, xf_old.e - sx) + w_old; + } + w = q - ldexp2_sat(z0f.m, z0f.e - sx); + xf = z0f; + xt = z0; + m = 0u; + } + } + + // Leave once both the delta and the derivative fit in f32 (an + // exactly-zero one always does); otherwise renormalize the + // mantissas when they drift (exact: powers of two only). + let kw = deep_log2(w, w_prev); + let kv = deep_log2(v, v_prev); + let w_ok = kw == DEEP_ZERO || sx + kw > DEEP_EXIT_LOG2; + let v_ok = !DE || kv == DEEP_ZERO || sv + kv > DEEP_EXIT_DZ_LOG2; + if w_ok && v_ok { + break; + } + if kw != DEEP_ZERO && abs(kw) > DEEP_RENORM_LOG2 { + w = ldexp2_sat(w, -kw); + w_prev = ldexp2_sat(w_prev, -kw); + sx = sx + kw; + sc = ldexp_sat(1.0, sx); + if !IS_JULIA { + d = ldexp2_sat(offset, scale_e - sx); + } + } + if DE && kv != DEEP_ZERO && abs(kv) > DEEP_RENORM_LOG2 { + v = ldexp2_sat(v, -kv); + v_prev = ldexp2_sat(v_prev, -kv); + sv = sv + kv; + if !IS_JULIA { + pd = ldexp_sat(px, scale_e - sv); + } + } + } + + // Hand over to the f32 loop. dc and px may flush to 0 here: they + // are below f32 rounding of the (now large enough) delta and + // derivative. + px_t = ldexp_sat(px, scale_e); + if !IS_JULIA { + step_add = ldexp2_sat(offset, scale_e); + } + e = ldexp2_sat(w, sx); + e_prev = ldexp2_sat(w_prev, sx); + dzs = ldexp2_sat(v, sv); + dzs_prev = ldexp2_sat(v_prev, sv); + xm = xt; + if rebase_exit { + // Rebase in plain f32: the new delta (y - X[0], and for Phoenix + // the previous full value) is no longer tiny. + if phoenix_w > 0.0 { + e_prev = xt_old + ldexp2_sat(w_old, sx); + } + e = (xt + e) - z0; + xm = z0; + m = 0u; + } + if cut_exit && e.y == 0.0 && w.y != 0.0 { + // Keep the side of the cut the pixel is on even if e.y flushed + // to 0: that decides the branch in the next (f32) step. + e.y = select(-1.17549435e-38, 1.17549435e-38, w.y > 0.0); + } + z = xm + e; + z2 = dot(z, z); + // Restart periodicity detection from here (check_at must stay ahead + // of n, or no window would ever close). Nothing is saved until the + // first window closes: `z` here is at an arbitrary phase, and the + // pixel still shadows the reference (exactly periodic when it's a + // minibrot nucleus), so comparing against it would flag exterior + // pixels as interior (see PERIOD_FIRST_CHECK). The sentinel is far + // outside the bailout radius, so no return can match it. + z_saved = vec2(1e18, 1e18); + z_cand = z; + while check_at <= n { + check_at = check_at * 2u; + } + } + loop { - if z2 > bailout_sq { + // (`escaped` may already be set by the deep phase.) + if escaped || z2 > bailout_sq { escaped = true; break; } @@ -428,7 +1006,7 @@ fn iterate_sample(offset: vec2, px: f32) -> Sample { if DE { var dzs_new = cmul(fp, dzs); if !IS_JULIA { - dzs_new.x = dzs_new.x + px; + dzs_new.x = dzs_new.x + px_t; } if phoenix_w > 0.0 { dzs_new = dzs_new + phoenix_w * cmul(u.phoenix_p, dzs_prev); @@ -457,7 +1035,7 @@ fn iterate_sample(offset: vec2, px: f32) -> Sample { escaped = true; break; } - xm = ref_orbit[m]; + xm = ref_at(m); z = xm + e; z2 = dot(z, z); if z2 < dot(e, e) { diff --git a/src/view.rs b/src/view.rs index af1a399..f2bc620 100644 --- a/src/view.rs +++ b/src/view.rs @@ -1,9 +1,11 @@ //! Camera / view state over the complex plane. //! //! The center is stored in arbitrary precision (`FBig`) — this is what lets us -//! zoom far past f64's ~1e13x limit. The pixel *scale* stays `f64`: even at -//! 10^30x zoom the scale is ~1e-33, comfortably inside f64's range. Only the -//! center needs the extra digits. +//! zoom far past f64's ~1e13x limit. The pixel *scale* stays `f64`, which +//! bounds zoom at ~10^300x (`MIN_HALF_HEIGHT`). Only the center needs the +//! extra digits. Once a pixel is smaller than `DEEP_PIXEL_SIZE` the GPU +//! switches to rescaled deltas (see `needs_deep`), since f32 alone bottoms +//! out near 1e-38. use core::str::FromStr; @@ -16,6 +18,34 @@ 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; +/// Smallest half-height the view can zoom to: f64's range (the pixel scale, +/// and the rescaled GPU uniforms, are computed in f64). +pub const MIN_HALF_HEIGHT: f64 = 1e-300; + +/// 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 +/// offsets (a quarter pixel, for the AA grid) are normal floats (>= 2^-126), +/// i.e. down to 2^-124 per pixel. Measured: pixel-identical to the deep path +/// down to 2^-124 (no AA), first errors at 2^-126, all black by 2^-136. The +/// deep path is slower, so the switch is as late as that allows, with two +/// binades of margin. +pub const DEEP_PIXEL_SIZE: f64 = 1.0 / (1u128 << 122) as f64; // 2^-122 + +/// Whether a view rendered `height_px` pixels tall needs the deep pipeline +/// (see `DEEP_PIXEL_SIZE`). +pub fn needs_deep(half_height: f64, height_px: f64) -> bool { + 2.0 * half_height / height_px.max(1.0) < DEEP_PIXEL_SIZE +} + +/// Binary exponent `E` of the deep view scale: `floor(log2(half_height))`, +/// so the rescaled span is in `[2, 4)`. Never 0, which means "not deep" +/// (see `Uniforms::scale_exp`). +pub fn deep_scale_exp(half_height: f64) -> i32 { + let e = half_height.max(MIN_HALF_HEIGHT).log2().floor() as i32; + if e == 0 { -1 } else { e } +} + /// Guard bits added on top of the zoom-dictated precision. const GUARD_BITS: usize = 48; /// Upper bound on center precision (f32 GPU perturbation degrades long before @@ -103,7 +133,7 @@ impl ViewState { let k = cpp * (1.0 - factor); self.center_re = &self.center_re + &big_from_f64(off_x * k, bits); self.center_im = &self.center_im + &big_from_f64(off_y * k, bits); - self.half_height *= factor; + self.half_height = (self.half_height * factor).max(MIN_HALF_HEIGHT); } /// Build a view from full-precision center coordinates and a half-height. @@ -111,7 +141,7 @@ impl ViewState { let mut v = Self { center_re, center_im, - half_height, + half_height: half_height.max(MIN_HALF_HEIGHT), }; v.sync_precision(); v @@ -138,6 +168,7 @@ pub fn parse_view_spec(spec: &str) -> Option<(ViewState, Option)> { if !(half_height > 0.0 && half_height.is_finite()) { return None; } + let half_height = half_height.max(MIN_HALF_HEIGHT); let bits = precision_for(half_height); let re = big_from_decimal_str(parts[0], bits)?; let im = big_from_decimal_str(parts[1], bits)?; @@ -153,7 +184,7 @@ pub fn parse_half_height_spec(spec: &str) -> Option { if !(half_height > 0.0 && half_height.is_finite()) { return None; } - Some(half_height) + Some(half_height.max(MIN_HALF_HEIGHT)) } /// Parse a "re,im" spec (re/im decimal, parsed at /// full precision) into a view. Shared by diff --git a/src/worker.rs b/src/worker.rs index 85dacf2..a427f1f 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::{FractalKind, RefOrbit, compute_reference, compute_set_reference}; use crate::view::{Big, big_from_f64}; pub struct RefRequest { @@ -36,7 +36,7 @@ pub struct RefResult { pub center_re: Big, pub center_im: Big, pub half_height: f64, - pub points: Vec<[f32; 2]>, + pub points: RefOrbit, /// The kind and morph `points` was computed with (echoed from the /// request). pub kind: FractalKind, @@ -102,7 +102,7 @@ impl RefWorker { } } -fn compute(req: &RefRequest) -> Vec<[f32; 2]> { +fn compute(req: &RefRequest) -> RefOrbit { if req.julia { let jr = big_from_f64(req.julia_c.0, req.precision); let ji = big_from_f64(req.julia_c.1, req.precision); diff --git a/tests/shader_valid.rs b/tests/shader_valid.rs index b73371d..51ca8a1 100644 --- a/tests/shader_valid.rs +++ b/tests/shader_valid.rs @@ -77,7 +77,7 @@ fn mandelbrot_shader_is_valid() { } /// Every specialization renderer.rs can build (`PipelineKey`: kind × Julia × -/// DE × morph), for every fragment entry point. +/// DE × morph × deep), for every fragment entry point. #[test] fn mandelbrot_shader_specializations_compile() { let (module, info) = validate("mandelbrot.wgsl", MANDELBROT_SRC); @@ -85,21 +85,24 @@ fn mandelbrot_shader_specializations_compile() { for julia in [0.0, 1.0] { for de in [0.0, 1.0] { for morph in [0.0, 1.0] { - let constants = [ - ("KIND", kind as f64), - ("IS_JULIA", julia), - ("DE", de), - ("MORPH", morph), - ]; - for entry in ["fs_data", "fs_refine", "fs_color"] { - specialize( - "mandelbrot.wgsl", - &module, - &info, - naga::ShaderStage::Fragment, - entry, - &constants, - ); + for deep in [0.0, 1.0] { + let constants = [ + ("KIND", kind as f64), + ("IS_JULIA", julia), + ("DE", de), + ("MORPH", morph), + ("DEEP", deep), + ]; + for entry in ["fs_data", "fs_refine", "fs_color"] { + specialize( + "mandelbrot.wgsl", + &module, + &info, + naga::ShaderStage::Fragment, + entry, + &constants, + ); + } } } }