From 4d04b61b69978935f51570791bf6a058cf77578a Mon Sep 17 00:00:00 2001 From: "Claude (numtower)" Date: Thu, 8 Oct 2026 21:53:37 -0700 Subject: [PATCH] feat(numtower): Num, contagion, exact comparison, number language and sums (improvise-59nw.1) WIP: the crate compiles and its unit tests pass; property and differential tests follow. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01VdfEGXA7aJ6ftMsCc9PJze --- Cargo.lock | 167 ++++++++++++- Cargo.toml | 2 +- crates/numtower/Cargo.toml | 25 ++ crates/numtower/examples/bench.rs | 1 + crates/numtower/src/arith.rs | 378 ++++++++++++++++++++++++++++++ crates/numtower/src/cmp.rs | 92 ++++++++ crates/numtower/src/fixed.rs | 232 ++++++++++++++++++ crates/numtower/src/lang.rs | 181 ++++++++++++++ crates/numtower/src/lib.rs | 20 ++ crates/numtower/src/num.rs | 202 ++++++++++++++++ crates/numtower/src/ratio.rs | 172 ++++++++++++++ crates/numtower/src/round.rs | 85 +++++++ crates/numtower/src/sum.rs | 258 ++++++++++++++++++++ 13 files changed, 1812 insertions(+), 3 deletions(-) create mode 100644 crates/numtower/Cargo.toml create mode 100644 crates/numtower/examples/bench.rs create mode 100644 crates/numtower/src/arith.rs create mode 100644 crates/numtower/src/cmp.rs create mode 100644 crates/numtower/src/fixed.rs create mode 100644 crates/numtower/src/lang.rs create mode 100644 crates/numtower/src/lib.rs create mode 100644 crates/numtower/src/num.rs create mode 100644 crates/numtower/src/ratio.rs create mode 100644 crates/numtower/src/round.rs create mode 100644 crates/numtower/src/sum.rs diff --git a/Cargo.lock b/Cargo.lock index 77a3c5b8..51a21dda 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -97,6 +97,12 @@ dependencies = [ "num-traits", ] +[[package]] +name = "arrayvec" +version = "0.7.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d3fb67a6e08acf24fdeccbac2cb6ac4305825bd1f117462e0e6f2f193345ad56" + [[package]] name = "atomic" version = "0.6.1" @@ -178,6 +184,30 @@ dependencies = [ "hybrid-array", ] +[[package]] +name = "borsh" +version = "1.8.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "553c5d846a6ba5150c65e3b1b8ec073bcf1abc20f9b7220de384a4443ea4e20a" +dependencies = [ + "borsh-derive", + "bytes", + "cfg_aliases", +] + +[[package]] +name = "borsh-derive" +version = "1.8.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "12cdfe656708a01f89b451a7d36466e6fe6c414de0aa18fc54f864f6f9ca9f56" +dependencies = [ + "once_cell", + "proc-macro-crate", + "proc-macro2", + "quote", + "syn 3.0.6", +] + [[package]] name = "bumpalo" version = "3.20.3" @@ -705,6 +735,17 @@ dependencies = [ "version_check", ] +[[package]] +name = "getrandom" +version = "0.2.17" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ff2abc00be7fca6ebc474524697ae276ad847ad0a6b3faa4bcb027e9a4614ad0" +dependencies = [ + "cfg-if", + "libc", + "wasi", +] + [[package]] name = "getrandom" version = "0.3.4" @@ -1216,6 +1257,16 @@ dependencies = [ "minimal-lexical", ] +[[package]] +name = "num-bigint" +version = "0.4.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c89e69e7e0f03bea5ef08013795c25018e101932225a656383bd384495ecc367" +dependencies = [ + "num-integer", + "num-traits", +] + [[package]] name = "num-conv" version = "0.2.2" @@ -1233,6 +1284,26 @@ dependencies = [ "syn 2.0.119", ] +[[package]] +name = "num-integer" +version = "0.1.47" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7ce2d95d4b3734dc35aa2f45e1aa22cd416814592a4f9d9205e11affd5b8e10b" +dependencies = [ + "num-traits", +] + +[[package]] +name = "num-rational" +version = "0.4.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f83d14da390562dca69fc84082e73e548e1ad308d24accdedd2720017cb37824" +dependencies = [ + "num-bigint", + "num-integer", + "num-traits", +] + [[package]] name = "num-traits" version = "0.2.19" @@ -1251,6 +1322,18 @@ dependencies = [ "libc", ] +[[package]] +name = "numtower" +version = "0.1.0" +dependencies = [ + "num-bigint", + "num-integer", + "num-rational", + "num-traits", + "proptest", + "rust_decimal", +] + [[package]] name = "once_cell" version = "1.21.4" @@ -1449,6 +1532,15 @@ dependencies = [ "zerocopy", ] +[[package]] +name = "proc-macro-crate" +version = "3.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e67ba7e9b2b56446f1d419b1d807906278ffa1a658a8a5d8a39dcb1f5a78614f" +dependencies = [ + "toml_edit", +] + [[package]] name = "proc-macro2" version = "1.0.107" @@ -1469,7 +1561,7 @@ dependencies = [ "bitflags 2.13.2", "num-traits", "rand 0.9.5", - "rand_chacha", + "rand_chacha 0.9.0", "rand_xorshift", "regex-syntax", "rusty-fork", @@ -1510,6 +1602,8 @@ version = "0.8.8" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "e058c7de0b26af77780c769414d6257830bb240f3c38477dbc2c16e5f54d6d4c" dependencies = [ + "libc", + "rand_chacha 0.3.1", "rand_core 0.6.4", ] @@ -1519,7 +1613,7 @@ version = "0.9.5" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "b9ef1d0d795eb7d84685bca4f72f3649f064e6641543d3a8c415898726a57b41" dependencies = [ - "rand_chacha", + "rand_chacha 0.9.0", "rand_core 0.9.5", ] @@ -1534,6 +1628,16 @@ dependencies = [ "rand_core 0.10.1", ] +[[package]] +name = "rand_chacha" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e6c10a63a0fa32252be49d21e7709d4d4baf8d231c2dbce1eaa8141b9b127d88" +dependencies = [ + "ppv-lite86", + "rand_core 0.6.4", +] + [[package]] name = "rand_chacha" version = "0.9.0" @@ -1549,6 +1653,9 @@ name = "rand_core" version = "0.6.4" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "ec0be4795e2f6a28069bec0b5ff3e2ac9bafc99e6a9a7dc3547996c5c816922c" +dependencies = [ + "getrandom 0.2.17", +] [[package]] name = "rand_core" @@ -1713,6 +1820,23 @@ version = "0.8.11" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "d6f6ff9a378485b298a5286656da665ba74413d36db0979633275d2e708145d4" +[[package]] +name = "rust_decimal" +version = "1.43.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7653272e75dcac41dc199fbea6f5797633994fafd339943c06c9af16bf29cd3a" +dependencies = [ + "arrayvec", + "borsh", + "bytes", + "num-traits", + "rand 0.8.8", + "rand 0.9.5", + "serde", + "serde_json", + "wasm-bindgen", +] + [[package]] name = "rustc_version" version = "0.4.1" @@ -2161,6 +2285,36 @@ dependencies = [ "tungstenite", ] +[[package]] +name = "toml_datetime" +version = "1.1.2+spec-1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2b86d767906c6c42421dcba507eb9d203e779497710a47782a224bb871653053" +dependencies = [ + "serde_core", +] + +[[package]] +name = "toml_edit" +version = "0.25.17+spec-1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e3641d5bbb5349a79e1020a242d251efbc546ad8048d133958323ce9c40a9c9c" +dependencies = [ + "indexmap", + "toml_datetime", + "toml_parser", + "winnow", +] + +[[package]] +name = "toml_parser" +version = "1.1.5+spec-1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "baa693a8032d7e1cada7d0041e96126df243179ff061456783ac7f12bda4744c" +dependencies = [ + "winnow", +] + [[package]] name = "tungstenite" version = "0.30.0" @@ -2525,6 +2679,15 @@ dependencies = [ "windows-link", ] +[[package]] +name = "winnow" +version = "1.0.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "23b97319f7b8343df12cc98938e5c3eb436064524c8d2b4e30a1d3a36eecdf81" +dependencies = [ + "memchr", +] + [[package]] name = "wit-bindgen" version = "0.57.1" diff --git a/Cargo.toml b/Cargo.toml index 1a5467a7..a07fe043 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,5 +1,5 @@ [workspace] -members = [".", "crates/improvise-formula", "crates/improvise-core", "crates/improvise-protocol", "crates/improvise-io", "crates/improvise-command", "crates/improvise-wasm-client", "crates/improvise-ws-server", "crates/improvise-worker-server", "crates/gp-tree", "crates/improvise-gp"] +members = [".", "crates/improvise-formula", "crates/improvise-core", "crates/improvise-protocol", "crates/improvise-io", "crates/improvise-command", "crates/improvise-wasm-client", "crates/improvise-ws-server", "crates/improvise-worker-server", "crates/gp-tree", "crates/improvise-gp", "crates/numtower"] [package] name = "improvise" diff --git a/crates/numtower/Cargo.toml b/crates/numtower/Cargo.toml new file mode 100644 index 00000000..54ca4e27 --- /dev/null +++ b/crates/numtower/Cargo.toml @@ -0,0 +1,25 @@ +[package] +name = "numtower" +version = "0.1.0" +edition = "2024" +license = "Apache-2.0" +description = "A three-level numeric tower: exact decimals, bounded exact fractions and inexact f64, with Scheme's contagion rules and exact, order-free sums." +keywords = ["decimal", "rational", "exact", "numeric-tower"] +categories = ["mathematics"] + +# ADR 005 clause 14: rust_decimal (std only), num-rational without its default +# features (for Ratio), and num-bigint always. No serde. +[dependencies] +rust_decimal = { version = "1.43", default-features = false, features = ["std"] } +num-rational = { version = "0.4.2", default-features = false } +num-bigint = { version = "0.4", default-features = false, features = ["std"] } +num-integer = { version = "0.1", default-features = false } +num-traits = { version = "0.2", default-features = false, features = ["std"] } + +[dev-dependencies] +num-rational = { version = "0.4.2", features = ["num-bigint-std"] } +proptest = "1" + +[[example]] +name = "bench" +test = false diff --git a/crates/numtower/examples/bench.rs b/crates/numtower/examples/bench.rs new file mode 100644 index 00000000..f328e4d9 --- /dev/null +++ b/crates/numtower/examples/bench.rs @@ -0,0 +1 @@ +fn main() {} diff --git a/crates/numtower/src/arith.rs b/crates/numtower/src/arith.rs new file mode 100644 index 00000000..31664681 --- /dev/null +++ b/crates/numtower/src/arith.rs @@ -0,0 +1,378 @@ +//! Arithmetic and contagion (ADR 005 clauses 4 and 5, question 13 (b)). +//! The reasoning for each rule is in `context/code/numtower.md`, "Arithmetic". + +use std::ops::{Add, Mul, Neg, Sub}; + +use num_bigint::{BigInt, BigUint}; +use num_integer::Integer; +use num_traits::{One, Signed, ToPrimitive, Zero}; +use rust_decimal::Decimal; + +use crate::fixed::I256; +use crate::num::{MANTISSA_LIMIT, MAX_SCALE, Num, NumRef, decimal, pow10, small_ratio}; +use crate::ratio::BigQ; + +/// The widest an exact intermediate may grow, in bits of numerator or denominator. +pub const CEILING_BITS: u64 = 8192; + +/// A value inside one evaluation: exact intermediates are unbounded up to the ceiling. +#[derive(Clone, Debug)] +pub struct Interim(Inner); + +#[derive(Clone, Debug)] +enum Inner { + Stored(Num), + Wide(Box), +} + +/// The divisor was exactly zero, exact or inexact (clause 5). +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub struct DivisionByZero; + +impl From for Interim { + fn from(n: Num) -> Interim { + Interim(Inner::Stored(n)) + } +} + +impl Interim { + fn wide(q: BigQ) -> Interim { + Interim(Inner::Wide(Box::new(q))) + } + + /// The stored bound of question 13 (b): a wide value is the nearest f64. + pub fn store(self) -> Num { + match self.0 { + Inner::Stored(n) => n, + Inner::Wide(q) => Num::inexact(q.to_f64()), + } + } + + pub fn is_exact(&self) -> bool { + match &self.0 { + Inner::Stored(n) => n.is_exact(), + Inner::Wide(_) => true, + } + } + + fn operand(&self) -> Operand<'_> { + match &self.0 { + Inner::Stored(n) => Operand::of(n), + Inner::Wide(q) => Operand::Big(q), + } + } + + pub fn checked_div(&self, o: &Interim) -> Result { + divide(self.operand(), o.operand()) + } + + pub fn pow(&self, o: &Interim) -> Result { + power(self.operand(), o.operand()) + } +} + +impl Num { + pub fn checked_div(&self, o: &Num) -> Result { + divide(Operand::of(self), Operand::of(o)).map(Interim::store) + } + + pub fn pow(&self, o: &Num) -> Result { + power(Operand::of(self), Operand::of(o)).map(Interim::store) + } +} + +macro_rules! binary_ops { + ($($trait:ident $method:ident $f:ident;)*) => {$( + impl $trait<&Interim> for &Interim { + type Output = Interim; + fn $method(self, o: &Interim) -> Interim { + $f(self.operand(), o.operand()) + } + } + impl $trait<&Num> for &Num { + type Output = Num; + fn $method(self, o: &Num) -> Num { + $f(Operand::of(self), Operand::of(o)).store() + } + } + )*}; +} + +binary_ops! { + Add add plus; + Sub sub minus; + Mul mul times; +} + +impl Neg for &Interim { + type Output = Interim; + fn neg(self) -> Interim { + negate(self.operand()) + } +} + +impl Neg for &Num { + type Output = Num; + fn neg(self) -> Num { + negate(Operand::of(self)).store() + } +} + +#[derive(Clone, Copy)] +enum Operand<'a> { + Dec(&'a Decimal), + Small(i128, i128), + Big(&'a BigQ), + Float(f64), +} + +impl<'a> Operand<'a> { + fn of(n: &'a Num) -> Operand<'a> { + match n.kind() { + NumRef::Exact(d) => Operand::Dec(d), + NumRef::Fraction(r) => Operand::Small(*r.numer(), *r.denom()), + NumRef::Inexact(x) => Operand::Float(x), + } + } + + fn parts(self) -> Option<(i128, i128)> { + match self { + Operand::Dec(d) => Some((d.mantissa(), pow10(d.scale()))), + Operand::Small(n, d) => Some((n, d)), + Operand::Big(_) | Operand::Float(_) => None, + } + } + + fn big(self) -> Option { + match self { + Operand::Big(q) => Some(q.clone()), + Operand::Float(_) => None, + exact => exact.parts().map(|(n, d)| BigQ::from_i128(n, d)), + } + } + + fn float(self) -> f64 { + match self { + Operand::Dec(d) => crate::num::decimal_to_f64(d), + Operand::Small(n, d) => crate::ratio::ratio_to_f64(&BigInt::from(n), &BigInt::from(d)), + Operand::Big(q) => q.to_f64(), + Operand::Float(x) => x, + } + } + + fn is_float(self) -> bool { + matches!(self, Operand::Float(_)) + } + + fn is_zero(self) -> bool { + match self { + Operand::Dec(d) => d.is_zero(), + Operand::Small(..) | Operand::Big(_) => false, + Operand::Float(x) => x == 0.0, + } + } +} + +fn inexact(x: f64) -> Interim { + Interim::from(Num::inexact(x)) +} + +/// An exact result in the first arm it fits; wide while under the ceiling, else demoted. +fn settle_wide(q: BigQ) -> Interim { + match Num::fit(&q) { + Some(n) => Interim::from(n), + None if q.bits() > CEILING_BITS => inexact(q.to_f64()), + None => Interim::wide(q), + } +} + +fn exact_or_wide(a: Operand, b: Operand, small: impl Fn(i128, i128, i128, i128) -> Option<(i128, i128)>, big: impl Fn(&BigQ, &BigQ) -> BigQ) -> Interim { + let fast = a.parts().zip(b.parts()).and_then(|((n1, d1), (n2, d2))| small(n1, d1, n2, d2)); + fast.and_then(|(n, d)| small_ratio(n, d)) + .map(Interim::from) + .unwrap_or_else(|| match a.big().zip(b.big()) { + Some((x, y)) => settle_wide(big(&x, &y)), + None => unreachable!("both operands are exact here"), + }) +} + +fn plus(a: Operand, b: Operand) -> Interim { + match (a, b) { + (x, y) if x.is_float() || y.is_float() => inexact(x.float() + y.float()), + (Operand::Dec(x), Operand::Dec(y)) => decimal_sum(x, y), + _ => exact_or_wide( + a, + b, + |n1, d1, n2, d2| { + let n = n1.checked_mul(d2)?.checked_add(n2.checked_mul(d1)?)?; + Some((n, d1.checked_mul(d2)?)) + }, + BigQ::add, + ), + } +} + +fn minus(a: Operand, b: Operand) -> Interim { + plus(a, negate(b).operand()) +} + +fn times(a: Operand, b: Operand) -> Interim { + match (a, b) { + (x, y) if x.is_float() || y.is_float() => inexact(x.float() * y.float()), + (Operand::Dec(x), Operand::Dec(y)) => decimal_product(x, y), + _ => exact_or_wide( + a, + b, + |n1, d1, n2, d2| { + let (g1, g2) = (n1.gcd(&d2).max(1), n2.gcd(&d1).max(1)); + Some(((n1 / g1).checked_mul(n2 / g2)?, (d1 / g2).checked_mul(d2 / g1)?)) + }, + BigQ::mul, + ), + } +} + +fn negate(a: Operand) -> Interim { + match a { + Operand::Dec(d) => Interim::from(Num::exact(-*d)), + Operand::Float(x) => inexact(-x), + Operand::Small(n, d) => small_ratio(-n, d) + .map(Interim::from) + .unwrap_or_else(|| settle_wide(BigQ::from_i128(n, d).neg())), + Operand::Big(q) => settle_wide(q.neg()), + } +} + +fn divide(a: Operand, b: Operand) -> Result { + match (a, b) { + (_, y) if y.is_zero() => Err(DivisionByZero), + (x, y) if x.is_float() || y.is_float() => Ok(inexact(x.float() / y.float())), + _ => Ok(exact_or_wide( + a, + b, + |n1, d1, n2, d2| { + let (g1, g2) = (n1.gcd(&n2).max(1), d1.gcd(&d2).max(1)); + Some(((n1 / g1).checked_mul(d2 / g2)?, (d1 / g2).checked_mul(n2 / g1)?)) + }, + |x, y| x.div(y).unwrap_or_else(|| unreachable!("the divisor is not zero")), + )), + } +} + +/// Clause 4's hint: a returned scale below the nominal one may mean rounding. +fn decimal_sum(x: &Decimal, y: &Decimal) -> Interim { + let nominal = x.scale().max(y.scale()); + match x.checked_add(*y).filter(|r| r.scale() == nominal) { + Some(r) => Interim::from(Num::exact(r)), + None => { + let align = |d: &Decimal| { + I256::from_i128(d.mantissa()) + .checked_mul_u128(pow10(nominal - d.scale()) as u128) + .expect("an aligned mantissa is below 2^190") + }; + let sum = align(x).checked_add(&align(y)).expect("below 2^191"); + settle_decimal_integer(sum, nominal, nominal) + } + } +} + +fn decimal_product(x: &Decimal, y: &Decimal) -> Interim { + let nominal = x.scale() + y.scale(); + match x.checked_mul(*y).filter(|r| r.scale() == nominal) { + Some(r) => Interim::from(Num::exact(r)), + None => settle_decimal_integer(I256::mul_i128(x.mantissa(), y.mantissa()), nominal, nominal), + } +} + +/// The recompute of clause 4: strip trailing zeros, fit 96 bits at scale <= 28, else a fraction. +fn settle_decimal_integer(m: I256, scale: u32, nominal: u32) -> Interim { + let (mut m, mut s) = (m, scale); + while s > 0 { + let (q, r) = m.div_rem_u64(10); + if r != 0 { + break; + } + m = q; + s -= 1; + } + let fits = m.to_i128().filter(|v| v.unsigned_abs() < MANTISSA_LIMIT as u128); + match fits { + Some(v) if s <= MAX_SCALE => Interim::from(padded(v, s, nominal.min(MAX_SCALE))), + _ => settle_wide(BigQ::new(m.to_bigint(), BigInt::from(10).pow(s))), + } +} + +/// The stored scale of clause 4: pad toward `target` while the mantissa fits 96 bits. +pub(crate) fn padded(mantissa: i128, scale: u32, target: u32) -> Num { + let (mut m, mut s) = (mantissa, scale); + while s < target && (m * 10).unsigned_abs() < MANTISSA_LIMIT as u128 { + m *= 10; + s += 1; + } + decimal(m, s).expect("a padded mantissa still fits") +} + +/// The exponent of an exact integer-valued operand. +fn integer_value(b: Operand) -> Option { + match b { + Operand::Dec(d) if d.fract().is_zero() => Some(BigInt::from(d.mantissa() / pow10(d.scale()))), + Operand::Small(n, 1) => Some(BigInt::from(n)), + Operand::Big(q) if q.den.is_one() => Some(q.num.clone()), + _ => None, + } +} + +fn power(a: Operand, b: Operand) -> Result { + match (a.is_float(), integer_value(b)) { + (false, Some(n)) => exact_power(a, &n), + _ => Ok(inexact(a.float().powf(b.float()))), + } +} + +/// Clause 5's `POW`: the size is checked before anything is computed. +fn exact_power(a: Operand, n: &BigInt) -> Result { + let q = a.big().expect("an exact base"); + let base_scale = match a { + Operand::Dec(d) => d.scale(), + _ => 0, + }; + let k = n.magnitude(); + let target = k + .to_u32() + .and_then(|k| k.checked_mul(base_scale)) + .map_or(MAX_SCALE, |s| s.min(MAX_SCALE)); + let unit = q.den.is_one() && q.num.magnitude().is_one(); + match (q.num.is_zero(), n.is_negative(), unit) { + (_, _, _) if n.is_zero() => Ok(Interim::from(padded(1, 0, 0))), + (true, true, _) => Err(DivisionByZero), + (true, false, _) => Ok(Interim::from(padded(0, 0, target))), + (false, _, true) => { + let negative = q.num.is_negative() && k.is_odd(); + Ok(Interim::from(padded(if negative { -1 } else { 1 }, 0, target))) + } + (false, negative, false) => { + let lower = |x: &BigInt| BigUint::from(x.bits().saturating_sub(1)) * k; + let too_big = lower(&q.num).max(lower(&q.den)) > BigUint::from(CEILING_BITS); + match (too_big, k.to_u32()) { + (false, Some(k)) => Ok(exact_integer_power(&q, k, negative, target)), + _ => Ok(inexact(a.float().powf(n.to_f64().unwrap_or(f64::INFINITY)))), + } + } + } +} + +fn exact_integer_power(q: &BigQ, k: u32, negative: bool, target: u32) -> Interim { + let (num, den) = (q.num.pow(k), q.den.pow(k)); + let value = match negative { + true => BigQ::new(den, num), + false => BigQ { num, den }, + }; + let settled = settle_wide(value); + match (&settled.0, negative) { + (Inner::Stored(n), false) => match n.kind() { + NumRef::Exact(d) => Interim::from(padded(d.mantissa(), d.scale(), target)), + _ => settled, + }, + _ => settled, + } +} diff --git a/crates/numtower/src/cmp.rs b/crates/numtower/src/cmp.rs new file mode 100644 index 00000000..ffa2e40a --- /dev/null +++ b/crates/numtower/src/cmp.rs @@ -0,0 +1,92 @@ +//! Exact comparison at every pairing of arms (clause 6), and the `MIN`/`MAX` +//! tie-breaks (clause 8). Why: `context/code/numtower.md`, "Comparison". + +use std::cmp::Ordering; + +use crate::fixed::I256; +use crate::num::{Num, NumRef}; +use crate::ratio::BigQ; + +impl Num { + /// The order of the two values, exactly; `None` only when either is NaN. + pub fn cmp_exact(&self, o: &Num) -> Option { + match (self.exact_parts(), o.exact_parts()) { + (Some(a), Some(b)) => Some(cross(a, b)), + (Some(_), None) => mixed(self, o.to_f64()), + (None, Some(_)) => mixed(o, self.to_f64()).map(Ordering::reverse), + (None, None) => self.to_f64().partial_cmp(&o.to_f64()), + } + } + + /// The smaller, NaN skipped; ties keep the most information (clause 8, question 11). + pub fn min_of<'a>(&'a self, o: &'a Num) -> &'a Num { + pick(self, o, Ordering::Less) + } + + pub fn max_of<'a>(&'a self, o: &'a Num) -> &'a Num { + pick(self, o, Ordering::Greater) + } +} + +impl PartialEq for Num { + fn eq(&self, o: &Num) -> bool { + self.cmp_exact(o) == Some(Ordering::Equal) + } +} + +impl PartialOrd for Num { + fn partial_cmp(&self, o: &Num) -> Option { + self.cmp_exact(o) + } +} + +/// `a.0/a.1` against `b.0/b.1`, denominators positive, in 256 bits. +fn cross(a: (i128, i128), b: (i128, i128)) -> Ordering { + I256::mul_i128(a.0, b.1).cmp(&I256::mul_i128(b.0, a.1)) +} + +/// An exact value against an f64: the rounded value decides unless it ties. +fn mixed(exact: &Num, x: f64) -> Option { + match (x.is_nan(), x.is_infinite()) { + (true, _) => None, + (false, true) => Some(if x > 0.0 { Ordering::Less } else { Ordering::Greater }), + (false, false) => match exact.to_f64().partial_cmp(&x)? { + Ordering::Equal => { + let q = exact.exact_q().expect("an exact operand"); + Some(q.cmp_value(&BigQ::from_f64(x))) + } + decided => Some(decided), + }, + } +} + +fn pick<'a>(a: &'a Num, b: &'a Num, want: Ordering) -> &'a Num { + match (a.cmp_exact(b), is_nan(a), is_nan(b)) { + (_, true, _) => b, + (_, false, true) => a, + (Some(o), _, _) if o == want => a, + (Some(Ordering::Equal), _, _) => tie(a, b, want), + _ => b, + } +} + +fn is_nan(n: &Num) -> bool { + matches!(n.kind(), NumRef::Inexact(x) if x.is_nan()) +} + +/// Equal values: exact over inexact, the larger scale, and IEEE's signed zero. +fn tie<'a>(a: &'a Num, b: &'a Num, want: Ordering) -> &'a Num { + match (a.kind(), b.kind()) { + (NumRef::Inexact(_), NumRef::Exact(_) | NumRef::Fraction(_)) => b, + (NumRef::Exact(_) | NumRef::Fraction(_), NumRef::Inexact(_)) => a, + (NumRef::Exact(x), NumRef::Exact(y)) if y.scale() > x.scale() => b, + (NumRef::Inexact(x), NumRef::Inexact(y)) => { + let negative_wanted = want == Ordering::Less; + match y.is_sign_negative() == negative_wanted && x.is_sign_negative() != negative_wanted { + true => b, + false => a, + } + } + _ => a, + } +} diff --git a/crates/numtower/src/fixed.rs b/crates/numtower/src/fixed.rs new file mode 100644 index 00000000..89175010 --- /dev/null +++ b/crates/numtower/src/fixed.rs @@ -0,0 +1,232 @@ +//! Fixed-width two's-complement integers for the hot paths: no allocation. +//! Why they are hand-rolled: `context/code/numtower.md`, "The wide integers". + +use std::cmp::Ordering; + +use num_bigint::{BigInt, Sign}; + +/// A signed integer of `64 * N` bits, little-endian limbs, two's complement. +#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)] +pub(crate) struct Fixed { + limbs: [u64; N], +} + +/// The 256-bit integer clauses 4–6 use, and the 320-bit decimal partition of clause 8. +pub(crate) type I256 = Fixed<4>; +pub(crate) type I320 = Fixed<5>; + +impl Fixed { + pub(crate) const ZERO: Self = Fixed { limbs: [0; N] }; + + pub(crate) fn from_i128(x: i128) -> Self { + let fill = if x < 0 { u64::MAX } else { 0 }; + let mut limbs = [fill; N]; + limbs[0] = x as u64; + limbs[1] = (x >> 64) as u64; + Fixed { limbs } + } + + pub(crate) fn is_negative(&self) -> bool { + self.limbs[N - 1] >> 63 == 1 + } + + pub(crate) fn is_zero(&self) -> bool { + self.limbs.iter().all(|&l| l == 0) + } + + fn wrapping_neg(&self) -> Self { + let mut out = [0u64; N]; + let mut carry = 1u64; + for (o, &l) in out.iter_mut().zip(&self.limbs) { + let (v, c) = (!l).overflowing_add(carry); + *o = v; + carry = c as u64; + } + Fixed { limbs: out } + } + + fn wrapping_add(&self, other: &Self) -> Self { + let mut out = [0u64; N]; + let mut carry = false; + for i in 0..N { + let (a, c1) = self.limbs[i].overflowing_add(other.limbs[i]); + let (b, c2) = a.overflowing_add(carry as u64); + out[i] = b; + carry = c1 || c2; + } + Fixed { limbs: out } + } + + pub(crate) fn checked_add(&self, other: &Self) -> Option { + let sum = self.wrapping_add(other); + let same = self.is_negative() == other.is_negative(); + (!same || sum.is_negative() == self.is_negative()).then_some(sum) + } + + pub(crate) fn checked_neg(&self) -> Option { + let neg = self.wrapping_neg(); + (self.is_zero() || neg.is_negative() != self.is_negative()).then_some(neg) + } + + /// The magnitude as unsigned limbs, with the sign split out. + fn magnitude(&self) -> (bool, [u64; N]) { + match self.is_negative() { + true => (true, self.wrapping_neg().limbs), + false => (false, self.limbs), + } + } + + fn from_magnitude(negative: bool, limbs: [u64; N]) -> Option { + let unsigned = Fixed { limbs }; + match (negative, unsigned.is_negative()) { + (false, false) => Some(unsigned), + (true, _) => { + let neg = unsigned.wrapping_neg(); + (neg.is_negative() || neg.is_zero()).then_some(neg) + } + (false, true) => None, + } + } + + pub(crate) fn checked_mul_u128(&self, k: u128) -> Option { + let (negative, mag) = self.magnitude(); + let mut out = [0u64; N]; + for (j, kk) in [k as u64, (k >> 64) as u64].into_iter().enumerate() { + if kk == 0 { + continue; + } + let mut carry = 0u128; + for i in 0..N { + let idx = i + j; + let prod = mag[i] as u128 * kk as u128 + carry; + match idx < N { + true => { + let s = out[idx] as u128 + (prod as u64) as u128; + out[idx] = s as u64; + carry = (prod >> 64) + (s >> 64); + } + false if prod != 0 => return None, + false => carry = 0, + } + } + if carry != 0 { + return None; + } + } + Self::from_magnitude(negative, out) + } + + /// Divides the magnitude by `k`, returning the quotient (sign kept) and the remainder. + pub(crate) fn div_rem_u64(&self, k: u64) -> (Self, u64) { + let (negative, mag) = self.magnitude(); + let mut out = [0u64; N]; + let mut rem = 0u128; + for i in (0..N).rev() { + let cur = (rem << 64) | mag[i] as u128; + out[i] = (cur / k as u128) as u64; + rem = cur % k as u128; + } + let q = Self::from_magnitude(negative, out).unwrap_or(Self::ZERO); + (q, rem as u64) + } + + pub(crate) fn to_i128(self) -> Option { + let fill = if self.is_negative() { u64::MAX } else { 0 }; + let hi_ok = self.limbs[2..].iter().all(|&l| l == fill); + let v = (self.limbs[0] as u128 | ((self.limbs[1] as u128) << 64)) as i128; + (hi_ok && (v < 0) == self.is_negative()).then_some(v) + } + + pub(crate) fn to_bigint(self) -> BigInt { + let (negative, mag) = self.magnitude(); + let digits: Vec = mag + .iter() + .flat_map(|&l| [l as u32, (l >> 32) as u32]) + .collect(); + let sign = if negative { Sign::Minus } else { Sign::Plus }; + BigInt::from_slice(sign, &digits) + } + + /// The product of two `i128`s; needs `N >= 4`. + pub(crate) fn mul_i128(a: i128, b: i128) -> Self { + let negative = (a < 0) != (b < 0); + let (x, y) = (a.unsigned_abs(), b.unsigned_abs()); + let (x0, x1) = (x as u64 as u128, x >> 64); + let (y0, y1) = (y as u64 as u128, y >> 64); + let p00 = x0 * y0; + let p01 = x0 * y1; + let p10 = x1 * y0; + let p11 = x1 * y1; + let mid = (p00 >> 64) + (p01 as u64 as u128) + (p10 as u64 as u128); + let mut limbs = [0u64; N]; + limbs[0] = p00 as u64; + limbs[1] = mid as u64; + let hi = (mid >> 64) + (p01 >> 64) + (p10 >> 64) + p11; + limbs[2] = hi as u64; + limbs[3] = (hi >> 64) as u64; + let unsigned = Fixed { limbs }; + match negative { + true => unsigned.wrapping_neg(), + false => unsigned, + } + } +} + +impl Default for Fixed { + fn default() -> Self { + Self::ZERO + } +} + +impl Ord for Fixed { + fn cmp(&self, other: &Self) -> Ordering { + match (self.is_negative(), other.is_negative()) { + (true, false) => Ordering::Less, + (false, true) => Ordering::Greater, + _ => self.limbs.iter().rev().cmp(other.limbs.iter().rev()), + } + } +} + +impl PartialOrd for Fixed { + fn partial_cmp(&self, other: &Self) -> Option { + Some(self.cmp(other)) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn the_256_bit_product_of_two_i128s_matches_bigint() { + let cases = [ + (i128::MAX, i128::MAX), + (i128::MAX, -i128::MAX), + (-1, 1), + (0, i128::MIN), + (123_456_789_012_345_678_901_234, -98_765_432_109_876_543_210), + ]; + for (a, b) in cases { + let want = BigInt::from(a) * BigInt::from(b); + assert_eq!(I256::mul_i128(a, b).to_bigint(), want, "{a} * {b}"); + } + } + + #[test] + fn checked_arithmetic_refuses_what_overflows_and_keeps_what_fits() { + let max = I256::from_i128(i128::MAX).checked_mul_u128(u128::MAX).unwrap(); + assert_eq!( + max.to_bigint(), + BigInt::from(i128::MAX) * BigInt::from(u128::MAX) + ); + assert!(max.checked_mul_u128(u128::MAX).is_none()); + let neg = max.checked_neg().unwrap(); + assert_eq!(neg.to_bigint(), -max.to_bigint()); + let (q, r) = neg.div_rem_u64(10); + assert_eq!(q.to_bigint() * 10 - BigInt::from(r), neg.to_bigint()); + assert_eq!(I256::from_i128(-5).to_i128(), Some(-5)); + assert_eq!(max.to_i128(), None); + assert!(I256::from_i128(-1) < I256::ZERO); + } +} diff --git a/crates/numtower/src/lang.rs b/crates/numtower/src/lang.rs new file mode 100644 index 00000000..7bb16c2b --- /dev/null +++ b/crates/numtower/src/lang.rs @@ -0,0 +1,181 @@ +//! Clause 3's one number language: the gate, the exact parse, the forced `n/d` arm, +//! and the canonical text. Why each rule: `context/code/numtower.md`, "The number language". + +use std::fmt; + +use num_bigint::BigInt; +use num_traits::{One, Zero}; + +use crate::arith::padded; +use crate::num::{MAX_SCALE, Num, NumRef}; +use crate::ratio::BigQ; + +/// Whether the type of the text is being guessed or was forced (question 12). +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum Typing { + Inferred, + Forced, +} + +/// The gate: the strings `f64::from_str` accepts and reads as a finite value. +pub fn is_number(text: &str, typing: Typing) -> bool { + gate(text) || (typing == Typing::Forced && fraction_lexeme(text).is_some()) +} + +/// A number if the text is one, read exactly; `None` means the text is text. +pub fn parse_number(text: &str, typing: Typing) -> Option { + match (gate(text), typing) { + (true, _) => Lexeme::read(text).map(|lx| lx.value(text)), + (false, Typing::Forced) => fraction_lexeme(text).map(|q| Num::settle(&q)), + (false, Typing::Inferred) => None, + } +} + +fn gate(text: &str) -> bool { + text.parse::().is_ok_and(f64::is_finite) +} + +/// `[+-]digits/digits` with a nonzero denominator, as an exact value. +fn fraction_lexeme(text: &str) -> Option { + let (n, d) = text.split_once('/')?; + let digits = |s: &str| !s.is_empty() && s.bytes().all(|b| b.is_ascii_digit()); + let unsigned = n.strip_prefix(['+', '-']).unwrap_or(n); + let ok = digits(unsigned) && digits(d); + let num: BigInt = ok.then(|| n.trim_start_matches('+').parse().ok())??; + let den: BigInt = d.parse().ok()?; + (!den.is_zero()).then(|| BigQ::new(num, den)) +} + +/// A lexeme of the gate's language: sign, significant digits, and the power of ten. +#[derive(Debug, PartialEq, Eq)] +pub(crate) struct Lexeme { + pub(crate) negative: bool, + pub(crate) digits: String, + pub(crate) exp10: i64, +} + +/// Exponents beyond this cannot reach either exact arm, and saturate. +const EXP_SATURATE: i64 = 1 << 40; + +impl Lexeme { + /// `[+-]? (d+ (. d*)? | . d+) ([eE] [+-]? d+)?`, and nothing else. + pub(crate) fn read(text: &str) -> Option { + let (negative, rest) = match text.as_bytes().first()? { + b'-' => (true, &text[1..]), + b'+' => (false, &text[1..]), + _ => (false, text), + }; + let (mantissa, exponent) = match rest.find(['e', 'E']) { + Some(i) => (&rest[..i], Some(&rest[i + 1..])), + None => (rest, None), + }; + let (int, frac) = mantissa.split_once('.').unwrap_or((mantissa, "")); + let all_digits = |s: &str| s.bytes().all(|b| b.is_ascii_digit()); + let mantissa_ok = all_digits(int) && all_digits(frac) && !(int.is_empty() && frac.is_empty()); + let exp = match exponent { + Some(e) => read_exponent(e)?, + None => 0, + }; + let digits = format!("{int}{frac}"); + mantissa_ok.then(|| Lexeme { + negative, + digits, + exp10: exp.saturating_sub(frac.len() as i64), + }) + } + + /// Clause 3's order: a decimal at the lexeme's scale, then a fraction, then the nearest f64. + fn value(&self, text: &str) -> Num { + let significant = self.digits.trim_start_matches('0'); + let stripped = significant.trim_end_matches('0'); + let implied_scale = (-self.exp10).clamp(0, MAX_SCALE as i64) as u32; + match stripped.is_empty() { + true => padded(0, 0, implied_scale), + false => { + let shift = self.exp10.saturating_add((significant.len() - stripped.len()) as i64); + self.exact(stripped, shift, implied_scale) + .unwrap_or_else(|| Num::inexact(text.parse().unwrap_or(f64::NAN))) + } + } + } + + fn exact(&self, stripped: &str, shift: i64, implied_scale: u32) -> Option { + let leading = stripped.len() as i64 + shift; + let reachable = stripped.len() <= 180 && (-130..=40).contains(&leading) && shift >= -400; + let q = reachable.then(|| self.rational(stripped, shift))?; + Num::fit(&q).map(|n| match n.kind() { + NumRef::Exact(d) => padded(d.mantissa(), d.scale(), implied_scale), + _ => n, + }) + } + + pub(crate) fn rational(&self, stripped: &str, shift: i64) -> BigQ { + let m: BigInt = stripped.parse().unwrap_or_default(); + let m = if self.negative { -m } else { m }; + let ten = BigInt::from(10); + match shift >= 0 { + true => BigQ::new(m * ten.pow(shift as u32), BigInt::one()), + false => BigQ::new(m, ten.pow((-shift) as u32)), + } + } + + /// The exact value, whatever its size: for display rounding. + pub(crate) fn exact_value(&self) -> BigQ { + let stripped = self.digits.trim_start_matches('0'); + match stripped.is_empty() { + true => BigQ::new(BigInt::zero(), BigInt::one()), + false => self.rational(stripped, self.exp10), + } + } +} + +fn read_exponent(e: &str) -> Option { + let (negative, digits) = match e.as_bytes().first()? { + b'-' => (true, &e[1..]), + b'+' => (false, &e[1..]), + _ => (false, e), + }; + let ok = !digits.is_empty() && digits.bytes().all(|b| b.is_ascii_digit()); + let significant = digits.trim_start_matches('0'); + let magnitude = match significant.len() > 12 { + true => EXP_SATURATE, + false => significant.parse::().unwrap_or(0).min(EXP_SATURATE), + }; + ok.then_some(if negative { -magnitude } else { magnitude }) +} + +/// The shortest text that reads back as `x`; `None` for a non-finite value. +pub fn shortest_round_trip(x: f64) -> Option { + x.is_finite().then(|| format!("{x:?}")) +} + +const INEXACT_MARK: &str = "#i"; + +impl Num { + /// Reads every spelling `Display` writes: a decimal, `n/d`, `#i…`, `inf`, `-inf`, `nan`. + pub fn from_canonical(text: &str) -> Option { + match text { + "inf" => Some(Num::inexact(f64::INFINITY)), + "-inf" => Some(Num::inexact(f64::NEG_INFINITY)), + "nan" => Some(Num::inexact(f64::NAN)), + _ => match text.strip_prefix(INEXACT_MARK) { + Some(rest) if gate(rest) => rest.parse().ok().map(Num::inexact), + Some(_) => None, + None => parse_number(text, Typing::Forced).filter(Num::is_exact), + }, + } + } +} + +/// The canonical text: a decimal with its scale, `n/d`, `#i` and the shortest spelling, or a token. +impl fmt::Display for Num { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self.kind() { + NumRef::Exact(d) => write!(f, "{d}"), + NumRef::Fraction(r) => write!(f, "{}/{}", r.numer(), r.denom()), + NumRef::Inexact(x) if x.is_nan() => f.write_str("nan"), + NumRef::Inexact(x) if x.is_infinite() => f.write_str(if x > 0.0 { "inf" } else { "-inf" }), + NumRef::Inexact(x) => write!(f, "{INEXACT_MARK}{x:?}"), + } + } +} diff --git a/crates/numtower/src/lib.rs b/crates/numtower/src/lib.rs new file mode 100644 index 00000000..458d3e9f --- /dev/null +++ b/crates/numtower/src/lib.rs @@ -0,0 +1,20 @@ +//! A three-level numeric tower: exact decimal, exact bounded fraction, inexact f64. +//! The decisions are ADR 005; the reasoning is `context/code/numtower.md`. + +mod arith; +mod cmp; +mod fixed; +mod lang; +mod num; +mod ratio; +mod round; +mod sum; + +pub use arith::{CEILING_BITS, DivisionByZero, Interim}; +pub use lang::{Typing, is_number, parse_number, shortest_round_trip}; +pub use num::{Num, NumRef}; +pub use round::Rounding; +pub use sum::{ExactSum, K, NumSum}; + +pub use num_rational::Ratio; +pub use rust_decimal::Decimal; diff --git a/crates/numtower/src/num.rs b/crates/numtower/src/num.rs new file mode 100644 index 00000000..50075814 --- /dev/null +++ b/crates/numtower/src/num.rs @@ -0,0 +1,202 @@ +//! The stored number: an exact decimal, an exact bounded fraction, or an f64. +//! The arms are private so a value is always normalized: `context/code/numtower.md`. + +use num_bigint::BigInt; +use num_integer::Integer; +use num_rational::Ratio; +use num_traits::{One, ToPrimitive}; +use rust_decimal::Decimal; + +use crate::ratio::{BigQ, ratio_to_f64}; + +/// The most significant digits a `Decimal` holds, and its mantissa bound. +pub(crate) const MAX_SCALE: u32 = 28; +pub(crate) const MANTISSA_LIMIT: i128 = 1 << 96; + +/// A number of the tower (ADR 005 clause 2), always in its one normal form. +#[derive(Clone, Debug)] +pub struct Num(Repr); + +#[derive(Clone, Debug)] +enum Repr { + Exact(Decimal), + Fraction(Box>), + Inexact(f64), +} + +/// A borrowed view of a `Num`'s arm, for exhaustive matching by callers. +#[derive(Clone, Copy, Debug, PartialEq)] +pub enum NumRef<'a> { + Exact(&'a Decimal), + Fraction(&'a Ratio), + Inexact(f64), +} + +impl Num { + /// An exact decimal. Zero loses its sign: the sign of zero is an f64 artefact. + pub fn exact(d: Decimal) -> Num { + let mut d = d; + if d.is_zero() { + d.set_sign_positive(true); + } + Num(Repr::Exact(d)) + } + + pub fn inexact(x: f64) -> Num { + Num(Repr::Inexact(x)) + } + + /// The exact value `num / den`, normalized; `None` when `den` is zero. + pub fn ratio(num: i128, den: i128) -> Option { + (den != 0).then(|| Num::settle(&BigQ::from_i128(num, den))) + } + + fn fraction(num: i128, den: i128) -> Num { + Num(Repr::Fraction(Box::new(Ratio::new_raw(num, den)))) + } + + pub fn kind(&self) -> NumRef<'_> { + match &self.0 { + Repr::Exact(d) => NumRef::Exact(d), + Repr::Fraction(r) => NumRef::Fraction(r), + Repr::Inexact(x) => NumRef::Inexact(*x), + } + } + + pub fn is_exact(&self) -> bool { + !matches!(self.0, Repr::Inexact(_)) + } + + /// Exactly zero, at either level; the one zero predicate of clause 6. + pub fn is_zero(&self) -> bool { + match &self.0 { + Repr::Exact(d) => d.is_zero(), + Repr::Fraction(_) => false, + Repr::Inexact(x) => *x == 0.0, + } + } + + /// The correctly rounded f64 (clause 4): a decimal through its text, a fraction exactly. + pub fn to_f64(&self) -> f64 { + match &self.0 { + Repr::Exact(d) => decimal_to_f64(d), + Repr::Fraction(r) => ratio_to_f64(&BigInt::from(*r.numer()), &BigInt::from(*r.denom())), + Repr::Inexact(x) => *x, + } + } + + /// The exact value as a rational; `None` for an inexact number. + pub(crate) fn exact_parts(&self) -> Option<(i128, i128)> { + match &self.0 { + Repr::Exact(d) => Some((d.mantissa(), pow10(d.scale()))), + Repr::Fraction(r) => Some((*r.numer(), *r.denom())), + Repr::Inexact(_) => None, + } + } + + pub(crate) fn exact_q(&self) -> Option { + self.exact_parts().map(|(n, d)| BigQ::from_i128(n, d)) + } + + /// The stored bound: a decimal if it fits, else an `i128` fraction, else the nearest f64. + pub(crate) fn settle(q: &BigQ) -> Num { + Num::fit(q).unwrap_or_else(|| Num::inexact(q.to_f64())) + } + + /// `q` in an exact arm, at its minimal scale, if it fits one. + pub(crate) fn fit(q: &BigQ) -> Option { + decimal_of(q).or_else(|| { + let num = to_i128(&q.num)?; + let den = to_i128(&q.den)?; + Some(Num::fraction(num, den)) + }) + } +} + +/// `num / den` normalized without allocating; `None` hands the case to the bigint path. +pub(crate) fn small_ratio(num: i128, den: i128) -> Option { + let representable = num != i128::MIN && den != i128::MIN && den != 0; + representable.then(|| { + let g = num.gcd(&den); + let (n, d) = match den < 0 { + true => (-num / g, -den / g), + false => (num / g, den / g), + }; + let p = d.trailing_zeros(); + let mut rest = d >> p; + let mut q5 = 0; + while rest % 5 == 0 { + rest /= 5; + q5 += 1; + } + let k = p.max(q5); + let as_decimal = (rest == 1 && k <= MAX_SCALE) + .then(|| 2i128.pow(k - p) * 5i128.pow(k - q5)) + .and_then(|pad| n.checked_mul(pad)) + .and_then(|m| decimal(m, k)); + as_decimal.unwrap_or_else(|| Num::fraction(n, d)) + }) +} + +fn to_i128(x: &BigInt) -> Option { + x.to_i128().filter(|v| *v != i128::MIN) +} + +pub(crate) fn pow10(k: u32) -> i128 { + 10i128.pow(k) +} + +pub(crate) fn decimal_to_f64(d: &Decimal) -> f64 { + d.to_string().parse().unwrap_or(f64::NAN) +} + +/// A decimal from a mantissa and scale, if the mantissa fits 96 bits. +pub(crate) fn decimal(mantissa: i128, scale: u32) -> Option { + let fits = mantissa.unsigned_abs() < MANTISSA_LIMIT as u128 && scale <= MAX_SCALE; + fits.then(|| Num::exact(Decimal::from_i128_with_scale(mantissa, scale))) +} + +/// The direct test of clause 5: a reduced `n / (2^p 5^q)` is a decimal iff it fits. +fn decimal_of(q: &BigQ) -> Option { + let p = q.den.trailing_zeros().unwrap_or(0); + let mut rest = &q.den >> p as usize; + let five = BigInt::from(5); + let mut q5 = 0u64; + while rest.is_multiple_of(&five) && q5 <= MAX_SCALE as u64 { + rest /= &five; + q5 += 1; + } + let k = p.max(q5); + if !rest.is_one() || k > MAX_SCALE as u64 { + return None; + } + let pad = (BigInt::one() << (k - p) as usize) * five.pow((k - q5) as u32); + let mantissa = (&q.num * pad).to_i128()?; + decimal(mantissa, k as u32) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn a_num_is_24_bytes() { + assert_eq!(std::mem::size_of::(), 24); + } + + #[test] + fn a_fraction_that_is_a_fitting_decimal_is_stored_as_the_decimal() { + let half = Num::ratio(1, 2).unwrap(); + assert_eq!(half.kind(), NumRef::Exact(&Decimal::new(5, 1))); + let third = Num::ratio(-10, 3).unwrap(); + assert_eq!(third.kind(), NumRef::Fraction(&Ratio::new_raw(-10, 3))); + let tiny = Num::ratio(1, 1 << 40).unwrap(); + assert!(matches!(tiny.kind(), NumRef::Fraction(_)), "{tiny:?}"); + let mut minus_zero = Decimal::new(0, 2); + minus_zero.set_sign_negative(true); + assert!(minus_zero.is_sign_negative()); + let zero = Num::exact(minus_zero); + assert!(!zero.to_f64().is_sign_negative()); + assert!(Num::ratio(1, 0).is_none()); + } +} diff --git a/crates/numtower/src/ratio.rs b/crates/numtower/src/ratio.rs new file mode 100644 index 00000000..0c1a2d5c --- /dev/null +++ b/crates/numtower/src/ratio.rs @@ -0,0 +1,172 @@ +//! Exact rationals in bigints: the widened intermediate of question 13 (b), and +//! the one correctly rounded conversion of a rational to f64 (clauses 4 and 5). + +use std::cmp::Ordering; + +use num_bigint::{BigInt, BigUint, Sign}; +use num_integer::Integer; +use num_traits::{One, Signed, ToPrimitive, Zero}; + +/// A reduced fraction with a positive denominator. +#[derive(Clone, Debug, PartialEq, Eq)] +pub(crate) struct BigQ { + pub(crate) num: BigInt, + pub(crate) den: BigInt, +} + +impl BigQ { + pub(crate) fn new(num: BigInt, den: BigInt) -> BigQ { + let (num, den) = match den.is_negative() { + true => (-num, -den), + false => (num, den), + }; + let g = num.gcd(&den); + match g.is_zero() || g.is_one() { + true => BigQ { num, den }, + false => BigQ { + num: num / &g, + den: den / &g, + }, + } + } + + pub(crate) fn from_i128(num: i128, den: i128) -> BigQ { + BigQ::new(BigInt::from(num), BigInt::from(den)) + } + + pub(crate) fn add(&self, o: &BigQ) -> BigQ { + BigQ::new(&self.num * &o.den + &o.num * &self.den, &self.den * &o.den) + } + + pub(crate) fn neg(&self) -> BigQ { + BigQ { + num: -&self.num, + den: self.den.clone(), + } + } + + pub(crate) fn mul(&self, o: &BigQ) -> BigQ { + BigQ::new(&self.num * &o.num, &self.den * &o.den) + } + + /// `None` when `o` is zero. + pub(crate) fn div(&self, o: &BigQ) -> Option { + (!o.num.is_zero()).then(|| BigQ::new(&self.num * &o.den, &self.den * &o.num)) + } + + pub(crate) fn cmp_value(&self, o: &BigQ) -> Ordering { + (&self.num * &o.den).cmp(&(&o.num * &self.den)) + } + + pub(crate) fn bits(&self) -> u64 { + self.num.bits().max(self.den.bits()) + } + + pub(crate) fn to_f64(&self) -> f64 { + ratio_to_f64(&self.num, &self.den) + } + + /// The exact value of a finite f64. + pub(crate) fn from_f64(x: f64) -> BigQ { + let bits = x.to_bits(); + let biased = ((bits >> 52) & 0x7ff) as i64; + let fraction = bits & ((1u64 << 52) - 1); + let (mant, exp) = match biased { + 0 => (fraction, -1074), + e => (fraction | (1 << 52), e - 1075), + }; + let signed = match x.is_sign_negative() { + true => -BigInt::from(mant), + false => BigInt::from(mant), + }; + match exp >= 0 { + true => BigQ::new(signed << exp as usize, BigInt::one()), + false => BigQ::new(signed, BigInt::one() << (-exp) as usize), + } + } +} + +/// `num / den` rounded to the nearest f64, ties to even, subnormals and overflow included. +pub(crate) fn ratio_to_f64(num: &BigInt, den: &BigInt) -> f64 { + let negative = (num.sign() == Sign::Minus) != (den.sign() == Sign::Minus); + let a = num.magnitude(); + let b = den.magnitude(); + let magnitude = match a.is_zero() { + true => 0.0, + false => magnitude_to_f64(a, b), + }; + match negative { + true => -magnitude, + false => magnitude, + } +} + +fn magnitude_to_f64(a: &BigUint, b: &BigUint) -> f64 { + let shift = 54 - (a.bits() as i64 - b.bits() as i64); + let (n, d) = match shift >= 0 { + true => (a << shift as usize, b.clone()), + false => (a.clone(), b << (-shift) as usize), + }; + let (q, r) = n.div_rem(&d); + let qbits = q.bits() as i64; + let drop = (qbits - 53).max(shift - 1074); + if drop > qbits + 1 { + return 0.0; + } + let kept = (&q >> drop as usize).to_u64().unwrap_or(0); + let round = q.bit(drop as u64 - 1); + let sticky = !r.is_zero() || q.trailing_zeros().is_some_and(|t| (t as i64) < drop - 1); + let up = round && (sticky || kept & 1 == 1); + ldexp(kept + up as u64, drop - shift) +} + +/// `m * 2^e` for `m <= 2^53`, exact unless it overflows to infinity. +fn ldexp(m: u64, e: i64) -> f64 { + let top = e + 63 - m.leading_zeros() as i64; + match (m == 0, top > 1023) { + (true, _) => 0.0, + (false, true) => f64::INFINITY, + (false, false) => { + let half = e / 2; + m as f64 * pow2(half) * pow2(e - half) + } + } +} + +fn pow2(e: i64) -> f64 { + match e { + e if e >= -1022 => f64::from_bits(((e + 1023) as u64) << 52), + e => f64::from_bits(1u64 << (e + 1074)), + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn rounding_to_f64_reaches_subnormals_ties_and_overflow() { + let one = BigInt::one(); + let min_sub = f64::from_bits(1); + assert_eq!(ratio_to_f64(&one, &(BigInt::one() << 1074)), min_sub); + assert_eq!(ratio_to_f64(&one, &(BigInt::one() << 1075)), 0.0); + let three = BigInt::from(3); + assert_eq!(ratio_to_f64(&three, &(BigInt::one() << 1076)), min_sub); + assert_eq!(ratio_to_f64(&one, &BigInt::from(3)), 1.0 / 3.0); + let tie = (BigInt::one() << 53) + 1; + assert_eq!(ratio_to_f64(&tie, &one), 9007199254740992.0); + let big = BigInt::one() << 1024; + assert_eq!(ratio_to_f64(&big, &one), f64::INFINITY); + assert_eq!(ratio_to_f64(&-big, &one), f64::NEG_INFINITY); + let max = BigQ::from_f64(f64::MAX); + assert_eq!(max.to_f64(), f64::MAX); + } + + #[test] + fn every_f64_converts_to_its_exact_value_and_back() { + for x in [0.1, -2.5e-310, 1e300, f64::MIN_POSITIVE, 123456.789, -7.0] { + assert_eq!(BigQ::from_f64(x).to_f64().to_bits(), x.to_bits(), "{x}"); + } + assert!(BigQ::from_f64(-0.0).num.is_zero()); + } +} diff --git a/crates/numtower/src/round.rs b/crates/numtower/src/round.rs new file mode 100644 index 00000000..ec6a060e --- /dev/null +++ b/crates/numtower/src/round.rs @@ -0,0 +1,85 @@ +//! Rounding any value to `places` for display, exactly (clause 12, "Format"). +//! An inexact value is rounded over its shortest round-trip spelling (improvise-felo.9). + +use num_bigint::BigInt; +use num_integer::Integer; +use num_traits::Signed; + +use crate::lang::{Lexeme, shortest_round_trip}; +use crate::num::Num; +use crate::ratio::BigQ; + +/// The rounding mode is the caller's policy, not the crate's. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum Rounding { + HalfAwayFromZero, + HalfEven, +} + +impl Num { + /// Plain decimal text with exactly `places` fractional digits; `None` for a non-finite value. + pub fn round_to_places(&self, places: u8, mode: Rounding) -> Option { + let q = match self.exact_q() { + Some(q) => q, + None => Lexeme::read(&shortest_round_trip(self.to_f64())?)?.exact_value(), + }; + Some(render(round(&q, places, mode), places)) + } +} + +/// `q * 10^places`, rounded to an integer by `mode`. +fn round(q: &BigQ, places: u8, mode: Rounding) -> BigInt { + let scaled = q.num.abs() * BigInt::from(10).pow(places as u32); + let (whole, rest) = scaled.div_rem(&q.den); + let twice: BigInt = rest * 2; + let up = match twice.cmp(&q.den) { + std::cmp::Ordering::Greater => true, + std::cmp::Ordering::Equal => mode == Rounding::HalfAwayFromZero || whole.is_odd(), + std::cmp::Ordering::Less => false, + }; + let magnitude = whole + BigInt::from(up as u8); + match q.num.is_negative() { + true => -magnitude, + false => magnitude, + } +} + +fn render(units: BigInt, places: u8) -> String { + let digits = units.abs().to_string(); + let width = places as usize + 1; + let padded = format!("{digits:0>width$}"); + let (int, frac) = padded.split_at(padded.len() - places as usize); + let sign = if units.is_negative() { "-" } else { "" }; + match places { + 0 => format!("{sign}{int}"), + _ => format!("{sign}{int}.{frac}"), + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::lang::{Typing, parse_number}; + + fn typed(s: &str) -> Num { + parse_number(s, Typing::Inferred).unwrap() + } + + #[test] + fn a_typed_value_rounds_as_it_was_typed() { + let at2 = |n: &Num| n.round_to_places(2, Rounding::HalfAwayFromZero).unwrap(); + for (text, want) in [("1.005", "1.01"), ("1.015", "1.02"), ("0.285", "0.29"), ("-2.675", "-2.68")] { + assert_eq!(at2(&typed(text)), want, "{text}"); + assert_eq!(at2(&Num::inexact(text.parse().unwrap())), want, "inexact {text}"); + } + assert_eq!(at2(&Num::ratio(10, 3).unwrap()), "3.33"); + assert_eq!(at2(&Num::ratio(2, 3).unwrap()), "0.67"); + assert_eq!(at2(&typed("-0.001")), "0.00"); + let even = typed("0.125").round_to_places(2, Rounding::HalfEven).unwrap(); + assert_eq!(even, "0.12"); + assert_eq!(typed("12.5").round_to_places(0, Rounding::HalfEven).unwrap(), "12"); + assert_eq!(Num::inexact(f64::NAN).round_to_places(2, Rounding::HalfEven), None); + let wide = Num::ratio(1, 3).unwrap().round_to_places(40, Rounding::HalfEven).unwrap(); + assert_eq!(wide, format!("0.{}", "3".repeat(40))); + } +} diff --git a/crates/numtower/src/sum.rs b/crates/numtower/src/sum.rs new file mode 100644 index 00000000..f9b27b11 --- /dev/null +++ b/crates/numtower/src/sum.rs @@ -0,0 +1,258 @@ +//! Order-free, retractable sums (clause 8): ADR 009's f64 `ExactSum`, and `NumSum`, +//! which adds the exact partitions. The design: `context/code/numtower.md`, "Sums". + +use std::collections::BTreeMap; + +use num_bigint::{BigInt, BigUint}; +use num_integer::Integer; +use num_traits::{One, Zero}; + +use crate::fixed::{I256, I320}; +use crate::num::{MAX_SCALE, Num, NumRef, pow10}; +use crate::ratio::{BigQ, ratio_to_f64}; + +/// The most non-decimal denominator groups a sum reads exactly (question 16: K = 8). +pub const K: usize = 8; + +const LIMB_BITS: u32 = 32; +/// A finite f64 is `mant * 2^(shift - 1074)`, `shift <= 2045`; weights add 60 bits. +const LIMBS: usize = 70; +const CHUNK: u64 = 1 << 30; + +/// An exact sum of weighted f64s, rounded once when read (ADR 009). +#[derive(Clone, Debug, PartialEq, Eq)] +pub struct ExactSum { + limbs: [i64; LIMBS], + pending: u64, + nan: i128, + pos_inf: i128, + neg_inf: i128, +} + +impl Default for ExactSum { + fn default() -> Self { + ExactSum { + limbs: [0; LIMBS], + pending: 0, + nan: 0, + pos_inf: 0, + neg_inf: 0, + } + } +} + +impl ExactSum { + pub fn new() -> ExactSum { + ExactSum::default() + } + + /// Adds `x` with a Z-set `weight`; the weight is split into base-2^30 digits. + pub fn add(&mut self, x: f64, weight: i64) { + match (x.is_nan(), x.is_infinite(), x > 0.0) { + (true, _, _) => self.nan += weight as i128, + (false, true, true) => self.pos_inf += weight as i128, + (false, true, false) => self.neg_inf += weight as i128, + (false, false, _) => { + let magnitude = weight.unsigned_abs(); + for j in 0..3 { + let digit = (magnitude >> (30 * j)) & (CHUNK - 1); + if digit != 0 { + let signed = if weight < 0 { -(digit as i64) } else { digit as i64 }; + self.add_digit(x, signed, 30 * j); + } + } + } + } + } + + fn add_digit(&mut self, x: f64, weight: i64, extra: u32) { + if self.pending + weight.unsigned_abs() > CHUNK { + propagate(&mut self.limbs); + self.pending = 0; + } + let bits = x.to_bits(); + let biased = ((bits >> 52) & 0x7ff) as u32; + let fraction = bits & ((1 << 52) - 1); + let (mant, shift) = match biased { + 0 => (fraction, 0), + e => (fraction | (1 << 52), e - 1), + }; + let sign = if x.is_sign_negative() { -weight } else { weight }; + let at = shift + extra; + let wide = u128::from(mant) << (at % LIMB_BITS); + let k = (at / LIMB_BITS) as usize; + for i in 0..3 { + let part = ((wide >> (LIMB_BITS * i as u32)) & 0xffff_ffff) as i64; + self.limbs[k + i] += sign * part; + } + self.pending += weight.unsigned_abs(); + } + + /// Adds another accumulator's contents. + pub fn absorb(&mut self, o: &ExactSum) { + let mut theirs = o.limbs; + propagate(&mut theirs); + propagate(&mut self.limbs); + for (mine, t) in self.limbs.iter_mut().zip(theirs) { + *mine += t; + } + propagate(&mut self.limbs); + self.pending = 0; + self.nan += o.nan; + self.pos_inf += o.pos_inf; + self.neg_inf += o.neg_inf; + } + + /// The exact sum rounded to the nearest f64, ties to even. + pub fn value(&self) -> f64 { + match (self.nan != 0, self.pos_inf != 0, self.neg_inf != 0) { + (true, _, _) | (false, true, true) => f64::NAN, + (false, true, false) => f64::INFINITY, + (false, false, true) => f64::NEG_INFINITY, + (false, false, false) => ratio_to_f64(&self.total(), &(BigInt::one() << 1074)), + } + } + + /// The finite part as an integer in units of 2^-1074. + fn total(&self) -> BigInt { + let mut limbs = self.limbs; + propagate(&mut limbs); + let digits: Vec = limbs[..LIMBS - 1].iter().map(|&l| l as u32).collect(); + let low = BigInt::from(BigUint::from_slice(&digits)); + (BigInt::from(limbs[LIMBS - 1]) << (LIMB_BITS as usize * (LIMBS - 1))) + low + } +} + +/// Carry each limb into the next, leaving every limb but the top in [0, 2^32). +fn propagate(limbs: &mut [i64; LIMBS]) { + for k in 0..LIMBS - 1 { + let carry = limbs[k] >> LIMB_BITS; + limbs[k] -= carry << LIMB_BITS; + limbs[k + 1] += carry; + } +} + +/// A `SUM` over any `Num`s: exact when every term is exact and at most K groups remain. +#[derive(Clone, Debug, Default)] +pub struct NumSum { + decimals: I320, + groups: BTreeMap, + group_floats: ExactSum, + inexact: ExactSum, + inexact_weight: i128, +} + +impl NumSum { + pub fn new() -> NumSum { + NumSum::default() + } + + pub fn insert(&mut self, x: &Num) { + self.add(x, 1); + } + + pub fn retract(&mut self, x: &Num) { + self.add(x, -1); + } + + /// Adds `x` with a Z-set `weight`. Exact and associative: the order cannot matter. + pub fn add(&mut self, x: &Num, weight: i64) { + match (x.kind(), x.exact_parts()) { + (NumRef::Inexact(f), _) => { + self.inexact.add(f, weight); + self.inexact_weight += weight as i128; + } + (_, Some((n, d))) => match at_scale_28(d) { + Some(unit) => self.add_scaled(n, unit, weight), + None => self.add_to_group(n, d, weight), + }, + (_, None) => unreachable!("an exact value has parts"), + } + } + + fn add_scaled(&mut self, n: i128, unit: u128, weight: i64) { + let term = I320::from_i128(n) + .checked_mul_u128(unit) + .and_then(|t| t.checked_mul_u128(weight.unsigned_abs() as u128)) + .and_then(|t| if weight < 0 { t.checked_neg() } else { Some(t) }) + .expect("a weighted term is below 2^284"); + self.decimals = self.decimals.checked_add(&term).expect("320 bits holds 2^36 terms"); + } + + fn add_to_group(&mut self, n: i128, d: i128, weight: i64) { + let old = self.groups.get(&d).copied().unwrap_or(I256::ZERO); + let new = old + .checked_add(&I256::mul_i128(n, weight as i128)) + .expect("256 bits holds 2^64 terms"); + let value = |num: I256| ratio_to_f64(&num.to_bigint(), &BigInt::from(d)); + if !old.is_zero() { + self.group_floats.add(value(old), -1); + } + match new.is_zero() { + true => { + self.groups.remove(&d); + } + false => { + self.group_floats.add(value(new), 1); + self.groups.insert(d, new); + } + } + } + + /// The current non-decimal denominator groups. + pub fn groups(&self) -> usize { + self.groups.len() + } + + /// One read: exact over at most K groups, else the f64 accumulators (clause 8). + pub fn read(&self) -> Num { + match (self.inexact_weight != 0, self.groups.len() > K) { + (false, false) => Num::settle(&self.exact_total()), + (_, true) => { + let mut acc = self.inexact.clone(); + acc.absorb(&self.group_floats); + acc.add(self.decimal_total().to_f64(), 1); + Num::inexact(acc.value()) + } + (true, false) => { + let mut acc = self.inexact.clone(); + acc.add(self.exact_total().to_f64(), 1); + Num::inexact(acc.value()) + } + } + } + + fn decimal_total(&self) -> BigQ { + BigQ::new(self.decimals.to_bigint(), BigInt::from(pow10(MAX_SCALE))) + } + + fn exact_total(&self) -> BigQ { + let decimal = (!self.decimals.is_zero()) + .then(|| (self.decimals.to_bigint(), BigInt::from(pow10(MAX_SCALE)))); + let terms: Vec<(BigInt, BigInt)> = decimal + .into_iter() + .chain(self.groups.iter().map(|(d, n)| (n.to_bigint(), BigInt::from(*d)))) + .collect(); + let (num, den) = combine_over_lcm(&terms); + BigQ::new(num, den) + } +} + +/// `Σ n/d` over the lcm of the denominators, not yet reduced: the read reduces once. +pub(crate) fn combine_over_lcm(terms: &[(BigInt, BigInt)]) -> (BigInt, BigInt) { + let lcm = terms.iter().fold(BigInt::one(), |l, (_, d)| l.lcm(d)); + let num = terms.iter().fold(BigInt::zero(), |s, (n, d)| s + n * (&lcm / d)); + (num, lcm) +} + +/// `10^28 / d` when `d = 2^p 5^q` with `max(p, q) <= 28`: the term joins the decimal partition. +fn at_scale_28(d: i128) -> Option { + let p = d.trailing_zeros(); + let mut rest = d >> p; + let mut q = 0; + while rest % 5 == 0 { + rest /= 5; + q += 1; + } + (rest == 1 && p.max(q) <= MAX_SCALE).then(|| (pow10(MAX_SCALE) / d) as u128) +}