Merge pull request #337 from fiddlerwoaroof/adr005/numtower
feat(numtower): ADR 005 step 1, the numeric tower crate (improvise-59nw.1)
This commit is contained in:
Ed Langley
GitHub
@@ -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<i128>), 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
|
||||
@@ -0,0 +1,187 @@
|
||||
//! ADR 005 migration step 1's benchmark: `cargo run --release -p numtower --example bench`.
|
||||
//! A microbenchmark, not the evaluator; what it does and does not show is in `context/code/numtower.md`.
|
||||
|
||||
use std::hint::black_box;
|
||||
use std::time::Instant;
|
||||
|
||||
use numtower::{Decimal, ExactSum, Num, NumSum, Rounding, Typing, parse_number};
|
||||
|
||||
struct SplitMix64(u64);
|
||||
|
||||
impl SplitMix64 {
|
||||
fn next(&mut self) -> u64 {
|
||||
self.0 = self.0.wrapping_add(0x9E37_79B9_7F4A_7C15);
|
||||
let mut z = self.0;
|
||||
z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
|
||||
z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
|
||||
z ^ (z >> 31)
|
||||
}
|
||||
}
|
||||
|
||||
/// Typed two-place values from 0.01 to 999,999.99, as the ADR's appendix draws them.
|
||||
fn typed_texts(seed: u64, n: usize) -> Vec<String> {
|
||||
let mut rng = SplitMix64(seed);
|
||||
(0..n)
|
||||
.map(|_| {
|
||||
let cents = 1 + rng.next() % 99_999_999;
|
||||
format!("{}.{:02}", cents / 100, cents % 100)
|
||||
})
|
||||
.collect()
|
||||
}
|
||||
|
||||
fn nums(texts: &[String]) -> Vec<Num> {
|
||||
texts
|
||||
.iter()
|
||||
.map(|t| parse_number(t, Typing::Inferred).unwrap())
|
||||
.collect()
|
||||
}
|
||||
|
||||
/// The best of `runs`, in nanoseconds per item.
|
||||
fn best(runs: usize, items: usize, mut f: impl FnMut()) -> f64 {
|
||||
(0..runs)
|
||||
.map(|_| {
|
||||
let t = Instant::now();
|
||||
f();
|
||||
t.elapsed().as_nanos() as f64 / items as f64
|
||||
})
|
||||
.fold(f64::INFINITY, f64::min)
|
||||
}
|
||||
|
||||
fn sums() {
|
||||
println!("SUM over typed two-place values, ns per element, best of 20");
|
||||
for n in [18_000, 90_000] {
|
||||
let texts = typed_texts(1, n);
|
||||
let floats: Vec<f64> = texts.iter().map(|t| t.parse().unwrap()).collect();
|
||||
let decimals: Vec<Decimal> = texts.iter().map(|t| t.parse().unwrap()).collect();
|
||||
let values = nums(&texts);
|
||||
let f64_fold = best(20, n, || {
|
||||
black_box(black_box(&floats).iter().sum::<f64>());
|
||||
});
|
||||
let decimal_fold = best(20, n, || {
|
||||
let s = black_box(&decimals)
|
||||
.iter()
|
||||
.fold(Decimal::ZERO, |a, b| a.checked_add(*b).unwrap());
|
||||
black_box(s);
|
||||
});
|
||||
let num_fold = best(20, n, || {
|
||||
let zero = Num::exact(Decimal::ZERO);
|
||||
black_box(black_box(&values).iter().fold(zero, |a, b| &a + b));
|
||||
});
|
||||
let accumulator = best(20, n, || {
|
||||
let mut s = NumSum::new();
|
||||
black_box(&values).iter().for_each(|v| s.insert(v));
|
||||
black_box(s.read());
|
||||
});
|
||||
println!(
|
||||
" {n:>6} cells: f64 {f64_fold:.2}, rust_decimal {decimal_fold:.2}, Num + fold {num_fold:.2}, NumSum {accumulator:.2}"
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
fn operations() {
|
||||
println!("Per operation over 18,000 typed pairs, ns, best of 20");
|
||||
let (a, b) = (
|
||||
nums(&typed_texts(21, 18_000)),
|
||||
nums(&typed_texts(22, 18_000)),
|
||||
);
|
||||
let decimal =
|
||||
|v: &[Num]| -> Vec<Decimal> { v.iter().map(|n| n.to_string().parse().unwrap()).collect() };
|
||||
let (da, db) = (decimal(&a), decimal(&b));
|
||||
let pairs = a.len();
|
||||
let rd_add = best(20, pairs, || {
|
||||
da.iter().zip(&db).for_each(|(x, y)| {
|
||||
black_box(x.checked_add(*y));
|
||||
})
|
||||
});
|
||||
let rd_mul = best(20, pairs, || {
|
||||
da.iter().zip(&db).for_each(|(x, y)| {
|
||||
black_box(x.checked_mul(*y));
|
||||
})
|
||||
});
|
||||
let rd_div = best(20, pairs, || {
|
||||
da.iter().zip(&db).for_each(|(x, y)| {
|
||||
black_box(x.checked_div(*y));
|
||||
})
|
||||
});
|
||||
let add = best(20, pairs, || {
|
||||
a.iter().zip(&b).for_each(|(x, y)| drop(black_box(x + y)))
|
||||
});
|
||||
let mul = best(20, pairs, || {
|
||||
a.iter().zip(&b).for_each(|(x, y)| drop(black_box(x * y)))
|
||||
});
|
||||
let div = best(20, pairs, || {
|
||||
a.iter()
|
||||
.zip(&b)
|
||||
.for_each(|(x, y)| drop(black_box(x.checked_div(y))))
|
||||
});
|
||||
let margin = best(20, pairs, || {
|
||||
a.iter()
|
||||
.zip(&b)
|
||||
.for_each(|(rev, cost)| drop(black_box((rev - cost).checked_div(rev))));
|
||||
});
|
||||
let to_f64 = best(20, pairs, || {
|
||||
a.iter().for_each(|x| {
|
||||
black_box(x.to_f64());
|
||||
})
|
||||
});
|
||||
let places = best(20, pairs, || {
|
||||
a.iter()
|
||||
.for_each(|x| drop(black_box(x.round_to_places(2, Rounding::HalfAwayFromZero))));
|
||||
});
|
||||
println!(" rust_decimal checked: add {rd_add:.1}, mul {rd_mul:.1}, div {rd_div:.1}");
|
||||
println!(" Num: add {add:.1}, mul {mul:.1}, div {div:.1}");
|
||||
println!(
|
||||
" Num (Rev - Cost) / Rev {margin:.1}, to_f64 {to_f64:.1}, round_to_places(2) {places:.1}"
|
||||
);
|
||||
}
|
||||
|
||||
/// K groups of near-`i128` denominators: the ADR's worst case for the exact read.
|
||||
fn grouped_read() {
|
||||
println!("NumSum read over K distinct 126-bit denominators, us per read, best of 200");
|
||||
for k in [2usize, 8] {
|
||||
let mut rng = SplitMix64(0x5EED);
|
||||
let mut sum = NumSum::new();
|
||||
for _ in 0..k {
|
||||
let den = ((rng.next() as i128) << 62 | rng.next() as i128 >> 2) | 1 | (1 << 125);
|
||||
let num = (rng.next() as i128) << 60 ^ rng.next() as i128;
|
||||
let term = Num::ratio(num, den).unwrap();
|
||||
sum.add(&term, i64::MAX);
|
||||
sum.add(&term, i64::MAX);
|
||||
}
|
||||
assert_eq!(sum.groups(), k);
|
||||
let read = best(200, 1, || drop(black_box(sum.read())));
|
||||
println!(
|
||||
" K = {k}: {:.1} us (exact: {})",
|
||||
read / 1000.0,
|
||||
sum.read().is_exact()
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
fn exact_sum() {
|
||||
println!("ExactSum (f64 superaccumulator), best of 20");
|
||||
let floats: Vec<f64> = typed_texts(1, 1_000_000)
|
||||
.iter()
|
||||
.map(|t| t.parse().unwrap())
|
||||
.collect();
|
||||
let mut acc = ExactSum::new();
|
||||
let insert = best(20, floats.len(), || {
|
||||
acc = ExactSum::new();
|
||||
floats.iter().for_each(|x| acc.add(*x, 1));
|
||||
});
|
||||
let read = best(20, 1, || {
|
||||
black_box(acc.value());
|
||||
});
|
||||
println!(
|
||||
" {insert:.2} ns per insert, {:.2} us per read",
|
||||
read / 1000.0
|
||||
);
|
||||
}
|
||||
|
||||
fn main() {
|
||||
println!("numtower bench (release); name the machine when you quote these");
|
||||
sums();
|
||||
operations();
|
||||
grouped_read();
|
||||
exact_sum();
|
||||
}
|
||||
@@ -0,0 +1,445 @@
|
||||
//! 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<BigQ>),
|
||||
}
|
||||
|
||||
/// The divisor was exactly zero, exact or inexact (clause 5).
|
||||
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
|
||||
pub struct DivisionByZero;
|
||||
|
||||
impl From<Num> 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<Interim, DivisionByZero> {
|
||||
divide(self.operand(), o.operand())
|
||||
}
|
||||
|
||||
pub fn pow(&self, o: &Interim) -> Result<Interim, DivisionByZero> {
|
||||
power(self.operand(), o.operand())
|
||||
}
|
||||
|
||||
/// Exactly zero; a wide value never is, because zero fits a stored arm.
|
||||
pub fn is_zero(&self) -> bool {
|
||||
match &self.0 {
|
||||
Inner::Stored(n) => n.is_zero(),
|
||||
Inner::Wide(_) => false,
|
||||
}
|
||||
}
|
||||
|
||||
/// Clause 6's exact order, before the stored bound; `None` only when either is NaN.
|
||||
pub fn cmp_exact(&self, o: &Interim) -> Option<std::cmp::Ordering> {
|
||||
match (&self.0, &o.0) {
|
||||
(Inner::Stored(a), Inner::Stored(b)) => a.cmp_exact(b),
|
||||
_ => match (self.operand().big(), o.operand().big()) {
|
||||
(Some(a), Some(b)) => Some(a.cmp_value(&b)),
|
||||
(Some(a), None) => crate::cmp::exact_against_f64(&a, o.operand().float()),
|
||||
(None, Some(b)) => {
|
||||
crate::cmp::exact_against_f64(&b, self.operand().float()).map(|c| c.reverse())
|
||||
}
|
||||
(None, None) => self.operand().float().partial_cmp(&o.operand().float()),
|
||||
},
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
impl Num {
|
||||
pub fn checked_div(&self, o: &Num) -> Result<Num, DivisionByZero> {
|
||||
divide(Operand::of(self), Operand::of(o)).map(Interim::store)
|
||||
}
|
||||
|
||||
pub fn pow(&self, o: &Num) -> Result<Num, DivisionByZero> {
|
||||
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),
|
||||
}
|
||||
}
|
||||
|
||||
/// The `i128` fast path's operands; `i128::MIN` cannot be negated, so it takes the bigint path.
|
||||
fn parts(self) -> Option<(i128, i128)> {
|
||||
match self {
|
||||
Operand::Dec(d) => Some((d.mantissa(), pow10(d.scale()))),
|
||||
Operand::Small(n, d) => (n != i128::MIN).then_some((n, d)),
|
||||
Operand::Big(_) | Operand::Float(_) => None,
|
||||
}
|
||||
}
|
||||
|
||||
fn big(self) -> Option<BigQ> {
|
||||
match self {
|
||||
Operand::Dec(d) => Some(BigQ::from_i128(d.mantissa(), pow10(d.scale()))),
|
||||
Operand::Small(n, d) => Some(BigQ::from_i128(n, d)),
|
||||
Operand::Big(q) => Some(q.clone()),
|
||||
Operand::Float(_) => None,
|
||||
}
|
||||
}
|
||||
|
||||
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) => n
|
||||
.checked_neg()
|
||||
.and_then(|m| small_ratio(m, 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<Interim, DivisionByZero> {
|
||||
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<BigInt> {
|
||||
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<Interim, DivisionByZero> {
|
||||
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<Interim, DivisionByZero> {
|
||||
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 = match n.is_negative() {
|
||||
true => 0,
|
||||
false => 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);
|
||||
let odd_negative = q.num.is_negative() && k.is_odd();
|
||||
match (too_big, k.to_u32()) {
|
||||
(false, Some(k)) => Ok(exact_integer_power(&q, k, negative, target)),
|
||||
_ => Ok(inexact(signed_powf(a.float(), n, odd_negative))),
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// `powf` on the magnitude, the sign taken from the exact exponent's parity.
|
||||
fn signed_powf(base: f64, n: &BigInt, odd_negative: bool) -> f64 {
|
||||
let magnitude = base.abs().powf(n.to_f64().unwrap_or(f64::INFINITY));
|
||||
match odd_negative {
|
||||
true => -magnitude,
|
||||
false => magnitude,
|
||||
}
|
||||
}
|
||||
|
||||
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,
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,107 @@
|
||||
//! 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<Ordering> {
|
||||
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<Ordering> {
|
||||
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<Ordering> {
|
||||
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),
|
||||
},
|
||||
}
|
||||
}
|
||||
|
||||
/// A wide exact value against an f64, with no fast path: the slow, total form of `mixed`.
|
||||
pub(crate) fn exact_against_f64(q: &BigQ, x: f64) -> Option<Ordering> {
|
||||
match (x.is_nan(), x.is_infinite(), x > 0.0) {
|
||||
(true, _, _) => None,
|
||||
(false, true, true) => Some(Ordering::Less),
|
||||
(false, true, false) => Some(Ordering::Greater),
|
||||
(false, false, _) => Some(q.cmp_value(&BigQ::from_f64(x))),
|
||||
}
|
||||
}
|
||||
|
||||
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,
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,234 @@
|
||||
//! 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<const N: usize> {
|
||||
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<const N: usize> Fixed<N> {
|
||||
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 (o, (a, b)) in out.iter_mut().zip(self.limbs.iter().zip(&other.limbs)) {
|
||||
let (s, c1) = a.overflowing_add(*b);
|
||||
let (s, c2) = s.overflowing_add(carry as u64);
|
||||
*o = s;
|
||||
carry = c1 || c2;
|
||||
}
|
||||
Fixed { limbs: out }
|
||||
}
|
||||
|
||||
pub(crate) fn checked_add(&self, other: &Self) -> Option<Self> {
|
||||
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<Self> {
|
||||
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<Self> {
|
||||
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<Self> {
|
||||
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, &limb) in mag.iter().enumerate() {
|
||||
let idx = i + j;
|
||||
let prod = limb 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<i128> {
|
||||
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<u32> = 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<const N: usize> Default for Fixed<N> {
|
||||
fn default() -> Self {
|
||||
Self::ZERO
|
||||
}
|
||||
}
|
||||
|
||||
impl<const N: usize> Ord for Fixed<N> {
|
||||
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<const N: usize> PartialOrd for Fixed<N> {
|
||||
fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
|
||||
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);
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,192 @@
|
||||
//! 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 && forced_fraction(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<Num> {
|
||||
match (gate(text), typing) {
|
||||
(true, _) => Lexeme::read(text).map(|lx| lx.value(text)),
|
||||
(false, Typing::Forced) => forced_fraction(text),
|
||||
(false, Typing::Inferred) => None,
|
||||
}
|
||||
}
|
||||
|
||||
/// A forced `n/d`; one that rounds to infinity is text, as a finite lexeme never reads as one.
|
||||
fn forced_fraction(text: &str) -> Option<Num> {
|
||||
let n = Num::settle(&fraction_lexeme(text)?);
|
||||
n.to_f64().is_finite().then_some(n)
|
||||
}
|
||||
|
||||
fn gate(text: &str) -> bool {
|
||||
text.parse::<f64>().is_ok_and(f64::is_finite)
|
||||
}
|
||||
|
||||
/// `[+-]digits/digits` with a nonzero denominator, as an exact value.
|
||||
fn fraction_lexeme(text: &str) -> Option<BigQ> {
|
||||
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<Lexeme> {
|
||||
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<Num> {
|
||||
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<i64> {
|
||||
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::<i64>().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<String> {
|
||||
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<Num> {
|
||||
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:?}"),
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -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;
|
||||
@@ -0,0 +1,198 @@
|
||||
//! 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<Ratio<i128>>),
|
||||
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<i128>),
|
||||
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<Num> {
|
||||
(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<BigQ> {
|
||||
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<Num> {
|
||||
decimal_of(q).or_else(|| {
|
||||
let num = q.num.to_i128()?;
|
||||
let den = q.den.to_i128()?;
|
||||
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<Num> {
|
||||
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))
|
||||
})
|
||||
}
|
||||
|
||||
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<Num> {
|
||||
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<Num> {
|
||||
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::<Num>(), 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());
|
||||
}
|
||||
}
|
||||
@@ -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<BigQ> {
|
||||
(!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());
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,107 @@
|
||||
//! 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<String> {
|
||||
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)));
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,396 @@
|
||||
//! 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)]
|
||||
pub struct ExactSum {
|
||||
limbs: [i64; LIMBS],
|
||||
pending: u64,
|
||||
nan: i128,
|
||||
pos_inf: i128,
|
||||
neg_inf: i128,
|
||||
}
|
||||
|
||||
/// Equal sums, not equal representations: carries propagated, `pending` ignored.
|
||||
impl PartialEq for ExactSum {
|
||||
fn eq(&self, o: &ExactSum) -> bool {
|
||||
let (mut mine, mut theirs) = (self.limbs, o.limbs);
|
||||
propagate(&mut mine);
|
||||
propagate(&mut theirs);
|
||||
let counts = |s: &ExactSum| (s.nan, s.pos_inf, s.neg_inf);
|
||||
mine == theirs && counts(self) == counts(o)
|
||||
}
|
||||
}
|
||||
|
||||
impl Eq for ExactSum {}
|
||||
|
||||
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<u32> = 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<i128, I256>,
|
||||
group_floats: ExactSum,
|
||||
inexact: ExactSum,
|
||||
/// Net weight per f64 bit pattern: a read is inexact while any term is present.
|
||||
inexact_terms: BTreeMap<u64, 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() {
|
||||
NumRef::Inexact(f) => {
|
||||
self.inexact.add(f, weight);
|
||||
let net = self.inexact_terms.entry(f.to_bits()).or_insert(0);
|
||||
*net += weight as i128;
|
||||
if *net == 0 {
|
||||
self.inexact_terms.remove(&f.to_bits());
|
||||
}
|
||||
}
|
||||
NumRef::Exact(d) => {
|
||||
let unit = pow10(MAX_SCALE - d.scale()) as u128;
|
||||
self.add_scaled(d.mantissa(), unit, weight);
|
||||
}
|
||||
NumRef::Fraction(r) => match at_scale_28(*r.denom()) {
|
||||
Some(unit) => self.add_scaled(*r.numer(), unit, weight),
|
||||
None => self.add_to_group(*r.numer(), *r.denom(), weight),
|
||||
},
|
||||
}
|
||||
}
|
||||
|
||||
fn add_scaled(&mut self, n: i128, unit: u128, weight: i64) {
|
||||
let magnitude = weight.unsigned_abs() as u128;
|
||||
let scaled = match unit.checked_mul(magnitude) {
|
||||
Some(both) => I320::from_i128(n).checked_mul_u128(both),
|
||||
None => I320::from_i128(n)
|
||||
.checked_mul_u128(unit)
|
||||
.and_then(|t| t.checked_mul_u128(magnitude)),
|
||||
};
|
||||
let term = scaled
|
||||
.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_terms.is_empty(), 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<u128> {
|
||||
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)
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
use num_rational::BigRational;
|
||||
use proptest::prelude::*;
|
||||
use proptest::test_runner::{Config, RngAlgorithm, TestRng, TestRunner};
|
||||
|
||||
type Combine = fn(&[(BigInt, BigInt)]) -> (BigInt, BigInt);
|
||||
|
||||
/// Up to K groups of wide numerators whose denominators share factors.
|
||||
fn groups() -> impl Strategy<Value = Vec<(BigInt, BigInt)>> {
|
||||
let shared = prop::sample::select(vec![1i128, 3, 7 * 11, 1 << 20, 3i128.pow(40)]);
|
||||
let group =
|
||||
(any::<i128>(), any::<i64>(), shared, 1i128..(1 << 60)).prop_map(|(n, w, s, d)| {
|
||||
(
|
||||
BigInt::from(n) * BigInt::from(w),
|
||||
BigInt::from(s) * BigInt::from(d),
|
||||
)
|
||||
});
|
||||
prop::collection::vec(group, 1..=K)
|
||||
}
|
||||
|
||||
fn check_grouped_read(combine: Combine, deterministic: bool) -> Result<(), String> {
|
||||
let config = Config {
|
||||
cases: 256,
|
||||
failure_persistence: None,
|
||||
..Config::default()
|
||||
};
|
||||
let mut runner = match deterministic {
|
||||
true => {
|
||||
TestRunner::new_with_rng(config, TestRng::deterministic_rng(RngAlgorithm::ChaCha))
|
||||
}
|
||||
false => TestRunner::new(config),
|
||||
};
|
||||
let result = runner.run(&groups(), |terms| {
|
||||
let (num, den) = combine(&terms);
|
||||
let lcm = terms.iter().fold(BigInt::one(), |l, (_, d)| l.lcm(d));
|
||||
prop_assert_eq!(&den, &lcm, "the read did not combine over the lcm");
|
||||
let oracle: BigRational = terms
|
||||
.iter()
|
||||
.map(|(n, d)| BigRational::new(n.clone(), d.clone()))
|
||||
.sum();
|
||||
let read = BigQ::new(num, den);
|
||||
let got = BigRational::new(read.num.clone(), read.den.clone());
|
||||
prop_assert_eq!(got, oracle, "the read is not the sum");
|
||||
prop_assert!(read.num.gcd(&read.den).is_one(), "the read is not reduced");
|
||||
Ok(())
|
||||
});
|
||||
result.map_err(|e| e.to_string())
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn the_grouped_read_combines_over_the_lcm_and_reduces_once() {
|
||||
check_grouped_read(combine_over_lcm, false).unwrap();
|
||||
}
|
||||
|
||||
fn over_the_product(terms: &[(BigInt, BigInt)]) -> (BigInt, BigInt) {
|
||||
let den = terms.iter().fold(BigInt::one(), |p, (_, d)| p * d);
|
||||
let num = terms
|
||||
.iter()
|
||||
.fold(BigInt::zero(), |s, (n, d)| s + n * (&den / d));
|
||||
(num, den)
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn the_grouped_read_combines_over_the_lcm_and_reduces_once_catches_the_product() {
|
||||
let err = check_grouped_read(over_the_product, true).expect_err("the rival passed");
|
||||
assert!(
|
||||
err.contains("the read did not combine over the lcm"),
|
||||
"{err}"
|
||||
);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn a_fraction_with_a_decimal_denominator_joins_the_decimal_partition() {
|
||||
let mut s = NumSum::new();
|
||||
s.insert(&Num::ratio(10i128.pow(35), 1).unwrap());
|
||||
s.insert(&Num::ratio(3, 1 << 20).unwrap());
|
||||
assert_eq!(s.groups(), 0);
|
||||
s.insert(&Num::ratio(1, 1i128 << 40).unwrap());
|
||||
s.insert(&Num::ratio(1, 3).unwrap());
|
||||
assert_eq!(
|
||||
s.groups(),
|
||||
2,
|
||||
"2^40 has scale 40, past 28: a group of its own"
|
||||
);
|
||||
s.retract(&Num::ratio(1, 3).unwrap());
|
||||
assert_eq!(
|
||||
s.groups(),
|
||||
1,
|
||||
"a group whose numerator returns to zero is dropped"
|
||||
);
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,189 @@
|
||||
//! The `BigRational` oracle and the generators the property tests share.
|
||||
//! Each property is a checker taking the implementation, so a rival can be run through it.
|
||||
|
||||
#![allow(dead_code)]
|
||||
|
||||
use std::cmp::Ordering;
|
||||
|
||||
use num_bigint::BigInt;
|
||||
use num_rational::BigRational;
|
||||
use num_traits::{One, Signed, ToPrimitive, Zero};
|
||||
use numtower::{Decimal, Num, NumRef};
|
||||
use proptest::prelude::*;
|
||||
use proptest::test_runner::{Config, RngAlgorithm, TestCaseError, TestRng, TestRunner};
|
||||
|
||||
pub fn q(n: &Num) -> Option<BigRational> {
|
||||
match n.kind() {
|
||||
NumRef::Exact(d) => Some(BigRational::new(
|
||||
d.mantissa().into(),
|
||||
BigInt::from(10).pow(d.scale()),
|
||||
)),
|
||||
NumRef::Fraction(r) => Some(BigRational::new((*r.numer()).into(), (*r.denom()).into())),
|
||||
NumRef::Inexact(x) => BigRational::from_float(x),
|
||||
}
|
||||
}
|
||||
|
||||
/// Whether `v` fits a stored exact arm: a decimal at scale <= 28, or an `i128` fraction.
|
||||
pub fn fits(v: &BigRational) -> bool {
|
||||
let limit = BigInt::one() << 96;
|
||||
let decimal = (0..=28u32).any(|s| {
|
||||
let scaled = v * BigRational::from_integer(BigInt::from(10).pow(s));
|
||||
scaled.is_integer() && scaled.to_integer().abs() < limit
|
||||
});
|
||||
let (min, max) = (BigInt::from(i128::MIN), BigInt::from(i128::MAX));
|
||||
decimal || (&min <= v.numer() && v.numer() <= &max && v.denom() <= &max)
|
||||
}
|
||||
|
||||
/// Whether `f` is the round-to-nearest-even f64 of `v`: the review's interval form.
|
||||
pub fn is_nearest(f: f64, v: &BigRational) -> bool {
|
||||
let exact = |x: f64| BigRational::from_float(x).expect("finite");
|
||||
let two = BigRational::from_integer(BigInt::from(2));
|
||||
let half_ulp_past = |x: f64, up: bool| {
|
||||
let (next, inward) = match up {
|
||||
true => (x.next_up(), x.next_down()),
|
||||
false => (x.next_down(), x.next_up()),
|
||||
};
|
||||
match next.is_finite() {
|
||||
true => (exact(x) + exact(next)) / &two,
|
||||
false => exact(x) + (exact(x) - exact(inward)) / &two,
|
||||
}
|
||||
};
|
||||
match f.is_finite() {
|
||||
false => {
|
||||
let edge = half_ulp_past(f64::MAX.copysign(f), f > 0.0);
|
||||
(f > 0.0 && v >= &edge) || (f < 0.0 && v <= &edge)
|
||||
}
|
||||
true => {
|
||||
let lo = half_ulp_past(f, false);
|
||||
let hi = half_ulp_past(f, true);
|
||||
let even = f.to_bits() & 1 == 0;
|
||||
match even {
|
||||
true => &lo <= v && v <= &hi,
|
||||
false => &lo < v && v < &hi,
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// A decimal's text built digit by digit, without `rust_decimal`'s `Display`.
|
||||
pub fn decimal_text(mantissa: i128, scale: u32) -> String {
|
||||
let digits = format!(
|
||||
"{:0>width$}",
|
||||
mantissa.unsigned_abs(),
|
||||
width = scale as usize + 1
|
||||
);
|
||||
let (int, frac) = digits.split_at(digits.len() - scale as usize);
|
||||
let sign = if mantissa < 0 { "-" } else { "" };
|
||||
format!("{sign}{int}.{frac}0")
|
||||
}
|
||||
|
||||
/// The correctly rounded f64 of an exact value, independent of `numtower`'s own.
|
||||
pub fn independent_f64(n: &Num) -> f64 {
|
||||
match n.kind() {
|
||||
NumRef::Exact(d) => decimal_text(d.mantissa(), d.scale()).parse().unwrap(),
|
||||
NumRef::Fraction(_) => q(n).unwrap().to_f64().unwrap(),
|
||||
NumRef::Inexact(x) => x,
|
||||
}
|
||||
}
|
||||
|
||||
/// Same arm, same value and same scale: what a fingerprint would see.
|
||||
pub fn same(a: &Num, b: &Num) -> bool {
|
||||
match (a.kind(), b.kind()) {
|
||||
(NumRef::Exact(x), NumRef::Exact(y)) => {
|
||||
x.mantissa() == y.mantissa() && x.scale() == y.scale()
|
||||
}
|
||||
(NumRef::Fraction(x), NumRef::Fraction(y)) => x == y,
|
||||
(NumRef::Inexact(x), NumRef::Inexact(y)) => x.to_bits() == y.to_bits(),
|
||||
_ => false,
|
||||
}
|
||||
}
|
||||
|
||||
pub fn arm(n: &Num) -> &'static str {
|
||||
match n.kind() {
|
||||
NumRef::Exact(_) => "exact",
|
||||
NumRef::Fraction(_) => "fraction",
|
||||
NumRef::Inexact(_) => "inexact",
|
||||
}
|
||||
}
|
||||
|
||||
pub fn mantissa96() -> impl Strategy<Value = i128> {
|
||||
let max = (1i128 << 96) - 1;
|
||||
prop_oneof![
|
||||
-max..=max,
|
||||
-1_000_000i128..1_000_000,
|
||||
-(10i128.pow(15))..10i128.pow(15),
|
||||
]
|
||||
}
|
||||
|
||||
pub fn decimal() -> impl Strategy<Value = Num> {
|
||||
(mantissa96(), 0u32..=28).prop_map(|(m, s)| Num::exact(Decimal::from_i128_with_scale(m, s)))
|
||||
}
|
||||
|
||||
fn of_width(bits: u32) -> impl Strategy<Value = i128> {
|
||||
(any::<u128>(), any::<bool>()).prop_map(move |(r, neg)| {
|
||||
let v = ((r >> (128 - bits)) | (1u128 << (bits - 1))) as i128;
|
||||
if neg { -v } else { v }
|
||||
})
|
||||
}
|
||||
|
||||
pub fn fraction() -> impl Strategy<Value = Num> {
|
||||
let widths = (1u32..=126, 1u32..=126)
|
||||
.prop_flat_map(|(a, b)| (of_width(a), of_width(b)))
|
||||
.prop_map(|(n, d)| Num::ratio(n, d).unwrap());
|
||||
let edge = (1i128..i128::MAX).prop_map(|d| Num::ratio(i128::MIN, d | 1).unwrap());
|
||||
prop_oneof![8 => widths, 1 => edge]
|
||||
}
|
||||
|
||||
pub fn finite_f64() -> impl Strategy<Value = f64> {
|
||||
use proptest::num::f64::{NEGATIVE, NORMAL, POSITIVE, SUBNORMAL, ZERO};
|
||||
prop_oneof![
|
||||
(POSITIVE | NEGATIVE | NORMAL | SUBNORMAL | ZERO),
|
||||
-1e6f64..1e6,
|
||||
(-1000i64..1000).prop_map(|k| k as f64 / 8.0),
|
||||
]
|
||||
}
|
||||
|
||||
pub fn exact() -> impl Strategy<Value = Num> {
|
||||
prop_oneof![decimal(), fraction()]
|
||||
}
|
||||
|
||||
pub fn any_num() -> impl Strategy<Value = Num> {
|
||||
prop_oneof![
|
||||
3 => exact(),
|
||||
1 => finite_f64().prop_map(Num::inexact),
|
||||
1 => prop_oneof![Just(f64::NAN), Just(f64::INFINITY), Just(f64::NEG_INFINITY)].prop_map(Num::inexact),
|
||||
]
|
||||
}
|
||||
|
||||
/// Runs a property; `deterministic` for a rival, so its failure is not luck.
|
||||
pub fn run<S: Strategy>(
|
||||
strategy: S,
|
||||
cases: u32,
|
||||
deterministic: bool,
|
||||
test: impl Fn(S::Value) -> Result<(), TestCaseError>,
|
||||
) -> Result<(), String> {
|
||||
let config = Config {
|
||||
cases,
|
||||
failure_persistence: None,
|
||||
..Config::default()
|
||||
};
|
||||
let mut runner = match deterministic {
|
||||
true => TestRunner::new_with_rng(config, TestRng::deterministic_rng(RngAlgorithm::ChaCha)),
|
||||
false => TestRunner::new(config),
|
||||
};
|
||||
runner.run(&strategy, test).map_err(|e| e.to_string())
|
||||
}
|
||||
|
||||
/// A rival must fail on the property's own assertion, never on something else.
|
||||
pub fn assert_caught(result: Result<(), String>, marker: &str) {
|
||||
let err = result.expect_err("the rival passed: the property cannot see it");
|
||||
assert!(err.contains(marker), "failed for another reason: {err}");
|
||||
}
|
||||
|
||||
pub fn oracle_cmp(a: &Num, b: &Num) -> Option<Ordering> {
|
||||
Some(q(a)?.cmp(&q(b)?))
|
||||
}
|
||||
|
||||
pub fn integer(v: &BigRational) -> bool {
|
||||
v.is_integer() && !v.numer().is_zero()
|
||||
}
|
||||
File diff suppressed because it is too large
Load Diff
Reference in New Issue
Block a user