diff --git a/CLAUDE.md b/CLAUDE.md index 01a20e3..66f8cd7 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -6,7 +6,8 @@ This file provides guidance to Claude Code (claude.ai/code) when working with co 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`), +reference orbit is computed on the CPU (arbitrary precision via `bignum::Big`: +`rug` natively, `malachite-float` on the web), 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 @@ -22,6 +23,7 @@ storage buffers, which the fragment shader needs for the reference orbit). ```sh cargo run --release # native, run (release matters: fractal math is hot) cargo test # reference-orbit math, share-link round-trip, WGSL validation +cargo test --features wasm # same, on the web build's malachite big-float backend cargo test --test shader_valid # just the WGSL parse/validate tests (naga, no GPU needed) cargo clippy cargo fmt # rustfmt.toml just pins edition = "2024" @@ -33,7 +35,7 @@ Web build (WebGPU): ```sh rustup target add wasm32-unknown-unknown cargo install wasm-bindgen-cli --version 0.2.128 # must match the wasm-bindgen crate version -./build-web.sh # -> ./dist +./build-web.sh # -> ./dist (builds with --features wasm) python3 -m http.server -d dist 8080 ``` @@ -105,8 +107,19 @@ runs out), rebase: `e ← y_n − X_0`, restart the reference index at 0. This i 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). The pixel scale (`half_height`) is a `Scale`, an f64 +- `src/bignum/` — `Big`, the arbitrary-precision binary float, with one + backend per library behind the same inherent API + operators: `rug` + (GMP/MPFR, default, fastest, can't target wasm32) and `malachite-float` + (pure Rust, the `wasm` feature, required for the web build; a + `compile_error!` enforces it). Cargo features are additive, so rug is a + non-wasm32 target dependency and the backend is picked by + `cfg(feature = "wasm")`. Both follow the precision rule: a result has the + larger operand precision, rounded to nearest; shifts are exact. Malachite's + zero has no precision, so its wrapper stores `prec` alongside. Any new + `Big` operation must be added to both backends (`cargo test` and + `cargo test --features wasm` run the same tests on each). +- `src/view.rs` — `ViewState`; center is arbitrary-precision `Big` + (re-exported as `view::Big`). 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 @@ -131,10 +144,10 @@ pixel is a handful of `f32` complex multiplies. `f32` pairs — that's the reference orbit the GPU perturbs from. At 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; + exists twice in this file (f64 + `Big`) and both must stay in sync; `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` + `Big` 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× @@ -190,7 +203,7 @@ pixel is a handful of `f32` complex multiplies. are `advance_delta_kind`/`fprime_kind`; `advance_delta`/`fprime` wrap them to blend two kinds during the kind-switch morph (`u.morph_from`, `u.morph_w`: each step is `(1-w)·f_kind + w·f_from`, mirrored on the CPU by - the `morph` argument of `compute_reference`, in both its f64 and `FBig` + the `morph` argument of `compute_reference`, in both its f64 and `Big` paths). The blend only exists in pipelines built with the `MORPH` override (part of `PipelineKey`, on while `morph_w > 0`); those also skip periodicity detection and the cardioid bypass. App side: `KindMorph` in diff --git a/Cargo.toml b/Cargo.toml index 39bb4b6..966727c 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -7,10 +7,12 @@ edition = "2024" default = ["gui"] # Windowed egui app. Without it only `--headless` rendering is built. gui = ["dep:eframe", "dep:egui"] +# Pure-Rust big floats for the web build (`rug` needs GMP/MPFR, which can't +# target wasm32). Required for wasm32, optional natively (to test that backend). +wasm = ["dep:malachite-float", "dep:malachite-base"] [dependencies] bytemuck = { version = "1.25.2", features = ["derive"] } -dashu-float = "0.6.0" eframe = { version = "0.36.2", default-features = false, features = ["wgpu", "default_fonts", "x11", "wayland", "accesskit"], optional = true } egui = { version = "0.36.2", optional = true } ecolor = { version = "0.36.2", features = ["bytemuck"] } @@ -18,11 +20,14 @@ wgpu = "30.0.1" glam = "0.33.8" log = "0.4.34" png = "0.18.1" +malachite-float = { version = "0.12", optional = true } +malachite-base = { version = "0.12", optional = true } [target.'cfg(not(target_arch = "wasm32"))'.dependencies] env_logger = "0.11.11" clap = { version = "4.5.51", features = ["derive"] } pollster = "1.0.1" +rug = { version = "1.30.0", default-features = false, features = ["float", "std"] } [target.'cfg(target_arch = "wasm32")'.dependencies] futures-channel = { version = "0.3.34", default-features = false, features = ["alloc", "std"] } @@ -39,7 +44,7 @@ opt-level = 3 # codegen-units = 1 debug = true -# Dev: keep our own crate debuggable, but optimize dependencies (dashu, wgpu, +# Dev: keep our own crate debuggable, but optimize dependencies (big floats, wgpu, # egui) so the explorer is actually interactive during development. [profile.dev] opt-level = 1 diff --git a/README.md b/README.md index d43d712..76a5288 100644 --- a/README.md +++ b/README.md @@ -3,7 +3,8 @@ A fast, interactive deep-zoom fractal explorer — Mandelbrot and Julia sets — built with **Rust + wgpu + egui + WGSL**. It zooms far 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 +computed on the CPU (arbitrary precision via `rug` natively, `malachite-float` +on the web), and every pixel is rendered on the GPU as a cheap `f32` delta from it, with **rebasing** to avoid glitches. Runs natively (Vulkan/Metal/DX12) and in the browser (WebGPU). @@ -27,6 +28,9 @@ The `f32` GPU tier reaches roughly **10³⁰× magnification** with sharp detail cargo run --release ``` +Native builds use `rug` (GMP/MPFR) for the reference orbit, which needs a C +toolchain and `m4` (on Windows, MSYS2). + ## Build & run — web (WebGPU) Requires the `wasm32-unknown-unknown` target and `wasm-bindgen-cli` (matching the @@ -37,6 +41,7 @@ rustup target add wasm32-unknown-unknown cargo install wasm-bindgen-cli --version 0.2.128 # once ./build-web.sh # outputs ./dist (index.html, .js, .wasm) + # (builds with --features wasm: pure-Rust big floats) python3 -m http.server -d dist 8080 # serve over http ``` @@ -60,7 +65,8 @@ qualifies. Deploy by serving the `dist/` directory as static files. ## How it works -- `src/view.rs` — view state. Center is arbitrary precision (`FBig`); the pixel +- `src/view.rs` — view state. Center is arbitrary precision (`Big`, from + `src/bignum/`); the pixel scale stays `f64` (even at 10³⁰× it is ~10⁻³³, within `f64` range). - `src/fractal/reference.rs` — high-precision reference orbit `Z_{n+1}=Z_n²+C`. - `src/shaders/mandelbrot.wgsl` — per-pixel perturbation `e_{n+1}=2·Z_n·e_n+e_n²+δc` @@ -84,6 +90,7 @@ reference computation to a Web Worker. ```sh cargo test +cargo test --features wasm # same suite on the web build's big-float backend ``` Covers the reference orbit (vs. a naive `f64` iteration, Mandelbrot and Julia) diff --git a/build-web.sh b/build-web.sh index 87e8f82..8ac13e3 100755 --- a/build-web.sh +++ b/build-web.sh @@ -8,7 +8,7 @@ export PATH="$HOME/.cargo/bin:$PATH" OUT="${1:-dist}" echo "==> cargo build (wasm32, release)" -cargo build --release --target wasm32-unknown-unknown +cargo build --release --target wasm32-unknown-unknown --features wasm echo "==> wasm-bindgen -> $OUT" mkdir -p "$OUT" diff --git a/src/app.rs b/src/app.rs index 604b4e1..7395477 100644 --- a/src/app.rs +++ b/src/app.rs @@ -1144,12 +1144,8 @@ impl FractalApp { fn drift_from(&self, key: &RequestKey) -> f64 { let hh = self.view.half_height; let k = -hh.exponent() as isize; - let dre = ((&self.view.center_re - &key.center_re) << k) - .to_f64() - .value(); - let dim = ((&self.view.center_im - &key.center_im) << k) - .to_f64() - .value(); + let dre = ((&self.view.center_re - &key.center_re) << k).to_f64(); + let dim = ((&self.view.center_im - &key.center_im) << k).to_f64(); (dre * dre + dim * dim).sqrt() / hh.scaled_f64(-hh.exponent()) } @@ -1198,12 +1194,8 @@ impl FractalApp { /// 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) << k) - .to_f64() - .value(); + let dre = ((&self.view.center_re - &self.ref_center_re) << k).to_f64(); + let dim = ((&self.view.center_im - &self.ref_center_im) << k).to_f64(); (dre, dim) } @@ -1494,8 +1486,8 @@ impl FractalApp { /// Buddhabrot mode doesn't support deep zoom (see `fractal::buddhabrot`). fn make_buddhabrot_uniforms(&self, aspect: f64) -> BuddhabrotUniforms { let center = [ - self.view.center_re.to_f64().value() as f32, - self.view.center_im.to_f64().value() as f32, + self.view.center_re.to_f64() as f32, + self.view.center_im.to_f64() as f32, ]; BuddhabrotUniforms { center, diff --git a/src/bignum/mod.rs b/src/bignum/mod.rs new file mode 100644 index 0000000..1ccc8ba --- /dev/null +++ b/src/bignum/mod.rs @@ -0,0 +1,225 @@ +//! Arbitrary-precision binary float [`Big`] (one coordinate of a deep-zoom +//! center, or of a reference orbit point), over one of two backends: +//! +//! - native (default): `rug` (GMP/MPFR), the fastest, but C code that can't +//! target `wasm32-unknown-unknown`; +//! - `wasm` feature: a pure-Rust library, required for the web build. +//! +//! Both expose the same inherent API plus the operators below, and follow +//! the same precision rule: an operation's result has the larger of its +//! operands' precisions (in bits), rounded to nearest. Shifts (`<<`/`>>` by +//! an `isize`) are exact multiplications by powers of two. + +#[cfg(not(feature = "wasm"))] +mod rug_backend; +#[cfg(not(feature = "wasm"))] +pub use rug_backend::Big; + +#[cfg(feature = "wasm")] +mod web_backend; +#[cfg(feature = "wasm")] +pub use web_backend::Big; +// Native `--features wasm` builds (to test the web backend) still link rug. +#[cfg(all(feature = "wasm", not(target_arch = "wasm32")))] +use rug as _; + +#[cfg(all(target_arch = "wasm32", not(feature = "wasm")))] +compile_error!("the web build needs `--features wasm` (rug can't target wasm32)"); + +use core::ops::{Add, Mul, Neg, Shl, Shr, Sub}; + +/// Implements `$Trait` for every owned/borrowed combination of `Big` +/// operands, through the backend's by-reference `$imp`. +macro_rules! forward_binop { + ($Trait:ident, $method:ident, $imp:ident) => { + impl $Trait<&Big> for &Big { + type Output = Big; + fn $method(self, rhs: &Big) -> Big { + self.$imp(rhs) + } + } + impl $Trait for &Big { + type Output = Big; + fn $method(self, rhs: Big) -> Big { + self.$imp(&rhs) + } + } + impl $Trait<&Big> for Big { + type Output = Big; + fn $method(self, rhs: &Big) -> Big { + (&self).$imp(rhs) + } + } + impl $Trait for Big { + type Output = Big; + fn $method(self, rhs: Big) -> Big { + (&self).$imp(&rhs) + } + } + }; +} + +forward_binop!(Add, add, add_ref); +forward_binop!(Sub, sub, sub_ref); +forward_binop!(Mul, mul, mul_ref); + +impl Neg for Big { + type Output = Big; + fn neg(self) -> Big { + self.negated() + } +} + +impl Neg for &Big { + type Output = Big; + fn neg(self) -> Big { + self.clone().negated() + } +} + +impl Shl for Big { + type Output = Big; + fn shl(self, k: isize) -> Big { + self.mul_pow2(k) + } +} + +impl Shl for &Big { + type Output = Big; + fn shl(self, k: isize) -> Big { + self.clone().mul_pow2(k) + } +} + +impl Shr for Big { + type Output = Big; + fn shr(self, k: isize) -> Big { + self.mul_pow2(-k) + } +} + +impl Shr for &Big { + type Output = Big; + fn shr(self, k: isize) -> Big { + self.clone().mul_pow2(-k) + } +} + +/// Significant decimal digits of a value, as a backend's +/// `to_decimal_parts` returns them: the value is `±0.digits × 10^exp10`, +/// `digits` has no trailing zeros, and it is empty for zero. +pub struct DecimalParts { + pub negative: bool, + pub digits: String, + pub exp10: isize, +} + +impl DecimalParts { + /// Build from a significand digit string (maybe with trailing zeros) + /// whose value is `0.digits × 10^exp10`. + fn new(negative: bool, digits: &str, exp10: isize) -> Self { + let digits = digits.trim_end_matches('0').to_string(); + Self { + negative: negative && !digits.is_empty(), + digits, + exp10, + } + } +} + +impl Big { + /// Decimal string with `sig_digits` significant digits, in plain + /// positional notation (`-0.000123`, `4500`), trailing zeros trimmed. + pub fn to_decimal_string(&self, sig_digits: usize) -> String { + let DecimalParts { + negative, + digits, + exp10, + } = self.to_decimal_parts(sig_digits.max(1)); + if digits.is_empty() { + return "0".to_string(); + } + let sign = if negative { "-" } else { "" }; + let len = digits.len() as isize; + if exp10 <= 0 { + let zeros = "0".repeat((-exp10) as usize); + format!("{sign}0.{zeros}{digits}") + } else if exp10 >= len { + let zeros = "0".repeat((exp10 - len) as usize); + format!("{sign}{digits}{zeros}") + } else { + let (int, frac) = digits.split_at(exp10 as usize); + format!("{sign}{int}.{frac}") + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + fn big(x: f64) -> Big { + Big::from_f64(x, 200) + } + + #[test] + fn decimal_string_is_positional() { + assert_eq!( + big(-0.7515).to_decimal_string(20), + "-0.75149999999999994582" + ); + assert_eq!(big(-0.7515).to_decimal_string(5), "-0.7515"); + assert_eq!(big(0.0).to_decimal_string(5), "0"); + assert_eq!(big(1.0).to_decimal_string(5), "1"); + assert_eq!(big(123.456).to_decimal_string(5), "123.46"); + assert_eq!(big(1e20).to_decimal_string(5), "100000000000000000000"); + assert_eq!(big(0.1).to_decimal_string(5), "0.1"); + let tiny = big(1.0) >> 200; + let s = tiny.to_decimal_string(5); + assert_eq!(s, format!("0.{}6223", "0".repeat(60))); + } + + #[test] + fn decimal_round_trip() { + let s = "-0.743643887037158704752191506114774"; + let x = Big::from_decimal_str(s, 200).unwrap(); + assert_eq!(x.to_decimal_string(33), s); + assert_eq!( + Big::from_decimal_str("1.5e-20", 64).unwrap().to_f64(), + 1.5e-20 + ); + assert!(Big::from_decimal_str("abc", 64).is_none()); + } + + #[test] + fn precision_is_max_of_operands() { + let a = Big::from_f64(1.0, 100); + let b = Big::from_f64(3.0, 200); + assert_eq!((&a + &b).precision(), 200); + assert_eq!((&a * &b).precision(), 200); + assert_eq!((&b - &a).precision(), 200); + assert_eq!(a.clone().with_precision(300).precision(), 300); + } + + #[test] + fn exact_shifts_and_log2() { + let x = Big::from_f64(0.75, 64) >> 5000; + assert_eq!(x.log2_floor(), Some(-5001)); + assert_eq!((x << 5000).to_f64(), 0.75); + assert_eq!(Big::zero(64).log2_floor(), None); + assert!(Big::from_f64(-1.0, 64).is_negative()); + assert!(!Big::zero(64).is_negative()); + } + + #[test] + fn transcendentals() { + let x = Big::from_f64(0.5, 128); + assert!((x.ln().to_f64() - 0.5f64.ln()).abs() < 1e-15); + assert!((x.exp().to_f64() - 0.5f64.exp()).abs() < 1e-15); + let (s, c) = x.sin_cos(); + assert!((s.to_f64() - 0.5f64.sin()).abs() < 1e-15); + assert!((c.to_f64() - 0.5f64.cos()).abs() < 1e-15); + let y = Big::from_f64(-0.3, 128); + assert!((y.atan2(&x).to_f64() - (-0.3f64).atan2(0.5)).abs() < 1e-15); + } +} diff --git a/src/bignum/rug_backend.rs b/src/bignum/rug_backend.rs new file mode 100644 index 0000000..bffe793 --- /dev/null +++ b/src/bignum/rug_backend.rs @@ -0,0 +1,128 @@ +//! [`Big`] over `rug::Float` (MPFR), for native builds. + +use core::cmp::Ordering; + +use rug::{Assign, Float}; + +use super::DecimalParts; + +/// Arbitrary-precision binary float. See the module docs. +#[derive(Clone, Debug, PartialEq)] +pub struct Big(Float); + +/// MPFR precision for `bits`, clamped to what it accepts. +fn prec(bits: usize) -> u32 { + (bits as u32).clamp(rug::float::prec_min(), rug::float::prec_max()) +} + +impl Big { + /// `x` at `bits` of precision (exact when `bits >= 53`). Non-finite `x` + /// reads as 0. + pub fn from_f64(x: f64, bits: usize) -> Self { + let x = if x.is_finite() { x } else { 0.0 }; + Self(Float::with_val(prec(bits), x)) + } + + pub fn zero(bits: usize) -> Self { + Self(Float::new(prec(bits))) + } + + /// Parse a decimal (`-0.75`, `1.5e-20`, any number of digits) at `bits` + /// (at least 53) of precision. + pub fn from_decimal_str(s: &str, bits: usize) -> Option { + let parsed = Float::parse(s.trim()).ok()?; + let x = Float::with_val(prec(bits.max(53)), parsed); + x.is_finite().then_some(Self(x)) + } + + pub fn precision(&self) -> usize { + self.0.prec() as usize + } + + /// Change the precision, rounding to nearest if it shrinks. + pub fn with_precision(mut self, bits: usize) -> Self { + self.0.set_prec(prec(bits)); + self + } + + pub fn to_f64(&self) -> f64 { + self.0.to_f64() + } + + pub fn is_zero(&self) -> bool { + self.0.is_zero() + } + + pub fn is_negative(&self) -> bool { + self.0.cmp0() == Some(Ordering::Less) + } + + pub fn abs(self) -> Self { + Self(self.0.abs()) + } + + pub fn sqr(&self) -> Self { + Self(Float::with_val(self.0.prec(), self.0.square_ref())) + } + + /// `floor(log2|x|)`, or `None` for zero. Exact, at any exponent. + pub fn log2_floor(&self) -> Option { + // MPFR normalizes the significand to [0.5, 1). + self.0.get_exp().map(|e| e as isize - 1) + } + + pub fn ln(&self) -> Self { + Self(Float::with_val(self.0.prec(), self.0.ln_ref())) + } + + pub fn exp(&self) -> Self { + Self(Float::with_val(self.0.prec(), self.0.exp_ref())) + } + + /// `atan2(self, x)`: the angle of `(x, self)`. + pub fn atan2(&self, x: &Big) -> Self { + let p = self.0.prec().max(x.0.prec()); + Self(Float::with_val(p, self.0.atan2_ref(&x.0))) + } + + pub fn sin_cos(&self) -> (Self, Self) { + let p = self.0.prec(); + let (mut s, mut c) = (Float::new(p), Float::new(p)); + (&mut s, &mut c).assign(self.0.sin_cos_ref()); + (Self(s), Self(c)) + } + + /// `sig` significant decimal digits (rounded to nearest). + pub fn to_decimal_parts(&self, sig: usize) -> DecimalParts { + if self.0.is_zero() { + return DecimalParts::new(false, "", 0); + } + let (negative, digits, exp) = self.0.to_sign_string_exp(10, Some(sig)); + DecimalParts::new(negative, &digits, exp.unwrap_or(0) as isize) + } + + pub(super) fn add_ref(&self, rhs: &Big) -> Big { + let p = self.0.prec().max(rhs.0.prec()); + Big(Float::with_val(p, &self.0 + &rhs.0)) + } + + pub(super) fn sub_ref(&self, rhs: &Big) -> Big { + let p = self.0.prec().max(rhs.0.prec()); + Big(Float::with_val(p, &self.0 - &rhs.0)) + } + + pub(super) fn mul_ref(&self, rhs: &Big) -> Big { + let p = self.0.prec().max(rhs.0.prec()); + Big(Float::with_val(p, &self.0 * &rhs.0)) + } + + pub(super) fn negated(self) -> Big { + Big(-self.0) + } + + /// `self · 2^k`, exact. `|k|` stays far below `i32::MAX` (precision and + /// zoom depth are capped around 2^20 bits). + pub(super) fn mul_pow2(self, k: isize) -> Big { + Big(self.0 << k as i32) + } +} diff --git a/src/bignum/web_backend.rs b/src/bignum/web_backend.rs new file mode 100644 index 0000000..b86b403 --- /dev/null +++ b/src/bignum/web_backend.rs @@ -0,0 +1,173 @@ +//! [`Big`] over `malachite-float`, pure Rust, for the web build. + +use malachite_base::num::basic::traits::Zero; +use malachite_base::num::conversion::string::options::ToSciOptions; +use malachite_base::num::conversion::traits::{RoundingFrom, ToSci}; +use malachite_base::rounding_modes::RoundingMode::Nearest; +use malachite_float::Float; + +use super::DecimalParts; + +/// Arbitrary-precision binary float. See the module docs. +/// +/// The precision is stored next to the value: malachite's zero has none, +/// but "zero at `p` bits" must still lift a sum with it to `p` bits. +#[derive(Clone, Debug)] +pub struct Big { + f: Float, + prec: u64, +} + +impl PartialEq for Big { + fn eq(&self, other: &Self) -> bool { + self.f == other.f + } +} + +fn prec(bits: usize) -> u64 { + (bits as u64).max(1) +} + +impl Big { + fn new(f: Float, prec: u64) -> Self { + // Only finite values reach here; anything else (overflow) is a bug. + debug_assert!(f.is_finite(), "non-finite Big: {f}"); + Self { f, prec } + } + + /// `x` at `bits` of precision (exact when `bits >= 53`). Non-finite `x` + /// reads as 0. + pub fn from_f64(x: f64, bits: usize) -> Self { + let p = prec(bits); + if !x.is_finite() { + return Self::zero(bits); + } + Self::new(Float::from_primitive_float_prec(x, p).0, p) + } + + pub fn zero(bits: usize) -> Self { + Self::new(Float::ZERO, prec(bits)) + } + + /// Parse a decimal (`-0.75`, `1.5e-20`, any number of digits) at `bits` + /// (at least 53) of precision. + pub fn from_decimal_str(s: &str, bits: usize) -> Option { + let p = prec(bits.max(53)); + let (f, _) = Float::from_sci_string_prec(s.trim(), p)?; + f.is_finite().then(|| Self::new(f, p)) + } + + pub fn precision(&self) -> usize { + self.prec as usize + } + + /// Change the precision, rounding to nearest if it shrinks. + pub fn with_precision(mut self, bits: usize) -> Self { + let p = prec(bits); + if self.f.get_prec().is_some() { + self.f.set_prec(p); + } + self.prec = p; + self + } + + pub fn to_f64(&self) -> f64 { + f64::rounding_from(&self.f, Nearest).0 + } + + pub fn is_zero(&self) -> bool { + self.f.is_zero() + } + + pub fn is_negative(&self) -> bool { + self.f.is_sign_negative() && !self.f.is_zero() + } + + pub fn abs(self) -> Self { + if self.f.is_sign_negative() { + self.negated() + } else { + self + } + } + + pub fn sqr(&self) -> Self { + Self::new(self.f.square_prec_ref(self.prec).0, self.prec) + } + + /// `floor(log2|x|)`, or `None` for zero. Exact, at any exponent. + pub fn log2_floor(&self) -> Option { + // The significand is normalized to [0.5, 1). + self.f.get_exponent().map(|e| e as isize - 1) + } + + pub fn ln(&self) -> Self { + Self::new(self.f.ln_prec_ref(self.prec).0, self.prec) + } + + pub fn exp(&self) -> Self { + Self::new(self.f.exp_prec_ref(self.prec).0, self.prec) + } + + /// `atan2(self, x)`: the angle of `(x, self)`. + pub fn atan2(&self, x: &Big) -> Self { + let p = self.prec.max(x.prec); + Self::new(self.f.atan2_prec_ref_ref(&x.f, p).0, p) + } + + pub fn sin_cos(&self) -> (Self, Self) { + let (s, c, _, _) = self.f.sin_cos_prec_ref(self.prec); + (Self::new(s, self.prec), Self::new(c, self.prec)) + } + + /// `sig` significant decimal digits (rounded to nearest). + pub fn to_decimal_parts(&self, sig: usize) -> DecimalParts { + if self.f.is_zero() { + return DecimalParts::new(false, "", 0); + } + let mut options = ToSciOptions::default(); + options.set_precision(sig as u64); + let s = self.f.to_sci_with_options(options).to_string(); + // `[-]int[.frac][e±N]`, positional or scientific depending on size. + let (negative, s) = match s.strip_prefix('-') { + Some(rest) => (true, rest), + None => (false, s.as_str()), + }; + let (mantissa, exp) = match s.split_once(['e', 'E']) { + Some((m, e)) => (m, e.trim_start_matches('+').parse::().unwrap_or(0)), + None => (s, 0), + }; + let (int, frac) = mantissa.split_once('.').unwrap_or((mantissa, "")); + let all = format!("{int}{frac}"); + let digits = all.trim_start_matches('0'); + let leading_zeros = (all.len() - digits.len()) as isize; + // 0.digits · 10^exp10 = int.frac · 10^exp. + let exp10 = int.len() as isize - leading_zeros + exp; + DecimalParts::new(negative, digits, exp10) + } + + pub(super) fn add_ref(&self, rhs: &Big) -> Big { + let p = self.prec.max(rhs.prec); + Big::new(self.f.add_prec_ref_ref(&rhs.f, p).0, p) + } + + pub(super) fn sub_ref(&self, rhs: &Big) -> Big { + let p = self.prec.max(rhs.prec); + Big::new(self.f.sub_prec_ref_ref(&rhs.f, p).0, p) + } + + pub(super) fn mul_ref(&self, rhs: &Big) -> Big { + let p = self.prec.max(rhs.prec); + Big::new(self.f.mul_prec_ref_ref(&rhs.f, p).0, p) + } + + pub(super) fn negated(self) -> Big { + Big::new(-self.f, self.prec) + } + + /// `self · 2^k`, exact (malachite's exponent range is ±2^30, far past + /// what zoom depth needs). + pub(super) fn mul_pow2(self, k: isize) -> Big { + Big::new(self.f << k as i64, self.prec) + } +} diff --git a/src/fractal/reference.rs b/src/fractal/reference.rs index 0bbe4fd..0d8be39 100644 --- a/src/fractal/reference.rs +++ b/src/fractal/reference.rs @@ -1,7 +1,7 @@ //! High-precision reference-orbit computation for perturbation rendering. //! //! We iterate the fractal's formula `Z_{n+1} = f(Z_n, C)` at high precision -//! (`dashu-float`), storing each `Z_n` as an `f32` pair. Every pixel is then +//! ([`Big`]: `rug` natively, pure Rust on the web), storing each `Z_n` as an `f32` pair. Every pixel is then //! rendered on the GPU as a small `f32` delta from this orbit — that is what //! makes deep zoom cheap. See `shaders/mandelbrot.wgsl` for the delta side; the //! delta formula there must match the orbit formula here. @@ -24,7 +24,7 @@ use crate::view::{Big, big_from_f64}; const REFERENCE_ESCAPE_SQ: f64 = 1.0e10; /// Up to this working precision (bits) the orbit is iterated in plain `f64` -/// instead of `FBig` — orders of magnitude faster, which matters most on the +/// instead of `Big` — orders of magnitude faster, which matters most on the /// web (where the reference is computed inline on the UI thread). /// /// `precision_for` asks for `zoom_bits + 48` guard bits, but the GPU only @@ -96,29 +96,21 @@ impl<'a> IntoIterator for &'a RefOrbit { } } -/// `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 (lr, li) = (zr.log2_floor(), zi.log2_floor()); 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; + let mr = (zr.clone() << -k).to_f64() as f32; + let mi = (zi.clone() << -k).to_f64() 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]), + _ => orbit.push([zr.to_f64() as f32, zi.to_f64() as f32]), } } @@ -147,23 +139,17 @@ pub fn compute_reference( let morph = morph.filter(|&(_, w)| w != 0.0); if precision <= F64_MAX_PRECISION { let k = StepConstsF64 { - c: (c_re.to_f64().value(), c_im.to_f64().value()), + c: (c_re.to_f64(), c_im.to_f64()), p: phoenix_p, l: lambda_l, cpow: complex_power, power, }; - return compute_reference_f64( - (z0_re.to_f64().value(), z0_im.to_f64().value()), - max_iter, - kind, - &k, - morph, - ); + return compute_reference_f64((z0_re.to_f64(), z0_im.to_f64()), max_iter, kind, &k, morph); } let k = StepConsts { - cr: c_re.clone().with_precision(precision).value(), - ci: c_im.clone().with_precision(precision).value(), + cr: c_re.clone().with_precision(precision), + ci: c_im.clone().with_precision(precision), pr: big_from_f64(phoenix_p.0, precision), pi: big_from_f64(phoenix_p.1, precision), lr: big_from_f64(lambda_l.0, precision), @@ -292,7 +278,7 @@ struct StepConsts { precision: usize, } -/// [`compute_reference`] at arbitrary precision (`FBig`), for deep views. +/// [`compute_reference`] at arbitrary precision ([`Big`]), for deep views. fn compute_reference_big( z0_re: &Big, z0_im: &Big, @@ -304,11 +290,11 @@ fn compute_reference_big( let precision = k.precision; let morph = morph.map(|(from, w)| (from, big_from_f64(w, precision))); - let mut zr = z0_re.clone().with_precision(precision).value(); - let mut zi = z0_im.clone().with_precision(precision).value(); + let mut zr = z0_re.clone().with_precision(precision); + let mut zi = z0_im.clone().with_precision(precision); // Previous iterate, for the Phoenix two-term recurrence (Y_{-1} = 0). - let mut zr_prev = big_zero(precision); - let mut zi_prev = big_zero(precision); + let mut zr_prev = Big::zero(precision); + let mut zi_prev = Big::zero(precision); let mut points = RefOrbit::with_capacity(max_iter as usize + 1); @@ -332,8 +318,8 @@ fn compute_reference_big( // Shift the previous iterate (only the Phoenix arm reads it). zr_prev = zr; zi_prev = zi; - zr = new_zr.with_precision(precision).value(); - zi = new_zi.with_precision(precision).value(); + zr = new_zr.with_precision(precision); + zi = new_zi.with_precision(precision); } points @@ -361,7 +347,7 @@ fn step( FractalKind::BurningShip => { // (|zr| + i|zi|)^2 = (zr^2 - zi^2) + 2|zr zi| i. let re = &zr.sqr() - &zi.sqr() + cr; - let im = big_abs((zr * zi) << 1) + ci; + let im = ((zr * zi) << 1).abs() + ci; (re, im) } FractalKind::Tricorn => { @@ -376,14 +362,14 @@ fn step( } FractalKind::Celtic => { // |Re(z^2)| + i·Im(z^2): abs the real output of the square. - let re = big_abs(&zr.sqr() - &zi.sqr()) + cr; + let re = (&zr.sqr() - &zi.sqr()).abs() + cr; let im = ((zr * zi) << 1) + ci; (re, im) } 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 < Big::ZERO { + let im = if zi.is_negative() { ci + ((zr * zi) << 1) } else { ci - ((zr * zi) << 1) @@ -392,8 +378,8 @@ fn step( } FractalKind::Buffalo => { // |Re(z^2)| - |Im(z^2)| i: abs both outputs. - let re = big_abs(&zr.sqr() - &zi.sqr()) + cr; - let im = ci - big_abs((zr * zi) << 1); + let re = (&zr.sqr() - &zi.sqr()).abs() + cr; + let im = ci - ((zr * zi) << 1).abs(); (re, im) } FractalKind::Phoenix => { @@ -406,7 +392,7 @@ fn step( } FractalKind::Lambda => { // λ·z(1 - z) + c: logistic map plus the usual additive `c`. - let re2 = 1 - zr; + let re2 = Big::from_f64(1.0, k.precision) - zr; let im2 = -zi; let lzr = &k.lr * zr - &k.li * zi; let lzi = &k.lr * zi + &k.li * zr; @@ -421,36 +407,20 @@ fn step( } } -fn big_zero(precision: usize) -> Big { - Big::from(0i32).with_precision(precision).value() -} - -/// 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 < Big::ZERO { -x } else { x } -} - /// `(zr + i zi)^power` by repeated complex multiply at `precision` bits. fn complex_pow(zr: &Big, zi: &Big, power: u32, precision: usize) -> (Big, Big) { - let mut rr = Big::from(1i32).with_precision(precision).value(); - let mut ri = big_zero(precision); + let mut rr = Big::from_f64(1.0, precision); + let mut ri = Big::zero(precision); for _ in 0..power { // (rr + i ri)(zr + i zi) = (rr zr - ri zi) + (rr zi + ri zr) i. - let nr = (&rr * zr - &ri * zi).with_precision(precision).value(); - let ni = (&rr * zi + &ri * zr).with_precision(precision).value(); + let nr = (&rr * zr - &ri * zi).with_precision(precision); + let ni = (&rr * zi + &ri * zr).with_precision(precision); rr = nr; ri = ni; } (rr, ri) } -/// `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.repr().significand().is_zero() -} - /// `(zr + i zi)^(pr + i pi)` for a complex exponent, via the principal branch /// `z^p = exp(p·ln z)` where `ln z = ln|z| + i·arg(z)`. Used by /// `ComplexMultibrot`; must be kept in sync with the shader's `cpow`. @@ -458,14 +428,14 @@ fn is_big_zero(x: &Big) -> bool { /// panic; this is the correct limit for the `Re(p) > 0` region the UI /// exposes). fn complex_pow_complex(zr: &Big, zi: &Big, pr: &Big, pi: &Big, precision: usize) -> (Big, Big) { - if is_big_zero(zr) && is_big_zero(zi) { - return (big_zero(precision), big_zero(precision)); + if zr.is_zero() && zi.is_zero() { + return (Big::zero(precision), Big::zero(precision)); } let r2 = &zr.sqr() + &zi.sqr(); let ln_r = r2.ln() >> 1; // 0.5 * ln(r2) = ln(sqrt(r2)); exact halving. let theta = zi.atan2(zr); - let exp_re = (pr * &ln_r - pi * &theta).with_precision(precision).value(); - let exp_im = (pr * &theta + pi * &ln_r).with_precision(precision).value(); + let exp_re = (pr * &ln_r - pi * &theta).with_precision(precision); + let exp_im = (pr * &theta + pi * &ln_r).with_precision(precision); let mag = exp_re.exp(); let (sin_a, cos_a) = exp_im.sin_cos(); (&mag * &cos_a, &mag * &sin_a) @@ -486,7 +456,7 @@ pub fn compute_set_reference( complex_power: (f64, f64), morph: Option<(FractalKind, f64)>, ) -> RefOrbit { - let zero = big_zero(precision); + let zero = Big::zero(precision); compute_reference( &zero, &zero, @@ -511,19 +481,19 @@ mod tests { /// `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); + let tiny = Big::from_f64(1.0, 64) >> 5000; + assert!(!tiny.is_zero()); + assert!(Big::zero(64).is_zero()); + assert_eq!((-tiny.clone()).abs(), tiny); + assert_eq!(tiny.clone().abs(), tiny); } /// The high-precision reference must agree with a plain f64 iteration for a /// shallow point (where f64 is accurate). #[test] fn reference_matches_naive_f64() { - let cr = Big::try_from(-0.75_f64).unwrap(); - let ci = Big::try_from(0.1_f64).unwrap(); + let cr = Big::from_f64(-0.75, 53); + let ci = Big::from_f64(0.1, 53); let points = compute_set_reference( &cr, &ci, @@ -656,8 +626,8 @@ mod tests { /// A point inside the main cardioid never escapes: full-length orbit. #[test] fn interior_orbit_runs_full_length() { - let cr = Big::try_from(-0.2_f64).unwrap(); - let ci = Big::try_from(0.0_f64).unwrap(); + let cr = Big::from_f64(-0.2, 53); + let ci = Big::from_f64(0.0, 53); let points = compute_set_reference( &cr, &ci, @@ -676,8 +646,8 @@ mod tests { /// Burning Ship reference matches a naive f64 iteration of the same formula. #[test] fn burning_ship_reference_matches_naive_f64() { - let cr = Big::try_from(-1.75_f64).unwrap(); - let ci = Big::try_from(-0.03_f64).unwrap(); + let cr = Big::from_f64(-1.75, 53); + let ci = Big::from_f64(-0.03, 53); let points = compute_set_reference( &cr, &ci, @@ -707,8 +677,8 @@ mod tests { /// Multibrot (power 3) reference matches a naive f64 cube iteration. #[test] fn multibrot3_reference_matches_naive_f64() { - let cr = Big::try_from(0.3_f64).unwrap(); - let ci = Big::try_from(0.2_f64).unwrap(); + let cr = Big::from_f64(0.3, 53); + let ci = Big::from_f64(0.2, 53); let points = compute_set_reference( &cr, &ci, @@ -740,10 +710,10 @@ mod tests { /// Julia orbit (fixed c, z0 = center) matches a naive f64 iteration. #[test] fn julia_reference_matches_naive_f64() { - let z0_re = Big::try_from(0.15_f64).unwrap(); - let z0_im = Big::try_from(-0.1_f64).unwrap(); - let c_re = Big::try_from(-0.8_f64).unwrap(); - let c_im = Big::try_from(0.156_f64).unwrap(); + let z0_re = Big::from_f64(0.15, 53); + let z0_im = Big::from_f64(-0.1, 53); + let c_re = Big::from_f64(-0.8, 53); + let c_im = Big::from_f64(0.156, 53); let points = compute_reference( &z0_re, &z0_im, @@ -778,10 +748,10 @@ mod tests { let (lr, li) = (-0.5_f64, 0.2_f64); let (cr, ci) = (0.1_f64, -0.3_f64); let points = compute_reference( - &Big::try_from(0.2_f64).unwrap(), - &Big::try_from(0.1_f64).unwrap(), - &Big::try_from(cr).unwrap(), - &Big::try_from(ci).unwrap(), + &Big::from_f64(0.2, 53), + &Big::from_f64(0.1, 53), + &Big::from_f64(cr, 53), + &Big::from_f64(ci, 53), 60, 200, FractalKind::Lambda, @@ -807,8 +777,8 @@ mod tests { /// Celtic reference matches a naive f64 iteration: real = |x^2 - y^2| + cr. #[test] fn celtic_reference_matches_naive_f64() { - let cr = Big::try_from(-0.6_f64).unwrap(); - let ci = Big::try_from(0.4_f64).unwrap(); + let cr = Big::from_f64(-0.6, 53); + let ci = Big::from_f64(0.4, 53); let points = compute_set_reference( &cr, &ci, @@ -839,8 +809,8 @@ mod tests { /// real = x^2 - y^2 + cr, imag = -2·x·|y| + ci. #[test] fn perpendicular_reference_matches_naive_f64() { - let cr = Big::try_from(-0.7_f64).unwrap(); - let ci = Big::try_from(-0.2_f64).unwrap(); + let cr = Big::from_f64(-0.7, 53); + let ci = Big::from_f64(-0.2, 53); let points = compute_set_reference( &cr, &ci, @@ -871,8 +841,8 @@ mod tests { /// real = |x^2 - y^2| + cr, imag = -|2·x·y| + ci. #[test] fn buffalo_reference_matches_naive_f64() { - let cr = Big::try_from(-1.2_f64).unwrap(); - let ci = Big::try_from(-0.35_f64).unwrap(); + let cr = Big::from_f64(-1.2, 53); + let ci = Big::from_f64(-0.35, 53); let points = compute_set_reference( &cr, &ci, @@ -903,8 +873,8 @@ mod tests { /// `z_{n+1} = z_n^2 + c + p·z_{n-1}` (z_0 = 0, z_{-1} = 0). #[test] fn phoenix_reference_matches_naive_f64() { - let cr = Big::try_from(0.5667_f64).unwrap(); - let ci = Big::try_from(0.0_f64).unwrap(); + let cr = Big::from_f64(0.5667, 53); + let ci = Big::from_f64(0.0, 53); let p = (-0.5_f64, 0.0_f64); let points = compute_set_reference( &cr, @@ -942,8 +912,8 @@ mod tests { /// iteration of `z^p = exp(p·ln z)`. #[test] fn complex_multibrot_reference_matches_naive_f64() { - let cr = Big::try_from(0.1_f64).unwrap(); - let ci = Big::try_from(-0.2_f64).unwrap(); + let cr = Big::from_f64(0.1, 53); + let ci = Big::from_f64(-0.2, 53); let power = (2.5_f64, 0.3_f64); let points = compute_set_reference( &cr, @@ -987,8 +957,8 @@ mod tests { 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(), + &Big::from_f64(cr, 53), + &Big::from_f64(ci, 53), 60, 200, kind, diff --git a/src/main.rs b/src/main.rs index e461e8c..24947e1 100644 --- a/src/main.rs +++ b/src/main.rs @@ -11,6 +11,7 @@ // and calls the wasm `main`, which boots eframe onto the page's . mod app; +mod bignum; mod camera; mod fractal; mod lights; diff --git a/src/view.rs b/src/view.rs index 6daae08..e5f0df6 100644 --- a/src/view.rs +++ b/src/view.rs @@ -1,6 +1,6 @@ //! Camera / view state over the complex plane. //! -//! The center is stored in arbitrary precision (`FBig`) — this is what lets us +//! The center is stored in arbitrary precision ([`Big`]) — this is what lets us //! 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 @@ -10,11 +10,7 @@ use core::str::FromStr; -use dashu_float::round::mode::HalfAway; -use dashu_float::{DBig, FBig}; - -/// Arbitrary-precision binary float (base 2, round-half-away). One coordinate. -pub type Big = FBig; +pub use crate::bignum::Big; /// Half-height (complex units) of the default view; also the zoom-1 reference. pub const DEFAULT_HALF_HEIGHT: f64 = 1.25; @@ -194,28 +190,37 @@ impl core::fmt::Display for Scale { }; } // Out of f64's range: round the exact decimal expansion instead. - let sig = f.precision().map_or(17, |p| p + 1); - let dec = self - .to_big() - .to_decimal() - .value() - .with_precision(sig) - .value(); - let 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(), + let big = self.to_big(); + let sci = |sig: usize, pad: Option| { + let parts = big.to_decimal_parts(sig); + let digits = if parts.digits.is_empty() { + "0" + } else { + &parts.digits + }; + // value = 0.digits · 10^exp10; move the point after the first digit. + let exp10 = parts.exp10 - 1; + let (head, tail) = digits.split_at(1); + let tail = match pad { + Some(p) => format!("{tail:0 tail.to_string(), + }; + if tail.is_empty() { + format!("{head}e{exp10}") + } else { + format!("{head}.{tail}e{exp10}") + } }; - if tail.is_empty() { - write!(f, "{head}e{exp10}") - } else { - write!(f, "{head}.{tail}e{exp10}") + match f.precision() { + Some(p) => f.write_str(&sci(p + 1, Some(p))), + None => { + // Like `{:e}` on f64: the shortest string that parses back. + let s = (1..17) + .map(|sig| sci(sig, None)) + .find(|s| s.parse::() == Ok(*self)) + .unwrap_or_else(|| sci(17, None)); + f.write_str(&s) + } } } } @@ -237,18 +242,12 @@ impl FromStr for Scale { }; } // 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 { + let bin = Big::from_decimal_str(s, 64).ok_or(())?; + if bin.is_negative() { 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 top = bin.log2_floor().ok_or(())?; + let m = (bin >> top).to_f64(); let e = top.clamp(i32::MIN as isize, i32::MAX as isize) as i32; Ok(Self::from_parts(m, e)) } @@ -312,10 +311,10 @@ impl ViewState { pub fn sync_precision(&mut self) { let bits = self.precision_bits(); if self.center_re.precision() < bits { - self.center_re = self.center_re.clone().with_precision(bits).value(); + self.center_re = self.center_re.clone().with_precision(bits); } if self.center_im.precision() < bits { - self.center_im = self.center_im.clone().with_precision(bits).value(); + self.center_im = self.center_im.clone().with_precision(bits); } } @@ -360,8 +359,7 @@ impl ViewState { /// Parse a decimal string (any number of digits) losslessly into a `Big` with at /// least `bits` of precision. Used for share links and debug view specs. pub fn big_from_decimal_str(s: &str, bits: usize) -> Option { - let dec = DBig::from_str(s.trim()).ok()?; - Some(dec.with_base_and_precision::<2>(bits.max(53)).value()) + Big::from_decimal_str(s, bits) } /// Parse a "re,im,half_height[,iterations]" spec (re/im decimal, parsed at @@ -440,10 +438,10 @@ pub fn interpolate_view(from: &ViewState, to: &ViewState, t: f64) -> ViewState { // 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 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(); - let im1 = to.center_im.clone().with_precision(bits).value(); + let re0 = from.center_re.clone().with_precision(bits); + let im0 = from.center_im.clone().with_precision(bits); + let re1 = to.center_re.clone().with_precision(bits); + let im1 = to.center_im.clone().with_precision(bits); let center_re = &re1 + &(&(&re0 - &re1) * &g_big); let center_im = &im1 + &(&(&im0 - &im1) * &g_big); ViewState::with_center(center_re, center_im, half_height) @@ -456,12 +454,7 @@ pub fn interpolate_f64(from: f64, to: f64, t: f64) -> f64 { /// Render a `Big` as a decimal string with `sig_digits` significant digits. pub fn big_to_decimal_str(x: &Big, sig_digits: usize) -> String { - let dec = x - .to_decimal() - .value() - .with_precision(sig_digits.max(1)) - .value(); - format!("{dec}") + x.to_decimal_string(sig_digits) } /// Precision (bits) needed to resolve the center at a given half-height. @@ -472,12 +465,9 @@ pub fn precision_for(half_height: Scale) -> usize { (zoom_bits + GUARD_BITS).clamp(53, MAX_PRECISION_BITS) } -/// Build an `FBig` from an f64 with an explicit precision context. +/// Build a `Big` from an f64 with an explicit precision. pub fn big_from_f64(x: f64, bits: usize) -> Big { - Big::try_from(x) - .unwrap_or_default() - .with_precision(bits) - .value() + Big::from_f64(x, bits) } #[cfg(test)] @@ -489,8 +479,8 @@ mod tests { } 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(); + let re: f64 = v.center_re.to_f64(); + let im: f64 = v.center_im.to_f64(); (re, im) } @@ -612,8 +602,8 @@ mod tests { 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 dre = ((&mid.center_re - &to.center_re) << k).to_f64(); + let dim = ((&mid.center_im - &to.center_im) << k).to_f64(); 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 d6cb4f5..c403f72 100644 --- a/src/worker.rs +++ b/src/worker.rs @@ -1,7 +1,7 @@ //! Native background worker for reference-orbit computation. //! //! At deep zoom the high-precision reference can take many milliseconds (tens of -//! thousands of `FBig` iterations), which would stutter the UI if done inline. +//! thousands of `Big` iterations), which would stutter the UI if done inline. //! This runs it on a thread and coalesces bursts of requests (e.g. during a //! drag) down to the most recent one. On the web we compute inline instead //! (browsers need a Web Worker for threads); see `app.rs`.