diff --git a/CLAUDE.md b/CLAUDE.md index 1cbe8f4..28317c4 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -11,8 +11,9 @@ and every pixel is rendered on the GPU as a cheap `f32` delta from it, with 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 +rescaled deltas (f32 mantissa × 2^i32). There's no practical depth limit: +`half_height` is a `view::Scale` (f64 mantissa × 2^i32), floored only at +`Scale::MIN` = 2^-(2^20) to keep shader exponent sums in i32. 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). @@ -94,8 +95,13 @@ 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` (so zoom is clamped at - `MIN_HALF_HEIGHT` = 1e-300). Precision (bits) scales with zoom depth + type alias). The pixel scale (`half_height`) is a `Scale`, an f64 + mantissa with its own i32 exponent, so it goes past f64's ~1e-308. Never + collapse it (or a center difference) to a plain `f64` on a path used at + depth. Rescale first: `Scale::scaled_f64(k)`, or shift the `Big` by + `-scale_exp` before `to_f64()`, as `dc_offset`/`drift_from` do. + `Display`/`FromStr` use scientific notation of any exponent (share links, + `--view`, the zoom field). 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 diff --git a/src/app.rs b/src/app.rs index d6e34f7..b40d121 100644 --- a/src/app.rs +++ b/src/app.rs @@ -19,7 +19,7 @@ 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, MIN_HALF_HEIGHT, ViewState, big_from_decimal_str, big_from_f64, + Big, DEFAULT_HALF_HEIGHT, MAX_PRECISION_BITS, Scale, ViewState, big_from_decimal_str, big_from_f64, big_to_decimal_str, deep_scale_exp, interpolate_view, needs_deep, parse_view_spec, precision_for, }; @@ -206,7 +206,7 @@ impl RefJob { struct RequestKey { center_re: Big, center_im: Big, - half_height: f64, + half_height: Scale, julia: bool, julia_c: (f64, f64), phoenix_p: (f64, f64), @@ -556,7 +556,7 @@ pub struct FractalApp { /// slightly from the live view; the shader compensates via `dc_offset`). ref_center_re: Big, ref_center_im: Big, - ref_half_height: f64, + ref_half_height: Scale, /// Kind and kind-switch morph the current `reference` was computed with. /// The shader iterates with these (not the live kind/morph) so its delta /// formula always matches the orbit, even while the worker lags a frame @@ -614,8 +614,8 @@ fn sig_digits_for(bits: usize) -> usize { } /// Format a magnification for the editable field (compact scientific). -fn format_zoom(m: f64) -> String { - format!("{m:.4e}") +fn format_zoom(m: Scale) -> String { + format!("{m:.4}") } /// Precision (bits) to parse a typed center at: at least what the current zoom @@ -624,7 +624,7 @@ fn format_zoom(m: f64) -> String { fn parse_bits_for(s: &str, min_bits: usize) -> usize { let digits = s.chars().filter(char::is_ascii_digit).count(); let from_input = (digits as f64 * std::f64::consts::LOG2_10).ceil() as usize + 16; - min_bits.max(from_input).min(2048) + min_bits.max(from_input).min(MAX_PRECISION_BITS) } impl FractalApp { @@ -746,7 +746,7 @@ impl FractalApp { if let Some(k) = cli.kind { self.kind = k.into(); if let Some(p) = cli.power { - self.power = p.clamp(2, 8); + self.power = p.clamp(2, 20); } self.view = Self::default_view_for(self.mode, self.kind); } @@ -979,6 +979,7 @@ impl FractalApp { /// Jump to a preset Mandelbrot location: decimal center (parsed at the /// precision the zoom needs), half-height, and a fitting iteration count. fn go_to_place(&mut self, re: &str, im: &str, half_height: f64, iterations: u32) { + let half_height = Scale::from_f64(half_height); let bits = precision_for(half_height); if let (Some(cre), Some(cim)) = ( big_from_decimal_str(re, bits), @@ -997,7 +998,7 @@ impl FractalApp { /// `auto_iterations` is on. Grows roughly linearly with zoom decades so deep /// zooms keep enough iterations to stay sharp instead of banding. fn auto_iteration_count(&self) -> u32 { - let decades = self.view.magnification().log10().max(0.0); + let decades = self.view.magnification_log10().max(0.0); let iters = 400.0 + 900.0 * decades; (iters.round() as u32).clamp(200, MAX_REF_POINTS as u32 - 1) } @@ -1033,7 +1034,7 @@ impl FractalApp { FractalMode::Mandelbrot }; self.kind = s.kind; - self.power = s.power.clamp(2, 8); + self.power = s.power.clamp(2, 200); self.julia_c = s.julia_c; self.phoenix_p = s.phoenix_p; self.lambda_l = s.lambda_l; @@ -1077,10 +1078,10 @@ impl FractalApp { /// each kind's interesting region. fn default_view_for(mode: FractalMode, kind: FractalKind) -> ViewState { if mode == FractalMode::Julia { - return ViewState::with_center(big_from_f64(0.0, 53), big_from_f64(0.0, 53), 1.5); + return ViewState::with_center(big_from_f64(0.0, 53), big_from_f64(0.0, 53), Scale::from_f64(1.5)); } let (cr, ci, hh) = kind.default_set_view(); - ViewState::with_center(big_from_f64(cr, 53), big_from_f64(ci, 53), hh) + ViewState::with_center(big_from_f64(cr, 53), big_from_f64(ci, 53), Scale::from_f64(hh)) } /// The request key for the current state. Its `iter` is the reference @@ -1104,10 +1105,14 @@ impl FractalApp { } /// Distance (complex units) the live view center has drifted from `key`. + /// Distance of the live center from `key`'s, in units of the live + /// half-height (measured at that scale, so it works past f64's range). fn drift_from(&self, key: &RequestKey) -> f64 { - let dre = (&self.view.center_re - &key.center_re).to_f64().value(); - let dim = (&self.view.center_im - &key.center_im).to_f64().value(); - (dre * dre + dim * dim).sqrt() + let hh = self.view.half_height; + let k = -hh.exponent() as isize; + let dre = ((&self.view.center_re - &key.center_re) << k).to_f64().value(); + let dim = ((&self.view.center_im - &key.center_im) << k).to_f64().value(); + (dre * dre + dim * dim).sqrt() / hh.scaled_f64(-hh.exponent()) } /// Whether the reference should be (re)computed: parameters changed, or the @@ -1143,19 +1148,22 @@ impl FractalApp { && self.morph.is_none() { // But still recompute on significant zoom changes for precision - let ratio = self.view.half_height / key.half_height; + let ratio = self.view.half_height.ratio(key.half_height); return !(0.5..=2.0).contains(&ratio); } - let ratio = self.view.half_height / key.half_height; - self.drift_from(key) > 0.5 * self.view.half_height || !(0.5..=2.0).contains(&ratio) + let ratio = self.view.half_height.ratio(key.half_height); + self.drift_from(key) > 0.5 || !(0.5..=2.0).contains(&ratio) } - /// 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) + /// Complex offset of the live view center from the reference center, in + /// units of `2^scale_exp` (shifted exactly in `Big`, so it doesn't + /// underflow f64 at deep zooms). + fn dc_offset(&self, scale_exp: i32) -> (f64, f64) { + let k = -scale_exp as isize; + let dre = ((&self.view.center_re - &self.ref_center_re) << k) .to_f64() .value(); - let dim = (&self.view.center_im - &self.ref_center_im) + let dim = ((&self.view.center_im - &self.ref_center_im) << k) .to_f64() .value(); (dre, dim) @@ -1178,7 +1186,7 @@ impl FractalApp { points: RefOrbit, cre: Big, cim: Big, - hh: f64, + hh: Scale, kind: FractalKind, morph: Option<(FractalKind, f32)>, ) { @@ -1404,10 +1412,12 @@ impl FractalApp { 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(); + let (dc_re, dc_im) = self.dc_offset(scale_exp); Uniforms { - span: [(span_x * inv_scale) as f32, (span_y * inv_scale) as f32], + span: [ + span_x.scaled_f64(-scale_exp) as f32, + span_y.scaled_f64(-scale_exp) 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, @@ -1420,7 +1430,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: [(dc_re * inv_scale) as f32, (dc_im * inv_scale) as f32], + dc_offset: [dc_re as f32, dc_im 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], @@ -1450,7 +1460,7 @@ impl FractalApp { ]; BuddhabrotUniforms { center, - half_height: self.view.half_height as f32, + half_height: self.view.half_height.to_f64() as f32, aspect: aspect 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], @@ -2016,11 +2026,10 @@ impl FractalApp { } } if self.anim.zoom && self.anim.zoom_speed != 0.0 { - let min_hh = DEFAULT_HALF_HEIGHT * 1.0e-26; // practical f32-perturbation depth - let max_hh = DEFAULT_HALF_HEIGHT * 4.0; + let max_hh = Scale::from_f64(DEFAULT_HALF_HEIGHT * 4.0); let factor = (-(self.anim.zoom_speed as f64) * dt).exp(); - let target = (self.view.half_height * factor).clamp(min_hh, max_hh); - let f = target / self.view.half_height; + let target = self.view.half_height.mul_f64(factor).clamp(Scale::MIN, max_hh); + let f = target.ratio(self.view.half_height); if (f - 1.0).abs() > 1.0e-9 { self.view .zoom_at_pixel(0.0, 0.0, self.last_size_px.y.max(1.0) as f64, f); @@ -2453,11 +2462,9 @@ impl FractalApp { } if zoom_resp.lost_focus() { if self.zoom_edited - && let Ok(hh) = self.zoom_edit.trim().parse::() - && hh > 0.0 - && hh.is_finite() + && let Ok(hh) = self.zoom_edit.parse::() { - self.view.half_height = hh.max(MIN_HALF_HEIGHT); + self.view.half_height = hh; self.view.sync_precision(); } self.zoom_edited = false; @@ -2560,7 +2567,7 @@ impl FractalApp { }); ui.checkbox(&mut self.buddha_accumulate, "Keep sampling") .on_hover_text("Dispatch a fresh batch of random samples every frame."); - if self.view.magnification() > 1.0e5 { + if self.view.magnification_log10() > 5.0 { ui.colored_label( egui::Color32::LIGHT_YELLOW, "deep zoom isn't supported here (f32 precision only)", diff --git a/src/fractal/reference.rs b/src/fractal/reference.rs index 9557b30..3b1afc0 100644 --- a/src/fractal/reference.rs +++ b/src/fractal/reference.rs @@ -380,7 +380,7 @@ fn step( FractalKind::Perpendicular => { // (x^2 - y^2) - 2·x·|y| i: abs the imaginary input. let re = &zr.sqr() - &zi.sqr() + cr; - let im = if zi.to_f64().value() < 0.0 { + let im = if *zi < Big::ZERO { ci + ((zr * zi) << 1) } else { ci - ((zr * zi) << 1) @@ -422,10 +422,10 @@ fn big_zero(precision: usize) -> Big { Big::from(0i32).with_precision(precision).value() } -/// Absolute value of a `Big`. The sign check via f64 is exact except for values -/// so tiny that |x| ≈ x either way — negligible against the f32 orbit storage. +/// Absolute value of a `Big`. The sign comes from the `Big` itself: through +/// f64, anything below ~1e-308 reads as ±0 and would keep its sign. fn big_abs(x: Big) -> Big { - if x.to_f64().value() < 0.0 { -x } else { x } + if x < Big::ZERO { -x } else { x } } /// `(zr + i zi)^power` by repeated complex multiply at `precision` bits. @@ -442,10 +442,10 @@ fn complex_pow(zr: &Big, zi: &Big, power: u32, precision: usize) -> (Big, Big) { (rr, ri) } -/// `true` if `x` is (numerically) zero. The f64 check is exact for a true -/// zero; only matters here to special-case `ln(0)`. +/// `true` if `x` is exactly zero (special-cases `ln(0)`). Not via f64, +/// which flushes values below ~1e-308 to zero. fn is_big_zero(x: &Big) -> bool { - x.to_f64().value() == 0.0 + x.repr().significand().is_zero() } /// `(zr + i zi)^(pr + i pi)` for a complex exponent, via the principal branch @@ -504,6 +504,17 @@ pub fn compute_set_reference( mod tests { use super::*; + /// Sign and zero tests must hold far below f64's range, where + /// `to_f64` reads as ±0. + #[test] + fn sign_and_zero_below_f64_range() { + let tiny = Big::try_from(1.0_f64).unwrap().with_precision(64).value() >> 5000; + assert!(!is_big_zero(&tiny)); + assert!(is_big_zero(&big_zero(64))); + assert_eq!(big_abs(-tiny.clone()), tiny); + assert_eq!(big_abs(tiny.clone()), tiny); + } + /// The high-precision reference must agree with a plain f64 iteration for a /// shallow point (where f64 is accurate). #[test] diff --git a/src/fractal/share.rs b/src/fractal/share.rs index c70574f..0f30ce4 100644 --- a/src/fractal/share.rs +++ b/src/fractal/share.rs @@ -2,13 +2,14 @@ //! iterations, Julia constant, coloring) as a compact URL fragment so deep-zoom //! locations can be shared or bookmarked. //! -//! Format: `m=m&f=&re=&im=&hh=&it=&cs=&co=` with +//! Format: `m=m&f=&re=&im=&hh=&it=&cs=&co=` with //! `m=j&jr=&ji=` added for Julia. `re`/`im` are full-precision decimal -//! strings. +//! strings; `hh` is a `Scale` in scientific notation (any exponent). use std::collections::HashMap; use crate::fractal::FractalKind; +use crate::view::Scale; #[derive(Clone, Debug)] pub struct ShareState { @@ -18,7 +19,7 @@ pub struct ShareState { pub power: u32, pub center_re: String, pub center_im: String, - pub half_height: f64, + pub half_height: Scale, pub iterations: u32, pub julia_c: (f64, f64), /// Distortion constant for the Phoenix kind (ignored by others). @@ -115,7 +116,7 @@ mod tests { power: 5, center_re: "-0.743643887037158704752191506114774".into(), center_im: "0.131825904205311970493132056385139".into(), - half_height: 1.5e-20, + half_height: Scale::from_f64(1.5e-20), iterations: 4000, julia_c: (-0.123, 0.745), phoenix_p: (-0.5, 0.1), @@ -146,5 +147,15 @@ mod tests { let d = ShareState::decode("#m=m&re=0.0&im=0.0&hh=1.25&it=256").unwrap(); assert!(!d.julia); assert_eq!(d.iterations, 256); + assert_eq!(d.half_height, Scale::from_f64(1.25)); + } + + /// Zooms past f64's range survive a round trip. + #[test] + fn round_trip_past_f64_range() { + let d = ShareState::decode("#m=m&re=0.0&im=0.0&hh=1.5e-1234&it=256").unwrap(); + assert_eq!(d.half_height, "1.5e-1234".parse().unwrap()); + let d2 = ShareState::decode(&d.encode()).unwrap(); + assert_eq!(d2.half_height, d.half_height); } } diff --git a/src/view.rs b/src/view.rs index f2bc620..6fff78e 100644 --- a/src/view.rs +++ b/src/view.rs @@ -1,11 +1,12 @@ //! 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`, 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. +//! zoom far past f64's ~1e13x limit. The pixel *scale* is a [`Scale`]: an +//! f64 mantissa with its own `i32` binary exponent, so it isn't bound by +//! f64's ~1e-308 range either (the floor, `Scale::MIN`, only keeps the GPU's +//! i32 exponent arithmetic far from overflow). 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; @@ -18,9 +19,6 @@ pub type Big = FBig; /// Half-height (complex units) of the default view; also the zoom-1 reference. pub const DEFAULT_HALF_HEIGHT: f64 = 1.25; -/// 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 @@ -34,62 +32,266 @@ 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 +pub fn needs_deep(half_height: Scale, height_px: f64) -> bool { + half_height.mul_f64(2.0 / height_px.max(1.0)) < Scale::from_f64(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; +pub fn deep_scale_exp(half_height: Scale) -> i32 { + let e = half_height.exponent(); if e == 0 { -1 } else { e } } +/// `x * 2^k` for any `k`, saturating to 0 / infinity like the true value +/// would (`powi` alone overflows at 2^±1024 even when the product fits). +fn ldexp(mut x: f64, mut k: i32) -> f64 { + while k > 1000 { + x *= 2f64.powi(1000); + k -= 1000; + if x.is_infinite() || x == 0.0 { + return x; + } + } + while k < -1000 { + x *= 2f64.powi(-1000); + k += 1000; + if x == 0.0 || x.is_infinite() { + return x; + } + } + x * 2f64.powi(k) +} + +/// A positive real with f64 precision and an `i32` binary exponent: +/// `m · 2^e`, `m` in `[1, 2)`. The view's half-height (and the pixel size +/// derived from it) is one of these, so zoom isn't bound by f64's range. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct Scale { + m: f64, + e: i32, +} + +impl Scale { + /// Deepest scale the view can reach: 2^-(2^20) (about 1e-315653). Far + /// past anything a reference orbit can practically be computed for; it + /// only keeps the shader's i32 exponent sums (scale × degree) from + /// overflowing. + pub const MIN: Scale = Scale { + m: 1.0, + e: -(1 << 20), + }; + + /// `m · 2^e`, normalized. Non-positive or NaN input gives `MIN`. + pub fn from_parts(m: f64, e: i32) -> Self { + if m.is_nan() || m <= 0.0 { + return Self::MIN; + } + if m.is_infinite() { + return Scale { m: 1.0, e: i32::MAX / 2 }; + } + // Bring m into [1, 2) through its own binary exponent (exact). + let k = m.log2().floor() as i32; + let mut m = ldexp(m, -k); + let mut e = e.saturating_add(k); + // log2 can round across a power of two. + if m >= 2.0 { + m /= 2.0; + e = e.saturating_add(1); + } else if m < 1.0 { + m *= 2.0; + e = e.saturating_sub(1); + } + Scale { m, e }.max(Self::MIN) + } + + pub fn from_f64(x: f64) -> Self { + Self::from_parts(x, 0) + } + + /// `2^l`. + pub fn from_log2(l: f64) -> Self { + let e = l.floor(); + Self::from_parts((l - e).exp2(), e as i32) + } + + /// The value as an f64 (0 or infinity outside its range). + pub fn to_f64(self) -> f64 { + ldexp(self.m, self.e) + } + + /// `self · 2^k` as an f64: the value in units of `2^-k`. + pub fn scaled_f64(self, k: i32) -> f64 { + ldexp(self.m, self.e.saturating_add(k)) + } + + /// `floor(log2(self))`. + pub fn exponent(self) -> i32 { + self.e + } + + pub fn log2(self) -> f64 { + self.m.log2() + self.e as f64 + } + + pub fn log10(self) -> f64 { + self.log2() * core::f64::consts::LOG10_2 + } + + /// `self · f` (`f > 0`). + pub fn mul_f64(self, f: f64) -> Self { + Self::from_parts(self.m * f, self.e) + } + + /// `self / other`, as an f64. + pub fn ratio(self, other: Scale) -> f64 { + ldexp(self.m / other.m, self.e.saturating_sub(other.e)) + } + + pub fn max(self, other: Scale) -> Self { + if other > self { other } else { self } + } + + pub fn min(self, other: Scale) -> Self { + if other < self { other } else { self } + } + + pub fn clamp(self, lo: Scale, hi: Scale) -> Self { + self.max(lo).min(hi) + } + + /// `f · self` as an exact `Big` at `bits` of precision (`f` any f64). + pub fn big_times(self, f: f64, bits: usize) -> Big { + big_from_f64(f * self.m, bits) << self.e as isize + } + + /// Exact binary value as a `Big`. + fn to_big(self) -> Big { + big_from_f64(self.m, 53) << self.e as isize + } +} + +impl PartialOrd for Scale { + fn partial_cmp(&self, other: &Self) -> Option { + // Normalized and positive: the exponent decides, then the mantissa. + Some( + self.e + .cmp(&other.e) + .then(self.m.partial_cmp(&other.m)?), + ) + } +} + +impl core::fmt::Display for Scale { + /// Scientific notation, `1.5e-20` / `3.7e-4000`. The precision flag + /// (`{:.4}`) sets mantissa digits after the point; without it, enough + /// digits to parse back to the same value. + fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result { + let x = self.to_f64(); + if x.is_normal() { + return match f.precision() { + Some(p) => write!(f, "{x:.p$e}"), + None => write!(f, "{x:e}"), + }; + } + // Out of f64's range: round the exact decimal expansion instead. + let sig = f.precision().map_or(17, |p| p + 1); + let dec = self.to_big().to_decimal().value().with_precision(sig).value(); + let repr = dec.repr(); + let digits = repr.significand().to_string(); + let digits = digits.trim_end_matches('0'); + let digits = if digits.is_empty() { "0" } else { digits }; + // value = significand · 10^exponent; move the point after the first digit. + let exp10 = repr.exponent() + repr.significand().to_string().len() as isize - 1; + let (head, tail) = digits.split_at(1); + let tail = match f.precision() { + Some(p) => format!("{tail:0 tail.to_string(), + }; + if tail.is_empty() { + write!(f, "{head}e{exp10}") + } else { + write!(f, "{head}.{tail}e{exp10}") + } + } +} + +impl FromStr for Scale { + type Err = (); + + /// Parses any positive decimal (`1.25`, `1.5e-20`, `3.7e-4000`), rounding + /// to the nearest f64 mantissa. Values past `MIN` clamp to it. + fn from_str(s: &str) -> Result { + let s = s.trim(); + if let Ok(x) = s.parse::() + && x.is_normal() + { + return if x > 0.0 { Ok(Self::from_f64(x)) } else { Err(()) }; + } + // Too small (or large) for f64: go through an exact decimal. + let dec = DBig::from_str(s).map_err(|_| ())?; + let bin: Big = dec.with_base_and_precision::<2>(64).value(); + if bin < Big::ZERO { + return Err(()); + } + let repr = bin.repr(); + let digits = repr.digits(); + if digits == 0 { + return Err(()); + } + let top = repr.exponent() + digits as isize - 1; + let m = (bin.clone() >> top).to_f64().value(); + let e = top.clamp(i32::MIN as isize, i32::MAX as isize) as i32; + Ok(Self::from_parts(m, 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 -/// this; the cap just prevents pathological allocation). -const MAX_PRECISION_BITS: usize = 2048; +/// Upper bound on center precision: what `Scale::MIN` needs. Only a guard +/// against pathological input; the reference orbit is impractically slow +/// long before this. +pub const MAX_PRECISION_BITS: usize = (1 << 20) + GUARD_BITS; #[derive(Clone, Debug)] pub struct ViewState { pub center_re: Big, pub center_im: Big, /// Half the view height in complex-plane units. Zooming in shrinks this. - pub half_height: f64, + pub half_height: Scale, } impl Default for ViewState { fn default() -> Self { - let bits = precision_for(DEFAULT_HALF_HEIGHT); + let bits = precision_for(Scale::from_f64(DEFAULT_HALF_HEIGHT)); Self { center_re: big_from_f64(-0.5, bits), center_im: big_from_f64(0.0, bits), - half_height: DEFAULT_HALF_HEIGHT, + half_height: Scale::from_f64(DEFAULT_HALF_HEIGHT), } } } impl ViewState { /// Complex-plane span (width, height) for the given pixel aspect ratio. - pub fn span(&self, aspect: f64) -> (f64, f64) { - let h = self.half_height * 2.0; - (h * aspect, h) + pub fn span(&self, aspect: f64) -> (Scale, Scale) { + let h = self.half_height.mul_f64(2.0); + (h.mul_f64(aspect), h) } /// Complex-plane units per pixel, given the viewport height in pixels. - pub fn complex_per_pixel(&self, height_px: f64) -> f64 { - (self.half_height * 2.0) / height_px + pub fn complex_per_pixel(&self, height_px: f64) -> Scale { + self.half_height.mul_f64(2.0 / height_px) } - /// Current magnification relative to the default view. - pub fn magnification(&self) -> f64 { - DEFAULT_HALF_HEIGHT / self.half_height + /// log10 of the current magnification relative to the default view. + pub fn magnification_log10(&self) -> f64 { + DEFAULT_HALF_HEIGHT.log10() - self.half_height.log10() } /// Current zoom level. - pub fn zoom(&self) -> f64 { + pub fn zoom(&self) -> Scale { self.half_height } @@ -116,8 +318,8 @@ impl ViewState { let cpp = self.complex_per_pixel(height_px); let bits = self.precision_bits(); // Grab-and-drag: moving the mouse right shows content to the left. - self.center_re = &self.center_re - &big_from_f64(dx * cpp, bits); - self.center_im = &self.center_im - &big_from_f64(dy * cpp, bits); + self.center_re = &self.center_re - &cpp.big_times(dx, bits); + self.center_im = &self.center_im - &cpp.big_times(dy, bits); } /// Zoom by `factor` (<1 zooms in) keeping the complex point currently under @@ -130,18 +332,18 @@ impl ViewState { // The cursor's complex offset from the center is (off * cpp). Keeping it // fixed while scaling the view by `factor` moves the center by // off * cpp * (1 - factor). (Derivation: new_c = fixed + (c-fixed)*f.) - 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 = (self.half_height * factor).max(MIN_HALF_HEIGHT); + let k = 1.0 - factor; + self.center_re = &self.center_re + &cpp.big_times(off_x * k, bits); + self.center_im = &self.center_im + &cpp.big_times(off_y * k, bits); + self.half_height = self.half_height.mul_f64(factor); } /// Build a view from full-precision center coordinates and a half-height. - pub fn with_center(center_re: Big, center_im: Big, half_height: f64) -> Self { + pub fn with_center(center_re: Big, center_im: Big, half_height: Scale) -> Self { let mut v = Self { center_re, center_im, - half_height: half_height.max(MIN_HALF_HEIGHT), + half_height, }; v.sync_precision(); v @@ -164,11 +366,7 @@ pub fn parse_view_spec(spec: &str) -> Option<(ViewState, Option)> { if parts.len() < 3 { return None; } - let half_height = parts[2].trim().parse::().ok()?; - if !(half_height > 0.0 && half_height.is_finite()) { - return None; - } - let half_height = half_height.max(MIN_HALF_HEIGHT); + let half_height = parse_half_height_spec(parts[2])?; let bits = precision_for(half_height); let re = big_from_decimal_str(parts[0], bits)?; let im = big_from_decimal_str(parts[1], bits)?; @@ -179,12 +377,8 @@ pub fn parse_view_spec(spec: &str) -> Option<(ViewState, Option)> { /// Parse a half_height spec. Shared by /// `FractalApp::apply_half_height_spec` (the `--zoom` CLI flag) and headless /// animation's `--to-zoom`. -pub fn parse_half_height_spec(spec: &str) -> Option { - let half_height = spec.trim().parse::().ok()?; - if !(half_height > 0.0 && half_height.is_finite()) { - return None; - } - Some(half_height.max(MIN_HALF_HEIGHT)) +pub fn parse_half_height_spec(spec: &str) -> Option { + spec.parse::().ok() } /// Parse a "re,im" spec (re/im decimal, parsed at /// full precision) into a view. Shared by @@ -213,16 +407,32 @@ pub fn parse_re_im_spec(spec: &str, bits: usize) -> Option<(Big, Big)> { /// offset/half_height ratio roughly constant, i.e. the target's on-screen /// position steady) but is shifted so it lands on exactly 1 at `t = 0` and /// exactly 0 at `t = 1`. +/// +/// With `d = log2(q)`, `g = q^t · (1 - q^(1-t)) / (1 - q)`, all in `Scale` +/// / `expm1` form: zooming in by more than f64's range, `q` (and `q^t`) +/// underflow, yet `g · (from - to)` must keep tracking the half-height. pub fn interpolate_view(from: &ViewState, to: &ViewState, t: f64) -> ViewState { - let q = to.half_height / from.half_height; - let half_height = from.half_height * q.powf(t); - let bits = precision_for(half_height); - let g = if (q - 1.0).abs() < 1e-12 { - 1.0 - t + let (l0, l1) = (from.half_height.log2(), to.half_height.log2()); + let d = l1 - l0; + let half_height = if t <= 0.0 { + from.half_height + } else if t >= 1.0 { + to.half_height } else { - (q.powf(t) - q) / (1.0 - q) + Scale::from_log2(l0 + t * d) + }; + let bits = precision_for(half_height); + let ln2 = core::f64::consts::LN_2; + let g_big = if d.abs() < 1e-12 { + big_from_f64(1.0 - t, bits) + } else if d < 0.0 { + // Zooming in: q^t may be far below f64's range, keep it as a Scale. + let f = ((1.0 - t) * d * ln2).exp_m1() / (d * ln2).exp_m1(); + Scale::from_log2(t * d).big_times(f, bits) + } else { + // Zooming out: g = (1 - q^(t-1)) / (1 - q^-1), every term bounded. + big_from_f64(((t - 1.0) * d * ln2).exp_m1() / (-d * ln2).exp_m1(), bits) }; - let g_big = big_from_f64(g, bits); let re0 = from.center_re.clone().with_precision(bits).value(); let im0 = from.center_im.clone().with_precision(bits).value(); let re1 = to.center_re.clone().with_precision(bits).value(); @@ -248,14 +458,10 @@ pub fn big_to_decimal_str(x: &Big, sig_digits: usize) -> String { } /// Precision (bits) needed to resolve the center at a given half-height. -pub fn precision_for(half_height: f64) -> usize { +pub fn precision_for(half_height: Scale) -> usize { // We need enough bits to distinguish points a pixel apart, i.e. roughly // log2(1 / half_height) significant bits, plus a guard margin. - let zoom_bits = if half_height > 0.0 && half_height.is_finite() { - (-half_height.log2()).ceil().max(0.0) as usize - } else { - 0 - }; + let zoom_bits = (-half_height.log2()).ceil().max(0.0) as usize; (zoom_bits + GUARD_BITS).clamp(53, MAX_PRECISION_BITS) } @@ -271,6 +477,10 @@ pub fn big_from_f64(x: f64, bits: usize) -> Big { mod tests { use super::*; + fn sc(x: f64) -> Scale { + Scale::from_f64(x) + } + fn re_im_f64(v: &ViewState) -> (f64, f64) { let re: f64 = v.center_re.to_decimal().value().to_f64().value(); let im: f64 = v.center_im.to_decimal().value().to_f64().value(); @@ -279,12 +489,12 @@ mod tests { #[test] fn interpolate_view_hits_exact_endpoints() { - let bits = precision_for(1.0); - let from = ViewState::with_center(big_from_f64(-0.5, bits), big_from_f64(0.0, bits), 1.5); + let bits = precision_for(sc(1.0)); + let from = ViewState::with_center(big_from_f64(-0.5, bits), big_from_f64(0.0, bits), sc(1.5)); let to = ViewState::with_center( - big_from_f64(-0.7515, precision_for(1e-20)), - big_from_f64(0.1013, precision_for(1e-20)), - 1e-20, + big_from_f64(-0.7515, precision_for(sc(1e-20))), + big_from_f64(0.1013, precision_for(sc(1e-20))), + sc(1e-20), ); let start = interpolate_view(&from, &to, 0.0); @@ -304,12 +514,12 @@ mod tests { /// instead stay roughly bounded throughout. #[test] fn interpolate_view_keeps_target_offset_bounded() { - let bits = precision_for(1.0); - let from = ViewState::with_center(big_from_f64(-0.5, bits), big_from_f64(0.0, bits), 1.5); + let bits = precision_for(sc(1.0)); + let from = ViewState::with_center(big_from_f64(-0.5, bits), big_from_f64(0.0, bits), sc(1.5)); let to = ViewState::with_center( - big_from_f64(-0.7515, precision_for(1e-20)), - big_from_f64(0.1013, precision_for(1e-20)), - 1e-20, + big_from_f64(-0.7515, precision_for(sc(1e-20))), + big_from_f64(0.1013, precision_for(sc(1e-20))), + sc(1e-20), ); let (to_re, to_im) = re_im_f64(&to); @@ -318,7 +528,7 @@ mod tests { let mid = interpolate_view(&from, &to, t); let (re, im) = re_im_f64(&mid); let offset = ((re - to_re).powi(2) + (im - to_im).powi(2)).sqrt(); - let ratio = offset / mid.half_height; + let ratio = offset / mid.half_height.to_f64(); assert!( ratio < 10.0, "t={t}: offset/half_height ratio {ratio} blew up (offset={offset}, half_height={})", @@ -326,4 +536,74 @@ mod tests { ); } } + + #[test] + fn scale_parse_display_round_trip() { + for s in ["1.25", "1e-20", "1.5e-20", "3.7e-4000", "1e-400", "9.99999e-310"] { + let a: Scale = s.parse().unwrap(); + let b: Scale = a.to_string().parse().unwrap(); + assert_eq!(a, b, "{s} -> {a}"); + } + assert_eq!("1.25".parse::().unwrap().to_f64(), 1.25); + assert_eq!(sc(1.5e-20).to_string(), "1.5e-20"); + let deep: Scale = "3.7e-4000".parse().unwrap(); + assert_eq!(deep.to_string(), "3.7e-4000"); + assert_eq!(format!("{deep:.2}"), "3.70e-4000"); + assert!((deep.log10() - (3.7f64.log10() - 4000.0)).abs() < 1e-9); + assert!("0".parse::().is_err()); + assert!("-1e-500".parse::().is_err()); + assert!("abc".parse::().is_err()); + } + + #[test] + fn scale_arithmetic() { + let a: Scale = "1e-1000".parse().unwrap(); + let b = a.mul_f64(0.25); + assert!((b.ratio(a) - 0.25).abs() < 1e-15); + assert!(b < a && a > b); + assert_eq!(a.mul_f64(3.0).mul_f64(1.0 / 3.0).exponent(), a.exponent()); + assert_eq!(sc(1.0).exponent(), 0); + assert_eq!(sc(0.75).exponent(), -1); + assert_eq!(sc(4.0).scaled_f64(-2), 1.0); + assert_eq!(a.scaled_f64(-a.exponent()), a.mul_f64(1.0).scaled_f64(-a.exponent())); + assert!((1.0..2.0).contains(&a.scaled_f64(-a.exponent()))); + assert_eq!(a.to_f64(), 0.0); + assert_eq!(Scale::MIN.mul_f64(0.5), Scale::MIN); + let p = precision_for(a); + assert!((3322 + 48..=3323 + 48).contains(&p), "{p}"); + } + + /// Past f64's range, the center must still land on the target at the + /// same geometric pace as the half-height. + #[test] + fn interpolate_view_past_f64_range() { + let from = ViewState::with_center( + big_from_f64(-0.5, 64), + big_from_f64(0.0, 64), + sc(1.5), + ); + let hh: Scale = "1e-1000".parse().unwrap(); + let bits = precision_for(hh); + let to = ViewState::with_center( + big_from_f64(-0.7515, bits), + big_from_f64(0.1013, bits), + hh, + ); + let end = interpolate_view(&from, &to, 1.0); + assert_eq!(end.half_height, hh); + assert_eq!(re_im_f64(&end), re_im_f64(&to)); + let mut prev = from.half_height; + for i in 1..20 { + let t = i as f64 / 20.0; + let mid = interpolate_view(&from, &to, t); + assert!(mid.half_height < prev); + prev = mid.half_height; + // Offset from the target, in units of the view's half-height. + let k = -mid.half_height.exponent() as isize; + let dre = ((&mid.center_re - &to.center_re) << k).to_f64().value(); + let dim = ((&mid.center_im - &to.center_im) << k).to_f64().value(); + let ratio = (dre * dre + dim * dim).sqrt() / mid.half_height.scaled_f64(k as i32); + assert!(ratio > 0.01 && ratio < 10.0, "t={t}: ratio {ratio}"); + } + } } diff --git a/src/worker.rs b/src/worker.rs index a427f1f..d6cb4f5 100644 --- a/src/worker.rs +++ b/src/worker.rs @@ -10,12 +10,12 @@ use std::sync::mpsc::{Receiver, Sender, TryRecvError, channel}; use std::thread; use crate::fractal::{FractalKind, RefOrbit, compute_reference, compute_set_reference}; -use crate::view::{Big, big_from_f64}; +use crate::view::{Big, Scale, big_from_f64}; pub struct RefRequest { pub center_re: Big, pub center_im: Big, - pub half_height: f64, + pub half_height: Scale, pub julia: bool, pub julia_c: (f64, f64), pub max_iter: u32, @@ -35,7 +35,7 @@ pub struct RefRequest { pub struct RefResult { pub center_re: Big, pub center_im: Big, - pub half_height: f64, + pub half_height: Scale, pub points: RefOrbit, /// The kind and morph `points` was computed with (echoed from the /// request).